#!/usr/bin/env python3
"""scan_C82.py -- independent direct evaluation (NOT a proof) of the auxiliary
function A(z) of the note on the two circles |z| = 1 and |z-1| = 1, from the
certificate file CERTIFICATE_C82_ZhangZagier_FINAL.json.  NumPy only.
Polynomials are evaluated from their integer coefficients in z (numpy.polyval);
no code is shared with prove_C82.py.  Prints the hypotheses of the certificate
and the minimum of A over 400000 equally spaced points per circle.
David Niedbala Giraudin, ORCID 0009-0009-1526-1178.
Technical assistance: Claude (Anthropic).
"""
import json, sys
import numpy as np

J = json.load(open(sys.argv[1] if len(sys.argv) > 1 else "CERTIFICATE_C82_ZhangZagier_FINAL.json"))
l0 = J["lambda0"]
C = np.array(J["eta"]["centres_re"]) + 1j * np.array(J["eta"]["centres_im"])
R = np.array(J["eta"]["radii"])
M = np.array(J["eta"]["masses"])
P = [(p["multiplier"], np.array(p["coeffs_desc"], float)) for p in J["polynomials"]]
print("budget", sum(c * (len(q) - 1) for c, q in P) + 2 * l0, " total mass", M.sum(),
      " min multiplier", min(c for c, _ in P))
D = np.abs(C[:, None] - C[None, :])
np.fill_diagonal(D, np.inf)
print("disc disjointness margin", (D - (R[:, None] + R[None, :])).min())
D2 = D.copy()
np.fill_diagonal(D2, 1.0)
K = np.log(D2)
np.fill_diagonal(K, np.log(R) - 0.25)
Ie = M @ K @ M
print("I(eta)", Ie, " (file:", J["I_eta"], ")")
print("density bound D = max m/(pi r^2) =", (M / (np.pi * R * R)).max())


def A(z):
    v = np.log(np.maximum(1, np.abs(z))) + np.log(np.maximum(1, np.abs(1 - z)))
    with np.errstate(divide="ignore"):
        for c, q in P:
            v = v - c * np.log(np.abs(np.polyval(q, z)))
        U = np.zeros(z.shape)
        for s in range(0, len(C), 500):
            d = np.abs(z[:, None] - C[None, s:s + 500])
            r = R[None, s:s + 500]
            U = U + np.where(d >= r, np.log(d), np.log(r) + (d * d / (r * r) - 1) / 2) @ M[s:s + 500]
    return v - 2 * l0 * U + l0 * Ie


best = np.inf
for shift in (0.0, 1.0):
    th = np.linspace(0, 2 * np.pi, 400001)[:-1]
    for s in range(0, len(th), 20000):
        a = A(shift + np.exp(1j * th[s:s + 20000]))
        best = min(best, float(np.nanmin(a)))
    print("circle |z -", shift, "| = 1 done")
print("minimum of A over the sample points:", best)