#!/usr/bin/env python3
"""
Finite identity checks for CSM_RH Paper 04.
This validates the exact centered-shift reconstruction on finite data.
It is NOT evidence for RH.
"""
import math
import csv
import numpy as np

def sieve_primes(limit):
    mark = np.ones(limit + 1, dtype=bool)
    mark[:2] = False
    for p in range(2, int(limit**0.5) + 1):
        if mark[p]:
            mark[p*p:limit+1:p] = False
    return np.flatnonzero(mark).tolist()

def von_mangoldt(limit):
    lam = np.zeros(limit + 1, dtype=float)
    for p in sieve_primes(limit):
        q = p
        lp = math.log(p)
        while q <= limit:
            lam[q] = lp
            if q > limit // p:
                break
            q *= p
    return lam

def twin_constant(prime_limit=50000):
    prod = 1.0
    for p in sieve_primes(prime_limit):
        if p > 2:
            prod *= 1.0 - 1.0/((p-1.0)**2)
    return prod

def singular_series_pair(h, C2):
    if h % 2:
        return 0.0
    n = h
    out = 2.0*C2
    p = 3
    while p*p <= n:
        if n % p == 0:
            out *= (p-1.0)/(p-2.0)
            while n % p == 0:
                n //= p
        p += 2
    if n > 2:
        out *= (n-1.0)/(n-2.0)
    return out

def endpoint_weight(N, n):
    if 1 <= n <= N:
        return float(N)
    if N < n < 2*N:
        return float(2*N-n)
    return 0.0

def W_mass(N, h):
    if h <= N:
        return N*(N-h) + N*(N-1)/2.0
    m = 2*N-h-1
    return m*(m+1)/2.0

def run(N, C2):
    limit = 2*N-1
    lam = von_mangoldt(limit)
    a = lam - 1.0
    a[0] = 0.0
    A = np.cumsum(a)
    js = np.arange(N, 2*N)
    J = float(np.sum(A[js]**2))
    linear = float(np.sum(A[js]))
    I_exact = J - linear + N/3.0

    w = np.array([endpoint_weight(N,n) for n in range(limit+1)])
    D = float(np.sum(w[1:]*a[1:]**2))

    M = 0.0
    Rsum = 0.0
    for h in range(1, 2*N-1):
        n = np.arange(h+1, 2*N)
        C = float(np.sum(w[n]*a[n]*a[n-h]))
        W = W_mass(N,h)
        S = singular_series_pair(h,C2)
        M += (S-1.0)*W
        Rsum += C - (S-1.0)*W

    reconstructed = D + 2*M + 2*Rsum
    signed_identity_residual = reconstructed - J
    halfJ_residual = Rsum - 0.5*J
    deterministic_correction = -0.5*(D + 2*M)
    halfJ_formula_residual = halfJ_residual - deterministic_correction

    return [
        N, J, I_exact, D, M, Rsum,
        signed_identity_residual,
        halfJ_formula_residual
    ]

if __name__ == "__main__":
    C2 = twin_constant()
    rows = [run(N,C2) for N in [80,120,200,320]]
    with open("campaign03_signed_gate_crosscheck.csv","w",encoding="utf-8",newline="") as f:
        w = csv.writer(f)
        w.writerow([
            "N","J_N","I_exact","D_N","M_N","A_N",
            "J_reconstruction_residual",
            "A_halfJ_correction_residual"
        ])
        w.writerows(rows)
    maxerr = max(max(abs(r[-1]), abs(r[-2])) for r in rows)
    print("PASS")
    print("max_identity_residual", maxerr)