"""Independent checks for Odd period polynomials -- numberdb.org/T696. Run: sage -python checks.py. Needs generate.py and optional T694-roster.json. The Sage exact generator is compared with PARI's independent modular-form construction and numerical mfperiodpol, and with explicit split Mellin sums. Numerical agreement is a check, not the source of the exact integer entries. """ import json,math,hashlib from pathlib import Path from functools import lru_cache import numberdb.sage as numberdb from sage.libs.pari import pari from sage.rings.rational_field import QQ from sage.rings.integer_ring import ZZ from sage.rings.real_mpfr import RealField from sage.rings.complex_mpfr import ComplexField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.qqbar import AA from sage.arith.misc import euler_phi from sage.env import SAGE_VERSION from generate import OddPeriodPolynomials,space_data,selected_spaces,exact_polynomial,QX,X from check import exactness,measure,names_its_rings BASE=Path(__file__).resolve().parent P=PolynomialRing(QQ,'x') def pol(c):return P(str(pari.lift(c)).replace('y','x')) @lru_cache(None) def dimensions(N,k): primes=ZZ(N).prime_divisors() mu=QQ(N) for p in primes:mu*=1+QQ(1)/p e2=0 if N%4==0 else math.prod(1+int(ZZ(-4).kronecker(p)) for p in primes) e3=0 if N%2==0 or N%9==0 else math.prod(1+int(ZZ(-3).kronecker(p)) for p in primes) cusps=sum(euler_phi(math.gcd(d,N//d)) for d in ZZ(N).divisors()) genus=1+mu/12-QQ(e2)/4-QQ(e3)/3-QQ(cusps)/2 full=(k-1)*(genus-1)+(k//2-1)*cusps+e2*(k//4)+e3*(k//3) new=full-sum(len(ZZ(N//d).divisors())*dimensions(int(d),k)[1] for d in ZZ(N).divisors() if d=1) (1-q^m)^24, independent integer convolution. a=[1]+[0]*(nmax-1) for m in range(1,nmax): b=[0]*nmax for j in range(0,min(24,(nmax-1)//m)+1): c=(-1)**j*math.comb(24,j) for n in range(nmax-j*m):b[n+j*m]+=c*a[n] a=b return [0]+a def upper_mellin_weights(t, T, k): """Integrals of y^(s-1) exp(-t*y), T <= y < infinity, 1 <= s < k.""" R=t.parent();z=t*T;decay=(-z).exp() out=[R(0)]*k;partial=R(1);term=R(1) for s in range(1,k): if s>1:term*=z/(s-1);partial+=term out[s]=decay*math.factorial(s-1)*partial/t**s return out def mellin_coefficients(coeff,N,k,epsilon,working,terms): R=RealField(numberdb.bits(working,losing=32));T=1/R(N).sqrt();twopi=2*R.pi() integrals=[R(0) for s in range(k)] for n in range(1,terms+1): a=R(coeff[n]);weights=upper_mellin_weights(twopi*n,T,k) for s in range(1,k):integrals[s]+=a*weights[s] moments={s:integrals[s]+epsilon*R(N)**(k//2-s)*integrals[k-s] for s in range(1,k)} w=k-2 c=[R(0)]*(w+1) for j in range(1,w+1,2):c[w-j]=math.comb(w,j)*(-1)**((j-1)//2)*moments[j+1] assert c[w-1]!=0 return [v/c[w-1] for v in c] def exact_identities(p,N,k,fricke): w=k-2 assert p.degree()==w-1 and p.leading_coefficient()>0 assert math.gcd(*(int(c) for c in p.list()))==1 assert p(-X)==-p(X) reciprocal=sum(QQ(p[w-j])*QQ(N)**(j-w//2)*X**j for j in range(w+1)) assert QX(p)==fricke*reciprocal if N==1: manin=p+sum(p[j]*(X-1)**j*X**(w-j) for j in range(w+1))+sum(p[j]*(-1)**j*(X-1)**(w-j) for j in range(w+1)) assert manin==0 and p(1)==p(-1)==0 def run(values=None): # Nonzero controls run first. assert eta_delta(5)==[0,1,-24,252,-1472,4830] control=list(pari.mfcoefs(pari.mfDelta(),100)) assert [int(c) for c in control]==eta_delta(100) delta=4*X**9-25*X**7+42*X**5-25*X**3+4*X assert exact_polynomial(1,12,1)==delta R=RealField(numberdb.bits(150)) controls=upper_mellin_weights(R(2),R(3),4) # Elementary antiderivatives give exp(-6)/2, 7*exp(-6)/4, 25*exp(-6)/4. for j,answer in enumerate([R(1)/2,R(7)/4,R(25)/4],1): expected=answer*R(-6).exp() assert controls[j]>0 and abs(controls[j]/expected-1)RR('1e-20') norm=[c/coef[w-1] for c in coef] err=max(abs(norm[j]-CC(QQ(p[j])/p[w-1])) for j in range(w+1)) assert not err.is_NaN() and not err.is_infinity() and err