"""Revision point 1 (2026-09-22): the order of the base rule's step, measured where it can be.

Near a simple root the base rule's step is Newton's step with the derivative of the parabola
through the last three points (tent similarity, step shorter than the distance to the far end,
projection inactive).  The error of that derivative is f'''/6 (x_n - x_{n-1})(x_n - x_{n-2}),
i.e. O(e_{n-1} e_{n-2}), so e_{n+1} ~ e_n e_{n-1} e_{n-2} and the order is the root of
p^3 = p^2 + p + 1, 1.839 - Muller's order.  On a quadratic f the parabola is exact and the
order is Newton's 2.  In double precision a run has 3-4 superlinear steps before the error
reaches the tolerance, too few to tell 1.84 from 2; here the iteration is carried to 10^-2900
in 3000-digit arithmetic and the computational order e_{n+1} ~ e_n^p is read per step.
The two-point (secant) scale is measured the same way for reference.

    python -u order_high_precision.py   -> results_order_high_precision.txt (stdout)
"""
from __future__ import annotations

import mpmath as mp

mp.mp.dps = 3000
FLOOR = mp.mpf(10) ** -2900


def iterate(f, x0, x1, x2, points, steps=60):
    xs = [mp.mpf(x0), mp.mpf(x1), mp.mpf(x2)]
    for _ in range(steps):
        a, b, c = xs[-3], xs[-2], xs[-1]
        fb, fc = f(b), f(c)
        d12 = (fc - fb) / (c - b)
        if points == 3:
            d01 = (fb - f(a)) / (b - a)
            slope = d12 + (d12 - d01) / (c - a) * (c - b)
        else:
            slope = d12
        xs.append(c - fc / slope)
        if xs[-1] == c:
            break
    return xs


def orders(xs, root):
    """Computational order per step, from errors above the precision floor only."""
    e = [abs(x - root) for x in xs]
    e = [v for v in e if v > FLOOR]
    return [mp.log(e[i + 1]) / mp.log(e[i]) for i in range(3, len(e) - 1)]


CASES = (
    ("cos x - x", lambda x: mp.cos(x) - x, (0.5, 0.9, 0.7), lambda: mp.findroot(lambda x: mp.cos(x) - x, 0.74)),
    ("exp(x) - 2", lambda x: mp.exp(x) - 2, (0.5, 1.0, 0.8), lambda: mp.log(2)),
    ("x^3 - x - 2", lambda x: x ** 3 - x - 2, (1.3, 1.7, 1.5), lambda: mp.findroot(lambda x: x ** 3 - x - 2, 1.52)),
    ("sin x - 1/2", lambda x: mp.sin(x) - mp.mpf(1) / 2, (0.3, 0.7, 0.5), lambda: mp.asin(mp.mpf(1) / 2)),
    ("x^2 - 2", lambda x: x ** 2 - 2, (1.2, 1.6, 1.3), lambda: mp.sqrt(2)),
)


def main():
    muller = mp.findroot(lambda p: p ** 3 - p ** 2 - p - 1, 1.8)
    golden = (1 + mp.sqrt(5)) / 2
    print(f"Muller order (root of p^3 = p^2 + p + 1): {float(muller):.4f}")
    print(f"secant order (golden ratio):              {float(golden):.4f}")
    print(f"{'function':<14}{'three-point: last orders':<40}{'two-point: last orders'}")
    for name, f, (a, b, c), root in CASES:
        r = root()
        o3 = orders(iterate(f, a, b, c, 3), r)
        o2 = orders(iterate(f, a, b, c, 2), r)
        print(f"{name:<14}{' '.join(f'{float(v):.3f}' for v in o3[-3:]):<40}"
              f"{' '.join(f'{float(v):.3f}' for v in o2[-3:])}")
    print("order high precision done")


if __name__ == "__main__":
    main()