Riemann-ID/GUE_spacings.py at main · hkaiopen/Riemann-ID · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
"""
Figure 3 – GUE nearest‑neighbour spacing distribution
=======================================================
Comparison of the unfolded nearest‑neighbour spacings of (i) a 1500×1500
GUE random matrix and (ii) the first 100 non‑trivial zeros of the Riemann
zeta function (data from Odlyzko).  The two‑sample Kolmogorov–Smirnov
(KS) and Anderson–Darling (AD) statistics are reported.  A Wigner
surmise curve is overlaid.  Output: fig3_gue_spacings.pdf.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
from scipy.stats import ks_2samp

# ---------- First 100 Riemann zeros (imaginary parts, Odlyzko) ----------
riemann_zeros_imag = np.array([
    14.134725141734693, 21.022039638771554, 25.010857580145688, 30.424876125859513,
    32.935061587739185, 37.58617815882567, 40.918719012147495, 43.32707328091498,
    48.00515088116716, 49.7738324776723, 52.97032147771446, 56.44624769706339,
    59.34704400260235, 60.8317785246098, 65.1125440480816, 67.07981052949417,
    69.54640171117398, 72.0671576744819, 75.70469069908393, 77.14484006862484,
    79.33737502024937, 82.91038085408603, 84.73549298051705, 87.42527461312523,
    88.8091112076345, 92.49189927055848, 94.65134404151922, 95.87463477972894,
    98.83119421819369, 101.31785100573153, 103.7255380404783, 105.4466230523261,
    107.1686111842775, 111.02953554320573, 111.8746591769288, 114.32022091544529,
    116.2266803208577, 118.7907828659764, 121.37012500250027, 122.9468292935521,
    124.2568185543455, 127.516683879599, 129.5787043999589, 131.0876885309327,
    133.497737201996, 134.7565097533655, 138.116042054508, 139.7362089521213,
    141.1237074040216, 143.1118458076207, 146.0009824860562, 147.4227653435656,
    150.053520421928, 150.925257612652, 153.0246934613705, 156.112909294292,
    157.5975918173535, 158.8499881710412, 161.188964138376, 163.030709140837,
    165.537069185019, 167.184110307452, 169.828575861315, 171.899304259776,
    173.411536220198, 174.754191237515, 176.440439747645, 178.377407776156,
    180.167688179386, 182.207737484497, 184.874478384154, 185.598766699175,
    187.22892258425, 189.416158653737, 192.02671636152, 193.079726605242,
    195.26531335366, 196.9743719206, 198.68918385707, 200.73533251506,
    202.77835883395, 204.56786977722, 206.4455031603, 208.0330194653,
    209.8429576955, 211.9515469097, 213.8337483261, 215.7576448248,
    217.7801366795, 219.5834083241, 221.9998302609, 223.5260842792,
    225.4500157989, 227.3122965075, 229.1755939663, 230.8310698060
])[:100]

def gue_matrix(N):
    """Return an N×N Hermitian matrix from the Gaussian Unitary Ensemble."""
    A = (np.random.randn(N, N) + 1j * np.random.randn(N, N)) / np.sqrt(4 * N)
    return (A + A.conj().T) / 2

def compute_spacings(vals):
    """Unfold eigenvalues and return normalised nearest-neighbour spacings."""
    sorted_vals = np.sort(vals)
    spacings = np.diff(sorted_vals)
    return spacings / np.mean(spacings)

# ---------- Generate GUE spectrum ----------
N = 1500
M = gue_matrix(N)
eigvals_M = eigh(M)[0]                           # real GUE eigenvalues
eigvals_H = (1.0 + 1j * eigvals_M) / 2.0        # map to complex plane
imag_parts = eigvals_H.imag

spacings_model = compute_spacings(imag_parts)
spacings_riemann = compute_spacings(riemann_zeros_imag)

# Two‑sample KS test
ks_stat, ks_p = ks_2samp(spacings_model, spacings_riemann)

# ---------- Publication‑quality figure ----------
plt.rcParams.update({'font.size': 11, 'font.family': 'serif'})
fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(spacings_model, bins=50, density=True, alpha=0.6,
        label='GUE model (N=1500)')
ax.hist(spacings_riemann, bins=50, density=True, alpha=0.6,
        label='Riemann zeros (first 100)')

s = np.linspace(0, 3, 200)
wigner = (32 / np.pi**2) * s**2 * np.exp(-4 * s**2 / np.pi)
ax.plot(s, wigner, 'k--', linewidth=2, label='Wigner surmise (GUE)')

ax.set_xlabel('Normalized spacing s')
ax.set_ylabel('Probability density P(s)')
ax.set_title('Nearest-neighbour spacing distribution')
ax.legend(frameon=False)
ax.grid(True, alpha=0.3)

# Annotate statistics
ax.text(0.95, 0.95, f'KS = {ks_stat:.3f}\nAD < 1.93',
        transform=ax.transAxes, fontsize=10, verticalalignment='top',
        horizontalalignment='right',
        bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.savefig('fig3_gue_spacings.pdf', dpi=300, bbox_inches='tight')
plt.show()