burnol_residue.py
13.5 kB · python · 390 lines
1import math2import sys3import time4from fractions import Fraction56import numpy as np7import mpmath8from mpmath import iv910iv.prec = 12811BOX = iv.mpc(iv.mpf([-1, 1]), iv.mpf([-1, 1]))1213# EXACT PRINTING1415def to_frac(t):16 sign, man, exp, bc = t17 if man == 0 and bc < 0:18 return None19 v = Fraction(man) * (Fraction(2) ** exp)20 return -v if sign else v2122def lo(x):23 return to_frac(x._mpi_[0])2425def hi(x):26 return to_frac(x._mpi_[1])2728def dec(n, d):29 neg = n < 030 n = abs(n)31 ip, fp = divmod(n, 10 ** d)32 return ("-" if neg else "") + f"{ip}.{fp:0{d}d}"3334def floor_str(fr, d):35 return dec((fr.numerator * 10 ** d) // fr.denominator, d)3637def ceil_str(fr, d):38 return dec(-((-fr.numerator * 10 ** d) // fr.denominator), d)3940def fmt_r(x, d):41 return f"[{floor_str(lo(x), d)}, {ceil_str(hi(x), d)}]"4243def fmt_c(z, d):44 return fmt_r(z.real, d) + " + i " + fmt_r(z.imag, d)4546def fmt_up(x, d=3):47 v = hi(x)48 if v <= 0:49 return f"{dec(0, d)}e+00"50 e = 051 while v >= 10:52 v /= 1053 e += 154 while v < 1:55 v *= 1056 e -= 157 m = -((-v.numerator * 10 ** d) // v.denominator)58 if m >= 10 ** (d + 1):59 m //= 1060 e += 161 return f"{dec(m, d)}e{e:+03d}"6263def dist_zero(z, d):64 dx = max(Fraction(0), lo(z.real), -hi(z.real))65 dy = max(Fraction(0), lo(z.imag), -hi(z.imag))66 r2 = dx * dx + dy * dy67 n = math.isqrt((r2.numerator * 10 ** (2 * d)) // r2.denominator)68 return dec(n, d), r2 > 06970def meets(z1, z2):71 ok_re = max(lo(z1.real), lo(z2.real)) <= min(hi(z1.real), hi(z2.real))72 ok_im = max(lo(z1.imag), lo(z2.imag)) <= min(hi(z1.imag), hi(z2.imag))73 return ok_re and ok_im7475def contains(z, cre, cim):76 return lo(z.real) <= cre <= hi(z.real) and lo(z.imag) <= cim <= hi(z.imag)7778def width(z):79 return max(float(hi(z.real) - lo(z.real)), float(hi(z.imag) - lo(z.imag)))8081def rat(p, q):82 return iv.mpf(p) / iv.mpf(q)8384# DESIGN8586class Design:87 def __init__(self, q, F):88 assert 0 in F and len(F) >= 289 self.q = q90 self.F = sorted(F)91 self.N = len(F)92 self.F1 = [a for a in self.F if a > 0]93 self.N1 = len(self.F1)94 self.amin = min(self.F1)95 self.amax = max(self.F1)96 self.logq = iv.log(q)97 self.s0 = iv.log(self.N) / self.logq98 self.name = f"base {q} digits {{{','.join(map(str, self.F))}}}"99100 def gamma(self, l):101 return sum(a ** l for a in self.F)102103 def s(self, k):104 return iv.mpc(self.s0, 2 * iv.pi * k / self.logq)105106 def lev0(self, w):107 acc = iv.mpc(0)108 for a in self.F1:109 acc = acc + (iv.mpc(1) if a == 1 else iv.exp(-w * iv.log(a)))110 return acc111112 def amin_pow(self, i):113 return (iv.mpf(self.amin) ** (-(self.s0 + i))).b114115 def j_start(self):116 j0 = 0117 while 3 * self.amax > self.amin * self.q ** (j0 + 1):118 j0 += 1119 return j0120121 def level(self, j):122 lev = list(self.F1)123 for _ in range(j):124 lev = [self.q * n + a for n in lev for a in self.F]125 return lev126127# TAILS128129def poch_row(a, T):130 pr = [iv.mpf(1)]131 for l in range(T + 1):132 pr.append(pr[-1] * (a + l) / (l + 1))133 return pr134135def hyp_tail(a, rho, T, pr=None):136 if pr is None:137 pr = poch_row(a, T)138 part = iv.mpf(0)139 rp = iv.mpf(1)140 for l in range(T + 1):141 part = part + pr[l] * rp142 rp = rp * rho143 return ((1 - rho) ** (-a) - part).b144145def choose_I(a, rho, eps):146 T = 8147 while True:148 if hi(hyp_tail(a, rho, T)) < eps:149 return T + 2150 T += 2151152# ENGINE ONE: THE LEVEL RECURSION153154def level_limit(D, k, L, I):155 s = D.s(k)156 q, N = D.q, D.N157 coeff, poch, absw = [], [], []158 for i in range(I + 1):159 w = s + i160 aw = abs(w).b161 T = I - i162 row, P = [], iv.mpc(1)163 for l in range(T + 1):164 row.append(P * rat(((-1) ** l) * D.gamma(l), q ** l))165 P = P * (w + l) / (l + 1)166 coeff.append(row)167 poch.append(poch_row(aw, T))168 absw.append(aw)169 j0 = D.j_start()170 lev = D.level(j0)171 z = [iv.exp(-s * iv.log(n)) if n > 1 else iv.mpc(1) for n in lev]172 pw = [iv.mpf(1) for n in lev]173 V = []174 for i in range(I + 1):175 acc = iv.mpc(0)176 for t, n in enumerate(lev):177 acc = acc + z[t] * pw[t]178 pw[t] = pw[t] / n179 V.append(acc)180 mult = [rat(1, N * q ** i) for i in range(I + 1)]181 apow = [D.amin_pow(i) for i in range(I + 1)]182 for j in range(j0, L):183 rho = rat(D.amax, D.amin * q ** (j + 1))184 newV = []185 for i in range(I + 1):186 T = I - i187 row = coeff[i]188 acc = iv.mpc(0)189 for l in range(T + 1):190 acc = acc + row[l] * V[i + l]191 tail = hyp_tail(absw[i], rho, T, poch[i])192 tau = (N * D.N1 * apow[i] * rat(1, q ** (j * i)) * tail).b193 newV.append(mult[i] * (acc + tau * BOX))194 V = newV195 rhoL = rat(D.amax, D.amin * q ** (L + 1))196 a = absw[0]197 E = (D.N1 * apow[0] * a * (1 - rhoL) ** (-(a + 1)) * rhoL * rat(q, q - 1)).b198 R = V[0] + E * BOX199 return R, R / D.logq, E200201# ENGINE TWO: THE FUNCTIONAL EQUATION WITH DIRECT SUMS202203def direct_enclosure(D, k, Lp, M):204 s = D.s(k)205 q, N = D.q, D.N206 nums = [n for j in range(Lp) for n in D.level(j)]207 z = [iv.exp(-s * iv.log(n)) if n > 1 else iv.mpc(1) for n in nums]208 inv = [rat(1, n) for n in nums]209 pw = [iv.mpf(1) for n in nums]210 a = abs(s).b211 R = D.lev0(s)212 P = iv.mpc(1)213 for m in range(1, M + 1):214 P = P * (s + m - 1) / m215 Km = iv.mpc(0)216 for t in range(len(nums)):217 pw[t] = pw[t] * inv[t]218 Km = Km + z[t] * pw[t]219 tailK = (D.N1 * D.amin_pow(m) * rat(1, q ** (Lp * m)) / (1 - rat(1, q ** m))).b220 Km = Km + tailK * BOX221 R = R + P * rat(((-1) ** m) * D.gamma(m), N * q ** m) * Km222 rho0 = rat(D.amax, D.amin * q)223 tailM = (D.N1 * D.amin_pow(0) / (1 - rat(1, q ** (M + 1))) * hyp_tail(a, rho0, M)).b224 R = R + tailM * BOX225 return R, R / D.logq226227# THE COLUMN228229def column(D, k, lam0, mmax):230 s = D.s(k)231 q, N = D.q, D.N232 lam = [lam0]233 for m in range(1, mmax + 1):234 acc = iv.mpc(0)235 fall = iv.mpc(1)236 for j in range(1, m + 1):237 fall = fall * (m - j + 1 - s)238 acc = acc + fall / math.factorial(j) * rat(D.gamma(j), q ** j) * lam[m - j]239 lam.append(-acc / N / (1 - rat(1, q ** m)))240 return lam241242def column_closed(D, k, lam0):243 s = D.s(k)244 mean = Fraction(sum(D.F), D.N)245 var = Fraction(sum(a * a for a in D.F), D.N) - mean * mean246 m1 = mean / (D.q - 1)247 m2 = var / (D.q * D.q - 1) + m1 * m1248 b1 = -m1249 b2 = m1 * m1 - m2 / 2250 c1 = -(s - 1) * lam0 * rat(b1.numerator, b1.denominator)251 c2 = (s - 1) * (s - 2) * lam0 * rat(b2.numerator, b2.denominator)252 return [c1, c2], [b1, b2]253254# CONTROLS IN FLOATS255256def dist_box(z, cre, cim):257 dx = max(Fraction(0), lo(z.real) - cre, cre - hi(z.real))258 dy = max(Fraction(0), lo(z.imag) - cim, cim - hi(z.imag))259 return math.hypot(float(dx), float(dy))260261def level_bound(D, k, J):262 s = math.log(D.N) / math.log(D.q) + 2j * math.pi * k / math.log(D.q)263 a = abs(s)264 rho = D.amax / (D.amin * D.q ** (J + 1))265 return D.N1 * D.amin ** (-math.log(D.N) / math.log(D.q)) * a * (1 - rho) ** (-a - 1) * rho * D.q / (D.q - 1) / math.log(D.q)266267def control_levels(D, k, J):268 s = math.log(D.N) / math.log(D.q) + 2j * math.pi * k / math.log(D.q)269 arr = np.array(D.F1, dtype=np.int64)270 out = []271 J = min(J, int(20 * math.log(2) / math.log(D.N)))272 for j in range(J + 1):273 if j > 0:274 arr = np.concatenate([D.q * arr + a for a in D.F])275 out.append(complex(np.sum(np.exp(-s * np.log(arr.astype(np.float64))))))276 return out277278def count_le(xs, D):279 q, F, N = D.q, D.F, D.N280 x = xs.copy()281 nd = 1 + int(math.log(float(xs.max())) / math.log(q)) + 1282 digits = np.zeros((len(xs), nd), dtype=np.int64)283 for p in range(nd):284 digits[:, p] = x % q285 x = x // q286 cnt = np.zeros(len(xs), dtype=np.int64)287 tight = np.ones(len(xs), dtype=bool)288 for p in reversed(range(nd)):289 dp = digits[:, p]290 less = sum((dp > a).astype(np.int64) for a in F)291 cnt += np.where(tight, less * (N ** p), 0)292 tight &= np.isin(dp, F)293 cnt += tight294 return cnt - 1295296def control_fourier(D, k, J, M):297 u = (np.arange(M) + 0.5) / M298 x = np.floor(float(D.q) ** (u + J)).astype(np.int64)299 A = count_le(x, D)300 Phi = A / (float(D.N) ** (u + J))301 ck = complex(np.mean(Phi * np.exp(-2j * np.pi * k * u)))302 s0 = math.log(D.N) / math.log(D.q)303 s = s0 + 2j * math.pi * k / math.log(D.q)304 top = u >= 1 - s0305 dev = float(np.abs(Phi[top] - D.N ** (1 - u[top])).max()) if D.F == [0, 1] and D.q == 3 else float("nan")306 return ck, s * ck * math.log(D.q), float(Phi.min()), float(Phi.max()), dev307308# REPORT309310def report(D, k, L, second=False, controls=False, col=0, d=15):311 s = D.s(k)312 a = abs(s).b313 j0 = D.j_start()314 I = choose_I(a, rat(D.amax, D.amin * D.q ** (j0 + 1)), Fraction(1, 10 ** 24))315 t0 = time.time()316 R, lam, E = level_limit(D, k, L, I)317 t1 = time.time()318 print(f"k = {k}: s = {fmt_c(s, 12)}, recursion from level {j0} ({len(D.level(j0))} terms summed directly) to L = {L}, I = {I}, {t1 - t0:.1f} s")319 print(f" R_k = lim_j Lev_j(s) = {fmt_c(R, d)}")320 print(f" lambda_(0,k) = R_k/log q = {fmt_c(lam, d)}")321 dz, excl = dist_zero(lam, 9)322 print(f" |lambda_(0,k)| >= {dz}, excludes zero: {excl}, width {width(lam):.1e}, level tail E_L = {fmt_up(E)}")323 out = {"R": R, "lam": lam, "I": I, "excl": excl, "dz": dz}324 if second:325 t0 = time.time()326 R2, lam2 = direct_enclosure(D, k, 13, 36)327 print(f" engine two (direct sums to level 13, M = 36, {time.time() - t0:.1f} s): lambda = {fmt_c(lam2, 8)}, width {width(lam2):.1e}, meets engine one: {meets(lam, lam2)}")328 if controls:329 lv = control_levels(D, k, 20)330 J = len(lv) - 1331 lamf = [v / math.log(D.q) for v in lv]332 dist = dist_box(lam, Fraction(lamf[J].real), Fraction(lamf[J].imag))333 print(f" control, level sums by enumeration (Burnol 5.1), lambda from Lev_{J - 2}, Lev_{J - 1}, Lev_{J}: "334 + ", ".join(f"{v.real:.9f}{v.imag:+.9f}i" for v in lamf[J - 2:J + 1]))335 print(f" successive differences |Lev_(j+1) - Lev_j|, j = {J - 4}..{J - 1}: "336 + ", ".join(f"{abs(lv[j + 1] - lv[j]):.2e}" for j in range(J - 4, J))337 + f"; distance of Lev_{J}/log q from engine one {dist:.1e}, its own bound E_{J}/log q = {level_bound(D, k, J):.1e}, within: {dist <= level_bound(D, k, J)}")338 for M in (2 ** 15, 2 ** 17):339 ck, Rf, pmin, pmax, dev = control_fourier(D, k, 30, M)340 lf = Rf / math.log(D.q)341 print(f" control, Fourier coefficient of Phi on {M} midpoints, J = 30: c_k = {ck.real:.9f}{ck.imag:+.9f}i, "342 f"lambda = s c_k = {lf.real:.9f}{lf.imag:+.9f}i, distance from engine one {dist_box(lam, Fraction(lf.real), Fraction(lf.imag)):.1e}; "343 f"Phi in [{pmin:.6f}, {pmax:.6f}], max deviation of Phi from N^(1-u) on [1 - s_0, 1]: {dev:.1e}")344 if col:345 lams = column(D, k, lam, col)346 closed, bs = column_closed(D, k, lam)347 for m in range(1, col + 1):348 line = f" lambda_({m},k) = {fmt_c(lams[m], 10)}"349 if m <= 2:350 line += f", closed form via Theorem 7.4 with [t^{m}](1/E) = {bs[m - 1]}: {fmt_c(closed[m - 1], 10)}, meets: {meets(lams[m], closed[m - 1])}"351 print(line)352 return out353354def main():355 band = "--band" in sys.argv356 t_all = time.time()357 print("BURNOL RESIDUE: certified enclosures of Res K at s_(0,k) = log_q N + 2 pi i k/log q, mpmath.iv at 128 bits")358 print("rests on: Burnol 2026 Proposition 5.1 (lambda_(0,k) log q = lim_j Lev_j(s_(0,k))), the level recursion, the two tail lemmas, outward-rounded interval arithmetic")359 D = Design(3, [0, 1])360 print(f"\n{D.name}: N = {D.N}, s_0 = {fmt_r(D.s0, 15)}, log 3 = {fmt_r(D.logq, 15)}, 2^(s_0) = {fmt_r(iv.mpf(2) ** D.s0, 9)}, 2^(1 - s_0) = {fmt_r(iv.mpf(2) ** (1 - D.s0), 9)}")361 L = 40362 res = {}363 res[0] = report(D, 0, L)364 res[1] = report(D, 1, L, second=True, controls=True, col=3)365 ks = list(range(2, 11)) if band else []366 for k in ks:367 res[k] = report(D, k, L, controls=True)368 print("\nBAND (printed by the generator; endpoints floored and ceiled at 12 decimals; |lambda| truncated at 9)")369 print("| `k` | `Re lambda_(0,k)` | `Im lambda_(0,k)` | `abs(lambda_(0,k)) >=` | zero excluded |")370 print("|---|---|---|---|---|")371 for k in sorted(res):372 z = res[k]["lam"]373 print(f"| {k} | `{fmt_r(z.real, 12)}` | `{fmt_r(z.imag, 12)}` | `{res[k]['dz']}` | {'yes' if res[k]['excl'] else 'NO'} |")374375 D2 = Design(3, [0, 2])376 print(f"\n{D2.name}: N = {D2.N}; every element is twice one of {{0,1}}, so lambda scales by 2^(-s_(0,k))")377 r2 = report(D2, 1, L, controls=True)378 scaled = res[1]["lam"] * iv.exp(-D.s(1) * iv.log(2))379 print(f" 2^(-s) times the {{0,1}} enclosure = {fmt_c(scaled, 15)}, meets the {{0,2}} enclosure: {meets(scaled, r2['lam'])}")380381 D3 = Design(3, [0, 1, 2])382 print(f"\n{D3.name}: K is the Riemann zeta function, s_0 = 1; the residue at 1 is 1 and every off-real lattice point is regular")383 r30 = report(D3, 0, L)384 print(f" contains 1: {contains(r30['lam'], Fraction(1), Fraction(0))}")385 r31 = report(D3, 1, L, controls=True)386 print(f" contains 0: {contains(r31['lam'], Fraction(0), Fraction(0))}")387 print(f"\ntotal {time.time() - t_all:.1f} s")388389if __name__ == "__main__":390 main()