import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import pearsonr

# Load zero heights
x = np.loadtxt("zeros4.txt")

# Spacings
d = np.diff(x)
d = d - d.mean()  # center

primes = [2,3,5,7,11,13,17,19,23,29,31,37]
vals = []

for p in primes:
    a = d[:-p]
    b = d[p:]
    c = np.mean(a*b)  # covariance at lag p
    vals.append(c)

vals = np.array(vals)
target = (np.log(primes)**2) / np.array(primes)

# BK amplitude law test
r, pv = pearsonr(vals, target)
print("correlation =", r)
print("p-value =", pv)

# Scatter plot
plt.scatter(target, vals)
plt.xlabel(r'$(\log p)^2/p$')
plt.ylabel('Amplitude')
plt.title('Prime-Locked Excess')
plt.savefig("scatter.png", dpi=200)

"""Тэгүүдийн power spectrum S(tau) тооцоолох"""
import numpy as np
import matplotlib.pyplot as plt

def power_spectrum(tau, gamma):
    """S(tau) = |Σ e^{i τ γ_n}|² / N"""
    N = len(gamma)
    Z = np.sum(np.exp(1j * tau * gamma))
    return np.abs(Z)**2 / N

# Prime frequencies τ_p = log(p) / 2π
primes = np.array([2,3,5,7,11,13,17,19,23,29,31,37])
tau_p = np.log(primes) / (2*np.pi)

# Full spectrum
tau_all = np.linspace(0, 0.8, 1000)
S_all = np.array([power_spectrum(t, gamma) for t in tau_all])
S_obs = np.array([power_spectrum(t, gamma) for t in tau_p])

## Hilbert space
H = L²((0,∞), dx/x)

## Mellin transform
(Mf)(τ) = (1/√(2π)) ∫₀^∞ f(x) x^{-iτ} dx/x

## Prime operator
(P_σ f)(x) = Σ_{n≥1} Λ(n)/n^σ · f(x/n)

## Multiplier theorem  ← KEY RESULT
M[P_σ f](τ) = -ζ'(σ+iτ)/ζ(σ+iτ) · Mf(τ)

## Energy identity
E(σ,ε,X) = ∫_R |ζ'/ζ(σ+iτ)|² |φ̂(τ)|² dτ

## ❌ WRONG normalization (low-T)
τ_p = log(p) / (2π)
# → r ≈ 0.4–0.6, BK not confirmed

## ✅ CORRECT normalization (high-T unfolded)
τ_p = log(p) / log(T / 2π)
# → r = 0.014

## Why this matters
# At height T, zeros have mean spacing ~2π/log(T/2π)
# The "unfolded" coordinate accounts for this
# Montgomery's pair correlation uses this normalization