"""sectors.py -- exact sector decomposition X = X31+X22+X211+X1111 over the lcm of each multiset,
with a valuation prefilter; reproduces the sector table and the exact jumps Delta X(21), Delta X(22).
"""
from fractions import Fraction as F
from math import lcm, gcd
from itertools import combinations

def mu(n):
    if n==1: return 1
    r=1; m=n; p=2
    while p*p<=m:
        if m%p==0:
            m//=p
            if m%p==0: return 0
            r=-r
        p+=1
    if m>1: r=-r
    return r

def mq(q,Y):
    return sum(F(mu(m),m) for m in range(1,int(Y)+1) if gcd(m,q)==1 and mu(m)!=0)

def sigvals(q):
    """values of sigma_q(x/q) on (j,j+1), j=0..q-1 (period q)"""
    out=[]
    for j in range(q):
        x=F(2*j+1,2); v=F(0)
        for t in range(1,q+1):
            if q%t==0 and mu(t)!=0:
                y=F(t)*x/q
                v+=F(mu(t),t)*(y-(y.numerator//y.denominator)-F(1,2))
        out.append(v)
    return out

def avg_prod(conds_powers):
    """exact long-run average of prod sigma_{q}(x/q)^e over its own period"""
    qs=[q for q,e in conds_powers]
    L=lcm(*qs)
    SV={q:sigvals(q) for q in qs}
    tot=F(0)
    for j in range(L):
        p=F(1)
        for q,e in conds_powers:
            p*=SV[q][j%q]**e
        tot+=p
    return tot/L

from fractions import Fraction as F
from math import gcd
from itertools import combinations

def pfac(n):
    s=set(); m=n; p=2
    while p*p<=m:
        if m%p==0:
            s.add(p)
            while m%p==0: m//=p
        p+=1
    if m>1: s.add(m)
    return s

def valuation_ok(ms):
    """every prime divides 0 or >=2 of the conductors (with multiplicity)"""
    from collections import Counter
    cnt=Counter()
    for q in ms:
        for p in pfac(q): cnt[p]+=1
    return all(c>=2 for c in cnt.values())

CACHE={}
def M(ms):
    key=tuple(sorted(ms))
    if key in CACHE: return CACHE[key]
    if not valuation_ok(ms):
        CACHE[key]=F(0); return F(0)
    from collections import Counter
    cp=[(q,e) for q,e in Counter(ms).items()]
    r=avg_prod(cp); CACHE[key]=r; return r

def X_sectors(D):
    qs=[q for q in range(2,D+1) if mu(q)!=0]
    c={q: mu(q)*mq(q,D//q) for q in qs}
    v={q: M([q,q]) for q in qs}
    X31=sum(4*c[q]**3*c[r]*M([q,q,q,r]) for q in qs for r in qs if r!=q)
    X22=sum(6*c[q]**2*c[r]**2*(M([q,q,r,r])-v[q]*v[r]) for q,r in combinations(qs,2))
    X211=sum(12*c[q]**2*c[r]*c[s]*M([q,q,r,s]) for q in qs for r,s in combinations([x for x in qs if x!=q],2))
    X1111=sum(24*c[a]*c[b]*c[cc]*c[d]*M([a,b,cc,d]) for a,b,cc,d in combinations(qs,4))
    return X31,X22,X211,X1111

res={}
for D in range(15,23):
    if mu(D)==0 and D not in (16,18,20):  pass
    a,b,cc,d=X_sectors(D)
    res[D]=(a,b,cc,d,a+b+cc+d)
    print(f"D={D}: X31={float(a):+.10f} X22={float(b):+.10f} X211={float(cc):+.10f} X1111={float(d):+.10f}  X={float(a+b+cc+d):+.10f}")

dx21=res[21][4]-res[20][4]; dx22=res[22][4]-res[21][4]
print("\nDelta X(21) =",dx21,"=",float(dx21))
print("Delta X(22) =",dx22,"=",float(dx22))
p=11
pure=F(-(p*p-1),32*p**4)
print("pure 2p law at p=11:",pure,"=",float(pure))
print("correction =",dx22-pure,"=",float(dx22-pure))
# sector split of the jump at 22
for i,name in enumerate(["X31","X22","X211","X1111"]):
    d=res[22][i]-res[21][i]
    print(f"  jump in {name}: {d} = {float(d):+.10f}")