1import json
 2import math
 3import numpy as np
 4
 5
 6def fold(p, j):
 7    """Apply the stated local rational fold to complex polygon coordinate j."""
 8    q = np.array(p, copy=True)
 9    n = len(q)
10    a, b, c = p[(j - 1) % n], p[j], p[(j + 1) % n]
11    den = (b - c) - (a - b)  # 2*b-a-c
12    q[j] = ((b - c) * a - (a - b) * c) / den
13    return q, abs(den)
14
15
16def apply_schedule(p, schedule):
17    q = np.array(p, copy=True)
18    denoms = []
19    for j in schedule:
20        q, d = fold(q, int(j))
21        denoms.append(d)
22    return q, np.asarray(denoms)
23
24
25def relerr(a, b):
26    return float(np.linalg.norm(a - b) / max(np.linalg.norm(b), 1e-30))
27
28
29def polygon(seed, n=8, dtype=np.complex128):
30    rng = np.random.default_rng(seed)
31    # A separated, mildly perturbed convex polygon avoids accidental singularities.
32    theta = np.arange(n) * (2 * np.pi / n)
33    radius = 1.0 + .04 * rng.normal(size=n)
34    p = radius * np.exp(1j * theta) + .02 * (rng.normal(size=n) + 1j*rng.normal(size=n))
35    return p.astype(dtype)
36
37
38def fold_math_check():
39    rng = np.random.default_rng(10)
40    errs = []
41    for _ in range(2000):
42        a, c = rng.normal(size=2) + 1j*rng.normal(size=2)
43        b = rng.normal() + 1j*rng.normal()
44        if abs(2*b-a-c) < .1:
45            continue
46        q, _ = fold(np.array([a, b, c], dtype=np.complex128), 1)
47        r, _ = fold(q, 1)
48        errs.append(abs(r[1]-b) / max(abs(b), 1e-12))
49    return {"samples": len(errs), "max_relative_error": max(errs), "median_relative_error": float(np.median(errs))}
50
51
52def schedule_check(dtype, schedule, repeats=1):
53    p = polygon(3, dtype=dtype)
54    sched = list(schedule) * repeats
55    y, den = apply_schedule(p, sched)
56    z, den2 = apply_schedule(y, list(reversed(sched)))
57    return {"dtype": str(dtype), "steps": len(sched), "reconstruction_relative_error": relerr(z, p),
58            "minimum_denominator": float(min(np.min(den), np.min(den2))),
59            "max_abs_state": float(max(np.max(np.abs(p)), np.max(np.abs(y)), np.max(np.abs(z))))}
60
61
62def denominator_sweep():
63    # Directly approach the singular locus 2b-a-c=0 while retaining nonzero a,b,c.
64    a, c = 0.2 + .4j, -0.7 + .1j
65    out = []
66    for eps in [1e-1, 1e-3, 1e-5, 1e-7, 1e-9]:
67        b = (a+c)/2 + eps*(1+0.3j)
68        p = np.array([a,b,c], dtype=np.complex128)
69        q, d = fold(p, 1)
70        r, _ = fold(q, 1)
71        out.append({"epsilon": eps, "denominator": float(d),
72                    "folded_abs": float(abs(q[1])), "involution_error": float(abs(r[1]-b))})
73    return out
74
75
76def main():
77    schedule = [0, 2, 4, 6, 1, 3, 5, 7]
78    results = {
79        "math_check": fold_math_check(),
80        "fixed_schedule": [schedule_check(np.complex128, schedule, k) for k in [1, 5, 20]],
81        "float32_schedule": [schedule_check(np.complex64, schedule, k) for k in [1, 5, 20]],
82        "singularity_sweep": denominator_sweep(),
83        "memory_note": {"reversible_saved_fold_states": 0,
84                        "ordinary_unroll_saved_fold_states_20_steps": 20,
85                        "state_size_points": 8},
86    }
87    with open("results.json", "w") as f:
88        json.dump(results, f, indent=2)
89    print(json.dumps(results, indent=2))
90
91if __name__ == "__main__":
92    main()