"""profile_limit.py -- the limit profile of Theorem thm:profile.

For squarefree m>=2 and primes p -> infinity, the family (r,s,E)=(p,mp,p)
samples xi = k/(mp), N = pk, and

    W_{p,mp}(k/(mp))  ->  W_inf(k) = pi*cot(pi k/m)/k - m/k^2   ((k,m)=1),

so Sigma*(p,mp,p) -> Sigma*_inf(m) = sum_{(k,m)=1} |W_inf(k)|/k.
Closed forms (the L(3,chi) cancellation, phi(m)=2):
    Sigma*_inf(3) = 6 - 2*pi/sqrt(3) = 2.3724012...   <- conjectured sup
    Sigma*_inf(6) = 12 - 2*sqrt(3)*pi
    Sigma*_inf(2) = (7/4) zeta(3)                      (cot vanishes)
and formally (m=4 not squarefree) 8 - 2*pi.

This script: (1) checks the closed forms against the direct series,
(2) computes the family Sigma*(p, 3p, p) through the collapsed cotangent
formula of Theorem thm:collapse and watches it climb toward 6-2*pi/sqrt(3),
(3) scans Sigma*_inf(m) for m <= 40 (maximum at m=3).
"""
import numpy as np
from math import gcd, pi, sqrt, tan, log


def cot(x):
    return 1.0 / tan(x)


def W_collapsed(r, s, kn, kd):
    """Exact collapsed cotangent evaluation of W_{r,s}(kn/kd)."""
    g = gcd(r, s); r1, s1 = r // g, s // g
    M = r * s
    N = M * kn // kd
    if gcd(N, r1 * s1) != 1:
        return 0.0
    tot = 0.0
    for a in range(1, r):
        if gcd(a, r) != 1: continue
        if (a * s - N) % r: continue
        v = (N - a * s) // r
        if g > 1 and gcd(v % g if v % g else g, g) != 1: continue
        tot += cot(pi * a / r)
    for b in range(1, s):
        if gcd(b, s) != 1: continue
        if (b * r - N) % s: continue
        u = (N - b * r) // s
        if g > 1 and gcd(u % g if u % g else g, g) != 1: continue
        tot += cot(pi * b / s)
    return (pi / N) * tot


def sigma_star(r, s, E, K=3000):
    g = gcd(r, s); r1, s1 = r // g, s // g
    den = r1 * s1 * E
    head = sum(abs(W_collapsed(r, s, k, den)) / k
               for k in range(1, K + 1) if (r * s * k) % den == 0)
    tail = 2 * (1 + log(r * s)) * (r + s) * den / (r * s) / K
    return head, tail


def sigma_inf(m, K=2_000_000):
    return sum(abs((pi / k) * cot(pi * k / m) - m / k**2) / k
               for k in range(1, K + 1) if gcd(k, m) == 1)


if __name__ == "__main__":
    SUP = 6 - 2 * pi / sqrt(3)
    z3 = sum(1.0 / n**3 for n in range(1, 2_000_000))
    print("== closed forms vs direct series ==")
    for m, cf, name in [(2, 1.75 * z3, "(7/4)zeta(3)"),
                        (3, SUP, "6-2pi/sqrt(3)"),
                        (4, 8 - 2 * pi, "8-2pi (formal)"),
                        (6, 12 - 2 * sqrt(3) * pi, "12-2sqrt(3)pi")]:
        v = sigma_inf(m)
        print(f"m={m}: series={v:.9f}  {name}={cf:.9f}  diff={abs(v-cf):.2e}")

    print("\n== the family (p, 3p, p): climb toward the sup ==")
    for p in [5, 11, 31, 101, 211]:
        h, t = sigma_star(p, 3 * p, p)
        print(f"p={p:4d}: Sigma* >= {h:.5f}  (tail bound {t:.4f});  sup-h = {SUP-h:+.5f}")

    print("\n== Sigma*_inf(m), m <= 40: max at m=3 ==")
    vals = {m: sigma_inf(m, K=400_000) for m in range(2, 41)}
    top = sorted(vals, key=vals.get, reverse=True)[:6]
    for m in top:
        print(f"m={m:2d}: {vals[m]:.6f}")