# -*- coding: utf-8 -*-
# Numerical value of the deficit c1(xi2) of the Harman-sieve minorant, with the formula of
# J. Stadlmann, arXiv:2608.31126, Proposition 2 (sum of two four-dimensional integrals I1 + I2).
# Method 1: nested Gauss-Legendre quadrature, each range split where its piecewise-linear limits change.
# Method 2 (independent check): Monte Carlo, with the regions written directly from the summation
#           conditions of the minorant, not from the integration limits.
# Not a certified computation. Requires numpy. Run: python3 deficit_c1.py
# Auteur : David Niedbala Giraudin (Fino), ORCID 0009-0009-1526-1178, avec assistance IA (Claude, Anthropic).
# Deficit c1(xi2), formula of Stadlmann arXiv:2608.31126 Prop. 2, by nested Gauss-Legendre with
# splitting of each range at the points where the piecewise-linear limits change formula.
import sys, numpy as np
from numpy.polynomial.legendre import leggauss
def gl(n):
    t, w = leggauss(n); return (t + 1)/2, w/2
def nested(f, limits, n, depth=0, outer=()):
    """limits[k](outer...) -> list of breakpoints (sorted, >=2 entries) for variable k given outer vars"""
    bps = limits[depth](*outer)
    t, w = gl(n); tot = 0.0
    for a, b in zip(bps[:-1], bps[1:]):
        if b <= a: continue
        xs = a + (b - a)*t
        if depth == len(limits) - 1:
            tot += (b - a)*np.dot(w, f(*outer, xs))
        else:
            tot += (b - a)*sum(wi*nested(f, limits, n, depth + 1, outer + (xi,)) for xi, wi in zip(xs, w))
    return tot
def clip(l, u, pts):
    pts = sorted(set([l, u] + [p for p in pts if l < p < u]))
    return pts if u > l else [l, l]
def c1(x2, n):
    lo, up = 1 - 2*x2, 8*x2 - 3
    # I1: a1, a2, a3, a4
    L1 = [lambda: clip(lo, x2, [x2/2, x2 - lo]),
          lambda a1: clip(lo, max(lo, min(a1, x2 - a1)), []),
          lambda a1, a2: clip(lo, max(lo, a2), [1 - x2 - a2 - lo, (1 - a1 - a2 - lo)/2, 1 - x2 - a2 - a2]),
          lambda a1, a2, a3: (lambda l: clip(l, max(l, min(a3, (1 - a1 - a2 - a3)/2)), []))(max(lo, 1 - x2 - a2 - a3))]
    f1 = lambda a1, a2, a3, a4: 1/(a1*a2*a3*a4*(1 - a1 - a2 - a3 - a4))
    I1 = nested(f1, L1, n)
    # I2: a2, a3, a4, a5
    L2 = [lambda: clip(lo, max(lo, up), [x2 - up, x2 - lo]),
          lambda a2: clip(lo, max(lo, min(up, a2)), [x2 - a2, 1 - x2 - up - (x2 - a2)]),
          lambda a2, a3: (lambda l: clip(l, max(l, min(up, x2 - a2)), [1 - x2 - a3 - up, 1 - x2 - a3 - lo]))(max(a3, lo)),
          lambda a2, a3, a4: (lambda l: clip(l, max(l, min(up, (1 - a2 - a3 - a4)/2)), []))(max(lo, 1 - x2 - a3 - a4))]
    f2 = lambda a2, a3, a4, a5: 1/(a2*a3*a4*a5*(1 - a2 - a3 - a4 - a5))
    I2 = nested(f2, L2, n)
    return I1, I2

def montecarlo(x2, N, seed=12345):
    lo, up = 1 - 2*x2, 8*x2 - 3; rng = np.random.default_rng(seed); B = 4_000_000
    # I1 in the box a1 in [lo, x2-lo], a2, a3, a4 in [lo, x2/2] (which contains its region)
    L = np.array([lo, lo, lo, lo]); U = np.array([x2 - lo, x2/2, x2/2, x2/2]); s = q = 0.0; t = 0
    while t < N:
        a = rng.uniform(L, U, size=(B, 4)); a1, a2, a3, a4 = a.T; a5 = 1 - a.sum(1)
        m = (a4 < a3) & (a3 < a2) & (a2 < a1) & (a1 + a2 < x2) & (a2 + a3 + a4 > 1 - x2) & (a5 > a4)
        v = np.where(m, 1/(a1*a2*a3*a4*np.where(m, a5, 1)), 0.0); s += v.sum(); q += (v*v).sum(); t += B
    V = np.prod(U - L); I1 = V*s/t; e1 = V*np.sqrt(max(q/t - (s/t)**2, 0)/t)
    # I2 in the box [lo, up]^4 (variables a2, a3, a4, a5; a6 = 1 - sum)
    s = q = 0.0; t = 0
    while t < N:
        b = rng.uniform(lo, up, size=(B, 4)); b2, b3, b4, b5 = b.T; b6 = 1 - b.sum(1)
        m = (b2 > b3) & (b4 > b3) & (b2 + b4 < x2) & (b3 + b4 + b5 > 1 - x2) & (b6 > b5) & (b6 < up)
        v = np.where(m, 1/(b2*b3*b4*b5*np.where(m, b6, 1)), 0.0); s += v.sum(); q += (v*v).sum(); t += B
    V = (up - lo)**4; I2 = V*s/t; e2 = V*np.sqrt(max(q/t - (s/t)**2, 0)/t)
    return I1, e1, I2, e2

if __name__ == "__main__":
    from fractions import Fraction as Fr
    for name, x2 in (("55/136", float(Fr(55, 136))), ("0.40481 (Stadlmann 2023, published bound c1 <= 2e-5)", 0.40481)):
        print("xi2 =", name)
        for n in (16, 24, 32):
            I1, I2 = c1(x2, n); print(f"   Gauss-Legendre n={n}: I1={I1:.6e}  I2={I2:.6e}  c1={I1 + I2:.6e}", flush=True)
        I1, e1, I2, e2 = montecarlo(x2, 2e8)
        print(f"   Monte Carlo 2e8 points: I1={I1:.4e} +- {e1:.1e}  I2={I2:.4e} +- {e2:.1e}  c1={I1 + I2:.4e} +- {np.hypot(e1, e2):.1e}", flush=True)