#!/usr/bin/env python3
"""
Structural identity checks for CSM_RH Paper 05.

The model coefficient m_h below is arbitrary. The identities checked here are
algebraic and do not depend on the Hardy-Littlewood singular series.

This is not evidence for RH.
"""
import numpy as np
import csv

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

def run(N, seed=0):
    rng = np.random.default_rng(seed + N)
    a = np.zeros(2*N)
    a[1:] = rng.normal(size=2*N-1)
    model = np.zeros(2*N-1)
    model[1:] = rng.normal(scale=0.2, size=2*N-2)

    A = np.cumsum(a)
    J = float(np.sum(A[N:2*N]**2))
    D = sum(wN(N,n)*a[n]**2 for n in range(1,2*N))

    C = {}
    W = {}
    R = {}
    for h in range(1,2*N-1):
        C[h] = sum(wN(N,n)*a[n]*a[n-h] for n in range(h+1,2*N))
        W[h] = sum(wN(N,n) for n in range(h+1,2*N))
        R[h] = C[h] - model[h]*W[h]

    M = sum(model[h]*W[h] for h in R)
    cssa = sum(R.values())

    # Paraproduct
    para = sum(wN(N,n)*a[n]*A[n-1] for n in range(2,2*N))

    # Endpoint reconstruction
    Bsum = 0.0
    Bend_resid = 0.0
    for j in range(N,2*N):
        Bj = 0.0
        for h in range(1,j):
            q = sum(a[n]*a[n-h] for n in range(h+1,j+1)) - model[h]*(j-h)
            Bj += q
        Lj = float(np.sum(a[1:j+1]**2))
        Mj = sum(model[h]*(j-h) for h in range(1,j))
        target = 0.5*A[j]**2 - 0.5*Lj - Mj
        Bend_resid = max(Bend_resid, abs(Bj-target))
        Bsum += Bj

    # Fourier zero-frequency identity is coefficient sum.
    H0 = sum(R.values())
    L2_coeff = np.sqrt(sum(abs(v)**2 for v in R.values()))
    Mcount = len(R)
    eval_ratio = abs(H0) / (np.sqrt(Mcount)*L2_coeff) if L2_coeff else 0.0

    return [
        N,
        J - (D + 2*sum(C.values())),
        para - 0.5*(J-D),
        cssa - (0.5*(J-D)-M),
        Bsum - cssa,
        Bend_resid,
        H0 - cssa,
        eval_ratio
    ]

rows = [run(N) for N in [16,24,40,64]]
with open("campaign04_zero_frequency_structural_crosscheck.csv","w",encoding="utf-8",newline="") as f:
    w = csv.writer(f)
    w.writerow([
        "N",
        "J_additive_residual",
        "paraproduct_energy_residual",
        "cssa_energy_model_residual",
        "endpoint_sum_residual",
        "max_endpoint_identity_residual",
        "zero_frequency_residual",
        "generic_eval_ratio_le_1"
    ])
    w.writerows(rows)

max_identity = max(max(abs(x) for x in r[1:7]) for r in rows)
max_ratio = max(r[7] for r in rows)
print("PASS")
print("max_identity_residual", max_identity)
print("max_generic_eval_ratio", max_ratio)