#!/usr/bin/env python3 # hexabox_vop_lib.py — exact closed-function algebra for the sec255 layer # recursion (support library for hexabox_vop.py), shipped with the BootLoops # hexabox page so the closure re-derives at runtime with no external imports. # Public deps: mpmath, python-flint. # # CHANGELOG 2026-07-05: certified-bound support for the refine loops — # tail_bound_geom (trailing-window geometric tail bound), Engine keeps the value-> # Chebyshev-coefficient matrix (self._CM) and exposes coeff_tail(); Engine.interp's # node-coincidence snap now tracks mp.mp.dps (was a fixed 1e-60). # """vop_lib.py — exact closed-function algebra for the sec255 VoP layer closure. v2 Function model: F(t) = sum_terms C_tag * sqrt(rad) * R(t) * G_word(t) - C_tag: named constant (registry CONSTS) - rad: tuple of letter indices, sqrt(rad) = prod_a sqrt(l_a(t)/l_a(0)) continued along the standard detour contour (=1 at t=0) - R: exact rational function (flint fmpq_poly pair) - G_word: iterated integral int_0^t K_1 (int K_2 (...)) with exact kernels; kernel keys: ('L', li)=dlog letter, ('R', ki)=rational remainder kernel, ('A', ai)=algebraic kernel Rat*sqrt(rad). G_()=1; G_word(0)=0. Closed under +, Rat*, sqrt-mul, d/dt (exact), int_0^t (exact Hermite for rational prefactors; algebraic prefactors absorbed into registered kernels). """ import time from fractions import Fraction import flint import mpmath as mp FQ = flint.fmpq FP = flint.fmpq_poly def fq(x): if isinstance(x, FQ): return x f = Fraction(x) return FQ(f.numerator, f.denominator) def fp(coeffs): # low-to-high return FP([fq(c) for c in coeffs]) ZERO_P = FP([0]); ONE_P = FP([1]) # ───────────────────────── Rat ───────────────────────── class Rat: __slots__ = ('n', 'd') def __init__(self, n, d=None): if d is None: d = ONE_P if not isinstance(n, FP): n = fp(n if isinstance(n, (list, tuple)) else [n]) if not isinstance(d, FP): d = fp(d if isinstance(d, (list, tuple)) else [d]) if d == ZERO_P: raise ZeroDivisionError if n == ZERO_P: self.n, self.d = ZERO_P, ONE_P; return g = n.gcd(d) if g.degree() > 0: n, d = n // g, d // g lc = d.coeffs()[-1] if lc != 1: n, d = n * FP([1/lc]), d * FP([1/lc]) self.n, self.d = n, d def is_zero(self): return self.n == ZERO_P def __add__(a, b): return Rat(a.n*b.d + b.n*a.d, a.d*b.d) def __sub__(a, b): return Rat(a.n*b.d - b.n*a.d, a.d*b.d) def __mul__(a, b): if isinstance(b, FQ): return Rat(a.n*FP([b]), a.d) return Rat(a.n*b.n, a.d*b.d) def __neg__(a): return Rat(-a.n, a.d) def inv(a): return Rat(a.d, a.n) def deriv(a): return Rat(a.n.derivative()*a.d - a.n*a.d.derivative(), a.d*a.d) def __eq__(a, b): return a.n == b.n and a.d == b.d def key(a): return (tuple(str(c) for c in a.n.coeffs()), tuple(str(c) for c in a.d.coeffs())) def __hash__(a): return hash(a.key()) def ev_fq(a, x): x = fq(x); return a.n(x) / a.d(x) def ev_mpc(a, z): num = mp.mpc(0) for c in reversed(a.n.coeffs()): num = num*z + mp.mpf(int(c.p))/int(c.q) den = mp.mpc(0) for c in reversed(a.d.coeffs()): den = den*z + mp.mpf(int(c.p))/int(c.q) return num/den def __repr__(a): return f"({a.n})/({a.d})" R_ZERO = Rat(ZERO_P); R_ONE = Rat(ONE_P) # ─────────────────── registries ─────────────────── LETTERS = []; _LK = {} def reg_letter(f): lc = f.coeffs()[-1] if lc != 1: f = f * FP([1/lc]) k = tuple(str(c) for c in f.coeffs()) if k in _LK: return _LK[k] LETTERS.append(f); _LK[k] = len(LETTERS)-1 return len(LETTERS)-1 RKERNS = []; _RK = {} def reg_rkern(r): k = r.key() if k in _RK: return _RK[k] RKERNS.append(r); _RK[k] = len(RKERNS)-1 return len(RKERNS)-1 AKERNS = []; _AK = {} # (Rat, rad) def reg_akern(r, rad): k = (r.key(), rad) if k in _AK: return _AK[k] AKERNS.append((r, rad)); _AK[k] = len(AKERNS)-1 return len(AKERNS)-1 CONSTS = {'1': mp.mpc(1)} def reg_const(tag, val): if tag in CONSTS: assert abs(CONSTS[tag]-val) < mp.mpf(10)**(-mp.mp.dps+12), f"const clash {tag}" CONSTS[tag] = val return tag def rad_mul(rad1, rad2): """sqrt(rad1)*sqrt(rad2) = Rfac * sqrt(rad). Rfac = prod_common l(t)/l(0).""" s1, s2 = set(rad1), set(rad2) common = s1 & s2 rad = tuple(sorted(s1 ^ s2)) R = R_ONE for li in common: f = LETTERS[li] R = R * Rat(f, FP([f(fq(0))])) return R, rad def rad_dlog_half(rad): """(1/2) sum dlog l_a as Rat.""" R = R_ZERO for li in rad: f = LETTERS[li] R = R + Rat(f.derivative(), f * FP([2])) return R def kern_eval_parts(kk): """kernel -> (Rat, rad)""" if kk[0] == 'L': f = LETTERS[kk[1]]; return Rat(f.derivative(), f), () if kk[0] == 'R': return RKERNS[kk[1]], () return AKERNS[kk[1]] # ─────────────────── Hermite split ─────────────────── _HCACHE = {} def hermite(r): """r = H' + sum c_a dlog(l_a) + sum c_k K_k (registered simple-pole kernels).""" ck = r.key() if ck in _HCACHE: return _HCACHE[ck] H = R_ZERO; logs = {}; rems = [] q, rem = divmod(r.n, r.d) if q != ZERO_P: qc = q.coeffs() H = H + Rat(FP([FQ(0)] + [qc[i]/(i+1) for i in range(len(qc))])) if rem != ZERO_P: c, facs = r.d.factor() num = rem * FP([1/c]) pieces = [] rest_facs = list(facs) while len(rest_facs) > 1: f, e = rest_facs.pop() fe = f**e D2 = ONE_P for (g, eg) in rest_facs: D2 = D2 * g**eg g0, a, b = fe.xgcd(D2) assert g0.degree() == 0 a, b = a * FP([1/g0.coeffs()[0]]), b * FP([1/g0.coeffs()[0]]) pieces.append(((num * b) % fe, f, e)) num = (num * a) % D2 f, e = rest_facs[0] pieces.append((num % (f**e), f, e)) for (nf, f, e) in pieces: fp_ = f.derivative() g0, u, v = fp_.xgcd(f) u = u * FP([1/g0.coeffs()[0]]) while e > 1: s = (nf * u) % f w = (nf - s * fp_) // f H = H + Rat(-s * FP([FQ(1, e-1)]), f**(e-1)) nf = w + s.derivative() * FP([FQ(1, e-1)]) e -= 1 if nf == ZERO_P: continue dn, df = nf.degree(), f.degree() cl = FQ(0) if dn == df - 1: cl = nf.coeffs()[-1] / (fq(df) * f.coeffs()[-1]) if cl != 0: li = reg_letter(f) logs[li] = logs.get(li, FQ(0)) + cl nf = nf - fp_ * FP([cl]) if nf != ZERO_P: rr = Rat(nf, f) lead = rr.n.coeffs()[-1] rk = reg_rkern(Rat(rr.n * FP([1/lead]), rr.d)) rems.append((lead, rk)) out = (H, logs, rems) _HCACHE[ck] = out return out # ─────────────────── containers: {(ctag, rad, word): Rat} ─────────────────── def cf_zero(): return {} def cf_const(ctag, r=None): return {(ctag, (), ()): (r if r is not None else R_ONE)} def cf_add(A, B): out = dict(A) for k, r in B.items(): if k in out: s = out[k] + r if s.is_zero(): del out[k] else: out[k] = s elif not r.is_zero(): out[k] = r return out def _neg_cf(A): return {k: -r for k, r in A.items()} def cf_scale(A, c): c = fq(c) if c == 0: return {} return {k: r * c for k, r in A.items()} def cf_mulrat(A, R): if R.is_zero(): return {} return {k: r * R for k, r in A.items()} def cf_mulalg(A, R, rad): """A * R*sqrt(rad)""" if R.is_zero(): return {} out = cf_zero() for (ctag, rd, word), r in A.items(): Rf, nrd = rad_mul(rd, rad) out = cf_add(out, {(ctag, nrd, word): r * R * Rf}) return out def cf_diff(A): out = cf_zero() for (ctag, rad, word), r in A.items(): rp = r.deriv() if rad: rp = rp + r * rad_dlog_half(rad) if not rp.is_zero(): out = cf_add(out, {(ctag, rad, word): rp}) if word: K, radk = kern_eval_parts(word[0]) Rf, nrd = rad_mul(rad, radk) out = cf_add(out, {(ctag, nrd, word[1:]): r * K * Rf}) return out def cf_int(A): out = cf_zero() for (ctag, rad, word), r in A.items(): out = cf_add(out, _int_term(ctag, rad, word, r)) return out def _int_term(ctag, rad, word, r): if rad: # absorb into an algebraic kernel: normalize scalar lead = r.n.coeffs()[-1] ai = reg_akern(Rat(r.n * FP([1/lead]), r.d), rad) return {(ctag, (), (('A', ai),) + word): Rat(FP([lead]))} H, logs, rems = hermite(r) out = cf_zero() if not H.is_zero(): out = cf_add(out, {(ctag, (), word): H}) if word == (): h0 = H.ev_fq(0) if h0 != 0: out = cf_add(out, {(ctag, (), ()): Rat(FP([-h0]))}) else: K, radk = kern_eval_parts(word[0]) out = cf_add(out, _neg_cf(_int_term(ctag, radk, word[1:], H * K))) for li, cl in logs.items(): out = cf_add(out, {(ctag, (), (('L', li),) + word): Rat(FP([cl]))}) for (cl, rk) in rems: out = cf_add(out, {(ctag, (), (('R', rk),) + word): Rat(FP([cl]))}) return out def cf_key_report(A): words = {} for (ctag, rad, word), r in A.items(): words.setdefault((rad, word), 0) words[(rad, word)] += 1 return words # ─────────────────── rational solutions of Y' = A Y ─────────────────── def solve_rat_solutions(Ain, extra_rad=None, dN_extra=30, eig_fn=None): """Exact rational solution vectors of Y' = (A - 1/2 dlog(g) I) Y. A: nb x nb Rat matrix. eig_fn(f) -> numeric residue-matrix eigenvalues at a root of irreducible factor f (used only for the denominator BOUND; result is exactly certified downstream).""" nb = len(Ain) A = [row[:] for row in Ain] if extra_rad: g = ONE_P for li in extra_rad: g = g * LETTERS[li] half = Rat(g.derivative(), g*FP([2])) for a in range(nb): A[a][a] = A[a][a] - half facs = {} for a in range(nb): for b in range(nb): if A[a][b].is_zero(): continue c, fl = A[a][b].d.factor() for f, e in fl: k2 = tuple(str(x) for x in f.coeffs()) if k2 not in facs: facs[k2] = [f, e] else: facs[k2][1] = max(facs[k2][1], e) D = ONE_P for k2, (f, emax) in facs.items(): ev = eig_fn(f) if (eig_fn is not None and emax == 1) else [] if ev != [] and emax == 1: lam = min(float(mp.re(x)) for x in ev) pf = max(0, int(mp.ceil(-lam - 0.25))) else: pf = emax + 1 if extra_rad: pf = pf + 1 if pf > 0: D = D * f**pf dN = D.degree() + dN_extra Vv = []; Wv = []; Uv = [] for a in range(nb): La = ONE_P for b in range(nb): if A[a][b].is_zero(): continue gg = La.gcd(A[a][b].d) La = La * (A[a][b].d // gg) Vv.append(La * D); Wv.append(-(La * D.derivative())) U = [] for b in range(nb): if A[a][b].is_zero(): U.append(None); continue U.append(-(D * A[a][b].n * (La // A[a][b].d))) Uv.append(U) maxdeg = max(Vv[a].degree() for a in range(nb)) + dN + 2 ncols = nb*(dN+1) rowsM = {} def addpoly(a, P, col): for kk, c in enumerate(P.coeffs()): if c == 0: continue key = a*(maxdeg+1)+kk rowsM.setdefault(key, {}) rowsM[key][col] = rowsM[key].get(col, FQ(0)) + c tpow = [FP([0]*k2+[1]) for k2 in range(dN+2)] for b in range(nb): for m2 in range(dN+1): col = b*(dN+1)+m2 for a in range(nb): acc = ZERO_P if Uv[a][b] is not None: acc = acc + Uv[a][b]*tpow[m2] if a == b: if m2 > 0: acc = acc + Vv[a]*FP([fq(m2)])*tpow[m2-1] acc = acc + Wv[a]*tpow[m2] if acc != ZERO_P: addpoly(a, acc, col) rowkeys = sorted(rowsM) Mz = flint.fmpq_mat(len(rowkeys), ncols) for ri, rk in enumerate(rowkeys): for col, c in rowsM[rk].items(): Mz[ri, col] = c R2, rank = Mz.rref() piv = []; ci = 0 for ri in range(rank): while ci < ncols and R2[ri, ci] == 0: ci += 1 piv.append(ci) free = [c2 for c2 in range(ncols) if c2 not in piv] sols = [] for fcol in free: vec = [FQ(0)]*ncols vec[fcol] = FQ(1) for ri in range(rank-1, -1, -1): pc = piv[ri] s = FQ(0) for c2 in range(pc+1, ncols): if R2[ri, c2] != 0 and vec[c2] != 0: s += R2[ri, c2]*vec[c2] vec[pc] = -s/R2[ri, pc] Nvecs = [Rat(fp(vec[b*(dN+1):(b+1)*(dN+1)]), D) for b in range(nb)] if all(x.is_zero() for x in Nvecs): continue sols.append(Nvecs) return sols # ─────────────────── numeric contour engine ─────────────────── def _mpf2arb(x): s, m, e, _ = mp.mpf(x)._mpf_ v = flint.arb(-int(m) if s else int(m)) return v*(flint.arb(2)**int(e)) if e else v def _mpc2acb(z): z = mp.mpc(z); return flint.acb(_mpf2arb(z.real), _mpf2arb(z.imag)) def _acb2mpc(z): def f(a): m, e = a.man_exp(); return mp.mpf(int(m))*mp.mpf(2)**int(e) return mp.mpc(f(z.real.mid()), f(z.imag.mid())) def tail_bound_geom(tail, prev, window=8): """Trailing-window geometric tail bound (2026-07-05): tail * r/(1-r) with the per-coefficient ratio r measured from the last two `window`-blocks of the Chebyshev spectrum, clipped to [1/4, 3/4]. Poor decay (r >= 3/4) gives the conservative factor 3, which fails the caller's tolerance test and escalates the refine loop.""" tail = mp.mpf(tail); prev = mp.mpf(prev) if tail == 0: return mp.mpf(0) if prev > 0 and tail < prev: r = (tail/prev)**(mp.mpf(1)/window) else: r = mp.mpf(1) r = min(max(r, mp.mpf(1)/4), mp.mpf(3)/4) return tail*r/(1-r) class Engine: def __init__(self, real_poles, NC=120, delta='0.08', pad='0.04'): self.NC = NC flint.ctx.prec = int(3.33*mp.mp.dps) + 20 delta = mp.mpf(delta); pad = mp.mpf(pad) rp = sorted(set(round(float(x), 10) for x in real_poles if 1e-8 < float(x) < 1-1e-8)) clusters = [] for x in rp: if clusters and x - clusters[-1][1] < 0.08: clusters[-1] = (clusters[-1][0], x) else: clusters.append((x, x)) verts = [mp.mpc(0)] for (a, b) in clusters: a_ = max(mp.mpf(a)-pad, mp.mpf(float(verts[-1].real))+mp.mpf('0.005')) b_ = min(mp.mpf(b)+pad, mp.mpf(1)-mp.mpf('0.005')) if a_ >= b_: continue verts += [mp.mpc(a_), mp.mpc(a_, delta), mp.mpc(b_, delta), mp.mpc(b_)] verts.append(mp.mpc(1)) self.segs = [(verts[k], verts[k+1]) for k in range(len(verts)-1)] self.uch = [(1-mp.cos(mp.pi*k/NC))/2 for k in range(NC+1)] self._COS = [[mp.cos(mp.pi*m*k/NC) for k in range(NC+1)] for m in range(NC+2)] self.zg = [[za+u*(zb-za) for u in self.uch] for (za, zb) in self.segs] self.zga = [[_mpc2acb(z) for z in row] for row in self.zg] self.dz = [zb-za for (za, zb) in self.segs] self.dza = [_mpc2acb(d) for d in self.dz] self._wcache = {} self._radcache = {} self._kcache = {} # spectral cumint as acb matrix product: out = W (C vals) N = NC Cm = flint.acb_mat(N+1, N+1) for m2 in range(N+1): for k in range(N+1): w = mp.mpf(1)/2 if k in (0, N) else mp.mpf(1) Cm[m2, k] = _mpc2acb(2*w*self._COS[m2][k]/N) Wm = flint.acb_mat(N+1, N+1) for k in range(N+1): Wm[k, 0] = _mpc2acb(self.uch[k]/2) Wm[k, 1] = _mpc2acb((1-self._COS[2][k])/8) for m2 in range(2, N+1): v = ((1-self._COS[m2+1][k])/(m2+1) - (1-self._COS[m2-1][k])/(m2-1))/4 if m2 == N: v = v/2 Wm[k, m2] = _mpc2acb(v) self._CUM = Wm*Cm self._CM = Cm # 2026-07-05: value -> Chebyshev-coefficient matrix (coeff_tail) def coeff_tail(self, vals_acb, window=8): """Trailing-window Chebyshev coefficient magnitudes of one node-table segment (2026-07-05, certified-bound support): returns (max |c_j| over the last `window` coefficients, max over the preceding `window`). Uses the unhalved DCT row for c_N (up to 2x conservative — a BOUND, never an underestimate).""" N = self.NC V = flint.acb_mat(N+1, 1) for k in range(N+1): V[k, 0] = vals_acb[k] C = self._CM*V lo = max(0, N+1-2*window) mags = [abs(_acb2mpc(C[j, 0])) for j in range(lo, N+1)] if len(mags) <= window: return max(mags), max(mags) return max(mags[-window:]), max(mags[:-window]) def _cumcheb_acb(self, vals_acb): N = self.NC V = flint.acb_mat(N+1, 1) for k in range(N+1): V[k, 0] = vals_acb[k] O = self._CUM*V return [O[k, 0] for k in range(N+1)] def _cumcheb(self, vals): return [_acb2mpc(x) for x in self._cumcheb_acb([_mpc2acb(v) for v in vals])] def kern_vals(self, kk): """acb node values of kernel kk (incl. continued radical) per segment.""" if kk in self._kcache: return self._kcache[kk] K, radk = kern_eval_parts(kk) rv = self.rad_vals(radk) if radk else None ncf = [_mpc2acb(mp.mpf(int(c.p))/int(c.q)) for c in reversed(K.n.coeffs())] dcf = [_mpc2acb(mp.mpf(int(c.p))/int(c.q)) for c in reversed(K.d.coeffs())] out = [] for s in range(len(self.segs)): row = [] for k in range(self.NC+1): z = self.zga[s][k] num = flint.acb(0) for c in ncf: num = num*z + c den = flint.acb(0) for c in dcf: den = den*z + c v = num/den if rv is not None: v = v*_mpc2acb(rv[s][k]) row.append(v) out.append(row) self._kcache[kk] = out return out def rad_vals(self, rad): """continued sqrt(prod l(z)/l(0)) at all nodes.""" if rad in self._radcache: return self._radcache[rad] w = [] for s in range(len(self.segs)): row = [] for k in range(self.NC+1): z = self.zg[s][k] v = mp.mpc(1) for li in rad: f = LETTERS[li] num = mp.mpc(0) for c in reversed(f.coeffs()): num = num*z + mp.mpf(int(c.p))/int(c.q) f0 = f(fq(0)) v *= num/(mp.mpf(int(f0.p))/int(f0.q)) row.append(v) w.append(row) out = []; prev = mp.mpc(1); prevw = mp.mpc(1) for s in range(len(self.segs)): row = [] for k in range(self.NC+1): cur = prev * mp.sqrt(w[s][k]/prevw) row.append(cur); prev = cur; prevw = w[s][k] out.append(row) self._radcache[rad] = out return out def word_vals_acb(self, word): if word in self._wcache: return self._wcache[word] if word == (): one = flint.acb(1) sv = [[one]*(self.NC+1) for _ in self.segs] self._wcache[word] = sv return sv inner = self.word_vals_acb(word[1:]) kv = self.kern_vals(word[0]) sv = []; cum = flint.acb(0) for s in range(len(self.segs)): dz = self.dza[s] g = [kv[s][k]*inner[s][k]*dz for k in range(self.NC+1)] segcum = self._cumcheb_acb(g) sv.append([cum+segcum[k] for k in range(self.NC+1)]) cum += segcum[-1] self._wcache[word] = sv return sv def word_vals(self, word): key = ('mpc', word) if key in self._wcache: return self._wcache[key] out = [[_acb2mpc(x) for x in row] for row in self.word_vals_acb(word)] self._wcache[key] = out return out def locate(self, z): for s, (za, zb) in enumerate(self.segs): d = self.dz[s] u = ((z-za)/d).real offax = abs(za + u*d - z) if -1e-12 <= u <= 1+1e-12 and offax < 1e-10: return s, mp.mpf(u) raise ValueError(f"point {z} not on contour") def interp(self, vals, u): num = mp.mpc(0); den = mp.mpc(0) for k in range(self.NC+1): d = u - self.uch[k] # node-coincidence snap tracks dps (2026-07-05; was a fixed 1e-60 cutoff) if abs(d) < mp.mpf(10)**(-mp.mp.dps-10): return vals[k] w = (-1)**k*(mp.mpf(1)/2 if k in (0, self.NC) else mp.mpf(1))/d num += w*vals[k]; den += w return num/den def eval_cf(self, A, z): z = mp.mpc(z) s, u = self.locate(z) tot = mp.mpc(0) for (ctag, rad, word), r in A.items(): gv = self.interp(self.word_vals(word)[s], u) if word else mp.mpc(1) rv = self.interp(self.rad_vals(rad)[s], u) if rad else mp.mpc(1) tot += CONSTS[ctag]*r.ev_mpc(z)*gv*rv return tot def eval_cf_end(self, A): return self.eval_cf(A, mp.mpc(1)) # ─────────────────── self-test ─────────────────── def _selftest(): mp.mp.dps = 60 r = Rat(fp([1]), fp([2, 1])) A = {('1', (), ()): r} I = cf_int(A) D = cf_diff(I) dr = cf_add(D, _neg_cf(A)) assert all(v.is_zero() for v in dr.values()) eng = Engine([], NC=96) assert abs(eng.eval_cf_end(I) - mp.log(mp.mpf(3)/2)) < mp.mpf(10)**-45 I2 = cf_int(cf_mulrat(I, Rat(fp([1]), fp([1, 1])))) ref2 = mp.quad(lambda s: mp.log((s+2)/2)/(s+1), [0, 1]) assert abs(eng.eval_cf_end(I2)-ref2) < mp.mpf(10)**-40 r3 = Rat(fp([1, 0, 3]), fp([2, 1])**2 * fp([5, 1, 1])) I3 = cf_int({('1', (), ()): r3}) dr3 = cf_add(cf_diff(I3), _neg_cf({('1', (), ()): r3})) assert all(v.is_zero() for v in dr3.values()) ref3 = mp.quad(lambda s: (3*s**2+1)/((s+2)**2*(s**2+s+5)), [0, 1]) assert abs(eng.eval_cf_end(I3)-ref3) < mp.mpf(10)**-40 # radical tests: y = sqrt((t+3)/3) satisfies y' = 1/2 dlog(t+3) y li = reg_letter(fp([3, 1])) Y = {('1', (li,), ()): R_ONE} D4 = cf_diff(Y) ref = cf_mulalg(Y, Rat(fp([1]), fp([6, 2])), ()) dr4 = cf_add(D4, _neg_cf(ref)) assert all(v.is_zero() for v in dr4.values()), dr4 # numeric radical continuation: value at 1 = sqrt(4/3) v4 = eng.eval_cf(Y, mp.mpc(1)) assert abs(v4 - mp.sqrt(mp.mpf(4)/3)) < mp.mpf(10)**-50 # int of algebraic: int_0^1 sqrt((s+3)/3)/(s+1) ds I5 = cf_int(cf_mulrat(Y, Rat(fp([1]), fp([1, 1])))) ref5 = mp.quad(lambda s: mp.sqrt((s+3)/3)/(s+1), [0, 1]) assert abs(eng.eval_cf_end(I5)-ref5) < mp.mpf(10)**-40 # negative-letter radical branch: l = 2/5 - t crosses zero; contour detours above. li2 = reg_letter(fp([Fraction(2,5), -1])) Y2 = {('1', (li2,), ()): R_ONE} eng2 = Engine([0.4], NC=96) v6 = eng2.eval_cf(Y2, mp.mpc(1)) # continued sqrt((2/5-t)/(2/5)) at t=1: ratio = -3/2, from above: exp(-i pi/2)*sqrt(3/2) ref6 = mp.sqrt(mp.mpf(3)/2)*mp.exp(-1j*mp.pi/2) assert abs(v6-ref6) < mp.mpf(10)**-45, (v6, ref6) print("vop_lib v2 selftest PASS", flush=True) if __name__ == "__main__": _selftest()