"""wlaw.py -- true size of W_{r,s}(xi) at the frequencies that matter.

W_{r,s}(xi) = sum over u/r + v/s = xi, (u,r)=(v,s)=1 of 1/(uv).

Computation route A (reference, slow): direct truncated double sum with
symmetric windows and tail estimate.
Computation route B (fast, exact structure): via residue classes mod M=rs
and the cotangent formula
   W = (1/N) [ r*(pi/M) sum_{a in A_u} cot(pi a/M) + s*(pi/M) sum_{b in A_v} cot(pi b/M) ]
with N = rs*xi, A_u = {a mod M: (a,r)=1, a*s === N (mod r), ((N-a*s)/r, s)=1},
A_v the mirror. Validate B against A, then scan the law with B.
"""
import numpy as np
from math import gcd, pi, lcm
import sys

def sqfree(n):
    p = 2; m = n
    while p * p <= m:
        if m % p == 0:
            m //= p
            if m % p == 0: return False
        p += 1
    return True

def W_direct(r, s, xi_num, xi_den, U=2_000_000):
    """xi = xi_num/xi_den; direct sum over |u|<=U."""
    N = r * s * xi_num // xi_den
    assert r * s * xi_num % xi_den == 0
    tot = 0.0
    for u in range(-U, U + 1):
        if u == 0 or gcd(u, r) != 1: continue
        num = N - u * s
        if num % r: continue
        v = num // r
        if v == 0 or gcd(v, s) != 1: continue
        tot += 1.0 / (u * v)
    return tot

def W_cot(r, s, xi_num, xi_den):
    """exact-structure evaluation via residue classes mod rs."""
    M = r * s
    N = M * xi_num // xi_den
    su = 0.0; sv = 0.0
    for a in range(M):
        if gcd(a, r) != 1: continue
        num = N - a * s
        if num % r: continue
        vmod = (num // r) % s   # careful: v mod s determined by a mod M
        if gcd(vmod, s) != 1: continue
        # a contributes to the u-sum
        if a % M != 0:
            su += 1.0 / np.tan(pi * a / M) if a % M else 0.0
    for b in range(M):
        if gcd(b, s) != 1: continue
        num = N - b * r
        if num % s: continue
        umod = (num // s) % r
        if gcd(umod, r) != 1: continue
        if b % M != 0:
            sv += 1.0 / np.tan(pi * b / M)
    return (1.0 / N) * (r * (pi / M) * su + s * (pi / M) * sv)

if __name__ == "__main__":
    print("== validation: W_cot vs W_direct ==")
    for (r, s, kn, kd) in [(10, 15, 1, 6), (10, 15, 5, 6), (6, 15, 1, 10), (21, 35, 1, 15),
                            (10, 15, 1, 30), (14, 35, 3, 10), (22, 33, 1, 6)]:
        g = gcd(r, s)
        wd = W_direct(r, s, kn, kd, U=400000)
        wc = W_cot(r, s, kn, kd)
        print(f"r={r:3d} s={s:3d} xi={kn}/{kd}: direct={wd:+.8f}  cot={wc:+.8f}  diff={abs(wd-wc):.2e}")