"""
Verification that the CHSH correlator for sign(phi + alpha) observables
respects |S| <= 2 when computed directly from the bivariate Gaussian
distribution of the vacuum state.

Reference: Azevedo et al., Eur. Phys. J. C 86, 230 (2026)
           Eq. (56)-(59)

Method: Monte Carlo sampling from the joint Gaussian distribution
of the smeared field operators on the vacuum, followed by direct
evaluation of sign(phi + alpha) for each realization.

Result: |S| = 2.000 (within statistical error), contradicting
the reported value |S| = 2.02 in eq. (60) of the paper.
"""

import numpy as np
from scipy.special import erf

# ============================================================
# Parameters from eq. (59) of the paper
# ============================================================
eta       = 0.00132432
eta_prime = 1.04948
sigma     = 0.0737113
sigma_prime = 4.8949
lam       = 0.811798
alpha     = -12.6591
alpha_prime = -8.04613
beta      = -6.00724
beta_prime = -6.05653

# ============================================================
# Covariance structure from eq. (48)
# ============================================================
# Variances: ||h||^2 = (norm parameter)^2 * (1 + lambda^2)
var_f  = eta**2       * (1 + lam**2)
var_fp = eta_prime**2 * (1 + lam**2)
var_g  = sigma**2     * (1 + lam**2)
var_gp = sigma_prime**2 * (1 + lam**2)

# Covariances: <f|g> = 2 * eta * sigma * lambda (real, since
# f and g=jf are spacelike separated, so PJ(f,g) = 0)
cov_fg   = 2 * eta * sigma * lam
cov_fpgp = 2 * eta_prime * sigma_prime * lam
cov_fpg  = 0.0   # from <f|jf'> = 0
cov_fgp  = 0.0   # from <f|jf'> = 0

# ============================================================
# Method 1: Monte Carlo (model-independent verification)
# ============================================================
def correlator_mc(a, b, v1, v2, cov12, N=200_000_000):
    """
    Compute E[sign(X+a) * sign(Y+b)] where (X,Y) ~ N(0, Sigma)
    with Sigma = [[v1, cov12], [cov12, v2]]
    """
    s1, s2 = np.sqrt(v1), np.sqrt(v2)
    rho = cov12 / (s1 * s2) if (s1 > 0 and s2 > 0) else 0.0
    
    Z1 = np.random.randn(N)
    Z2 = rho * Z1 + np.sqrt(1 - rho**2) * np.random.randn(N)
    
    X = s1 * Z1
    Y = s2 * Z2
    
    return np.mean(np.sign(X + a) * np.sign(Y + b))

# ============================================================
# Method 2: Exact analytical (bivariate Gaussian + erf)
# ============================================================
def correlator_exact(a, b, v1, v2, cov12):
    """
    For the case cov12 = 0, the correlator factors:
    E[sign(X+a)*sign(Y+b)] = erf(a/sqrt(2*v1)) * erf(b/sqrt(2*v2))
    
    For cov12 != 0 with extreme arguments (|a/sqrt(v1)| >> 1),
    the result approaches sign(a)*sign(b) = +1 or -1.
    """
    s1, s2 = np.sqrt(v1), np.sqrt(v2)
    u = a / (s1 * np.sqrt(2))
    v = b / (s2 * np.sqrt(2))
    
    if cov12 == 0:
        return erf(u) * erf(v)
    else:
        # For extreme arguments, erf -> +/-1
        # General case would need bivariate normal CDF
        # Here we use MC for the correlated cases
        return None  # Use MC for this case

# ============================================================
# Compute all four correlators
# ============================================================
print("=" * 65)
print("CHSH VERIFICATION: sign(phi + offset) on vacuum Gaussian state")
print("=" * 65)
print()
print("Parameters from Azevedo et al., Eur. Phys. J. C 86, 230 (2026)")
print(f"  eta={eta}, eta'={eta_prime}, sigma={sigma}, sigma'={sigma_prime}")
print(f"  lambda={lam}")
print(f"  alpha={alpha}, alpha'={alpha_prime}, beta={beta}, beta'={beta_prime}")
print()

print("Covariance structure (from eq. 48):")
print(f"  Var(phi(f))  = {var_f:.4e}   [std = {np.sqrt(var_f):.4e}]")
print(f"  Var(phi(f')) = {var_fp:.4e}   [std = {np.sqrt(var_fp):.4f}]")
print(f"  Var(phi(g))  = {var_g:.4e}   [std = {np.sqrt(var_g):.4f}]")
print(f"  Var(phi(g')) = {var_gp:.4e}  [std = {np.sqrt(var_gp):.4f}]")
print(f"  Cov(f,g)     = {cov_fg:.4e}  [rho = {cov_fg/np.sqrt(var_f*var_g):.6f}]")
print(f"  Cov(f',g')   = {cov_fpgp:.4e}  [rho = {cov_fpgp/np.sqrt(var_fp*var_gp):.6f}]")
print(f"  Cov(f',g)    = 0")
print(f"  Cov(f,g')    = 0")
print()

print("Note: offset/std ratios (how many std devs the offset is from zero):")
print(f"  alpha  / std(phi(f))  = {alpha/np.sqrt(var_f):.0f}")
print(f"  alpha' / std(phi(f')) = {alpha_prime/np.sqrt(var_fp):.2f}")
print(f"  beta   / std(phi(g))  = {beta/np.sqrt(var_g):.1f}")
print(f"  beta'  / std(phi(g')) = {beta_prime/np.sqrt(var_gp):.4f}")
print()

# Analytical results for zero-correlation pairs
E_fpg_exact  = erf(alpha_prime/np.sqrt(2*var_fp)) * erf(beta/np.sqrt(2*var_g))
E_fgp_exact  = erf(alpha/np.sqrt(2*var_f)) * erf(beta_prime/np.sqrt(2*var_gp))

print("-" * 65)
print("RESULTS")
print("-" * 65)
print()
print("Analytical (exact for uncorrelated pairs):")
print(f"  E(f',g)  = erf({alpha_prime/np.sqrt(2*var_fp):.3f}) * erf({beta/np.sqrt(2*var_g):.1f})")
print(f"           = {erf(alpha_prime/np.sqrt(2*var_fp)):.6f} * {erf(beta/np.sqrt(2*var_g)):.6f}")
print(f"           = {E_fpg_exact:.6f}")
print(f"  E(f,g')  = erf({alpha/np.sqrt(2*var_f):.1f}) * erf({beta_prime/np.sqrt(2*var_gp):.4f})")
print(f"           = {erf(alpha/np.sqrt(2*var_f)):.6f} * {erf(beta_prime/np.sqrt(2*var_gp)):.6f}")
print(f"           = {E_fgp_exact:.6f}")
print()

# Monte Carlo for all four (including correlated pairs)
np.random.seed(42)
N = 200_000_000
print(f"Monte Carlo (N = {N:,} samples):")

E_fg  = correlator_mc(alpha, beta, var_f, var_g, cov_fg, N)
E_fpg = correlator_mc(alpha_prime, beta, var_fp, var_g, cov_fpg, N)
E_fgp = correlator_mc(alpha, beta_prime, var_f, var_gp, cov_fgp, N)
E_fpgp = correlator_mc(alpha_prime, beta_prime, var_fp, var_gp, cov_fpgp, N)

print(f"  E(f,g)   = {E_fg:.6f}")
print(f"  E(f',g)  = {E_fpg:.6f}")
print(f"  E(f,g')  = {E_fgp:.6f}")
print(f"  E(f',g') = {E_fpgp:.6f}")
print()

# Statistical error estimate
err = 1.0 / np.sqrt(N)
print(f"  Statistical error: ~{err:.6f}")
print()

# CHSH
S = E_fg + E_fpg + E_fgp - E_fpgp
print("=" * 65)
print(f"CHSH = E(f,g) + E(f',g) + E(f,g') - E(f',g')")
print(f"     = {E_fg:.6f} + {E_fpg:.6f} + {E_fgp:.6f} - {E_fpgp:.6f}")
print(f"     = {S:.6f}")
print(f"|S|  = {abs(S):.6f}")
print()
print(f"Classical bound:   |S| <= 2")
print(f"Paper's claim:     |S| = 2.02034")
print(f"This computation:  |S| = {abs(S):.4f} (= 2 within statistical error)")
print()

if abs(S) <= 2.0 + 3*err:
    print("CONCLUSION: The CHSH bound is RESPECTED.")
    print("The value 2.02 in the paper is a numerical artifact.")
else:
    print("CONCLUSION: VIOLATION detected (unexpected).")

print()
print("=" * 65)
print("WHY the paper's integral (eq. 58) gives a wrong result:")
print("=" * 65)
print()
print("The integrand in eq. (58) for the E(f,g) term has:")
print(f"  - Gaussian envelope width in k: ~1/sqrt(var_f) = {1/np.sqrt(var_f):.0f}")
print(f"  - cos(alpha*k) oscillation period: 2*pi/|alpha| = {2*np.pi/abs(alpha):.3f}")
print(f"  - Number of oscillations within envelope: ~{1/np.sqrt(var_f)*abs(alpha)/(2*np.pi):.0f}")
print()
print("With ~1200 oscillations in one integration variable, the")
print("QuasiMonteCarlo method (even with 10^9 points) cannot resolve")
print("the integral accurately. The 0.02 excess is numerical noise.")