(* ::Package:: *)

(*
Consolidated exact hypergeometric-closure certificate for Ramanujan
Challenge 2.8.

Index convention:

  n = N+1 >= 1,                    u = 2 n + 1,
  x = -z/(1-z),
  R_N(x) = C(x) M(0,x)...M(N-1,x),
  P_n(z) = (-z)^n R_N(-z/(1-z)).

Because x=-z/(1-z), this is exactly the qhat transform used in the proof:

  P_n(z) = (1-z)^n x^n R_N(x),  n=N+1.

The certificate proves that the first component is

  P_n(z)[[1]] = a_n 4F3(
      -n,-n-1/6,-n-1/2,-n-5/6;
      1-2n,1-2n,1-2n; z),

where

  a_1 = 18,
  a_(n+1)/a_n =
    576 n^2 (2n+1)^2 / ((6n+1)(6n+5)).

Everything below is an identity over QQ(n,z) or QQ(n,x).  No numerical
sampling is used.  In one stateless run it certifies:

  * the nonterminating 4F3 kernel and M_N K_(N+1)=K_N;
  * the terminating denominator formula and its normalization;
  * the initial CM-error annihilation;
  * the rank-one/Smith-valuation facts behind the error lattice; and
  * the one-step lattice cancellation which gains one power of x;
  * the exact characteristic quartic, spectral separation, and cyclic frame
    needed to transfer the first-column limit to all four official columns.

The proof uses the exact differential gauge behind the authoritative matrix
and one coefficientwise rational contiguity identity.
*)

ClearAll["Global`*"];

check[label_, condition_] := If[
  TrueQ[condition],
  Print["PASS: " <> label],
  Print["FAIL: " <> label]; Abort[]
];


(* The authoritative matrix, with x=1/R and
   236337691420383=(14R-567)/9 already substituted. *)
mm[u_, x_] := Module[
  {r = 1/x, w, a1, a2, a3, a4, b1, b2, b3, c1, c2, c3, c4},

  w = u (3u-2) (3u+2);

  a1 = r (144u^5-288u^4+144u^3)
       + (-99u^5+333u^4-229u^3-114u^2+40u+64);
  a2 = r (432u^4-864u^3+432u^2)
       + (-243u^4+909u^3-868u^2-80u+272);
  a3 = r (432u^3-864u^2+432u)
       + (-153u^3+648u^2-860u+360);
  a4 = r 144 (u-1)^2;

  b1 = r (-144u^3) + (9u^4+63u^3+158u^2+168u+64);
  b2 = r (216u^2) + (36u^3-189u^2-316u-168);
  b3 = r (108u) + (54u^2-189u-158);

  c1 = r^2 (-288u^3)
       + r (54u^4+378u^3+948u^2+1008u+384)
       + (18u^5+45u^4-251u^3-1086u^2-1384u-576);
  c2 = r^2 (-432u^2)
       + r (153u^4-657u^3+1292u^2+2064u+1072)
       + (-72u^4+702u^3-1069u^2-2508u-1512);
  c3 = r^2 (-216u)
       + r (180u^3-891u^2+1450u+1116)
       + (-108u^3+864u^2-1385u-1422);
  c4 = r^2 (-4)
       + r (6u^2-33u+58+14/9)
       + (-4u^2+32u-63);

  {
    {a1/w, a2/w, a3/w, a4/w},
    {-u^3, -3u^2, -3u, -1},
    {x b1/144, -x b2/72, -x b3/36, x (-2r-(2u-7))/2},
    {x^2 c1/288, x^2 c2/144, x^2 c3/72, x^2 c4/4}
  }
];


(* H_n is the nonterminating tail

     4F3(n,n+1/6,n+1/2,n+5/6;2n,2n,2n;z).

   Its Euler equation is lTail(theta) H_n=0. *)
lTail[n_, z_, t_] :=
  Expand[
    t (t+2n-1)^3
    - z (t+n) (t+n+1/6) (t+n+1/2) (t+n+5/6)
  ];

companion[n_, z_] := Module[{t, cc},
  cc = Table[Coefficient[lTail[n,z,t],t,j],{j,0,4}];
  {
    {0,1,0,0},
    {0,0,1,0},
    {0,0,0,1},
    -Take[cc,4]/cc[[5]]
  }
];

thetaMatrix[a_] := Map[z D[#,z]&,a,{2}];


(* In z-coordinates one recurrence step is G_n=(-z)M(2n+1).
   This is an exact differential gauge from the n+1 tail system to the
   n tail system. *)
gauge = Together[-z mm[2n+1,-z/(1-z)]];
gaugeResidual = Map[
  Factor,
  Together[
    companion[n,z].gauge
    - thetaMatrix[gauge]
    - gauge.companion[n+1,z]
  ],
  {2}
];
check[
  "authoritative matrix is the exact tail differential gauge",
  gaugeResidual === ConstantArray[0,{4,4}]
];


(* Initial horizontal row. *)
cRow[x_] := {18/x+159/4,54/x+131/2,54/x+27,18/x};
pBase = Together[-z cRow[-z/(1-z)]];
baseResidual = Map[
  Factor,
  Together[
    Map[z D[#,z]&,pBase]
    + pBase.companion[1,z]
  ]
];
check[
  "the compact denominator seed is a horizontal adjoint row",
  baseResidual === ConstantArray[0,4]
];
check[
  "base terminating polynomial",
  Factor[pBase[[1]]-18 (1-77z/24)] === 0
];


(* Recover every horizontal-row component from its first component p.
   Here t acts as theta=z d/dz on p. *)
tailCoefficients = Table[
  Coefficient[lTail[n,z,t],t,j],
  {j,0,4}
];

thetaOperator[poly_] := Together[
  Sum[
    z D[Coefficient[poly,t,j],z] t^j
    + Coefficient[poly,t,j] t^(j+1),
    {j,0,Exponent[poly,t]}
  ]
];

pOp3 = Together[tailCoefficients[[5]]/tailCoefficients[[1]] t];
pOp2 = Together[
  tailCoefficients[[4]]/tailCoefficients[[5]] pOp3
  - thetaOperator[pOp3]
];
pOp1 = Together[
  tailCoefficients[[3]]/tailCoefficients[[5]] pOp3
  - thetaOperator[pOp2]
];
pOps = {1,pOp1,pOp2,pOp3};

(* The last horizontal equation is precisely the terminating 4F3 equation

     theta(theta-2n)^3 p
       - z(theta-n)(theta-n-1/6)(theta-n-1/2)(theta-n-5/6)p=0.
*)
adjointResidualOperator = Factor[
  Together[
    thetaOperator[pOp1] + 1
    - tailCoefficients[[2]]/tailCoefficients[[5]] pOp3
  ]
];
terminatingOperator =
  Expand[
    t (t-2n)^3
    - z (t-n) (t-n-1/6) (t-n-1/2) (t-n-5/6)
  ];
adjointFactor = -72/(n(2n+1)(6n+1)(6n+5)z);
check[
  "horizontal elimination gives the terminating 4F3 operator",
  Factor[
    Together[
      adjointResidualOperator-adjointFactor terminatingOperator
    ]
  ] === 0
];


(* The first component after one gauge step is D_n(theta) p. *)
stepOperator = Factor[Together[pOps.gauge[[All,1]]]];
check[
  "the scalar step operator has only z-degrees zero and one",
  Together[
    stepOperator
    - Coefficient[stepOperator,z,0]
    - z Coefficient[stepOperator,z,1]
  ] === 0
];

d0[q_] := Together[Coefficient[stepOperator,z,0] /. t->q];
d1[q_] := Together[Coefficient[stepOperator,z,1] /. t->q];


(* Coefficients h_(n,k) of the normalized terminating 4F3.

   Instead of expanding Pochhammer symbols, the proof needs only:
     withinRatio = h_(n,k)/h_(n,k-1),
     crossRatio  = h_(n+1,k)/h_(n,k).
*)
withinRatio =
  ((k-1-n) (k-1-n-1/6) (k-1-n-1/2) (k-1-n-5/6)) /
  ((k-2n)^3 k);

crossRatio =
  Product[
    (n+1+a)/(n+1+a-k),
    {a,{0,1/6,1/2,5/6}}
  ] * (((k-2n-1)(k-2n))/(2n(2n+1)))^3;

normalizationRatio =
  576 n^2 (2n+1)^2 / ((6n+1)(6n+5));

(* k=0. *)
check[
  "constant-term normalization recurrence",
  Factor[Together[d0[0]-normalizationRatio]] === 0
];

(* Generic coefficient 1<=k<=n:

   d0(k) h_(n,k) + d1(k-1) h_(n,k-1)
      = normalizationRatio h_(n+1,k).
*)
genericCoefficientResidual = Factor[
  Together[
    d0[k]
    + d1[k-1]/withinRatio
    - normalizationRatio crossRatio
  ]
];
check[
  "generic all-n terminating contiguity coefficient",
  genericCoefficientResidual === 0
];

(* The new top coefficient k=n+1 is a boundary case because h_(n,n+1)=0. *)
topCoefficientRatio =
  -(n+7/6)(n+3/2)(n+11/6)/(8(2n+1)^3);
check[
  "terminating top-coefficient boundary",
  Factor[
    Together[
      d1[n]-normalizationRatio topCoefficientRatio
    ]
  ] === 0
];


(* ---------------------------------------------------------------------- *)
(* Nonterminating kernel contiguity.                                       *)
(* ---------------------------------------------------------------------- *)

(* The analytic tail solution H_(n+1) has initial Euler jet e_1 at z=0.
   The gauge sends that vector to normalizationRatio e_1.  Uniqueness of the
   exponent-zero Frobenius solution therefore gives

      gauge J_(n+1) = normalizationRatio J_n.

   Since gauge=(-z)M and K_n=kappa_n z^n J_n, the following reciprocal
   normalization is exactly M_N K_(N+1)=K_N. *)
tailColumnAtOrigin = Map[
  Factor,
  Limit[gauge[[All,1]],z->0]
];
check[
  "tail gauge selects the exponent-zero analytic solution",
  tailColumnAtOrigin === {normalizationRatio,0,0,0}
];

kappaRatio =
  -(6n+1)(6n+5)/(576 n^2 (2n+1)^2);
check[
  "nonterminating kernel normalization and exact contiguity",
  Factor[Together[-kappaRatio normalizationRatio-1]] === 0
];


(* ---------------------------------------------------------------------- *)
(* Exact initial annihilation of the CM error rows.                        *)
(* ---------------------------------------------------------------------- *)

binomialRows = {
  {1,0,0,0},
  {1,1,0,0},
  {1,2,1,0},
  {1,3,3,1}
};

(* K_0=((theta-1)^j(y-1))_(j=0)^3. *)
kFormal = Table[(t-1)^j y-(-1)^j,{j,0,3}];
fFormal = {y-1,t y,t^2 y,t^3 y};
check[
  "binomial jet converts K_0 to the CM first jet",
  Expand[binomialRows.kFormal-fFormal] === ConstantArray[0,4]
];

(* The 3F2 equation in x coordinates is

     72 theta^3 y + 108 x theta^2 y + 46 x theta y + 5 x y = 0.

   This makes C K_0=-5/4, hence
   ((5/4)binomialRows+fFormal C) K_0=0. *)
cOperator = Factor[
  Sum[cRow[x][[j+1]] (t-1)^j,{j,0,3}]
];
cConstant = Factor[
  -Sum[cRow[x][[j+1]] (-1)^j,{j,0,3}]
];
cmDifferentialOperator =
  72t^3+108x t^2+46x t+5x;
check[
  "compact denominator row reduces to the 3F2 operator",
  Factor[cOperator/cmDifferentialOperator] === 1/(4x)
];
check[
  "exact CM-error seed annihilation constant",
  cConstant === -5/4
];


(* ---------------------------------------------------------------------- *)
(* Rank-one and regular-annihilator-lattice algebra.                       *)
(* ---------------------------------------------------------------------- *)

(* Multiplying the pole row by x makes the transfer analytic at x=0. *)
mHat = Map[
  Cancel,
  Together[
    DiagonalMatrix[{x,1,1,1}].mm[2n+1,x]
  ],
  {2}
];
mHat0 = Map[Factor[(#/.x->0)]&,mHat,{2}];
check[
  "regularized transfer has rank one at x=0",
  MatrixRank[mHat0] === 1
];

detExpected =
  x^3 (2n-1)^3 (2n)^2 (6n-5)^3 (6n-1)^3 /
  (20736 (2n+1)(6n+1)(6n+5));
check[
  "regularized determinant has exact x-adic order three",
  Factor[Together[Det[mHat]-detExpected]] === 0
];
check[
  "Smith valuations are exactly (0,1,1,1)",
  Factor[Cancel[Det[mHat]/x^3]/.x->0] =!= 0
];

(* For H_n=1+h_n z+O(z^2), every nonzero Euler derivative divided by
   H_n starts as h_n z=-h_n x+O(x^2).  Thus the normalized syzygy basis

      W_n = [ -theta^i H_n/H_n | e_i ], i=1,2,3,

   has the following first-order truncation.  The exact zero below proves
   that (W_n M_N)[:,2:4] is divisible by x.  This is the algebraic step
   which gains one valuation at every recurrence step. *)
hTail = Factor[
  n(n+1/6)(n+1/2)(n+5/6)/(2n)^3
];
wLeading = {
  {hTail x,1,0,0},
  {hTail x,0,1,0},
  {hTail x,0,0,1}
};
latticeStepLeading = Map[
  Cancel,
  Together[
    (wLeading.mm[2n+1,x])[[All,2;;4]]
  ],
  {2}
];
check[
  "regular annihilator transfer gains one power of x",
  Map[Factor[(#/.x->0)]&,latticeStepLeading,{2}]
    === ConstantArray[0,{3,3}]
];


(* ---------------------------------------------------------------------- *)
(* Exact spectral and cyclic-frame closure for all four official columns. *)
(* ---------------------------------------------------------------------- *)

rOfficial = 151931373056001;
sLimit = {
  {64rOfficial-44,96rOfficial-54,48rOfficial-17,8rOfficial},
  {-8,-12,-6,-1},
  {1/rOfficial,-4/rOfficial,-6/rOfficial,-2/rOfficial},
  {
    2/rOfficial^2,
    (17rOfficial-8)/rOfficial^2,
    4(5rOfficial-3)/rOfficial^2,
    (6rOfficial-4)/rOfficial^2
  }
};
qQuartic[lam_] := (
  rOfficial^2 lam^4
  - (64rOfficial^3-56rOfficial^2-4) lam^3
  + (48rOfficial^2-262rOfficial+220) lam^2
  - (12rOfficial-8) lam+1
);

check[
  "balanced characteristic polynomial is the authoritative quartic",
  Factor[
    CharacteristicPolynomial[sLimit,lam]
    -qQuartic[lam]/rOfficial^2
  ] === 0
];
check[
  "characteristic quartic is irreducible",
  TrueQ[IrreduciblePolynomialQ[qQuartic[lam]]]
];

roucheMargin =
  64rOfficial^3-105rOfficial^2+250rOfficial-217;
check[
  "Rouche separation has three roots in the unit disk",
  roucheMargin > 0
];

eFirst = {1,0,0,0};
cyclicFrame = Transpose[{
  eFirst,
  sLimit.eFirst,
  MatrixPower[sLimit,2].eFirst,
  MatrixPower[sLimit,3].eFirst
}];
cyclicDetExpected =
  -4(27rOfficial-11)(128rOfficial^2-149rOfficial-43)/
  rOfficial^6;
check[
  "limiting first-coordinate cyclic frame is invertible",
  Factor[Together[Det[cyclicFrame]-cyclicDetExpected]] === 0
  && cyclicDetExpected != 0
];

wLeft = {
  10+(44-7rOfficial)lam+(4+12rOfficial^2)lam^2
    +rOfficial^2 lam^3,
  2((-23+40rOfficial)
    +(-108+194rOfficial-28rOfficial^2)lam
    +(-27rOfficial^2+48rOfficial^3)lam^2),
  (-32+71rOfficial)
    +(-68+198rOfficial-8rOfficial^2)lam
    +(-17rOfficial^2+48rOfficial^3)lam^2,
  2rOfficial(8+(17+3rOfficial)lam+4rOfficial^2 lam^2)
};
leftResidual = Map[
  Factor,
  Together[
    wLeft.(lam IdentityMatrix[4]-sLimit)
    -{qQuartic[lam],0,0,0}
  ]
];
check[
  "dominant left eigenvector has four nonzero coordinates",
  leftResidual === ConstantArray[0,4]
  && And@@Map[
    Exponent[PolynomialGCD[qQuartic[lam],#],lam] === 0&,
    wLeft
  ]
];


Print["PASS: consolidated exact hypergeometric-closure certificate"];
Print[
  "M_N K_(N+1)=K_N with ",
  "kappa_(n+1)/kappa_n=-(6n+1)(6n+5)/",
  "(576n^2(2n+1)^2)"
];
Print[
  "P_n(z)[[1]]/a_n = ",
  "4F3(-n,-n-1/6,-n-1/2,-n-5/6;",
  "1-2n,1-2n,1-2n;z)"
];
Print[
  "a_(n+1)/a_n = ",
  "576 n^2 (2n+1)^2/((6n+1)(6n+5)), a_1=18"
];