diff --git a/.github/workflows/rh-verifier.yml b/.github/workflows/rh-verifier.yml new file mode 100644 index 0000000..c23f7dd --- /dev/null +++ b/.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 diff --git a/rh_verifier/verify_candidate.py b/rh_verifier/verify_candidate.py new file mode 100644 index 0000000..129d262 --- /dev/null +++ b/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.')