# Python 3.x + numpy
import numpy as np, math

def is_prime(n):
    if n < 2: return False
    if n % 2 == 0: return n == 2
    r = int(n**0.5); f = 3
    while f <= r:
        if n % f == 0: return False
        f += 2
    return True

def residues_mod_n(n):
    R=set()
    for x in range(1,n):
        r=(x*x) % n
        if r != 0: R.add(r)
    return R

def residue_first_row(n):
    R = residues_mod_n(n)
    a = np.zeros(n, float)
    for k in range(1, n):
        if (k in R) or ((n-k) in R): a[k] = 1.0
    return a

def lambda2_residue_fft(n):
    # Laplacian L = D - A for the residue-circulant
    a = residue_first_row(n)
    d = float(a.sum())
    eig_adj = np.fft.fft(a).real
    mu = d - eig_adj    # Laplacian eigenvalues of circulant
    mu.sort()
    # lambda_2 := first positive Laplacian eigenvalue
    for v in mu:
        if v > 1e-12: return float(v)
    return float(mu[1])

def paley_formula(n):
    return 0.5*(n - math.sqrt(n))

def paley_gap(n):  # only meaningful for n == 1 (mod 4)
    if n % 4 != 1: return float('inf')
    return abs(lambda2_residue_fft(n) - paley_formula(n))

def verify(limit=2601, eps=1e-12):
    primes14 = [m for m in range(5, limit+1, 4) if is_prime(m)]
    zeros = [m for m in range(5, limit+1, 4) if paley_gap(m) <= eps]
    extra = sorted(set(zeros) - set(primes14))
    miss  = sorted(set(primes14) - set(zeros))
    print(f"tested up to m={limit}: zeros={len(zeros)}, primes(1 mod 4)={len(primes14)}")
    print("composites flagged as zero:", extra)
    print("primes(1 mod 4) missed:", miss)

if __name__ == "__main__":
    verify(2601, 1e-12)