
(* double integration of rational functions via single integration of algebraic functions *)

If[ !MemberQ[$ContextPath, "HolonomicFunctions`"], << HolonomicFunctions.m ] 

(* Several methods for computing double integrals of rational functions in three variables. 
   They all take as input the integrand rat, and the three variables t, x, y. The last two
   variables are integrated out. Output is a pair {T, P} where T is the timing and P is the 
   telescoper. The methods are timeconstrained to 1000h. *)

(* 1st way: Chyzak + Chyzak *)
ctRat1[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {ct1, ct2},
    ct1 = CreativeTelescoping[rat, Der[x], {Der[y], Der[t]}];
    ct2 = CreativeTelescoping[ct1[[1]], Der[y], {Der[t]}];
    ct2[[1]]       
  ], 3600000, $Failed]]  

(* 2nd way: Chyzak + Koutschan *)
ctRat2[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {ct1, ct2},
    ct1 = CreativeTelescoping[rat, Der[x], {Der[y], Der[t]}];
    ct2 = FindCreativeTelescoping[ct1[[1]], Der[y], {Der[t]}];
    ct2[[1]]       
  ], 3600000, $Failed]]  

(* 3rd way: Koutschan *)
ctRat3[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {ct1},
    ct1 = FindCreativeTelescoping[rat, {Der[x], Der[y]}, {Der[t]}];
    ct1[[1]]
  ], 3600000, $Failed]]

(* 4th way: Equivalence + AE2DE + Chyzak *)
ctRat4[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {num, den, res, z, r, f, ops},
    {num, den} = {Numerator[#], Denominator[#]}&[Last[HermiteReduce[rat, y]]]; 
    res = FactorList[Resultant[num - z D[den, y], den, y]];
    ops = {};
    Do[
      If[ Exponent[res[[i, 1]], z] > 0, 
        AppendTo[ops, 
          CreativeTelescoping[r[f[z, res[[i, 1]]], 1] /. {f -> Function, r -> Root}, Der[x], {Der[t]}][[1]]];
      ]
    , {i, 1, Length[res]}];
    DFinitePlus@@ops
  ], 3600000, $Failed]]

(* 5th way: Equivalence + Residues + AE2DE + Chyzak *)
ctRat5[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {num, den, res, z, r, f, ops, kill, ann},
    {num, den} = {Numerator[#], Denominator[#]}&[Last[HermiteReduce[rat, y]]]; 
    res = FactorList[Resultant[num - z D[den, y], den, y]];
    ops = {};
    Do[
      If[ Exponent[res[[i, 1]], z] > 0, 
        kill = ResidueKiller[res[[i, 1]], t, x, z];
        ann = r[f[z, res[[i, 1]]], 1] /. {f -> Function, r -> Root};
        ann = Annihilator[ann, {Der[t], Der[x]}];
        ann = DFiniteOreAction[ann, ChangeOreAlgebra[kill, OreAlgebra[ann[[1]]]]];
        AppendTo[ops, (#**kill)& /@ CreativeTelescoping[ann, Der[x], {Der[t]}][[1]]];
      ]
    , {i, 1, Length[res]}];
    DFinitePlus@@ops
  ], 3600000, $Failed]]

(* 6th way: Equivalence + Residues + Chyzak-like *)
ctRat6[rat_, t_, x_, y_] := Timing[TimeConstrained[Module[ {num, den, res, z, r, f, op, out, kill, alg, one},
    {num, den} = {Numerator[#], Denominator[#]}&[Last[HermiteReduce[rat, y]]];
    res = FactorList[Resultant[num - z D[den, y], den, y]];
    out = one = Together[ToOrePolynomial[1, OreAlgebra[Der[t]]]];
    Do[
      If[ Exponent[res[[i, 1]], z] > 0, 
        kill = ResidueKiller[res[[i, 1]], t, x, z];
        alg = ApplyOperatorToAlgebraicFunction[kill, z, res[[i, 1]], t, x, z];
        op = If[ alg === 0, {one, 0}, ctAlgNew[alg, res[[i, 1]], t, x, z]];
        out = DFinitePlus[{out}, {op[[1]]**kill}][[1]];
      ]
    , {i, 1, Length[res]}];
    {out}
  ], 3600000, $Failed]]

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

ResidueKiller[minpoly_, t_, x_, y_] := Module[ {m, d, b, lc, kill, z, r, f, split, u, a, n, h, derY, derA, A},
   m = minpoly; d = Exponent[minpoly, y]; b = 1;
   (* make generator regular at finite places *)
   lc = Coefficient[m, y, d];
   kill = Together[lc/PolynomialGCD[lc, D[lc, x]]];
   While[ Exponent[kill, x] > 0,
      m = Numerator[Together[m /. y -> y/kill]];
      b = b*kill;
      kill = PolynomialGCD[kill, Coefficient[m, y, d]];
   ];
   (* determine the residues *)
   If[ Exponent[PolynomialGCD[b, D[b, x]], x] > 0, 
     split = Apart[y/b, x];
     If[ Head[split] === Plus, split = List@@split, split = {split} ];
     res = 1; derY = NormalizeAlgebraicFunction[-D[m, x]/D[m, y], m, t, x, y];
     Do[
       {a, b, n} = {Numerator[split[[i]]], Denominator[split[[i]]], 1};
       If[ Head[b] === Power, n = b[[2]]; b = First[b] ];
       (* handle a/b^n *) 
       (* begin untested code *)
       Print["handling multiplicity >1"];
       derA = {NormalizeAlgebraicFunction[D[a, x] + D[a, y]*derY, m, t, x, y]};
       While[Length[derA] < n - 1, 
         AppendTo[derA, NormalizeAlgebraicFunction[D[Last[derA], x] + D[Last[derA], y]*derY, m, t, x, y]]];
       h = Together[D[A[x]/u[x]^(n-1), {x, n - 1}]/(n - 1)!] 
             /. Derivative[i_][A][x] :> derA[[i]] /. A[x] -> a
             /. Derivative[i_][u][x] :> D[b, {x, i + 1}]/(i + 1) /. u[x] -> D[b, x]
             // Together;
       res = PolynomialLCM[res, Resultant[Resultant[Numerator[h] - z*Denominator[h], m, y], b, x]];
       (* end untested code *)
     , {i, 1, Length[split]}];
     res = Resultant[Resultant[y - z*D[b, x], m, y], b, x];
   ];
   res = Together[res / PolynomialGCD[res, D[res, z]]];
   If[FreeQ[res, z], Return[Annihilator[0, Der[t]][[1]]]];
   (* ae2de *)
   res = r[f[z, res], 1] /. {f -> Function, r -> Root};
   Annihilator[res, Der[t]][[1]]   
];

(* write rat = D[u, y] + v where v has a denominator which is squarefree wrt y *)
HermiteReduce[rat_, y_] := 
  Module[ {xgcd, quorem, dMinus, dStar, dMinus2, dMinusStar, gcd, g, a, b, c, d, q, r, u, v},

    xgcd = HolonomicFunctions`Private`MyPolynomialExtendedGCD;
    quorem = HolonomicFunctions`Private`MyPolynomialQuotientRemainder;
    
    {a, d} = {Numerator[#], Denominator[#]}&[Together[rat]];
    {g, a} = quorem[a, d, y]; g = Integrate[g, y];
    dMinus = xgcd[d, D[d, y], y][[1]]; dStar = Together[d/dMinus];

    While[ Exponent[dMinus, y] > 0, 

        dMinus2 = PolynomialGCD[dMinus, D[dMinus, y]];
        dMinusStar = Together[dMinus / dMinus2];

        {gcd, {b, c}} = xgcd[u = Together[-dStar * D[dMinus, y] / dMinus], v = dMinusStar, y];
        q = Together[a/gcd];
        {b, c} = Together[q*{b, c}];
        If[ Together[b] =!= 0 && Exponent[b, y] >= Exponent[v, y],
           {q, r} = quorem[b, v, y];
           {b, c} = {r, Together[c + q*u]};
        ];

        a = Together[c - D[b, y]*dStar/dMinusStar];
        g = Together[g + b/dMinus];
        dMinus = dMinus2
    ];

    Return[ {g, Together[a/dStar]} ];
];

(* INPUT
   op: an element of OreAlgebra[Der[t]]
   alg: an element of k(t,x)[y]/<m>
   m:  a polynomial in t, x, y
   t, x, y: variables
   OUTPUT
      a polynomial in y with rational function coefficients in t, x representing
      the algebraic function obtained by applying op to the algebraic function alg
      defined by m.
*)
ApplyOperatorToAlgebraicFunction[op_, alg_, m_, t_, x_, y_] :=
   Module[ {yt, coeffs, der, out},
      yt = NormalizeAlgebraicFunction[-D[m, t]/D[m, y], m, t, x, y];
      coeffs = CoefficientList[Normal[op], Der[t]];
      der = alg; out = coeffs[[1]]*alg;
      Do[
        der = NormalizeAlgebraicFunction[D[der, t] + D[der, y]*yt, m, t, x, y];
        out = NormalizeAlgebraicFunction[out + coeffs[[i]]*der, m, t, x, y];
      , {i, 2, Length[coeffs]}];
      Return[out];
   ];

(* INPUT
   alg: an element of k(t,x)[y]/<m> which has no residues
   m: minpoly
   t, x, y: variables
   OUTPUT:
   a pair (P, Q) in OreAlgebra[Der[t]] times k(t,x)[y]/<m> such that
   P(alg) = Der[x](Q)
*)
ctAlgNew[alg_, m_, t_, x_, y_] := Module[ {d, yt, yx, lhs, r, P, Q, p, q, sys, sol, vars},
   d = Exponent[m, y];
   yx = NormalizeAlgebraicFunction[-D[m, x]/D[m, y], m, t, x, y];
   Q = Sum[ D[q[i][x], x]y^i + q[i][x]*NormalizeAlgebraicFunction[i*y^(i-1)*yx, m, t, x, y], {i, 0, d - 1} ];
   P = {}; sol = {};
   yt = NormalizeAlgebraicFunction[-D[m, t]/D[m, y], m, t, x, y];
   While[ Length[sol] == 0,
     If[ Length[P] == 0, 
         P = {alg}, 
         AppendTo[P, NormalizeAlgebraicFunction[D[Last[P], t] + D[Last[P], y]*yt, m, t, x, y]] ];
     sys = Numerator[Together[CoefficientList[Q - Sum[p[i]P[[i]], {i, 1, Length[P]}], y]]];
     sol = SolveCoupledSystem[sys, Table[q[i], {i, 0, d-1}], {x}, ExtraParameters -> Array[p, Length[P]]];
     sol = {ToOrePolynomial[Sum[p[i]Der[t]^(i-1), {i, 1, Length[P]}], OreAlgebra[Der[t]]],
            Sum[q[i][x]y^i, {i, 0, d - 1}]} /. First[sol];
     vars = Union[Cases[sol, C[_], Infinity]];
     sol = sol /. Table[{vars[[i]] -> 1}, {i, 1, Length[vars]}] /. C[_] -> 0 // Together;
     sol = Select[sol, (First[#]=!=0*First[#])&];
   ];
   {p, q} = First[sol];
   {c, p} = NormalizeCoefficients[p, Extended -> True];
   {p, q/c}
];

NormalizeAlgebraicFunction[alg_, m_, t_, x_, y_] :=
  Module[ {num, den, g, u, v, out},
    {num, den} = {Numerator[#], Denominator[#]}&[Together[alg]];
    {g, {u, v}} = HolonomicFunctions`Private`MyPolynomialExtendedGCD[den, m, y];
    {num, den} = Together[{num*u, g}];
    {num, den} = {Numerator[num], den*Denominator[num]};
    num = Last[HolonomicFunctions`Private`MyPolynomialQuotientRemainder[num, m, y]];
    {num, den} = {Numerator[num], den*Denominator[num]};
    Return[num/den];
  ];

