"""ch11_false_position.py

Find one root of f(x) = 0 by the false position (regula falsi) method.

Origin:
    Python version of ch11_false_position.f90 by Alexander Godunov.
    The interpolation formula and convergence criterion match that file.

Example:
    f(x) = sqrt(x) - cos(x) on [0, 1].
    Reference root = 0.64171437087288... .

The program reports the iteration count, root estimate, and f(root).
The requested tolerance applies to the residual abs(f(root)). Unlike
bisection, ordinary false position may leave one endpoint fixed for
many iterations, so the bracket width is not the stopping test.
"""

import math


def main():
    """Set up the example and print its result."""
    # Initial bracket and requested residual tolerance
    x1 = 0.0
    x2 = 1.0
    eps = 1.0e-6

    print("A single root of f(x)")
    print("Method - False position")

    root, flag = false_position(f, x1, x2, eps)

    # Report a missing bracket or invalid input before printing a root
    if flag == 0:
        print("invalid interval or no initial sign change")
        return

    print(f"iterations         = {abs(flag)}")
    print(f"residual tolerance = {eps:.5e}")
    if flag < 0:
        print("iteration limit reached before the requested residual")

    print(f"root               = {root:.5e}")
    # Check the residual of the reported estimate
    print(f"f(root)            = {f(root):.5e}")


def false_position(fun, x1, x2, eps):
    """Locate one root of fun(x) on [x1, x2] by ordinary false position.

    Input:
        fun  - function to evaluate, continuous on the initial bracket
        x1   - left endpoint of the initial bracket
        x2   - right endpoint of the initial bracket
        eps  - requested maximum residual abs(fun(root)) (positive)

    Output:
        root - most recent false-position estimate when flag is nonzero
        flag - convergence indicator and iteration count
               flag > 0: abs(fun(root)) <= eps; flag is the count
               flag = 0: invalid interval/tolerance or no sign change
               flag < 0: iteration limit reached; abs(flag) is the count

    Notes:
        A sign change guarantees a root only when fun is continuous.
        The maximum number of iterations is 200. One bracket endpoint
        can remain nearly fixed, so a small interval width is not assumed.
        An endpoint root is returned immediately with a positive flag.
    """
    max_iterations = 200

    # Check that the inputs describe a usable initial bracket
    if x1 >= x2 or eps <= 0.0:
        return math.nan, 0

    a = x1
    b = x2
    fa = fun(a)
    fb = fun(b)
    if fa * fb > 0.0:
        return math.nan, 0

    # Handle an exact zero at an endpoint before forming a secant line.
    # Flag 1 follows the existing convention that positive means success.
    if fa == 0.0:
        return a, 1
    if fb == 0.0:
        return b, 1

    # Interpolate between the current endpoints and retain the sign change.
    for i in range(1, max_iterations + 1):
        c = b - fb * (b - a) / (fb - fa)
        fc = fun(c)

        # False position uses the residual, not the bracket width, to stop.
        if abs(fc) <= eps:
            return c, i

        if fa * fc <= 0.0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc

    return c, -max_iterations


def f(x):
    """Example function; choose a bracket in its domain x >= 0."""
    # To try the alternative from the original programs, use [-1, 0].
    # return x + math.cos(x)
    return math.sqrt(x) - math.cos(x)


if __name__ == "__main__":
    main()