#!/usr/bin/env python3
"""
test_352 -- CHARACTERISATIONS OF THE ROUND SPHERE: PARALLEL HYPERSURFACES, STEINER'S FORMULA, ALEXANDROV-FENCHEL IN
            SCHNEIDER'S NORMALISATION, MINKOWSKI, THE CUBIC BOUND, AND THE CORRUGATED CHAIN

For MB §8 (the companion of MA). Conventions as in MA. W_k = int sigma_k(K) dV for a closed hypersurface of E^4
(k = 0..3), so the sphere of radius a has W_k = C(3,k) 2 pi^2 a^{3-k}. Schneider's quermassintegrals V_j are defined by
Vol(K + delta B) = sum_{j=0}^{4} C(4,j) V_j delta^j.

PRE-REGISTERED:
  P1 (parallel hypersurfaces) For x_a = x + a N the principal curvatures are k_i/(1 + a k_i), so
     H(a) = sum k_i/(1 + a k_i) has H(0) = p1, H'(0) = -p2, H''(0) = 2 p3 (p_j = sum k_i^j), and Newton's identities
     recover sigma_1, sigma_2, sigma_3 from p1, p2, p3 (symbolic). So constant mean curvature on every parallel
     hypersurface near x makes the principal curvatures constant.
  P2 (Steiner) sigma_k(k_i/(1 + d k_i)) prod(1 + d k_i) = sum_{j>=k} C(j,k) d^{j-k} sigma_j(k) for k = 0..3 (symbolic), so
     dW_k/dd = (k+1) W_{k+1} along any parallel family, and the round values satisfy it.
  P3 (normalisation and the capsule) Vol(K + d B) = Vol + W0 d + W1 d^2/2 + W2 d^3/3 + W3 d^4/4, so W0 = 4V1, W1 = 12V2,
     W2 = 12V3, W3 = 4V4 (on the ball, symbolically); Alexandrov-Fenchel's V3^2 >= V2 V4 is W2^2 >= 3 W1 W3. For the
     capsule B_r + [-t, t]v, Vol(d) = pi^2 (r+d)^4/2 + (8 pi t/3)(r+d)^3 gives W0 = 2pi^2 r^3 + 8pi t r^2,
     W1 = 6pi^2 r^2 + 16 pi t r, W2 = 6 pi^2 r + 8 pi t, W3 = 2 pi^2, and W2^2 - 3 W1 W3 = 64 pi^2 t^2 exactly.
  P4 (Minkowski) W1 >= 3 W3^{1/3} W0^{2/3}: equality on the ball (symbolic), strict on the capsule for t/r in
     {1/10, 1, 10}.
  P5 (the cubic bound) sum k_i^2 >= 3 |k1 k2 k3|^{2/3} (AM-GM; 10^4 random triples, with equality checked at
     (1,1,1)); hence int |A|^3 >= 3 sqrt3 int |det K| >= 6 sqrt3 pi^2. Numerically on the radial graphs
     r = 1 + eps cos 2chi: int |A|^3 = 6 sqrt3 pi^2 at eps = 0 (to 1e-25), and strictly larger at eps = 1/10 (convex)
     and 3/10 (not convex).
  P6 (the corrugated chain) For the hypersurface of revolution with tangent angle psi = sigma + A sin 2k sigma on
     [0, pi] (rho' = cos psi): int sigma_2 dV = 4 pi int (1 + cos^2 psi) dsigma (by parts, symbolic for arbitrary psi),
     and = 6 pi^2 exactly for integer k >= 2 at every A (to 1e-25 at k = 2, 3 and A = 1/10, 3/10, 1/2, including
     2kA > 1); and W1 = 6 pi^2 - 4 pi^2 k^2 (4k^2 - 3) A^2/(4k^2 - 1)^2 + O(A^3) at k = 2, 3 (exact series).
"""
import sys
import random
import sympy as sp
import mpmath as mp

FAILURES = []


def check(name, ok, detail=""):
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}" + (f"  -- {detail}" if detail else ""))
    if not ok:
        FAILURES.append(name)


k1, k2, k3, a, d = sp.symbols("k1 k2 k3 a d", real=True)
ks = [k1, k2, k3]


def sig(j, v):
    if j == 0:
        return sp.Integer(1)
    from itertools import combinations
    return sum(sp.prod(c) for c in combinations(v, j))


# ------------------------------------------------------------------------------------------------------------- P1
print("P1  parallel hypersurfaces")
Ha = sum(k / (1 + a * k) for k in ks)
p = [sum(k**j for k in ks) for j in (1, 2, 3)]
ok1 = sp.simplify(Ha.subs(a, 0) - p[0]) == 0
ok1 &= sp.simplify(sp.diff(Ha, a).subs(a, 0) + p[1]) == 0
ok1 &= sp.simplify(sp.diff(Ha, a, 2).subs(a, 0) - 2 * p[2]) == 0
e1 = p[0]
e2 = (p[0] * e1 - p[1]) / 2
e3 = (p[2] - p[1] * e1 + p[0] * e2) / 3
ok1 &= all(sp.expand(ei - sig(j, ks)) == 0 for j, ei in zip((1, 2, 3), (e1, e2, e3)))
check("H(0) = p1, H'(0) = -p2, H''(0) = 2p3; Newton gives sigma_1..3", ok1)

# ------------------------------------------------------------------------------------------------------------- P2
print("P2  Steiner's formula")
kd = [k / (1 + d * k) for k in ks]
J = sp.prod([1 + d * k for k in ks])
ok2 = True
for k in range(4):
    lhs = sp.simplify(sig(k, kd) * J)
    rhs = sum(sp.binomial(j, k) * d**(j - k) * sig(j, ks) for j in range(k, 4))
    ok2 &= sp.expand(lhs - rhs) == 0
    if k < 3:
        ok2 &= sp.expand(sp.diff(rhs, d).subs(d, 0) - (k + 1) * sig(k + 1, ks)) == 0
ar = sp.Symbol("ar", positive=True)
Wr = [sp.binomial(3, k) * 2 * sp.pi**2 * ar**(3 - k) for k in range(4)]
ok2 &= all(sp.simplify(sp.diff(Wr[k], ar) - (k + 1) * Wr[k + 1]) == 0 for k in range(3))
check("sigma_k of the parallel surface times the Jacobian = sum C(j,k) d^(j-k) sigma_j; dW_k = (k+1) W_{k+1}", ok2)

# ------------------------------------------------------------------------------------------------------------- P3
print("P3  Schneider's normalisation; the capsule")
r, tt = sp.symbols("r t", positive=True)


def Ws_from(vol):
    c = [sp.expand(sp.diff(vol, d, j).subs(d, 0) / sp.factorial(j)) for j in range(5)]
    W = [c[1], 2 * c[2], 3 * c[3], 4 * c[4]]
    V = [c[j] / sp.binomial(4, j) for j in range(5)]
    return W, V


Wb, Vb = Ws_from(sp.pi**2 * (r + d)**4 / 2)
ok3 = all(sp.simplify(Wb[k] - sp.binomial(3, k) * 2 * sp.pi**2 * r**(3 - k)) == 0 for k in range(4))
ok3 &= all(sp.simplify(x1 - x2) == 0 for x1, x2 in
           ((Wb[0], 4 * Vb[1]), (Wb[1], 12 * Vb[2]), (Wb[2], 12 * Vb[3]), (Wb[3], 4 * Vb[4])))
Wc, Vc = Ws_from(sp.pi**2 * (r + d)**4 / 2 + 8 * sp.pi * tt / 3 * (r + d)**3)
want = [2 * sp.pi**2 * r**3 + 8 * sp.pi * tt * r**2, 6 * sp.pi**2 * r**2 + 16 * sp.pi * tt * r,
        6 * sp.pi**2 * r + 8 * sp.pi * tt, 2 * sp.pi**2]
ok3 &= all(sp.expand(Wc[k] - want[k]) == 0 for k in range(4))
ok3 &= sp.expand(Wc[2]**2 - 3 * Wc[1] * Wc[3] - 64 * sp.pi**2 * tt**2) == 0
ok3 &= sp.expand((Vc[3]**2 - Vc[2] * Vc[4]) * 144 - (Wc[2]**2 - 3 * Wc[1] * Wc[3])) == 0
check("W0 = 4V1, W1 = 12V2, W2 = 12V3, W3 = 4V4; capsule W_k; W2^2 - 3W1W3 = 64 pi^2 t^2 = 144(V3^2 - V2V4)", ok3)

# ------------------------------------------------------------------------------------------------------------- P4
print("P4  Minkowski's inequality")
ok4 = sp.simplify(Wb[1] - 3 * Wb[3]**sp.Rational(1, 3) * Wb[0]**sp.Rational(2, 3)) == 0
vals = []
for q in (sp.Rational(1, 10), 1, 10):
    v = sp.N((Wc[1] - 3 * Wc[3]**sp.Rational(1, 3) * Wc[0]**sp.Rational(2, 3)).subs({r: 1, tt: q}), 30)
    ok4 &= v > 0
    vals.append(f"t/r={q}: {sp.N(v, 6)}")
check("W1 = 3 W3^(1/3) W0^(2/3) on the ball, strictly larger on the capsule", ok4, "; ".join(vals))

# ------------------------------------------------------------------------------------------------------------- P5
print("P5  the cubic bound")
random.seed(352)
ok5 = True
for _ in range(10000):
    kk = [random.uniform(-5, 5) for _ in range(3)]
    ok5 &= sum(v * v for v in kk) >= 3 * abs(kk[0] * kk[1] * kk[2])**(2 / 3) - 1e-12
ok5 &= sp.simplify(3 - 3 * 1) == 0
mp.mp.dps = 40


def cubic_int(ev):
    rr = lambda c: 1 + ev * mp.cos(2 * c)
    r1 = lambda c: -2 * ev * mp.sin(2 * c)
    r2 = lambda c: -4 * ev * mp.cos(2 * c)
    A = lambda c: mp.sqrt(rr(c)**2 + r1(c)**2)
    km = lambda c: (rr(c)**2 + 2 * r1(c)**2 - rr(c) * r2(c)) / A(c)**3
    kp = lambda c: (rr(c) + 4 * ev * mp.cos(c)**2) / (rr(c) * A(c))      # r' cot chi = -4 eps cos^2 chi
    f = lambda c: (km(c)**2 + 2 * kp(c)**2)**mp.mpf(1.5) * rr(c)**2 * mp.sin(c)**2 * A(c) * 4 * mp.pi
    return mp.quad(f, [0, mp.pi / 4, mp.pi / 2, 3 * mp.pi / 4, mp.pi])


c0 = 6 * mp.sqrt(3) * mp.pi**2
v0, v1, v3 = cubic_int(mp.mpf(0)), cubic_int(mp.mpf(1) / 10), cubic_int(mp.mpf(3) / 10)
ok5 &= abs(v0 - c0) < mp.mpf("1e-25") and v1 > c0 + mp.mpf("1e-3") and v3 > c0 + mp.mpf("1e-3")
check("AM-GM on 10^4 triples; int |A|^3 = 6 sqrt3 pi^2 on the sphere, larger at eps = 1/10 and 3/10", ok5,
      f"ratios to 6sqrt3 pi^2: {mp.nstr(v1 / c0, 8)}, {mp.nstr(v3 / c0, 8)}")

# ------------------------------------------------------------------------------------------------------------- P6
print("P6  the corrugated chain")
s = sp.Symbol("sigma", real=True)
psi = sp.Function("psi")(s)
rho = sp.Function("rho")(s)
# sigma_2 dV over the orbit spheres: (2 k_mer k_sph + k_sph^2) rho^2 = 2 psi' rho sin psi + sin^2 psi; integrating
# 2 psi' rho sin psi = -2 rho (cos psi)' by parts with rho' = cos psi and rho = 0 at both ends gives 2 cos^2 psi
integrand = 2 * sp.diff(psi, s) * sp.sin(psi) / rho * rho**2 + sp.sin(psi)**2
byparts = sp.diff(-2 * rho * sp.cos(psi), s).subs(sp.diff(rho, s), sp.cos(psi)) + 2 * sp.cos(psi)**2 + sp.sin(psi)**2
ok6 = sp.simplify(integrand - byparts) == 0
ok6 &= sp.simplify(2 * sp.cos(psi)**2 + sp.sin(psi)**2 - (1 + sp.cos(psi)**2)) == 0
rows = []
for k in (2, 3):
    for A_ in (mp.mpf(1) / 10, mp.mpf(3) / 10, mp.mpf(1) / 2):
        val = 4 * mp.pi * mp.quad(lambda u: 1 + mp.cos(u + A_ * mp.sin(2 * k * u))**2, mp.linspace(0, mp.pi, 4 * k + 1))
        ok6 &= abs(val - 6 * mp.pi**2) < mp.mpf("1e-25")
A, tau = sp.symbols("A tau", real=True)
for k in (2, 3):
    ps = s + A * sp.sin(2 * k * s)
    cosp = sp.series(sp.cos(ps), A, 0, 3).removeO()
    sinp = sp.series(sp.sin(ps), A, 0, 3).removeO()
    rh = sp.expand(sp.integrate(sp.expand(sp.expand_trig(cosp.subs(s, tau))), (tau, 0, s)))
    W1i = sp.series(sp.diff(ps, s) * rh**2 + 2 * rh * sinp, A, 0, 3).removeO()
    W1 = sp.expand(4 * sp.pi * sp.integrate(sp.expand(W1i), (s, 0, sp.pi)))
    P = sp.Poly(W1, A)
    c0_, c1_, c2_ = [sp.simplify(P.coeff_monomial(A**j)) for j in range(3)]
    form = -4 * sp.pi**2 * k**2 * (4 * k**2 - 3) / sp.Integer(4 * k**2 - 1)**2
    ok6 &= sp.simplify(c0_ - 6 * sp.pi**2) == 0 and c1_ == 0 and sp.simplify(c2_ - form) == 0
    rows.append(f"k={k}: A^2 coefficient {c2_}")
check("int sigma_2 = 4pi int (1 + cos^2 psi) = 6 pi^2 at every amplitude; W1's A^2 coefficient as stated", ok6,
      "; ".join(rows))

print()
if FAILURES:
    print(f"test_352: {len(FAILURES)} FAILED")
    sys.exit(1)
print("test_352: ALL PASS")