2607.20408v1 / Repro/parity_check.py

all files

#!/usr/bin/env python3
"""Numerically check the parity defect in Lemma 2.1 using only Python's stdlib.

For d = -4 the character is chi_4 and the exact approximate functional equation
uses Q(3/4, pi (n/2)^2), while the paper prints Q(1/4, pi (n/2)^2).
Here Q is the regularized upper incomplete gamma function.
"""

import math

EPS = 3e-15
FPMIN = 1e-300
ITMAX = 10_000


def gammaincc(a: float, x: float) -> float:
    """Regularized upper incomplete gamma Q(a,x), Numerical Recipes algorithm."""
    gln = math.lgamma(a)
    if x == 0:
        return 1.0
    if x < a + 1:
        ap = a
        total = term = 1.0 / a
        for _ in range(ITMAX):
            ap += 1
            term *= x / ap
            total += term
            if abs(term) < abs(total) * EPS:
                break
        return 1 - total * math.exp(-x + a * math.log(x) - gln)
    b = x + 1 - a
    c = 1 / FPMIN
    d = 1 / b
    h = d
    for i in range(1, ITMAX + 1):
        an = -i * (i - a)
        b += 2
        d = an * d + b
        if abs(d) < FPMIN:
            d = FPMIN
        c = b + an / c
        if abs(c) < FPMIN:
            c = FPMIN
        d = 1 / d
        delta = d * c
        h *= delta
        if abs(delta - 1) < EPS:
            break
    return math.exp(-x + a * math.log(x) - gln) * h


def chi4(n: int) -> int:
    if n % 2 == 0:
        return 0
    return 1 if n % 4 == 1 else -1


def afe_rhs(gamma_parameter: float, terms: int = 100) -> float:
    return 2 * sum(
        chi4(n) / math.sqrt(n)
        * gammaincc(gamma_parameter, math.pi * (n / 2) ** 2)
        for n in range(1, terms + 1)
    )


if __name__ == "__main__":
    printed = afe_rhs(0.25)
    correct = afe_rhs(0.75)
    print(f"paper-even-weight RHS: {printed:.17g}")
    print(f"correct odd-weight RHS: {correct:.17g}")
    print(f"difference:             {printed - correct:.17g}")