"""eval_closed_form_V_all.py -- V(X;P) (conformally coupled one-loop triangle, flat-space wavefunction coefficient, loop-integrated) on the WHOLE region {X_v > 0, P a non-degenerate triangle}: both sides of X_v = P_v, all radicand signs, boundaries included. V = (sqrt2 pi/8) [ -1/4 LS_0 F_0 - 1/2 sum_e LS_e F_e - sum_(ij) LS^-_ij F^-_ij + sum_(ij) LS^+_ij F^+_ij ] (compact_representation.json) Per sector s with radicand rho_s and root w = sqrt(|rho_s|): rho_s < 0 : tube F = 2i sum_k c_k D(z_k), D(z) = Im Li2(z) + arg(1-z) log|z| (Bloch-Wigner, single valued), z_k = (p_a - i w)/(p_b - i w) vertex F = 2i sum_k c_k Cl2(theta_k), theta_k = arg x_k, |x_k| = 1 (Clausen) LS F is real; no chamber constants. rho_s > 0 : F = S_real - J. S_real = the same sum with real arguments, Re Li2 and log|.| (continuous except for jumps below). J = chamber constant, an integer multiple of pi^2/2, obtained by ANALYTIC CONTINUATION along the ray P -> tP, 0 < t <= 1: start: t_0 = last zero of rho_s on the ray (there F = S_real = 0, J = 0) or t_0 = 0+ (F -> 0; J_0 = S_real(0+), asserted to be in (pi^2/2) Z); walls: rational t where a linear letter of the sector vanishes; S_real jumps only where an argument passes through infinity (tube: Re Li2(+inf) - Re Li2(-inf) = pi^2/2) or through 0 / infinity (vertex: D(0+) - D(0-) = -pi^2/2, D(+inf) - D(-inf) = +pi^2/2); F itself is real-analytic there. All signs are decided in exact rational arithmetic. If two different letters vanish at the same t (or other coincidences) the constant is computed at a perturbed point of the SAME chamber. Removable loci at the evaluation point (rho_s = 0, a linear letter = 0 incl. X_i = P_i and P_k = |X_i - X_j|, poles of LS: X_k = X_i + X_j, any codimension): V is real-analytic there; the value is obtained by even-order polynomial interpolation of V along a fixed generic line through the point (nodes +-k delta), at raised working precision. No data enter any constant. Input must be exact (integers, 'p/q', decimal strings, Fractions; an mpf is converted to its exact binary value): all sign decisions are rational. Tested (continuation_results.json): 1010 + 372 data points incl. 61 + 328 with some X_v < P_v, exact boundaries (X_i = P_i in 1, 2, 3 sites, P_k = |X_i - X_j|, X_k = X_i + X_j), coincident walls on the ray (constant taken from a perturbed point of the same chamber). NOT tested: exact rational points on a radicand zero (none constructed; points within 1e-30 go through the interpolation branch), degenerate triangles, X_v -> 0. usage: python eval_closed_form_V_all.py X1 X2 X3 P1 P2 P3 [dps]""" import os, sys, json, itertools from fractions import Fraction as Fr import mpmath as mp _here = os.path.dirname(os.path.abspath(__file__)) SPEC = json.load(open(os.path.join(_here, 'compact_representation.json'))) PI2H = lambda: mp.pi**2/2 OPTIONS = {'drop_J': False, 'flip_jump_sign': False, 'naive_bw': False} # controls only class Degenerate(Exception): pass class Coincidence(Exception): pass def toFr(v): if isinstance(v, Fr): return v if isinstance(v, int): return Fr(v) if isinstance(v, str): return Fr(v) if isinstance(v, mp.mpf): s, m, e, _ = v._mpf_; r = Fr(int(m))*(Fr(2)**e); return -r if s else r return Fr(v) def _mp(v): return mp.mpf(v.numerator)/v.denominator def sgn(v): return (v > 0) - (v < 0) def _fr(s): n, d = (str(s).split('/') + ['1'])[:2]; return mp.mpf(int(n))/int(d) def det(M): n = len(M) if n == 1: return M[0][0] if n == 2: return M[0][0]*M[1][1] - M[0][1]*M[1][0] return sum((-1)**j*M[0][j]*det([r[:j] + r[j+1:] for r in M[1:]]) for j in range(n) if M[0][j] != 0) def interp_poly(f, deg): """exact coefficients (low -> high) of a polynomial of degree <= deg from values at t = 0..deg""" xs = [Fr(k) for k in range(deg + 1)]; ys = [f(x) for x in xs]; coef = [Fr(0)]*(deg + 1) for i, xi in enumerate(xs): num = [Fr(1)]; den = Fr(1) for j, xj in enumerate(xs): if j == i: continue num = [Fr(0)] + num; for k in range(len(num) - 1): num[k] -= xj*num[k+1] den *= (xi - xj) for k in range(len(num)): coef[k] += ys[i]*num[k]/den return coef def last_root_below_one(coef, fexact): """largest real root in (0,1) of the polynomial (exact coefs) as an mp number, or None; sign changes confirmed with exact evaluations""" c = list(coef) while c and c[-1] == 0: c.pop() lo = 0 while lo < len(c) and c[lo] == 0: lo += 1 c = c[lo:] if len(c) <= 1: return None with mp.workdps(max(60, mp.mp.dps)): rts = mp.polyroots([_mp(x) for x in reversed(c)], maxsteps=3000, extraprec=2000, error=False) real = sorted(mp.re(r) for r in rts if abs(mp.im(r)) < mp.mpf(10)**(-30) and 0 < mp.re(r) < 1) return real[-1] if real else None # ---------------------------------------------------------------- sector geometry (exact) class Tube: def __init__(s, typ, X, P): s.typ = typ; s.X = X; s.P = P; s.d = SPEC['tube'][typ]; s.sg = 1 if typ == 'p' else -1 s.codes = {a: compile(e, 'p', 'eval') for a, e in s.d['p'].items()}; s.terms = [(str(a), str(b), c) for a, b, c in s.d['terms']] def env(s, t): P = [p*t for p in s.P]; return {'X1': s.X[0], 'X2': s.X[1], 'X3': s.X[2], 'P1': P[0], 'P2': P[1], 'P3': P[2]} def rho(s, t): X, P = s.X, [p*t for p in s.P]; m = P[0]**2 + P[1]**2 - P[2]**2 return m*m - 4*P[0]**2*P[1]**2 + 4*P[0]**2*X[1]**2 + 4*P[1]**2*X[0]**2 + 4*s.sg*X[0]*X[1]*m rho_deg = 4 def ls_den(s): X = s.X; return (X[0] - X[1] - X[2])*(X[0] - X[1] + X[2]) if s.typ == 'm' else (X[0] + X[1] - X[2]) def ls_num(s): X = s.X; return X[2]*(X[0] + X[1] + X[2]) if s.typ == 'm' else X[2] def linear_letters(s): out = {} for rec in list(s.d['N'].values()) + list(s.d['D'].values()): for cv, ex in rec['letters']: al = sum(Fr(c)*x for c, x in zip(cv[:3], s.X)); be = sum(Fr(c)*p for c, p in zip(cv[3:], s.P)); out[tuple(cv)] = (al, be) return out def wall_letters(s): out = {} for a in s.d['N']: for cv, ex in s.d['N'][a]['letters']: al = sum(Fr(c)*x for c, x in zip(cv[:3], s.X)); be = sum(Fr(c)*p for c, p in zip(cv[3:], s.P)); out[tuple(cv)] = (al, be) return out def signs(s, t): """exact signs of the atoms p_a -+ w at ray parameter t (rho > 0)""" e = s.env(t); T = s.rho(t); assert T > 0; out = {} for a, code in s.codes.items(): p = eval(code, {}, e); N = p*p - T out[a] = (-1 if p <= 0 else sgn(N), 1 if p >= 0 else -sgn(N)) return out def args_signs(s, t): S = s.signs(t); return [((S[a][0], S[b][0]), (S[a][1], S[b][1])) for a, b, c in s.terms] def jump(s, A0, A1): RL = lambda sg_: mp.pi**2/3 if sg_ > 0 else -mp.pi**2/6; J = mp.mpf(0) for (a, b, c), (z0, zb0), (z1, zb1) in zip(s.terms, A0, A1): for (n0, d0), (n1, d1), sign in ((z0, z1, 1), (zb0, zb1, -1)): if 0 in (n0, d0, n1, d1): raise Coincidence('atom exactly zero at a sample') if d0 != d1 and n0 == n1: J += sign*c*(RL(n1*d1) - RL(n0*d0)) return J def S_real(s, t): e = s.env(t); T = _mp(s.rho(t)); w = mp.sqrt(T); p = {a: _mp(eval(code, {}, e)) for a, code in s.codes.items()}; tot = mp.mpf(0) for a, b, c in s.terms: pa, pb = p[a], p[b]; z = (pa - w)/(pb - w); zb = (pa + w)/(pb + w) tot += c*(mp.re(mp.polylog(2, z)) - mp.re(mp.polylog(2, zb)) + mp.log(abs(z*zb))*mp.log(abs((pb + w)/(pb - w)))/2) return tot def F_over_w_negative(s): """rho < 0: returns F/(i|w|) * ... i.e. the real number (sum_k c_k 2 D(z_k))/|w|""" e = s.env(Fr(1)); T = _mp(s.rho(Fr(1))); w = mp.sqrt(-T); p = {a: _mp(eval(code, {}, e)) for a, code in s.codes.items()}; tot = mp.mpf(0) for a, b, c in s.terms: z = mp.mpc(p[a], -w)/mp.mpc(p[b], -w) if OPTIONS['naive_bw']: zb = mp.conj(z); v = mp.polylog(2, z) - mp.polylog(2, zb) + mp.log(z*zb)*mp.log((1 - z)/(1 - zb))/2; tot += c*mp.im(v) else: tot += c*2*(mp.im(mp.polylog(2, z)) + mp.arg(1 - z)*mp.log(abs(z))) return tot/w class Vertex: rho_deg = 6 def __init__(s, which, X, P): s.which = which; s.X = X; s.P = P; X1, X2, X3 = X; Q = X1 + X2 + X3; h = Fr(1, 2); s.Q = Q s.y = [(X3 - X1 - X2)*h, (X1 - X2 - X3)*h, (X2 - X1 - X3)*h] if which == 'F0' else [-Q*h, Q*h - X2, Q*h - X1] s.terms = [(_fr(c), [tuple(e) for e in cyc], list(eps), (list(cyc[0]) if (len(cyc) == 2 and cyc[0] == cyc[1]) else sorted({i for e in cyc for i in e}))) for c, cyc, eps in SPEC['vertex'][which]] s.edges = list(itertools.combinations(range(1, 5), 2)) def dist(s, t): P1, P2, P3 = [p*t for p in s.P]; return {(1, 2): P2, (1, 3): P1, (2, 3): P3, (1, 4): s.y[0], (2, 4): s.y[1], (3, 4): s.y[2]} def M(s, t): d = s.dist(t); return [[0, 1, 1, 1, 1]] + [[1] + [0 if a == b else d[tuple(sorted((a, b)))]**2 for b in range(1, 5)] for a in range(1, 5)] def rho(s, t): return -2*det(s.M(t)) def cof(s, M, i, j): return det([[M[r][c] for c in range(5) if c != j] for r in range(5) if r != i])*(-1)**(i + j) def ls_den(s): return Fr(1) def ls_num(s): return s.Q def face_forms(s, t=None): """linear forms (alpha + beta t) whose product is -C_vv""" out = {} dX = {(1, 4): s.y[0], (2, 4): s.y[1], (3, 4): s.y[2]}; dP = {(1, 2): s.P[1], (1, 3): s.P[0], (2, 3): s.P[2]} for v in range(1, 5): ks = [k for k in sorted(list(dX) + list(dP)) if v not in k] for sg_ in [(1, 1, 1), (-1, 1, 1), (1, -1, 1), (1, 1, -1)]: al = sum(g*dX[k] for g, k in zip(sg_, ks) if k in dX); be = sum(g*dP[k] for g, k in zip(sg_, ks) if k in dP); out[(v, sg_)] = (Fr(al), Fr(be)) return out wall_letters = face_forms linear_letters = face_forms def atoms(s, t): M = s.M(t); R = -2*det(M); assert R > 0; d = s.dist(t); C = {}; U = {} for v in range(1, 5): C[v] = s.cof(M, v, v) for e in s.edges: rem = tuple(k for k in range(1, 5) if k not in e); c = s.cof(M, *e); l = d[rem]; nrm = C[e[0]]*C[e[1]] sp_ = sgn(c) if sgn(c)*sgn(l) >= 0 and (c != 0 or l != 0) else sgn(c)*sgn(nrm) # sign(c + l w), w > 0 if c == 0: sp_ = sgn(l) sm_ = sgn(c) if sgn(c)*sgn(-l) >= 0 else sgn(c)*sgn(nrm) # sign(c - l w) if c == 0: sm_ = -sgn(l) U[e] = {1: sp_, -1: sm_} return C, U def args_signs(s, t): C, U = s.atoms(t); return [([U[e][ep] for e, ep in zip(cyc, eps)], [sgn(C[v]) for v in vs]) for c, cyc, eps, vs in s.terms] def jump(s, A0, A1): D0 = lambda sg_: -mp.pi**2/3 if sg_ > 0 else mp.pi**2/6; Dinf = lambda sg_: mp.pi**2/3 if sg_ > 0 else -mp.pi**2/6; J = mp.mpf(0) for (c, cyc, eps, vs), (n0, d0), (n1, d1) in zip(s.terms, A0, A1): if 0 in n0 + d0 + n1 + d1: raise Coincidence('atom exactly zero at a sample') order = sum(a != b for a, b in zip(n0, n1)) - sum(a != b for a, b in zip(d0, d1)) s0 = 1; s1 = 1 for a in n0 + d0: s0 *= a for a in n1 + d1: s1 *= a if order < 0: J += c*(Dinf(s1) - Dinf(s0)) elif order > 0: J += c*(D0(s1) - D0(s0)) elif s0 != s1: raise Coincidence('sign change without net order') return J def _xs(s, t, w, negative=False): """arguments x_k; cofactors are computed EXACTLY (rational) and converted: floating determinants cancel completely near t = 0 (C_44 ~ t^4)""" Me = s.M(t); d = {k: _mp(Fr(v)) for k, v in s.dist(t).items()} C = {v: _mp(s.cof(Me, v, v)) for v in range(1, 5)}; U = {} for e in s.edges: rem = tuple(k for k in range(1, 5) if k not in e); c = _mp(s.cof(Me, *e)); lw = d[rem]*w if negative: U[e] = {1: mp.mpc(c, lw), -1: mp.mpc(c, -lw)} else: nrm = C[e[0]]*C[e[1]] if c*lw >= 0: big = c + lw; U[e] = {1: big, -1: nrm/big} else: big = c - lw; U[e] = {-1: big, 1: nrm/big} xs = [] for c, cyc, eps, vs in s.terms: x = mp.mpf(1) for e, ep in zip(cyc, eps): x = x*U[e][ep] for v in vs: x = x/C[v] xs.append(x) return xs def S_real(s, t): w = mp.sqrt(_mp(s.rho(t))); tot = mp.mpf(0) for (c, cyc, eps, vs), x in zip(s.terms, s._xs(t, w)): if abs(x) > 1: x = 1/x; sg_ = -1 else: sg_ = 1 L = mp.log(abs(x)); tot += sg_*c*(2*mp.re(mp.polylog(2, x)) + L*L/2 + (-mp.pi**2/3 if x > 0 else mp.pi**2/6)) return tot def F_over_w_negative(s): w = mp.sqrt(-_mp(s.rho(Fr(1)))); tot = mp.mpf(0) for (c, cyc, eps, vs), x in zip(s.terms, s._xs(Fr(1), w, negative=True)): if OPTIONS['naive_bw']: tot += c*mp.im(mp.polylog(2, x) - mp.polylog(2, 1/x)) else: tot += c*2*mp.clsin(2, mp.arg(x)) return tot/w # ---------------------------------------------------------------- chamber constant by continuation along the ray def chamber_constant(sec, info=None): """J such that F = S_real - J at t = 1 (rho(1) > 0); exact sign bookkeeping""" coef = interp_poly(sec.rho, sec.rho_deg); t0 = last_root_below_one(coef, sec.rho) walls = {} for key, (al, be) in sec.wall_letters().items(): if be == 0: continue tw = -al/be if tw == 1: raise Degenerate('letter vanishes at the evaluation point') if 0 < tw < 1 and (t0 is None or _mp(tw) > t0): walls.setdefault(tw, []).append(key) ts = sorted(walls) if any(len(walls[t]) > 1 for t in ts): raise Coincidence('two different letters vanish at the same t on the ray') lo = Fr(0) if t0 is None else toFr(mp.mpf(t0)) if t0 is not None: if ts and _mp(ts[0]) - t0 < mp.mpf(10)**-12: raise Coincidence('wall too close to radicand zero') if 1 - t0 < mp.mpf(10)**-25: raise Degenerate('radicand zero at the evaluation point (numerically)') samples = []; prev = lo for tw in ts + [Fr(1)]: samples.append((prev + tw)/2 if tw != 1 or True else Fr(1)); prev = tw samples[-1:] = [samples[-1], Fr(1)] if ts or True else samples if t0 is not None: k = 0 while sec.rho(samples[0]) <= 0 and k < 60: samples[0] = (samples[0] + (ts[0] if ts else Fr(1)))/2; k += 1 for sm in samples: if sec.rho(sm) <= 0: raise Coincidence('radicand not positive at a sample (unresolved radicand zeros)') A = [sec.args_signs(sm) for sm in samples]; J = mp.mpf(0) for A0, A1 in zip(A[:-1], A[1:]): J += sec.jump(A0, A1) if OPTIONS['flip_jump_sign']: J = -J J0 = mp.mpf(0) if t0 is None: eta = Fr(1, 10**30)*min([t for t in ts] + [Fr(1)]) if sec.rho(eta) <= 0: raise Coincidence('radicand not positive at 0+') A00 = sec.args_signs(eta) if A00 != A[0]: raise Coincidence('sign pattern changes between 0+ and the first sample') with mp.workdps(mp.mp.dps + 40): S0 = sec.S_real(eta); n0 = S0/PI2H(); n0r = mp.nint(n0) if abs(n0 - n0r) > mp.mpf(10)**-12: raise AssertionError('S_real(0+) is not a multiple of pi^2/2: %s' % mp.nstr(n0, 20)) J0 = n0r*PI2H() if info is not None: info.update({'t0': None if t0 is None else float(t0), 'n_walls': len(ts), 'n_jump': int(mp.nint(J/PI2H())), 'n_start': int(mp.nint(J0/PI2H()))}) return J + J0 SECTORS = [('S0', 'F0', (0, 1, 2), '-1/4')] + [(nm, 'Fe', tuple(pm), '-1/2') for nm, pm in SPEC['orbits']['vertex_e'].items()] + \ [('t%s%s' % (ij, typ), typ, tuple(pm), '-1' if typ == 'm' else '1') for ij, pm in SPEC['orbits']['tube'].items() for typ in 'mp'] def make_sector(kind, pm, X, P): Xr = [X[i] for i in pm]; Pr = [P[i] for i in pm] return Vertex(kind, Xr, Pr) if kind in ('F0', 'Fe') else Tube(kind, Xr, Pr) def _perturb(X, P, k): import random rnd = random.Random(1000 + k); e = Fr(1, 10**(14 - 2*min(k, 4))) return [x*(1 + e*Fr(rnd.randint(-1000, 1000), 1000)) for x in X], [p*(1 + e*Fr(rnd.randint(-1000, 1000), 1000)) for p in P] def _endpoint_pattern(sec): return (sgn(sec.rho(Fr(1))),) + tuple(sgn(al + be) for al, be in sec.linear_letters().values()) def sector_value(name, kind, pm, X, P, info=None): """LS_s F_s (real) at a point where nothing degenerates""" sec = make_sector(kind, pm, X, P); r1 = sec.rho(Fr(1)) if r1 == 0 or sec.ls_den() == 0: raise Degenerate(name) if any(al + be == 0 for al, be in sec.linear_letters().values()): raise Degenerate(name + ': linear letter zero') LSn = _mp(Fr(sec.ls_num())/Fr(sec.ls_den())); inf = {} if r1 < 0: val = LSn*sec.F_over_w_negative(); inf = {'rho': -1} else: J = None if OPTIONS['drop_J']: J = mp.mpf(0) else: try: J = chamber_constant(sec, inf) except Coincidence as ex: pat = _endpoint_pattern(sec) for k in range(8): Xp, Pp = _perturb(X, P, k); sp = make_sector(kind, pm, Xp, Pp) if _endpoint_pattern(sp) != pat: continue try: J = chamber_constant(sp, inf); inf['perturbed_for_constant'] = k + 1; break except Coincidence: continue if J is None: raise val = LSn*(sec.S_real(Fr(1)) - J)/mp.sqrt(_mp(r1)); inf['rho'] = 1 if info is not None: info[name] = inf return val def V_generic(X, P, info=None): tot = mp.mpf(0); terms = {} for name, kind, pm, cf in SECTORS: v = _fr(cf)*sector_value(name, kind, pm, X, P, info); terms[name] = v; tot += v if info is not None: info['terms'] = {k: mp.nstr(v, 30) for k, v in terms.items()}; info['cond'] = float(max(abs(v) for v in terms.values())/abs(tot)) if tot != 0 else None return mp.sqrt(2)*mp.pi/8*tot def smallness(X, P): """min over all removable loci of |quantity|/scale, exact; 0 if the point lies exactly on one""" sc = sum(X) + sum(P); m = Fr(1) for name, kind, pm, cf in SECTORS: sec = make_sector(kind, pm, X, P); deg = sec.rho_deg m = min(m, abs(sec.rho(Fr(1)))/sc**deg, abs(Fr(sec.ls_den()))/sc**(2 if kind == 'm' else 1)) for al, be in sec.linear_letters().values(): m = min(m, abs(al + be)/sc) return m _DIR = ([Fr(3, 7), Fr(-5, 11), Fr(2, 13)], [Fr(-4, 9), Fr(1, 3), Fr(6, 17)]) def V_all(X, P, info=None, K=5): X = [toFr(x) for x in X]; P = [toFr(p) for p in P] if min(X) <= 0 or min(P) <= 0 or min(sum(P) - 2*p for p in P) <= 0: raise ValueError('outside the region X_v > 0, P a non-degenerate triangle') dps = mp.mp.dps; sm = smallness(X, P); lost = 0 if sm == 0 else max(0, int(-mp.log10(_mp(sm)))) if sm != 0 and lost <= dps//3: with mp.workdps(dps + 25 + 3*lost): v = V_generic(X, P, info) if info is not None: info['mode'] = 'direct'; info['digits_margin'] = 25 + 3*lost return +v # removable locus (exact or numerically indistinguishable): interpolate V along a generic line through the point dlog = max(6, dps//(2*K) + 2); delta = Fr(1, 10**dlog) with mp.workdps(dps + 40 + 4*dlog): nodes = [k*sg_ for k in range(1, K + 1) for sg_ in (1, -1)]; vals = {} for k in nodes: Xk = [x*(1 + k*delta*a) for x, a in zip(X, _DIR[0])]; Pk = [p*(1 + k*delta*b) for p, b in zip(P, _DIR[1])] if smallness(Xk, Pk) == 0: raise Degenerate('interpolation node on a removable locus') vals[k] = V_generic(Xk, Pk) tot = mp.mpf(0) for k in nodes: wgt = mp.mpf(1) for j in nodes: if j != k: wgt *= mp.mpf(-j)/(k - j) tot += wgt*vals[k] if info is not None: info['mode'] = 'interpolated'; info['delta'] = float(delta); info['K'] = K return +tot if __name__ == '__main__': mp.mp.dps = int(sys.argv[7]) if len(sys.argv) > 7 else 40 print(mp.nstr(V_all(sys.argv[1:4], sys.argv[4:7]), mp.mp.dps - 5))