# SCRIPT: BH-KNOT-14
# The (3,4) torus knot 8_19 beside the trefoil T(2,3): quantum invariant at q = exp(2 pi i/N), its growth
# and its mirror. Also the classical facts quoted for 8_19.
# Method (Rosso-Jones, derived from the Adams operation, no external table):
#   For T(m,p), colour = the N-dimensional irreducible representation, quantum integer [d] = (q^{d/2}-q^{-d/2})/(q^{1/2}-q^{-1/2}),
#   twist theta_d = q^{(d^2-1)/4}. The Adams operation psi^m sends the character sum_k x^{N-1-2k} to
#   sum_k x^{a_k}, a_k = m(N-1-2k). This sum is symmetric under x -> 1/x, and chi_{a+1} - chi_{a-1} = x^a + x^{-a}
#   (chi_{-d} = -chi_d), so sum_k x^{a_k} = (1/2) sum_k (chi_{a_k+1} - chi_{a_k-1}). Then
#     J_N(T(m,p); q) = theta_N^{-mp} [N]^{-1} (1/2) sum_k sum_{e=+-1} e theta_{a_k+e}^{p/m} [a_k+e],
#   normalised so the unknot is 1. Why this can be trusted here: it is checked below against three
#   independent facts before it is used: (i) m = 1 gives the unknot, J = 1; (ii) N = 2 gives the Jones
#   polynomials of the trefoil (q + q^3 - q^4) and of 8_19 (q^3 + q^5 - q^8), up to the mirror q -> 1/q;
#   (iii) for the trefoil, |J_N| at q = exp(2 pi i/N) equals |F| of the Kontsevich-Zagier series used in
#   BH-KNOT-08. At a root of unity [N] = 0, so J is evaluated as the limit q = exp(2 pi i/N + h), h -> 0,
#   at high precision, and the limit is checked to be stable in h.
# Kill conditions, written before the run:
#  K1  checks (i)-(iii) pass to 1e-20 (exact cases) and 1e-12 (root of unity). Kill otherwise (formula wrong).
#  K2  8_19 grows like a power, not exponentially: |J_N|/N^{3/2} changes by less than 2% from N = 400 to 800,
#      and (2 pi/N) log|J_N| -> 0 (volume zero): it stays below (2 pi/N)(1.5 log N + 2). Kill otherwise.
#  K5  (added after the first run showed |J_N|/N^{3/2} differs between N = 200 and N = 800): the constant depends
#      on N mod 16. Test: at N in [640, 688), |c(N) - c(N+16)| < 1e-3 for all N, while |c(N) - c(N+8)| and
#      |c(N) - c(N+12)| exceed 0.1 somewhere (period exactly 16 among 8, 12, 16). Kill otherwise.
#  K3  mirror: q -> 1/q conjugates J_N(8_19) at the root of unity, and the invariant is not real (chiral). Kill otherwise.
#  K4  classical facts: Alexander polynomial (t^12-1)(t-1)/((t^3-1)(t^4-1)) = t^6 - t^5 + t^3 - t + 1, genus 3
#      (half the span); Moser: 1/n surgery on T(p,q) is the Brieskorn sphere Sigma(p,q,|pqn-1|), giving
#      Sigma(2,3,5), Sigma(2,3,7) for the trefoil and Sigma(3,4,11), Sigma(3,4,13) for 8_19. Kill if any differs.
import os as _os
_HERE=_os.path.dirname(_os.path.abspath(__file__))
_FIG=_os.path.join(_HERE,'..','paper','fig') if _os.path.isdir(_os.path.join(_HERE,'..','paper','fig')) else _os.path.join(_HERE,'..','fig')
import mpmath as mp, sympy as sp
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
def qint(q,d): return (q**(mp.mpf(d)/2)-q**(-mp.mpf(d)/2))/(q**mp.mpf(0.5)-q**(-mp.mpf(0.5)))
def J(m,p,N,q):
    th=lambda d,e=1: q**(mp.mpf(e)*(d*d-1)/4)
    s=0
    for k in range(N):
        a=m*(N-1-2*k)
        for e in (1,-1):
            d=a+e
            if d==0: continue
            s+=e*q**(mp.mpf(p)*(d*d-1)/(4*m))*qint(q,d)
    return s/(2*th(N)**(m*p)*qint(q,N))
mp.mp.dps=50
q0=mp.mpf('1.3')+mp.mpf('0.2')*1j
e_unknot=max(abs(J(1,p,N,q0)-1) for p in (2,3,5) for N in (2,3,6))
j31=J(2,3,2,q0); j819=J(3,4,2,q0)
jones31=[q0+q0**3-q0**4, 1/q0+1/q0**3-1/q0**4]; jones819=[q0**3+q0**5-q0**8, q0**-3+q0**-5-q0**-8]
e31=min(abs(j31-x) for x in jones31); e819=min(abs(j819-x) for x in jones819)
mirror_tag='q' if abs(j819-jones819[0])<abs(j819-jones819[1]) else '1/q'
def kz(N):
    mp.mp.dps=int(40+0.1*N); q=mp.exp(2j*mp.pi/N); s=0; p=mp.mpc(1)
    for n in range(N):
        s+=p; p*=(1-q**(n+1))
    return s
def Jroot(m,p,N,mirror=False,h=None):
    mp.mp.dps=int(60+0.25*N)
    h=h if h is not None else mp.mpf(10)**(-(mp.mp.dps//2))
    sgn=-1 if mirror else 1
    return J(m,p,N,mp.exp(sgn*2j*mp.pi/N+h))
etr=0
for N in (20,50,101):
    a=abs(Jroot(2,3,N)); b=abs(kz(N)); etr=max(etr,abs(a-b)/b)
st=abs(Jroot(3,4,60)-Jroot(3,4,60,h=mp.mpf(10)**-20))/abs(Jroot(3,4,60))
report('K1 formula checks: unknot, Jones N=2, trefoil vs Kontsevich-Zagier, h-stability',
       e_unknot<1e-20 and e31<1e-20 and e819<1e-20 and etr<1e-12 and st<1e-12,
       f'unknot {float(e_unknot):.1e}; Jones 3_1 {float(e31):.1e}, 8_19 {float(e819):.1e} (formula gives 8_19 in variable {mirror_tag}); |J|/|F| trefoil {float(etr):.1e}; h-stability {float(st):.1e}')
res={}
for N in (100,200,400,800):
    v=Jroot(3,4,N); res[N]=v
    print(f'   N={N:4d}: |J_N(8_19)|/N^1.5 = {float(abs(v)/N**1.5):.6f}   (2pi/N) log|J_N| = {float(2*mp.pi/N*mp.log(abs(v))):.6f}   arg = {float(mp.arg(v)/mp.pi):+.6f} pi')
r4=abs(res[400])/400**1.5; r8=abs(res[800])/800**1.5
report('K2 8_19 grows like N^{3/2}, no volume', abs(r8-r4)/r8<0.02 and float(2*mp.pi/800*mp.log(abs(res[800])))<float(2*mp.pi/800*(1.5*mp.log(800)+2)), f'{float(r4):.6f} -> {float(r8):.6f}')
N=200; v=Jroot(3,4,N); vm=Jroot(3,4,N,mirror=True)
report('K3 mirror conjugates 8_19; invariant not real', abs(vm-mp.conj(v))<1e-12*abs(v) and abs(mp.im(v))>1e-3*abs(v), f'J = {complex(v):.4f}, mirror = {complex(vm):.4f}')
t=sp.symbols('t')
A=sp.expand(sp.cancel((t**12-1)*(t-1)/((t**3-1)*(t**4-1))))
genus=sp.degree(A,t)//2
mos=lambda p,q,n: abs(p*q*n-1)
report('K4 Alexander, genus, Moser surgeries', A==t**6-t**5+t**3-t+1 and genus==3 and (mos(2,3,1),mos(2,3,-1),mos(3,4,1),mos(3,4,-1))==(5,7,11,13),
       f'Delta = {A}, genus {genus}; Sigma(2,3,{mos(2,3,1)}), Sigma(2,3,{mos(2,3,-1)}), Sigma(3,4,{mos(3,4,1)}), Sigma(3,4,{mos(3,4,-1)})')
def c(N):
    mp.mp.dps=int(80+0.3*N); h=mp.mpf(10)**(-int(40+0.15*N)); return float(abs(J(3,4,N,mp.exp(2j*mp.pi/N+h)))/N**1.5)
A={N:c(N) for N in range(640,688)}
d={per:max(abs(A[N]-A[N+per]) for N in range(640,688-per)) for per in (8,12,16)}
cls=[A[640+r] for r in range(16)]
print('   period test: max|c(N)-c(N+P)| =', {k:round(v,6) for k,v in d.items()})
print('   constants by N mod 16 (N = 640..655):', [round(x,4) for x in cls])
report('K5 the N^{3/2} constant of 8_19 has period 16 in N', d[16]<1e-3 and d[8]>0.1 and d[12]>0.1, f'range {min(cls):.4f} to {max(cls):.4f}')
import json
json.dump({'res':{str(N):[float(abs(res[N])),float(mp.arg(res[N]))] for N in res},'c640':A},open(_os.path.join(_HERE,'bh_knot_14.json'),'w'))
print('ALL PASS' if ok else 'SOME FAILED')