"""Evaluator for the anisotropic simple-cubic Watson integral. The Watson integral of the paper "Exact results for the anisotropic Watson integral" is W_S(alpha, beta, gamma; w) = (1/pi^3) * int_0^pi int_0^pi int_0^pi dk1 dk2 dk3 / (w - alpha*cos k1 - beta*cos k2 - gamma*cos k3), the screened lattice Green function at the origin of the simple-cubic lattice with hopping rates (alpha, beta, gamma), defined for positive rates and w >= alpha + beta + gamma (the band edge, or threshold, is w* = alpha + beta + gamma). watson(w, a, b, c, dps) evaluates W_S by the paper's closed-form procedure ("The explicit evaluation", seven steps): (i) elementary symmetric functions e1, e2, e3 of the squared rates (A, B, C) = (a^2, b^2, c^2) and the spectral point x = 1/w^2; (ii) the threshold quartic Q4(x) = prod_{+/-}(w +/- a +/- b +/- c)/w^8 and the marked branch R = -sqrt(Q4); (iii) the four Igusa-Clebsch invariants I2, I4, I6, I10 of the marked branch, polynomials in e1, e2, e3, x, R; (iv) a real genus-two sextic model y^2 = f(t) with these invariants from Mestre's conic-and-cubic construction (a real point on Mestre's conic substituted into his cubic covariant, the point chosen so that all six roots of f are real); (v) the marked partition f = G1*G2*G3: among the fifteen pairings of the six branch points into three quadratics, the one whose multiset {disc(Gi)*Res(Gj,Gk)} is proportional to x*{A,B,C} (the disc-resultant normalization test); Delta is the determinant of the 3x3 coefficient matrix of (G1, G2, G3); (vi) the two integration segments fixed by cyclic adjacency (the second and fourth of the five root gaps for the cut pairing, the first and third for the gap pairing), the period integrals p_g(J) = 2*int_J t^g dt/sqrt(|f(t)|) for g = 0, 1, and the minor M = p_0(J1)*p_1(J2) - p_1(J1)*p_0(J2); (vii) W_S = |Delta * M| / (4 pi^2 w). The evaluation is performed at increasing working precision until two consecutive precision levels agree to the requested number of digits, so the returned value is verified rather than assumed. The period quadratures substitute t = endpoint +/- s^2, which removes the inverse-square-root endpoint singularities exactly (|f| is evaluated in product form over the computed roots), leaving analytic integrands for tanh-sinh quadrature. The closed form lives off the rate-coincidence loci: when two or three rates are exactly equal the genus-two model degenerates (Mestre's conic drops rank). watson() then evaluates the limit instead: the coincident rates are split symmetrically by exact binary epsilons, W_S is evaluated at a ladder of split tuples, and Richardson extrapolation in epsilon^2 recovers the coincident value, again with an internal agreement check. watson_bessel(w, a, b, c, dps) is an independent check on all of the above: the Bessel-Laplace form of the defining integral, W_S = int_0^inf exp(-w t) I0(a t) I0(b t) I0(c t) dt (w > a+b+c), evaluated by direct quadrature. It shares no code with the closed form. Running the file as a script (python3 watson.py) computes a validation table: closed form against the Bessel-Laplace oracle at tuples spanning the chamber, and against two classical gamma-function values computed in place -- Watson's 1939 constant W_S(1,1,1;3) = sqrt(6)/(96 pi^3) * Gamma(1/24) Gamma(5/24) Gamma(7/24) Gamma(11/24) and the elliptic-CM value W_S(1,1,1; 3 sqrt(6)/2) = 3/(2^(31/6) pi^4) * Gamma(1/3)^6. Every number is computed at run time; nothing is compared against stored digits. Requires python3 and mpmath. No other dependencies. """ import itertools import sys from mpmath import mp, mpf, mpmathify # ---------------------------------------------------------------------- # polynomial helpers: coefficient lists, ascending degree # ---------------------------------------------------------------------- def _padd(p, q): n = max(len(p), len(q)) return [(p[i] if i < len(p) else mp.zero) + (q[i] if i < len(q) else mp.zero) for i in range(n)] def _pmul(p, q): r = [mp.zero] * (len(p) + len(q) - 1) for i, pi in enumerate(p): if pi: for j, qj in enumerate(q): r[i + j] += pi * qj return r def _pscale(p, s): return [s * pi for pi in p] def _ptrim(p, tol): top = max(abs(pi) for pi in p) n = len(p) while n > 1 and abs(p[n - 1]) <= tol * top: n -= 1 return p[:n] # ---------------------------------------------------------------------- # step (iv): Mestre's construction # ---------------------------------------------------------------------- def _mestre_sextic(i2, i4, i6, i10, sgn): """Sextic f(t) (ascending coefficient list) whose genus-two curve y^2 = f(t) has Igusa-Clebsch invariants (i2 : i4 : i6 : i10). Mestre's conic-and-cubic construction: L is the matrix of the conic, P a real point on it built from the eigenvector split (sgn selects between the two real points), and f is the cubic covariant evaluated along the line V = (0, 1, t) through P.""" x = 8 * (1 + 20 * i4 / i2**2) / 225 y = 16 * (1 + 80 * i4 / i2**2 - 600 * i6 / i2**3) / 3375 z = -64 * (-10800000 * i10 / i2**5 - 9 - 700 * i4 / i2**2 + 3600 * i6 / i2**3 + 12400 * i4**2 / i2**4 - 48000 * i4 * i6 / i2**5) / 253125 L = mp.matrix([ [x + 6 * y, 6 * x**2 + 2 * y, 2 * z], [6 * x**2 + 2 * y, 2 * z, 9 * x**3 + 4 * x * y + 6 * y**2], [2 * z, 9 * x**3 + 4 * x * y + 6 * y**2, 6 * x**2 * y + 2 * y**2 + 3 * x * z]]) E, Q = mp.eigsy(L) # eigenvalues ascending lam_min, lam_max = E[0], E[2] if not (lam_max > 0 and lam_min < 0): raise _ModelNotFound("Mestre conic matrix is not indefinite") P = [Q[k, 2] / mp.sqrt(lam_max) + sgn * Q[k, 0] / mp.sqrt(-lam_min) for k in range(3)] # V = (0, 1, s): each component of U = (V.L.V) P - 2 (P.L.V) V # is a quadratic in s (coefficient lists, ascending) VLV = [L[1, 1], 2 * L[1, 2], L[2, 2]] PL = [sum(P[k] * L[k, j] for k in range(3)) for j in range(3)] PLV = [PL[1], PL[2]] U = [] for i in range(3): Ui = _pscale(VLV, P[i]) if i == 1: Ui = _padd(Ui, _pscale(PLV + [mp.zero], -2)) if i == 2: Ui = _padd(Ui, _pmul([mp.zero, mp.one], _pscale(PLV, -2))) U.append(Ui) coeffs = { (0, 0, 0): 12*x*y - 2*y/3 - 4*z, (0, 0, 1): -18*x**3 - 12*x*y - 36*y**2 - 2*z, (0, 0, 2): -9*x**3 - 36*x**2*y - 4*x*y - 6*x*z - 18*y**2, (0, 1, 1): -9*x**3 - 36*x**2*y - 4*x*y - 6*x*z - 18*y**2, (0, 1, 2): -54*x**4 - 36*x**2*y - 36*x*y**2 - 6*x*z - 4*y**2 - 24*y*z, (0, 2, 2): (-27*x**4/2 - 72*x**3*y - 6*x**2*y - 9*x**2*z - 39*x*y**2 - 36*y**3 - 2*y*z), (1, 1, 1): -27*x**4 - 18*x**2*y - 6*x*y**2 - 8*y**2/3 + 2*y*z, (1, 1, 2): 9*x**3*y - 27*x**2*z + 6*x*y**2 + 18*y**3 - 8*y*z, (1, 2, 2): (-81*x**5/2 - 27*x**3*y - 9*x**2*y**2 - 4*x*y**2 + 3*x*y*z - 6*z**2), (2, 2, 2): (27*x**4*y/2 - 27*x**3*z/2 + 9*x**2*y**2 + 3*x*y**3 - 6*x*y*z + 4*y**3/3 - 10*y**2*z), } f = [mp.zero] * 7 for (i, j, k), cval in coeffs.items(): f = _padd(f, _pscale(_pmul(_pmul(U[i], U[j]), U[k]), cval)) return f # ---------------------------------------------------------------------- # steps (i)-(vii): one evaluation at fixed working precision # ---------------------------------------------------------------------- class _ModelNotFound(ValueError): """No real sextic model resolved at the current working precision.""" # the fifteen pairings of six ordered branch points into three pairs _PAIRINGS = [(1, 2, 3, 4, 5, 6), (1, 3, 2, 4, 5, 6), (1, 4, 2, 3, 5, 6), (1, 2, 3, 5, 4, 6), (1, 2, 3, 6, 4, 5), (1, 3, 2, 5, 4, 6), (1, 3, 2, 6, 4, 5), (1, 4, 2, 5, 3, 6), (1, 4, 2, 6, 3, 5), (1, 5, 2, 3, 4, 6), (1, 5, 2, 4, 3, 6), (1, 5, 2, 6, 3, 4), (1, 6, 2, 3, 4, 5), (1, 6, 2, 4, 3, 5), (1, 6, 2, 5, 3, 4)] def _core(w, a, b, c, wp): """One pass of the seven-step evaluation at working precision wp. Raises _ModelNotFound when the sextic model does not resolve.""" with mp.workdps(wp): w_, a_, b_, c_ = mpf(w), mpf(a), mpf(b), mpf(c) A, B, C = a_**2, b_**2, c_**2 e1, e2, e3 = A + B + C, A*B + A*C + B*C, A*B*C x = 1 / w_**2 # (ii) threshold quartic in eight-plane product form (exact at # the band edge, no cancellation) and the marked branch q4 = mp.one for s1, s2, s3 in itertools.product((1, -1), repeat=3): q4 *= (w_ + s1*a_ + s2*b_ + s3*c_) q4 /= w_**8 if q4 < 0: raise ValueError("threshold quartic negative: " "the tuple lies below the band edge") rr = -mp.sqrt(q4) # (iii) Igusa-Clebsch invariants of the marked branch i2 = 48*e1**2*x**2 - 32*e1*x - 192*e2*x**2 + 80 - 48*rr i4 = 544*e1**2*x**2 - 1088*e1*x - 1152*e2*x**2 + 544 - 480*rr i6 = (16384*e1**4*x**4 - 48384*e1**3*x**3 - 114688*e1**2*e2*x**4 + 64256*e1**2*x**2 + 177152*e1*e2*x**3 - 48896*e1*x + 196608*e2**2*x**4 - 99328*e2*x**2 - 417792*e3*x**3 + 16640 + (-16384*e1**2*x**2 + 17152*e1*x + 49152*e2*x**2 - 16128)*rr) i10 = 32768*e3*x**3*(e1**2*x**2 - 2*e1*x - 4*e2*x**2 + 1 + rr) # (iv) real sextic model; try both real points on Mestre's conic f = roots = None for sgn in (1, -1): try: fs = _mestre_sextic(i2, i4, i6, i10, sgn) except _ModelNotFound: continue fs = _ptrim(fs, mpf(10)**(-(wp - 5))) if len(fs) != 7: continue try: rts = mp.polyroots(list(reversed(fs)), maxsteps=200 + 5*wp, extraprec=2*wp) except mp.NoConvergence: continue if max(abs(mp.im(rt)) for rt in rts) <= mpf(10)**-10: f = fs roots = sorted(mp.re(rt) for rt in rts) break if f is None: raise _ModelNotFound("no real sextic model at this precision") lc = f[-1] r = roots # (v) marked partition by the disc-resultant normalization test; # for monic quadratics with known roots, disc and Res are the # closed-form root-difference products def _disc(i, j): return (r[i-1] - r[j-1])**2 def _res(p, q): (i, j), (k, l) = p, q return ((r[i-1] - r[k-1]) * (r[i-1] - r[l-1]) * (r[j-1] - r[k-1]) * (r[j-1] - r[l-1])) tgt = sorted([x*A, x*B, x*C]) best = best_score = None for pr in _PAIRINGS: p1, p2, p3 = (pr[0], pr[1]), (pr[2], pr[3]), (pr[4], pr[5]) v = sorted([abs(_disc(*p1) * _res(p2, p3)), abs(_disc(*p2) * _res(p3, p1)), abs(_disc(*p3) * _res(p1, p2))]) ratios = [v[k] / tgt[k] for k in range(3)] score = max(ratios) / min(ratios) if best_score is None or score < best_score: best, best_score = pr, score # Delta: leading coefficient times the determinant of the # coefficient matrix of the three monic quadratic factors rows = [[mp.one, -(r[i-1] + r[j-1]), r[i-1] * r[j-1]] for (i, j) in ((best[0], best[1]), (best[2], best[3]), (best[4], best[5]))] delta = lc * mp.det(mp.matrix(rows)) # (vi) segments by cyclic adjacency (cut vs gap pairing) cut = (best[1] == best[0] + 1) and (best[3] == best[2] + 1) seg = ((1, 2), (3, 4)) if cut else ((0, 1), (2, 3)) # period integrals: t = endpoint +/- s^2 removes the endpoint # singularities; |f| in product form over the roots alc = abs(lc) maxdeg = 6 + wp // 100 def _period(g, iu, iv): u, v = r[iu], r[iv] m = (u + v) / 2 other_u = [r[j] for j in range(6) if j != iu] other_v = [r[j] for j in range(6) if j != iv] def _half(endpoint, sign, others): def h(s): t = endpoint + sign * s * s prod = alc for rj in others: prod *= abs(t - rj) return 2 * t**g / mp.sqrt(prod) return h h1 = mp.quad(_half(u, 1, other_u), [0, mp.sqrt(m - u)], method='tanh-sinh', maxdegree=maxdeg) h2 = mp.quad(_half(v, -1, other_v), [0, mp.sqrt(v - m)], method='tanh-sinh', maxdegree=maxdeg) return 2 * (h1 + h2) p0a, p1a = _period(0, *seg[0]), _period(1, *seg[0]) p0b, p1b = _period(0, *seg[1]), _period(1, *seg[1]) minor = p0a * p1b - p1a * p0b # (vii) return +(abs(delta * minor) / (4 * mp.pi**2 * w_)) def _evaluate_distinct(w, a, b, c, dps): """Closed-form value at a distinct-rate tuple, verified: the working precision is raised until two consecutive levels agree to dps.""" wp = max(40, dps + 15) prev = None tol = mpf(10)**(-(dps + 2)) for _ in range(8): try: val = _core(w, a, b, c, wp) except _ModelNotFound: val = None if val is not None and prev is not None: if abs(val - prev) <= tol * abs(val): return val prev = val wp = int(wp * 1.6) + 20 raise ValueError( "closed-form evaluation did not stabilize by working precision " "%d; the genus-two model is too degenerate at this tuple " "(rates nearly coincident, or w very far above the band edge). " "watson_bessel() evaluates the defining integral directly." % wp) def _evaluate_coincident(w, a, b, c, dps): """Limit of the closed form onto a rate-coincidence locus, where the genus-two model degenerates: split the coincident rates by exact binary epsilons e_j = 2^-(7+j) (W_S is even in the split), evaluate at each split tuple, and Richardson-extrapolate in e^2 to zero. Two extrapolation windows must agree to the requested digits.""" node_dps = dps + 12 max_nodes = 10 win = 6 nodes = [] def _node(j): e = mpf(2)**(-(7 + j)) if a == b == c: aa, bb, cc = a*(1 + e), a, a*(1 - e) elif a == b: aa, bb, cc = a*(1 + e), a*(1 - e), c elif b == c: aa, bb, cc = a, b*(1 + e), b*(1 - e) else: # a == c aa, bb, cc = a*(1 + e), b, a*(1 - e) return (e*e, _evaluate_distinct(w, aa, bb, cc, node_dps)) def _extrapolate(nds): w0 = mp.zero for i, (ui, wi) in enumerate(nds): li = mp.one for j, (uj, _) in enumerate(nds): if j != i: li *= uj / (uj - ui) w0 += li * wi return w0 tol = mpf(10)**(-(dps + 2)) with mp.workdps(node_dps + 10): for j in range(win): nodes.append(_node(j)) while True: full = _extrapolate(nodes[-win:]) part = _extrapolate(nodes[-(win - 1):]) if abs(full - part) <= tol * abs(full): return full if len(nodes) >= max_nodes: raise ValueError( "rate-coincidence extrapolation did not stabilize " "to %d digits; try a lower dps" % dps) nodes.append(_node(len(nodes))) def watson(w, a=1, b=1, c=1, dps=30): """The anisotropic simple-cubic Watson integral W_S(a, b, c; w) by the closed-form evaluation (threshold quartic, Igusa-Clebsch invariants, Mestre sextic model, marked partition, period minor). w, a, b, c : the spectral parameter and the three hopping rates, with a, b, c > 0 and w >= a + b + c. Any mpmath-readable number; pass exact values (ints, fractions of ints, strings) when working at the band edge w = a + b + c. dps : decimal digits of the result (default 30). Returns an mpf correct to about dps digits (internally verified by agreement across working precisions). Raises ValueError below the band edge or when the evaluation cannot be stabilized.""" if dps < 5: raise ValueError("dps must be at least 5") with mp.workdps(dps + 25): w_, a_, b_, c_ = (mpmathify(v) for v in (w, a, b, c)) if not (a_ > 0 and b_ > 0 and c_ > 0): raise ValueError("the rates a, b, c must be positive") if w_ < a_ + b_ + c_: raise ValueError("w must satisfy w >= a + b + c " "(at or above the band edge)") if a_ == b_ or b_ == c_ or a_ == c_: val = _evaluate_coincident(w_, a_, b_, c_, dps) else: val = _evaluate_distinct(w_, a_, b_, c_, dps) with mp.workdps(dps): return +val # ---------------------------------------------------------------------- # independent check: the Bessel-Laplace form of the defining integral # ---------------------------------------------------------------------- def watson_bessel(w, a=1, b=1, c=1, dps=30): """W_S(a, b, c; w) from the Bessel-Laplace form of the defining integral, W_S = int_0^inf exp(-w t) I0(a t) I0(b t) I0(c t) dt, by direct quadrature (valid for w > a + b + c; the integrand decays like exp(-(w-a-b-c) t) t^(-3/2)). Independent of the closed form: no code is shared with watson(). Returns an mpf correct to about dps digits.""" wp = dps + 10 with mp.workdps(wp): w_, a_, b_, c_ = (mpmathify(v) for v in (w, a, b, c)) if not (a_ > 0 and b_ > 0 and c_ > 0): raise ValueError("the rates a, b, c must be positive") decay = w_ - (a_ + b_ + c_) if decay <= 0: raise ValueError("the Bessel-Laplace representation needs " "w > a + b + c") def h(t): return (mp.exp(-w_*t) * mp.besseli(0, a_*t) * mp.besseli(0, b_*t) * mp.besseli(0, c_*t)) # beyond T the integrand is below the target precision T = (wp + 5) * mp.log(10) / decay val = mp.quad(h, [0, 1/w_, T, mp.inf], method='tanh-sinh') with mp.workdps(dps): return +val # ---------------------------------------------------------------------- # self-test: closed form vs the oracle and vs classical gamma values # ---------------------------------------------------------------------- def _digits(u, v): """Number of decimal digits to which u and v agree.""" d = abs(u - v) if d == 0: return mp.inf return -mp.log10(d / abs(v)) def _selftest(): import time t_start = time.time() mp.dps = 40 failures = [] print("Watson integral evaluator self-test", flush=True) print("closed form (watson) vs Bessel-Laplace oracle (watson_bessel)," " dps=30, gate >= 20 digits", flush=True) print("-" * 78, flush=True) print("%-28s %-26s %9s %8s" % ("(w, a, b, c)", "W_S", "digits", "time"), flush=True) points = [ ((17, 1, 5, 7), "paper tuple"), ((22, 5, 6, 9), "paper tuple"), ((16, 1, 5, 8), "paper tuple"), ((25, 4, 7, 8), "paper tuple"), ((14, 1, 5, 7), "near band edge"), (('27/2', 1, 5, 7), "near band edge, non-integer"), ((100, 1, 2, 3), "large w"), ((5, 1, 1, 1), "isotropic (coincident limit)"), ((10, 2, 3, 3), "tetragonal (coincident limit)"), ] for args, label in points: t0 = time.time() cf = watson(*args, dps=30) orc = watson_bessel(*args, dps=32) dg = _digits(cf, orc) ok = dg >= 20 if not ok: failures.append("%s: only %.1f digits vs oracle" % (args, dg)) print("%-28s %-26s %9.1f %7.1fs %s [%s]" % (str(args), mp.nstr(cf, 20), float(dg), time.time() - t0, "ok" if ok else "FAIL", label), flush=True) print("-" * 78, flush=True) # Watson's 1939 constant at the isotropic band edge, computed from # its gamma-function expression t0 = time.time() cf = watson(3, 1, 1, 1, dps=30) ref = (mp.sqrt(6) / (96 * mp.pi**3) * mp.gamma(mpf(1)/24) * mp.gamma(mpf(5)/24) * mp.gamma(mpf(7)/24) * mp.gamma(mpf(11)/24)) dg = _digits(cf, ref) ok = dg >= 20 if not ok: failures.append("Watson constant: only %.1f digits" % dg) print("W_S(1,1,1;3) closed form %s" % mp.nstr(cf, 25), flush=True) print(" sqrt(6)/(96 pi^3) Gamma(1/24) Gamma(5/24) Gamma(7/24)" " Gamma(11/24)\n" " gamma value %s\n" " agreement %.1f digits, %.1fs %s" % (mp.nstr(ref, 25), float(dg), time.time() - t0, "ok" if ok else "FAIL"), flush=True) # the elliptic-CM point on the isotropic line, non-integer w t0 = time.time() with mp.workdps(60): w_cm = 3 * mp.sqrt(6) / 2 cf = watson(w_cm, 1, 1, 1, dps=30) ref = 3 / (mpf(2)**(mpf(31)/6) * mp.pi**4) * mp.gamma(mpf(1)/3)**6 dg = _digits(cf, ref) ok = dg >= 20 if not ok: failures.append("CM point: only %.1f digits" % dg) print("W_S(1,1,1;3 sqrt(6)/2) closed form %s" % mp.nstr(cf, 25), flush=True) print(" 3/(2^(31/6) pi^4) Gamma(1/3)^6\n" " gamma value %s\n" " agreement %.1f digits, %.1fs %s" % (mp.nstr(ref, 25), float(dg), time.time() - t0, "ok" if ok else "FAIL"), flush=True) print("-" * 78, flush=True) print("total runtime %.1fs" % (time.time() - t_start), flush=True) if failures: print("SELF-TEST FAILED:", flush=True) for msg in failures: print(" " + msg, flush=True) return 1 print("SELF-TEST PASSED: all points >= 20 digits", flush=True) return 0 if __name__ == '__main__': sys.exit(_selftest())