"""Independent character-trace continuations for numberdb.org/T633.
Adapted from the attached checks of T631 and T632, read from the live database.
This module imports no production generator or multiplicity functions.
"""
from functools import lru_cache
import mpmath as mp
mp.mp.dps=175

def mpq(a,b=1):
    return mp.mpf(int(a))/int(b)


@lru_cache(None)
def hz(j,n,p):
    a=mpq(n,p)
    return mp.zeta(-j,a),mp.zeta(-j,a,derivative=1)


@lru_cache(None)
def scalar_traces(p,q):
    # Coefficients of m^2,m,1 on m = r mod p. The central antipodal
    # element contributes (-1)^(m-1)*m^2. Noncentral equal-angle
    # elements contribute m*sin(m*alpha)/sin(alpha).
    out=[]
    for j in range(p):
        rows=[]
        v=(q*j)%p
        alpha=2*mp.pi*j/p; beta=2*mp.pi*v/p
        for r in range(2,p+2):
            if j==0:
                row=(mp.mpf(1),mp.mpf(0),mp.mpf(0))
            elif 2*j==p:
                row=(mp.mpf((-1)**(r-1)),mp.mpf(0),mp.mpf(0))
            elif v==j or v==p-j:
                row=(mp.mpf(0),mp.sin(r*alpha)/mp.sin(alpha),mp.mpf(0))
            else:
                row=(mp.mpf(0),mp.mpf(0),(mp.cos(r*beta)-mp.cos(r*alpha))/(mp.cos(beta)-mp.cos(alpha)))
            rows.append(row)
        out.append(rows)
    return out


@lru_cache(None)
def scalar_logdet(p,q,u):
    traces=scalar_traces(p,q)
    total=mp.mpf(0); lp=mp.log(p); mean_linear=mp.mpf(0)
    for ri,r in enumerate(range(2,p+2)):
        A,B,C=[mp.fsum(mp.cos(2*mp.pi*u*j/p)*traces[j][ri][k] for j in range(p))/p for k in range(3)]
        mean_linear+=B
        for c in (-1,1):
            shifted=(C-c*B+c*c*A,B-2*c*A,A)
            for k in range(3):
                z,dz=hz(k,r+c,p)
                total+=shifted[k]*p**k*(lp*z-dz)
    assert abs(mean_linear)<mp.mpf('1e-165')
    return total


def torsion_log(p,q,u):
    if u==0:
        return mp.log(2*mp.pi**2/p**2)
    v=u*pow(q,-1,p)%p
    return mp.log(4*mp.sin(mp.pi*u/p)*mp.sin(mp.pi*v/p))


def rational_text(text):
    if '/' in text:
        a,b=text.split('/'); return mpq(a,b)
    return mp.mpf(text)


def scalar_multiplicity_direct(p,q,u,m):
    # Harmonic homogeneous polynomials: monomials of degree m-1 in
    # z1,z1bar,z2,z2bar, minus |z|^2 times degree m-3.
    def monomials(n):
        if n<0:return 0
        return sum((a-b+q*(c-(n-a-b-c))-u)%p==0
                   for a in range(n+1) for b in range(n-a+1)
                   for c in range(n-a-b+1))
    return monomials(m-1)-monomials(m-3)


@lru_cache(None)
def hurwitz(j,r,period):
    a = mp.mpf(r)/period
    return mp.zeta(-j,a),mp.zeta(-j,a,derivative=1)


@lru_cache(None)
def coexact_traces(p,q):
    """Trace as C+B*m+A*m^2 on m=r mod 2p, for each group element."""
    period = 2*p
    def su2(h,r,shift):
        # Character of the (m+shift)-dimensional SU(2) representation.
        if h % p == 0:
            sign = (-1)**((r+shift-1)*(h//p))
            return mp.mpf(sign*shift),mp.mpf(sign)
        theta = mp.pi*h/p
        return mp.sin((r+shift)*theta)/mp.sin(theta),mp.mpf(0)
    out = []
    for j in range(p):
        rows = []
        for r in range(2,period+2):
            row = [mp.mpf(0) for _ in range(3)]
            for shift in (-1,1):
                left = su2((1+q)*j,r,shift)
                right = su2((1-q)*j,r,-shift)
                for a in range(2):
                    for b in range(2):
                        row[a+b] += left[a]*right[b]
            rows.append(tuple(row))
        out.append(rows)
    return out


@lru_cache(None)
def coexact_logdet(p,q,u):
    period = 2*p
    traces = coexact_traces(p,q)
    total = mp.mpf(0)
    lp = mp.log(period)
    for ri,r in enumerate(range(2,period+2)):
        coeff = [mp.fsum(mp.cos(2*mp.pi*u*j/p)*traces[j][ri][k]
                         for j in range(p))/p for k in range(3)]
        for k in range(3):
            z,dz = hurwitz(k,r,period)
            total += 2*coeff[k]*period**k*(lp*z-dz)
    return total


def logtorsion(p,q,u):
    if u == 0:
        return mp.log(2*mp.pi**2/p**2)
    v = u*pow(q,-1,p) % p
    return mp.log(4*mp.sin(mp.pi*u/p)*mp.sin(mp.pi*v/p))