# pft-riemann-demo.py — minimal reproducible demo
# This notice must remain-2025-10-02
# Author: James Johan Sebastian Allen
# PatternFieldTheory.com
import numpy as np, math
import matplotlib.pyplot as plt
from numpy.linalg import eigvalsh

def hex_nodes(R):
    nodes=[(0,0,1)]; idx=2
    for r in range(1,R+1):
        for k in range(6*r):
            nodes.append((r,k,idx)); idx+=1
    return nodes

def neighbor_map(R):
    nodes = hex_nodes(R)
    coord2idx = {(r,k):gi for (r,k,gi) in nodes}
    def wrap(k,m): return k % m if m>0 else None
    nbrs = {gi:set() for (_,_,gi) in nodes}
    if R>=1:
        g0 = coord2idx[(0,0)]
        for k in range(6):
            g1 = coord2idx[(1,k)]
            nbrs[g0].add(g1); nbrs[g1].add(g0)
    for r in range(1,R+1):
        m=6*r
        for k in range(m):
            gi = coord2idx[(r,k)]
            for dk in (-1,1):
                gj = coord2idx[(r, wrap(k+dk,m))]
                nbrs[gi].add(gj); nbrs[gj].add(gi)
            mi = 6*(r-1)
            if mi==0:
                gj = coord2idx[(0,0)]
                nbrs[gi].add(gj); nbrs[gj].add(gi)
            else:
                a=(k*mi)//m
                for off in (0,1):
                    gj = coord2idx[(r-1, wrap(a+off,mi))]
                    nbrs[gi].add(gj); nbrs[gj].add(gi)
            if r+1<=R:
                mo=6*(r+1)
                a=(k*mo)//m
                for off in (0,1):
                    gj = coord2idx[(r+1, wrap(a+off,mo))]
                    nbrs[gi].add(gj); nbrs[gj].add(gi)
    return nodes, nbrs

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

def lam(n):
    if n<2: return 0.0
    # simplified Mangoldt: log p if n is prime else 0
    # (good enough for demo; replace with true prime-power test if desired)
    for p in range(2, int(n**0.5)+1):
        if n%p==0:
            return 0.0
    return math.log(n)

def phase_rule(i, j, alpha=1/1.618, use_lambda=True):
    Phi = lam if use_lambda else (lambda n: 1.0 if n>1 else 0.0)
    theta = 2*math.pi*alpha*(Phi(i)+Phi(j))
    return theta

def build_operator(R=12, t=1.0, V0=2.5, alpha=1/1.618, use_lambda=True):
    nodes, nbrs = neighbor_map(R)
    max_idx = max(gi for _,_,gi in nodes)
    is_prime = sieve_primes(max_idx)
    N = max_idx+1
    H = np.zeros((N,N), dtype=complex)
    for _,_,gi in nodes:
        H[gi,gi] = V0 if (gi<=max_idx and is_prime[gi]) else 0.0
    for _,_,i in nodes:
        for j in nbrs[i]:
            if i<j:
                th = phase_rule(i,j,alpha=alpha,use_lambda=use_lambda)
                val = t*(math.cos(th)+1j*math.sin(th))
                H[i,j]=val; H[j,i]=np.conj(val)
    return H

if __name__ == "__main__":
    R=12
    H = build_operator(R=R, V0=2.5, alpha=1/1.618, use_lambda=True)
    nodes,_ = neighbor_map(R)
    idxs = [gi for _,_,gi in nodes]
    Hs = H[np.ix_(idxs,idxs)]
    E = np.sort(eigvalsh(Hs))
    n=len(E); lo=int(0.2*n); hi=int(0.8*n)
    Eb = E[lo:hi]
    s = np.diff(Eb); s = s/np.mean(s)
    x = np.linspace(0,3,400)
    P = (32/np.pi**2)*(x**2)*np.exp(-4*x**2/np.pi)
    import matplotlib.pyplot as plt
    plt.hist(s, bins=60, density=True, alpha=0.6, label="AOL operator (demo)")
    plt.plot(x, P, lw=2, label="GUE Wigner surmise")
    plt.xlabel("Unfolded spacing s"); plt.ylabel("Density")
    plt.legend(); plt.tight_layout(); plt.show()