
(* MK 2020-12-20 *)

MyNullSpace = If[ MemberQ[$ContextPath, "HolonomicFunctions`"], HolonomicFunctions`Private`EANullSpaceNew, NullSpace ];

MyLinearSolve[A_, b_] := Module[ {B},
  B = DeleteCases[MyNullSpace[Transpose[Append[Transpose[A], -b]]], {___, 0}];
  If[ Length[B] == 0, Throw["linear system has no solution"] ];
  B = First[B];
  Return[Most[B]/Last[B]];
];

MyRowReduce = RowReduce

ClearAll[MyPolynomialRemainder];
If[ MemberQ[$ContextPath, "HolonomicFunctions`"],
  MyPolynomialRemainder[p_, q_, x_] :=
    Last[HolonomicFunctions`Private`MyPolynomialQuotientRemainder[Together[p], Together[q], x]];
, MyPolynomialRemainder = PolynomialRemainder
];

xGCD = If[ MemberQ[$ContextPath, "HolonomicFunctions`"], HolonomicFunctions`Private`MyPolynomialExtendedGCD, PolynomialExtendedGCD ];

(**
 * Modular Inverse p^(-1) mod q
 **)
ModInv[p_, q_, y_] := xGCD[p, q, y][[2, 1]];

(**
 * square free part
 **)
SQFP[p_, x_] := Together[p/PolynomialGCD[D[p,x],p]];

(**
 * Normal form of an algebraic function given as rational function in generators
 **)
NF[fun_, minpoly_, x_, y_] :=
  Module[ {num, den},
    {num, den} = {Numerator[#], Denominator[#]}&[Together[fun]];
    If[ Exponent[den, y]==0 && Exponent[num, y] < Exponent[minpoly, y], Return[num/den] ];
    Return[Together[MyPolynomialRemainder[num*ModInv[den, minpoly, y], minpoly, y]]];
  ];

(**
 * Normal form of a vector wrt a subspace.
 **)
NF[v_, U_] :=
  Module[ {eval, red, supp, M, ker},
    M = Append[Append[#, 0]& /@ U, Append[v, 1]];
    eval = Thread[# -> (Prime /@ (10 + Range[Length[#]]))]&[Variables[{v, U}]];
    red = Last[DeleteCases[RowReduce[M /. eval, Modulus -> 45007], {___, 0}]];
    supp = Select[Range[Length[red]], red[[#]]===0&];
    If[ Length[supp]==0, Return[v] ];
    ker = DeleteCases[MyNullSpace[Table[M[[i, j]], {j, supp}, {i, 1, Length[M]}]], {___, 0}];
    If[ Length[ker] == 0, Throw["normal form computation failed; unlucky prime?"] ];
    red = Last[Together[{Last[ker]}.M]];
    Return[Together[Most[red]/Last[red]]];
  ];

(**
 * Write the given function as a linear combination of the given basis elements. assuming fun is in normal form.
 **)
Coeffs[fun_, basis_, y_] := Together[MyLinearSolve[Table[Coefficient[b, y, i], {i, 0, Length[basis] - 1}, {b, basis}], Table[Coefficient[fun, y, i], {i, 0, Length[basis] - 1}]]]
 
(**
 * Derivative of an algebraic function
 **)
Diff[fun_, minpoly_, x_, y_] := NF[D[#, x] - D[#, y]*D[minpoly, x]*DiffGen[minpoly, y], minpoly, x, y]&[NF[fun, minpoly, x, y]];
DiffGen[minpoly_, y_] := DiffGen[minpoly, y] = ModInv[D[minpoly, y], minpoly, y] (* cached *)

(** 
 * Derivation operator of an algebraic function field wrt a given basis. 
 **)
Basis2Matrix[basis_, minpoly_, x_, y_] :=
   Basis2Matrix[basis, minpoly, x, y] = Table[Coeffs[Diff[b, minpoly, x, y], basis, y], {b, basis}];

(**
 * Compute a basis for a k[x]-module that contains fun as well as the submodule generated by the basis elements. Assuming w is in normal form. 
 **)
ExtendModule[fun_, basis_, x_, y_] := Module[ {c, d, b, g, n, t},
  Print["extending module"];
  n = Length[basis];
  c = Together[Coeffs[fun, basis, y]];
  d = PolynomialLCM@@(Denominator/@c); c = Together[d*c];
  If[ Exponent[d, x] == 0, Return[basis] ]; (* fun is already in the module *)
  b = GroebnerBasis[Join[d Array[g, n], {c . Array[g, n]}, Flatten[Table[g[i]g[j],{i,1,n},{j,1,n}]]], Join[{x}, Reverse[Array[g, n]]], CoefficientDomain -> RationalFunctions]; (* hnf *)
  b = Select[b /. g[i_] -> t g[i], Exponent[#, t]==1&] /. t -> 1;
  Return[Together[Table[Sum[Coefficient[p, g[i]]basis[[i]], {i, 1, n}] /. g[_] -> 0, {p, b}]/d]];
]

BasisOfIntegralElements[minpoly_, x_, y_] := Module[ {n, m, p},
  m = CoefficientList[Numerator[Together[minpoly]], y]; p = 1; n = Exponent[minpoly, y];
  While[ Exponent[Last[m], x] > 0, q = SQFP[Last[m], x]; p *= q; m = m * (q^(n+1 - Range[n+1])); m = Together[m/PolynomialGCD@@m] ];
  Return[(p y)^(Range[n]-1)];   
]

(**
 * Construct a suitable basis extending the module generated by the given basis. also returns the corresponding derivation matrix. 
 * the given basis must consist of integral elements; specify basis==Null in order to get any suitable basis for the field k(x,y).
 **)
MakeSuitable[basis_, minpoly_, x_, y_] :=
  Module[ {n, m, p, q, d, g},
     n = Exponent[minpoly, y];
     b = If[ ListQ[basis] && Length[basis] === n, basis, BasisOfIntegralElements[minpoly, x, y] ];
     m = Basis2Matrix[b, minpoly, x, y];
     d = Extract[m, Position[(Exponent[PolynomialGCD[D[#, x], #], x]& /@ PolynomialLCM@@@Map[Denominator, m, {2}]) /. 0 -> Null, _Integer]];
     While[ Length[d] > 0,
       fun = First[d]; d = PolynomialLCM@@(Denominator /@ fun);
       fun = Together[(SQFP[d, x]*fun).b];
       b = ExtendModule[fun, b, x, y];
       m = Basis2Matrix[b, minpoly, x, y]; (* recomputation avoidable? *)
       d = Extract[m, Position[(Exponent[PolynomialGCD[D[#, x], #], x]& /@ PolynomialLCM@@@Map[Denominator, m, {2}]) /. 0 -> Null, _Integer]];
     ];
     Return[b];
  ];

MakeSuitableAtInfinity[basis_, minpoly_, x_, y_] :=
  Module[ {deg, rat, V, m, d, e, i, j, k, v, failed, min, xbar},
    deg[rat_] := (Exponent[Numerator[#], x] - Exponent[Denominator[#], x])&[Together[rat]];
    n = Exponent[minpoly, y]; V = basis; m = Basis2Matrix[V, minpoly, x, y];
    failed = False;
    Do[ (* fast but naive *)
      e = -Max[Table[deg[m[[i,j]]], {j, 1, n}]]-1;
      If[ e < 0,
        If[ Max[Table[deg[m[[k,i]]], {k, 1, i-1}]] > e || e === -deg[m[[i,i]]], failed=True; Break[] ];
        V[[i]] = V[[i]] x^e; v = Table[1, {n}]; v[[i]] = x^e;
        m = Together[(D[DiagonalMatrix[v], x] + DiagonalMatrix[v].m).Inverse[DiagonalMatrix[v]]];
      ];
    , {i, 1, n}];
    If[ !failed, Return[{V, m}] ];
    (* slow but complete fall back code *)
    d = 0; min = Numerator[Together[minpoly /. x -> xbar]];
    While[Exponent[Coefficient[min, y, n], xbar, Min]>0, d+=1; min = Numerator[Together[min /. y -> y/xbar]] ];
    V = Together[Table[v x^(-Max[Table[i d + Exponent[Coefficient[v, y, i], x], {i, 0, n-1}]]), {v, V}] /. x -> 1/xbar];
    d = Extract[m, Position[Max /@ Map[deg, m, {2}], _?(#>0&)]] /. x -> 1/xbar;
    While[ Length[d] > 0,
      fun = Together[1/xbar*(First[d].V)]; (* v integral at infinity => x*v' integral at infinity *)
      V = ExtendModule[fun, V, xbar, y];
      m = Basis2Matrix[V /. xbar -> 1/x, minpoly, x, y]; (* recomputation avoidable? *)
      d = Extract[m, Position[Max /@ Map[deg, m, {2}], _?(#>0&)]] /. x -> 1/xbar;
    ];
    Return[{V /. xbar -> 1/x, m}];
  ];

(**
 * Perform one step of lazy hermite reduction. if basis or mat are Null, they are chosen by the method. 
 **)
LHRstep[fun_, basis_, mat_, minpoly_, x_, y_] :=
   Module[ {n, b, m, c, d, e, u, v, w, f, g, h, A},
     If[ !ListQ[fun], Return[LHRstep[Coeffs[fun, basis, y], basis, mat, minpoly, x, y]] ];
     n = Exponent[minpoly, y];
     b = If[ ListQ[basis] && Length[basis] === n, basis, MakeSuitable[Null, minpoly, x, y] ];
     m = Together[If[ Dimensions[mat] === {n, n} && b === basis, mat, Basis2Matrix[b, minpoly, x, y] ]];
     e = PolynomialLCM@@(Denominator/@Flatten[m]); m = Together[e*m];
     If[ Exponent[PolynomialGCD[D[e,x], e], x] > 0, Throw["The basis is not suitable"] ];
     c = Together[If[ ListQ[fun], If[b === basis, fun, Throw["Integrand given by coefficients, but no basis given"]], Coeffs[NF[fun, minpoly, x, y], b, y] ]];
     d = PolynomialLCM@@Append[Denominator/@c, e]; c = d*c;
     If[ Exponent[d, x] <= 1, Return[ {0, fun, basis, mat} ] ]; (* already reduced *)
     f = Select[FactorList[d], (True && Exponent[First[#], x]>0)&]; k = Max[Last /@ f];
     If[ k == 1, Return[ {0, fun, basis, mat} ] ]; (* already reduced *)
     v = Times@@Cases[f, {u_, k} -> u]; u = Together[d/v^k];
     A = MyNullSpace[Join[PolynomialQuotient[u*v, e, x]Transpose[m] - (k-1) u D[v, x] IdentityMatrix[n], List /@ c, 2]];
     A = Together[((PolynomialLCM@@(Denominator/@#))& /@ A)*A];
     Which[
       Length[A] > 1 || A[[1,-1]] === 0, (* homogeneous system has a solution over k(x) ==> extend module by Thm 2 *)
         Print["Updating basis according to Thm 2"];
         f = (c.b)/d;
         b = ExtendModule[u/v (Reverse[Rest[First[Cases[MyRowReduce[Reverse /@ A], {0, __}]]]].b), b, x, y];
         m = Basis2Matrix[b, minpoly, x, y]; (* recomputation avoidable? *)
         Return[ {0, Coeffs[f, b, y], b, m} ],
       Exponent[PolynomialGCD[A[[1, -1]], v], x] === 0, (* inhomogeneous system has a solution over k[x]/<v> ==> hermite reduction succeeds *)
         Print["Reduction step succeeded, v=", InputForm[v]];
         g = (MyPolynomialRemainder[#, v, x]/v^(k-1)&) /@ (-Most[First[A]]ModInv[A[[1,-1]], v, x]);
         h = Together[c/d - (D[g, x] + (Transpose[m].g)/e)];
         Return[ {g.b, h, b, m/e} ],
       True, (* inhomogeneous system has a unique solution over k(x) but not over k[x]/<v> ==> extend module by Thm 3 *)
         Print["Updating basis according to Thm 3"];
         f = (c.b)/d;
         b = ExtendModule[Together[u*v/PolynomialGCD[A[[1, -1]], v]^2] (Most[First[A]].b), b, x, y];
         m = Basis2Matrix[b, minpoly, x, y]; (* recomputation avoidable? *)
         Return[ {0, Coeffs[f, b, y], b, m} ];
     ];
   ];

(**
 * Perform the full lazy hermite reduction. 
 **)
LHR[fun_, basis_, minpoly_, x_, y_] := Module[ {g, h, b, m},
  b = MakeSuitable[basis, minpoly, x, y];
  {g, h, b, m} = FixedPoint[({#[[1]], 0, 0, 0} + LHRstep[#[[2]], #[[3]], #[[4]], minpoly, x, y])&, {0, NF[fun, minpoly, x, y], b, Basis2Matrix[b, minpoly, x, y]}];
  {g, h, b, m}
];

CheckLHR[fun_, minpoly_, x_, y_, opts___] := Module[ {g, h, b, m},
  b = Basis /. {opts} /. {Basis -> Null};
  {g, h, b, m} = LHR[fun, b, minpoly, x, y];
  If[ ListQ[h], h = h.b ]; 
  NF[fun - (Diff[g, minpoly, x, y] + h), minpoly, x, y]===0
];

(* =================================================================== *)

(**
 * Polynomial reduction: write fun in the form (p1/u V)' + p2/a V with p2 minimal.
 **)
PR[u_, fun_, mat_] := Module[ {lambda, n, B, a, P, P1, P2, d, degB, dega, lca, imphi, as, ker, basis, cap, mod2vec, vec2mod},
  n = Length[fun];
  If[ !ListQ[fun] || Dimensions[mat] =!= {n, n}, Throw["matrix format mismatch."] ];
  a = PolynomialLCM@@(Denominator/@Flatten[mat]); B = Together[a*mat]; dega = Exponent[a, x]; lambda = Exponent[a, x, Min];
  If[ Exponent[PolynomialGCD[a/x^lambda, D[a/x^lambda, x]], x] > 0, Throw["denominator of matrix not square free."] ];
  If[ Max[Exponent[B, x]] >= dega, Throw["matrix numerator degree exceeds denominator degree"] ];
  B = Together[B - D[u, x] a/u IdentityMatrix[n]]; (* Lixin's modification *)
  degB = Max[dega, Exponent[B, x]];
  P = Together[u*a*fun];
  If[ Exponent[PolynomialLCM@@(Denominator/@P), x] > 0 || MyPolynomialRemainder[PolynomialLCM@@P, u, x] =!= 0,
    Throw["denominator of integrand does not divide denominator of matrix"] ];
  d = Max[Exponent[P, x]];
  If[ Max[Exponent[B, x]] == dega - 1,
    lca = Coefficient[a, x, dega];
    d = Max[d, Select[-Together[Eigenvalues[Coefficient[B, x, Max[Exponent[B, x]]]]/lca], IntegerQ]];
  ];
  If[ d < 0, Return[ {fun, 0 fun} ] ];
  mod2vec = Reverse[Flatten[Table[Coefficient[#[[j]], x, i], {j, 1, n}, {i, 0, d + degB}]]]&; (* checked *)
  vec2mod = Table[Sum[#[[-((d+degB+1)(j-1) + i + 1)]]x^i, {i, 0, d + degB}], {j, 1, n}]&; (* checked *)
  basis = Flatten[Table[x^i IdentityMatrix[n], {i, 0, d}], 1];
  If[ False, 
    imphi = mod2vec[a D[#, x] + #.B]& /@ basis;
    as = mod2vec /@ Expand[u Flatten[Table[x^i IdentityMatrix[n], {i, 0, d + degB - Exponent[u, x]}], 1]];
    ker = MyNullSpace[ Transpose[ Join[imphi, as] ] ];
    ker = DeleteCases[MyRowReduce[Take[#, Length[imphi]]& /@ ker], {0..}];
    cap = ker.Join[imphi, Most[IdentityMatrix[Length[imphi] + 1]], 2]; (* im(phi) cap u k[x]^n with bookkeeping columns *)
    red = Last[DeleteCases[MyRowReduce[Append[cap, Join[mod2vec[P], Last[IdentityMatrix[Length[imphi] + 1]]]]], {___, 0}]];
    red = -Most[red]/Last[red];
    P2 = -Together[vec2mod[Take[red, (d+degB+1)*n]]/u/a];
    P1 = Together[(Drop[red, (d+degB+1)*n].basis)/u];
  , (* alternative code: intersection using polynomial remainder and omit computation of certificate *)
    imphi = (a D[#, x] + #.B)& /@ basis;
    ker = MyNullSpace[ Transpose[mod2vec /@ Map[MyPolynomialRemainder[#, u, x]&, imphi, {2}]] ];
    cap = ker.(mod2vec /@ imphi);
    (*
    red = Last[DeleteCases[MyRowReduce[Append[Append[#, 0]& /@ cap, Append[mod2vec[P], -1]]], {___, 0}]];
    P2 = -Together[vec2mod[Most[red]/Last[red]]/u/a];
    Print[MatchQ[Together[P2 - vec2mod[NF[mod2vec[P], cap]]/u/a], {0..}]];
    *)
    P2 = Together[vec2mod[NF[mod2vec[P], cap]]/u/a];
    P1 = 0*P2; (* certificate not computed *)
  ];
  Return[{P1, P2}];
];

(*** only checks whether fun=D(P1)+P2, not whether P2 is minimal ***)
CheckPR[u_, fun_, minpoly_, x_, y_, opts___] := Module[ {B},
  B = Basis /. {opts} /. {Basis :> BasisOfIntegralElements[minpoly, x, y]};
  mat = Basis2Matrix[B, minpoly, x, y];
  {P1, P2} = PR[u, Coeffs[fun, B, y], mat];
  P1 = P1.B; P2 = P2.B;
  NF[fun - (Diff[P1, minpoly, x, y] + P2), minpoly, x, y]===0
];

(**
 * Given an algebraic function
 **)
Int[fun_, minpoly_, x_, y_, opts___] := Module[ {W, V, d, e, u, t, a, g, gg, p, P, Q, Q1, Q2},
  W = Basis /. {opts} /. {Basis -> BasisOfIntegralElements[minpoly, x, y]};
  (* lazy hermite reduction *)
  {g, h, W, m} = LHR[fun, W, minpoly, x, y];
  (* split h *)
  e = PolynomialLCM@@(Denominator /@ Flatten[Together[m]]);
  d = PolynomialLCM@@(Denominator /@ Together[e h]);
  {gg, {s, t}} = xGCD[d, e, x];
  If[ Exponent[gg, x] =!= 0, Throw["unexpected common divisor encountered."] ];
  {s, t} = Together[{s, t}/gg];
  Q = ((# s + e PolynomialQuotient[# t, d, x])& /@ Together[d e h])/e;
  P = Together[h - Q];
  (* make basis suitable at infinity *)
  V = W; 
  a = PolynomialLCM@@(Denominator /@ Flatten[m]);
  If[ Max[Exponent[a*m, x]] >= Exponent[a, x],
    {V, m} = MakeSuitableAtInfinity[V, minpoly, x, y];
    Q = Coeffs[Q.W, V, y];
  ];
  (* denominator bound for certificate part *)
  a = PolynomialLCM@@(Denominator /@ Flatten[m]);
  u = Product[If[Exponent[First[p], x]>0 && MyPolynomialRemainder[a, First[p], x]===0, First[p]^Floor[Last[p]/2], 1],
              {p, FactorList[Discriminant[minpoly, y]]}];
  (* polynomial reduction *)
  {Q1, Q2} = PR[u, Q, m];
  (* done *)
  Return[{g + Q1.V, {P, W}, {Q2, V}, h}];
];

CheckInt[fun_, minpoly_, x_, y_, opts___] := Module[ {g,h,p,q},
  {g,p,q,h} = Int[fun, minpoly, x, y, opts];
  NF[fun - (Diff[g, minpoly, x, y] + p[[1]].p[[2]] + q[[1]].q[[2]]), minpoly, x, y] === 0
];

(**
 * Creative telescoping: given fun, find L such that there is a g such that L(fun)=D(g) 
 **)
CT[fun_, minpoly_, x_, y_, t_] := Module[ {B, rem, L, q, m, p, n},
  If[ Exponent[Denominator[Together[minpoly]], x] > 0,
    Return[CT[fun, Numerator[Together[minpoly]], x, y, t]];
  ];
  If[ Exponent[Last[CoefficientList[minpoly, y]], x] > 0,
    (* preprocessing: replace y by an integral generator *)
    m = CoefficientList[Numerator[Together[minpoly]], y]; p = 1; n = Exponent[minpoly, y];
    While[ Exponent[Last[m], x] > 0,
       q = SQFP[Last[m], x]; p *= q; m = m * (q^(n+1 - Range[n+1])); m = Together[m/PolynomialGCD@@m]
    ];
    Return[CT[fun /. y -> y/p, m.(y^(Range[n + 1] - 1)), x, y, t]];
  ];
  rem = {Rest[Int[fun, minpoly, x, y]]}; (* list of remainders as produced by Int *)
  If[ MatchQ[Join[rem[[1,1,1]], rem[[1,2,1]]], {0..}], Return[1] ];
  While[True,
    Print["ct order ", Length[rem]-1];
    (* check if remainders are linearly dependent *)
    L = Numerator[Together[Sum[q[i]Join[rem[[i, 1, 1]], rem[[i, 2, 1]]], {i, 1, Length[rem]}]]];
    L = Flatten[CoefficientList[#, x]& /@ L];
    L = Table[Coefficient[#, q[i]], {i, 1, Length[rem]}]& /@ L;
    L = MyNullSpace[L];
    (* if so, done *)
    If[ Length[L] > 0, Return[Collect[Numerator[Together[First[L].Table[Der[t]^i, {i, 0, Length[rem]-1}]]], Der[t]]]];
    (* otherwise, compute a new remainder *)
    (* f = Diff[rem[[-1,1,1]].rem[[-1,1,2]] + rem[[-1,2,1]].rem[[-1,2,2]], minpoly, t, y]; *)
    f = Diff[rem[[-1,-1]].rem[[-1,1,2]], minpoly, t, y];
    AppendTo[rem, Rest[Int[f, minpoly, x, y, Basis -> rem[[-1,1,2]]]]];
    (* if basis has changed, adjust earlier remainders to new basis *)
    If[ rem[[-1,1,2]]=!=rem[[-2,1,2]], Throw["NYI: basis adaption"] ];    
  ]
];