Audit finite Weil-form candidate with Arb

.github/workflows/rh-verifier.yml

@@ -0,0 +1,20 @@
+name: RH Arb verifier
+
+on:
+  push:
+    branches: [rh-verifier]
+  pull_request:
+  workflow_dispatch:
+
+jobs:
+  verify:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.12'
+      - name: Install python-flint
+        run: python -m pip install --upgrade pip python-flint
+      - name: Run interval verifier
+        run: python rh_verifier/verify_candidate.py

rh_verifier/verify_candidate.py

@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+"""Independent Arb audit of a finite Weil-form candidate.
+
+This script uses python-flint real balls.  It does not import spreadsheet data.
+The tested primitive is
+
+    F(t) = sin(pi*t/5) - sin(3*pi*t/5)/3,   0 <= t <= 5,
+
+and the zero-mean test function is phi = -F'.
+"""
+from flint import arb, ctx, fmpq
+
+ctx.prec = 256
+L = arb(5)
+pi = arb.pi()
+k = pi / L
+N_LERCH = 5001  # j = 0,...,5000
+
+
+def H(r, m, n):
+    """Integral of exp(r*|t-u|) sin(mkt) sin(nku)."""
+    wm = arb(m) * k
+    wn = arb(n) * k
+    sm = -1 if m % 2 else 1
+    sn = -1 if n % 2 else 1
+    dm = wm*wm + r*r
+    dn = wn*wn + r*r
+    ans = arb(0)
+    if m == n:
+        ans += -r * L / dn
+    ans += wm * wn * (arb(1 + sm*sn) - (r*L).exp() * arb(sm + sn)) / (dm * dn)
+    return ans
+
+
+def V(r, m, n):
+    """Contribution of exp(r|h|) after two integrations by parts."""
+    ans = -(r*r) * H(r, m, n)
+    if m == n:
+        ans -= r * L
+    return ans
+
+
+def shifted_overlap(m, n, a):
+    """Integral_0^(L-a) sin(mk(t+a)) sin(nkt) dt."""
+    b = L - a
+    wm = arb(m) * k
+    wn = arb(n) * k
+    if m == n:
+        first = b * (wm*a).cos()
+    else:
+        first = (((wm-wn)*b + wm*a).sin() - (wm*a).sin()) / (wm-wn)
+    second = (((wm+wn)*b + wm*a).sin() - (wm*a).sin()) / (wm+wn)
+    return (first - second) / 2
+
+
+def is_prime(n):
+    if n < 2:
+        return False
+    p = 2
+    while p*p <= n:
+        if n % p == 0:
+            return False
+        p += 1
+    return True
+
+
+def prime_power_base(n):
+    """Return p if n is p^a for a prime p, else None."""
+    for p in range(2, n+1):
+        if not is_prime(p):
+            continue
+        q = p
+        while q < n:
+            q *= p
+        if q == n:
+            return p
+    return None
+
+
+def matrix_entry(m, n):
+    g1 = arb(fmpq(1,4)).digamma() - pi.log()
+    arch = -4*V(arb(fmpq(1,2)), m, n) - 4*V(arb(fmpq(-1,2)), m, n)
+    if m == n:
+        arch += g1 * L / 2
+    for j in range(N_LERCH):
+        q = arb(fmpq(4*j+1, 4))
+        arch += V(-2*q, m, n) / (4*q*q)
+
+    prime = arb(0)
+    # exp(5) is between 148 and 149, so n <= 148 exactly.
+    for value in range(2, 149):
+        p = prime_power_base(value)
+        if p is None:
+            continue
+        a = arb(value).log()
+        weight = arb(p).log() / arb(value).sqrt()
+        prime -= weight * (shifted_overlap(m, n, a) + shifted_overlap(n, m, a))
+    return arch, prime
+
+
+entries = {}
+for m, n in [(1,1), (1,3), (3,3)]:
+    entries[(m,n)] = matrix_entry(m,n)
+
+c1 = arb(1)
+c3 = arb(fmpq(-1,3))
+arch = c1*c1*entries[(1,1)][0] + 2*c1*c3*entries[(1,3)][0] + c3*c3*entries[(3,3)][0]
+prime = c1*c1*entries[(1,1)][1] + 2*c1*c3*entries[(1,3)][1] + c3*c3*entries[(3,3)][1]
+truncated = arch + prime
+
+# The omitted Lerch contribution is nonnegative.  By Fourier/Plancherel,
+# each omitted q contributes at most ||F'||_2^2/(4q^3), and
+# ||F'||_2^2 = pi^2/L.  For q_j=j+1/4 and j >= M,
+# sum q_j^-3 <= q0^-3 + integral_M^infty (x+1/4)^-3 dx.
+M = N_LERCH
+q0 = arb(fmpq(4*M+1,4))
+tail_upper = (pi*pi/(4*L)) * (1/(q0*q0*q0) + 1/(2*q0*q0))
+certified_upper = truncated.upper() + tail_upper.upper()
+
+print('Q11 arch, prime:', entries[(1,1)][0], entries[(1,1)][1])
+print('Q13 arch, prime:', entries[(1,3)][0], entries[(1,3)][1])
+print('Q33 arch, prime:', entries[(3,3)][0], entries[(3,3)][1])
+print('candidate arch    =', arch)
+print('candidate prime   =', prime)
+print('truncated total   =', truncated)
+print('tail upper bound  =', tail_upper)
+print('certified upper   =', certified_upper)
+print('prime-power count =', sum(prime_power_base(n) is not None for n in range(2,149)))
+
+if certified_upper < 0:
+    print('CERTIFIED_NEGATIVE')
+else:
+    raise SystemExit('No certified negative sign at this precision.')