back to table · edit · history · where entries came from · files · download
5459 bytes, as of the version from 2026-11-02 15:29 (current). Recorded here, not run.
"""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