"""
Feigenbaum's delta for the one-parameter family

        f_{r,p}(x) = r * (1 - |2x - 1|^p)      on x in [0,1]

- p = 2 reduces EXACTLY to the logistic map:  f = r(1-(2x-1)^2) = 4 r x(1-x).
  Known first Feigenbaum constant (quadratic maximum): delta = 4.66920160910299...
- p = 4 (quartic maximum) is a DIFFERENT universality class:  delta ~ 7.2846862...
- general p: the map's maximum at x=1/2 is of order p; delta depends only on that order
  and varies smoothly with p.

We locate the SUPERSTABLE parameters R_n -- the r at which the critical point x=1/2
lies on a period-2^n cycle, i.e.  f^{2^n}(1/2; r) = 1/2 -- and estimate

        delta_n = (R_n - R_{n-1}) / (R_{n+1} - R_n)  ->  delta.

Root-finding is bracketed bisection (robust), in mpmath at high precision.

Outputs:
  delta_of_p.json   -- {p: delta} on a fine grid (the reference curve), plus the
                       high-precision headline values for p=2 and p=4 with the
                       literature comparison.

Run:  python3 compute.py
"""
import json, os, mpmath as mp

def make_f(p):
    p = mp.mpf(p)
    def f(x, r):
        return r * (1 - abs(2*x - 1)**p)
    return f

def residual(r, f, n):
    x = mp.mpf('0.5')
    for _ in range(2**n):
        x = f(x, r)
    return x - mp.mpf('0.5')

def find_superstable(f, n, lo, hi, scan, tol_exp):
    """First sign change of the residual in (lo, hi], then bisect."""
    step = (hi - lo) / scan
    prev = residual(lo, f, n)
    r = lo
    for i in range(1, scan + 1):
        cur_r = lo + step * i
        cur = residual(cur_r, f, n)
        if prev == 0:
            return r
        if (prev < 0) != (cur < 0):
            a, b, fa = r, cur_r, prev
            for _ in range(300):
                m = (a + b) / 2
                fm = residual(m, f, n)
                if fm == 0:
                    return m
                if (fa < 0) != (fm < 0):
                    b = m
                else:
                    a, fa = m, fm
                if b - a < mp.mpf(10) ** (-tol_exp):
                    break
            return (a + b) / 2
        r, prev = cur_r, cur
    raise ValueError(f"no sign change for n={n} in ({lo},{hi})")

def deltas_for_p(p, nmax, scan=300, tol_exp=40):
    f = make_f(p)
    R = {0: mp.mpf('0.5')}
    R[1] = find_superstable(f, 1, mp.mpf('0.55'), mp.mpf('0.97'), scan, tol_exp)
    for n in range(2, nmax + 1):
        step = R[n-1] - R[n-2]
        R[n] = find_superstable(f, n, R[n-1] + step*mp.mpf('0.02'),
                                R[n-1] + step*mp.mpf('0.97'), scan, tol_exp)
    ds = []
    for n in range(1, nmax):
        ds.append((R[n+1] - R[n-1] - (R[n+1]-R[n]),  # = R[n]-R[n-1]
                   (R[n] - R[n-1]) / (R[n+1] - R[n])))
    return [d for _, d in ds]

# Literature values (see README for sources)
LIT = {
    2: mp.mpf('4.66920160910299067185320382046620161725'),
    4: mp.mpf('7.284686217656525'),
    6: mp.mpf('9.296246158651231'),
}

def main():
    out = {"family": "f(x) = r*(1 - |2x-1|^p) on [0,1]; p=2 is the logistic map",
           "method": "superstable cascade R_n: f^{2^n}(1/2)=1/2; delta=(R_n-R_{n-1})/(R_{n+1}-R_n)",
           "headline": {}, "curve": []}

    # Headline high-precision values
    mp.mp.dps = 55
    for p in (2, 4, 6):
        ds = deltas_for_p(p, 13, scan=400, tol_exp=48)
        best = ds[-1]
        lit = LIT[p]
        out["headline"][str(p)] = {
            "computed": mp.nstr(best, 13),
            "literature": mp.nstr(lit, 16),
            "abs_diff": mp.nstr(abs(best - lit), 3),
            "n": 13,
        }
        print(f"p={p}: computed={mp.nstr(best,12)} lit={mp.nstr(lit,12)} diff={mp.nstr(abs(best-lit),3)}")

    # Reference curve delta(p) on a fine grid (4-digit-ish; for plotting)
    mp.mp.dps = 35
    grid = [round(1.5 + 0.1*i, 2) for i in range(0, 36)]  # 1.5 .. 5.0
    for p in grid:
        try:
            ds = deltas_for_p(p, 11, scan=260, tol_exp=30)
            d = float(ds[-1])
            out["curve"].append({"p": p, "delta": round(d, 6)})
            print(f"  p={p:>4}  delta={d:.6f}")
        except Exception as e:
            print(f"  p={p:>4}  FAILED: {e}")

    here = os.path.dirname(os.path.abspath(__file__))
    with open(os.path.join(here, "delta_of_p.json"), "w") as fh:
        json.dump(out, fh, indent=2)
    print("\nwrote delta_of_p.json")

if __name__ == "__main__":
    main()