"""Taylor and conifold Frobenius continuation for the T680 canonical coordinates.

Pure mpmath arithmetic; heuristic truncation and rounding. No certified bounds.
The Taylor coefficient construction adapts the public checks attached to T679.
The logarithmic solution and conifold matching are implemented here.
"""
import math
from functools import lru_cache
from mpmath import mp


def falling(n,j):
    return math.prod(range(n-j+1,n+1))


def falling_derivative(n,j):
    return sum(math.prod(n-k for k in range(j) if k!=l) for l in range(j))


def ode_polynomials(alpha):
    r=len(alpha); sym=[mp.mpf(1)]
    for a in alpha:
        nxt=[mp.mpf(0)]*(len(sym)+1)
        for j,c in enumerate(sym): nxt[j]+=a*c; nxt[j+1]+=c
        sym=nxt
    S=[[0]*(r+1) for _ in range(r+1)]; S[0][0]=1
    for n in range(1,r+1):
        for j in range(1,n+1): S[n][j]=S[n-1][j-1]+j*S[n-1][j]
    ans=[]
    for j in range(r+1):
        p=[mp.mpf(0)]*(j+2);p[j]=S[r][j]
        p[j+1]=-sum(sym[k]*S[k][j] for k in range(j,r+1))
        ans.append(p)
    return ans


def initial_jets(alpha,x,terms):
    r=len(alpha); F=[mp.mpf(0)]*r;G=F.copy();c=mp.mpf(1);h=mp.mpf(0)
    power=mp.mpf(1)
    for n in range(terms):
        for j in range(min(r,n+1)):
            term=c*math.comb(n,j)*power/x**j
            F[j]+=term;G[j]+=term*h
        h+=sum(1/(a+n) for a in alpha)-mp.mpf(r)/(n+1)
        c*=mp.fprod(a+n for a in alpha)/mp.mpf(n+1)**r
        power*=x
    logjet=[mp.log(abs(x))]+[(-1)**(k-1)/(mp.mpf(k)*x**k) for k in range(1,r)]
    Y=[G[j]+sum(F[j-k]*logjet[k] for k in range(j+1)) for j in range(r)]
    return [F,Y]


def shift_polynomials(p,x):
    return [[sum(v[k]*math.comb(k,l)*x**(k-l) for k in range(l,len(v)))
             for l in range(len(v))] for v in p]


def evaluate_jet(c,x,r=4):
    # Horner on each normalized derivative, avoiding repeated power evaluation.
    jet=[]
    for j in range(r):
        val=mp.mpf(0)
        for n in range(len(c)-1,j-1,-1): val=val*x+math.comb(n,j)*c[n]
        jet.append(val)
    return jet


def step(p,jets,x,h,terms):
    r=len(p)-1;sh=shift_polynomials(p,x);coeff=[list(jet) for jet in jets]
    sparse=[(j,l,v) for j,poly in enumerate(sh) for l,v in enumerate(poly)
            if v and (j,l)!=(r,0)]
    for n in range(terms-r):
        row=[(n-l+j,v*falling(n-l+j,j)) for j,l,v in sparse if n>=l]
        den=sh[r][0]*falling(n+r,r)
        for c in coeff: c.append(-sum(v*c[k] for k,v in row)/den)
    return [evaluate_jet(c,h,r) for c in coeff]


def continue_jets(alpha,z,working):
    terms=4*working+80
    if z<0:
        x=mp.mpf('-0.25'); targets=[mp.mpf(t) for t in ['-0.375','-0.5','-0.75','-1']]
    else:
        x=mp.mpf('0.25')
        if z<=1:
            targets=[mp.mpf(t) for t in ['0.375','0.5']]
            if z==1: targets += [mp.mpf('0.625'),mp.mpf('0.75')]
        else:
            targets=[mp.mpc(a,b) for a,b in [('0.25','0.125'),('0.25','0.25'),('0.25','0.375'),('0.25','0.5'),('0.5','0.5'),('0.75','0.5'),('1','0.5'),('1.25','0.5'),('1.5','0.5'),('1.75','0.5'),('2','0.5'),('2','0.25'),('2','0')]]
    p=ode_polynomials(alpha);jets=initial_jets(alpha,x,terms)
    for y in targets:jets=step(p,jets,x,y-x,terms);x=y
    return jets


def conifold_basis(alpha,working):
    # w=1-z; local exponents are 0,1,1,2. Three analytic solutions and
    # B(w)log(w)+A(w) span the full solution space. B(0)=0, so all limits exist.
    p=ode_polynomials(alpha)
    p=[[(-1)**(j+l)*sum(v[k]*math.comb(k,l) for k in range(l,len(v)))
        for l in range(len(v))] for j,v in enumerate(p)]
    terms=4*working+80
    sparse=[(j,l,v) for j,poly in enumerate(p) for l,v in enumerate(poly) if v]
    def row(k,derivative=False):
        f=falling_derivative if derivative else falling
        return [(k-l+j,v*f(k-l+j,j)) for j,l,v in sparse if k-l+j>=0]
    def analytic(seed,forcing=None):
        c=list(seed)
        for n in range(3,terms):
            rr=row(n-3);den=sum(v for k,v in rr if k==n)
            rhs=sum(v*c[k] for k,v in rr if k<n)
            if forcing is not None:rhs+=sum(v*forcing[k] for k,v in row(n-3,True) if k<len(forcing))
            c.append(-rhs/den)
        return c
    bases=[analytic(seed) for seed in ([1,0,0],[0,1,0],[0,0,1])]
    # Cancel the coefficient of w^-1 in L(B log w).
    rr=row(-1,True)
    b2=-sum(v for k,v in rr if k==1)/sum(v for k,v in rr if k==2)
    B=[bases[1][n]+b2*bases[2][n] for n in range(terms)]
    A=analytic([0,0,0],B)
    w=mp.mpf('0.25');jets=[evaluate_jet(c,w) for c in bases]
    bj=evaluate_jet(B,w);aj=evaluate_jet(A,w)
    lj=[mp.log(w)]+[(-1)**(k-1)/(mp.mpf(k)*w**k) for k in range(1,4)]
    jets.append([aj[j]+sum(bj[j-k]*lj[k] for k in range(j+1)) for j in range(4)])
    # All jets currently in w; target jets are in z.
    return mp.matrix([[(-1)**j*jet[j] for jet in jets] for j in range(4)])


def fg(alpha,z,working):
    with mp.workdps(working):
        jets=continue_jets(alpha,z,working)
        if z==1:
            basis=conifold_basis(alpha,working)
            F=mp.lu_solve(basis,mp.matrix(jets[0]))[0]
            G=mp.lu_solve(basis,mp.matrix(jets[1]))[0]
        else:
            F,Y=jets[0][0],jets[1][0]
            G=Y-F*mp.log(abs(z))
        if z<=1:
            assert mp.im(F)==0 and mp.im(G)==0
        if not mp.isfinite(F) or not mp.isfinite(G) or abs(F)<mp.mpf('0.1'):
            raise ArithmeticError('Nonfinite or numerically small holomorphic denominator')
        return F,G