(* by M. Kauers *)

<< Singular.m;
<< NumberTheory`AlgebraicNumberFields`;

(**
 * Polynomial manipulation tools
 * by Manuel Kauers
 **)

Canonic[rat_] := Cancel[Together[rat, Extension -> Automatic], Extension -> Automatic];

(**
 * Quotient and remainder of univariat polynomials over rational functions
 **)
QuoRem[p_, q_, x_] :=
  Module[ {quo, rem},
    {quo, rem} = Developer`PolynomialDivision[p, q, x]; 
    Return[{quo, rem}];
  ];
NumDen[rat_, x_] :=
  Module[ {a, d}, 
    a = Canonic[rat];
    {a, d} = {Numerator[a], Denominator[a]};
    Return[Cancel[{a, d}/LeadingCoeff[d, x]]];
  ];

(**
 * returns the primitive part of a polynomial over a rational function field 
 **)
ClearContent[p_, x_] := 
  Module[ {p0}, 
    p0 = Numerator[Canonic[p]];
    If[ p0 === 0, Return[0] ];
    p0 = Cancel[p0/PolynomialGCD@@Append[CoefficientList[p0, x], Extension -> Automatic]];
    Return[p0]
  ];  

(**
 * Polynomial gcd with rational function coefficients 
 **)
UnivariatePolynomialGCD[p_, q_, x_] :=
  Module[ {a, b, g},
    a = Numerator[Canonic[p]];
    If[ a === 0, Return[ q ]];
    a = ClearContent[a, x];
    b = Numerator[Canonic[q]];
    If[ b === 0, Return[ a ]];
    b = ClearContent[b, x];
    g = PolynomialGCD[a, b, Extension -> Automatic];
    Return[g];
  ]; 

UnivariatePolynomialGCD[p_List, x_] := 
     Which[
       Length[p] == 1, 
         Return[First[p]],
       Length[p] == 2,
         Return[UnivariatePolynomialGCD[First[p], Last[p], x]],
       Length[p] >= 3,
         Return[UnivariatePolynomialGCD[First[p], UnivariatePolynomialGCD[Rest[p], x], x]]
       ];

UnivariatePolynomialLCM[p_, q_, x_] :=
   PolynomialQuotient[p * q, UnivariatePolynomialGCD[p, q, x], x];

UnivariatePolynomialLCM[p_List, x_] :=
     Which[
       Length[p] == 1, 
         Return[First[p]],
       Length[p] == 2,
         Return[UnivariatePolynomialLCM[First[p], Last[p], x]],
       Length[p] >= 3,
         Return[UnivariatePolynomialLCM[First[p], UnivariatePolynomialLCM[Rest[p], x], x]]
       ];

(** 
 * Extraction of leading coefficient 
 **)
LeadingCoefficient = LeadingCoeff;
LeadingCoeff[p_, x_] := 
  Module[ {p0}, 
    p0 = Canonic[p];
    Return[If[ p0 === 0, 0, Coefficient[p0, x, Exponent[p0, x]] ]]
  ];

(**
 * Cofactor computation (Bronstein p. 14)
 *
 * INPUT:
 * polynomials a, b, c in x
 *
 * OUTPUT:
 * s, t such that s a + t b = c
 **)
ExtendedEuclidean[a_, b_, c_, x_] :=
  Module[ {s, t, r},

    s = HalfExtendedEuclidean[a, b, c, x];
    {t, r} = QuoRem[Canonic[c - s*a], b, x];
    If[ Canonic[r] =!= 0, Throw["ExtendedEuclidean: r=!=0"]];
    Return[Canonic[{s, t}]];
  ];
HalfExtendedEuclidean[a_, b_, c_, x_] :=
  Module[ {s, g, q, r},
    {s, g} = HalfExtendedEuclidean[a, b, x];
    {q, r} = QuoRem[c, g, x];
    If[ Canonic[r] =!= 0, Throw["HalfExtendedEuclidean: r=!=0"]];
    s = q s;
    If[ Canonic[s] =!= 0 && Exponent[s, x] >= Exponent[b, x],
      s = PolynomialRemainder[s, b, x, Extension -> Automatic];
    ];
    Return[s];
  ];
HalfExtendedEuclidean[a_, b_, x_] :=
  Module[ {a0, b0, a1, b1, q, r, g},
    {a0, b0} = {a, b};
    {a1, b1} = {1, 0};
    While[ b0 =!= 0, 
      {q, r} = QuoRem[a0, b0, x];
      {a0, b0} = {b0, r};
      {a1, b1} = {b1, Canonic[a1 - q b1]};
    ];
    (* clear 'common content' of a1 and a0 for efficiency *)
    {a1, a0} = Canonic[{a1, a0}];
    {a1, a0} = {a1, a0} * PolynomialLCM[Denominator[a1], Denominator[a0]];
    g = PolynomialGCD@@Union[CoefficientList[a1, x], CoefficientList[a0, x]];
    {a1, a0} = Cancel[{a1, a0}/g];
    (* done *)
    Return[{a1, a0}];
  ];

(**
 *
 **)
RationalOrder[f_, p_, x_] :=
  Module[ {num, den, n, d, r},

    {num, den} = NumDen[f, x];

    If[ p === Infinity, 
      Return[If[ num === 0, Infinity, Exponent[den, x] - Exponent[num, x] ]]
    ];

    If[ num === 0, n = Infinity,
      n = 0; {num, r} = QuoRem[num, p, x];
      While[ Canonic[r] === 0, {num, r} = QuoRem[num, p, x]; n++]
    ];    

    If[ den === 0, d = Infinity,
      d = 0; {den, r} = QuoRem[den, p, x];
      While[ Canonic[r] === 0, {den, r} = QuoRem[den, p, x]; d++]
    ]; 

    Return[ n - d ];
  ];

(**
 * INPUT: f: in k(x,y)
 *        x, y: field generators
 *        p: minimal polynomial of y over k(x)
 *
 * OUTPUT: Derivative of f in k(x,y) with respect to x.
 **)
ClearAll[Diff];
Diff::usage = "Diff[f,x,y,p] differentiates an algebraic function f in k(x)[y]/<p> with respect to x. Example: Diff[x^2+y,x,y,y^2-(x^2+1)] gives 2x + x/y.";
Diff[f_, x_, y_, _] /; FreeQ[f, x]&&FreeQ[f, y] := 0;
Diff[x_, x_, _, _] = 1;
Diff[y_, x_, y_, p_] := -D[p, x]/D[p, y];
Diff[a_Plus, x_, y_, p_] := Diff[#, x, y, p]& /@ a;
Diff[a_List, x_, y_, p_] := Diff[#, x, y, p]& /@ a;
Diff[a_Times, x_, y_, p_] := Canonic[a*Plus@@((Diff[#, x, y, p]/#)& /@ List@@a)];
Diff[a_^i_Integer, x_, y_, p_] := i*a^(i-1)*Diff[a, x, y, p];
Diff[Log[a_], x_, y_, p_] := Diff[a, x, y, p]/a;
Diff[RootSum[u_, Function[v_]], x_, y_, p_] := 
  Module[ {r, f}, r[u, f[Diff[v, x, y, p]]] /. {r -> RootSum, f -> Function}];

(**
 * INPUT: f: in k(x,y),
 *        x, y: field generators,
 *        p: minimal polynomial of y over k(x),
 *        w: a basis for k(x,y) over k(x)
 * 
 * OUTPUT: {{a[1],...,a[n]}, d} such that f == (a[1]w[1]+...+a[n]w[n])/d
 **)
CanonicAlg[f_, x_, y_, p_, w_] :=
  Module[ {num, den, a, sys, gcd, lcm},

    {num, den} = Last /@ (QuoRem[#, p, y]&) /@ NumDen[f, y];
    num *= First[ExtendedEuclidean[den, p, 1, y]];
    {num, den} = NumDen[Last[QuoRem[num, p, y]], x];
    sys = CoefficientList[num - w . Array[a, Length[w]], y];
    a = Array[a, Length[w]]/.First[Solve[Thread[sys==0],Array[a,Length[w]]]];
    a = Canonic /@ a;
    lcm = UnivariatePolynomialLCM[Denominator /@ a, x];
    {a, den} = {a, den}*lcm;
    gcd = UnivariatePolynomialGCD[Append[a, den], x];
    a = PolynomialQuotient[#,gcd,x]& /@ a;
    Return[ {a, PolynomialQuotient[den,gcd,x]} ];
  ];

CanonicAlg[f_, x_, y_, p_] := 
  Factor[First[#].y^Range[0,Exponent[p,y]-1]/Last[#]]&[CanonicAlg[f,x,y,p,y^Range[0,Exponent[p,y]-1]]];
(**
 * INPUT: f in k(x,y)
 *
 * OUTPUT: p in k[x,t] such that p(x,f) = 0
 **)
MinPoly[f_, x_, y_, p_] :=
  Module[ {a, d, t, m, w},

    w = y^Range[0, Exponent[p, y] - 1];
    {a, d} = CanonicAlg[f, x, y, p, w];
    m = Resultant[a.w - t*d, Numerator[Together[p]], y];
    m = Select[FactorList[m], !FreeQ[#,t]&];
    m = Times@@(m /. {a_, b_Integer} -> a);
    Return[m /. t -> y];    
  ];

(**
 * INPUT: f in k(x,y)
 *        x0 in k or x0 == Infinity
 * OUTPUT: nu_x0(f): the order of f at all places over x0
 **)
Orders[f_, x0_, x_, y_, p_] :=
  Module[ {g, f0, p0, p1, alpha},
    Which[ 
      x0 === Infinity, Return[ Orders[f /. x -> 1/x, 0, x, y, p /. x -> 1/x] ],
      x0 =!= 0, Return[ Orders[f /. x -> x + x0, 0, x, y, p /. x -> x + x0 ] ],
      FreeQ[f, y] && PolynomialQ[f, x],
         f0 = Expand[f];
         If[ f0 === 0, Return[{Infinity}] ];
         f0 = CoefficientList[f0, x];
         Return[{Min[Cancel[f0/(f0 /. 0 -> 1)]*Range[Length[f0]] /. 0 -> Infinity] - 1}], 
      True,
         p0 = First[Orders[#, 0, x, y, p]]& /@ CoefficientList[MinPoly[f, x, y, p], y];
         p1 = Select[p0 + alpha*Range[0, Length[p0] - 1], FreeQ[#, Infinity]&];
         p0 = Apply[Equal, Subsets[p1, {2}], {1}];
         p0 = alpha /. (Solve[#, alpha]& /@ p0) // Flatten;
         p0 = Select[p0, (Count[p1 /. alpha->#, Min[p1 /. alpha -> #]] > 1)&];
         Return[Union[p0]];
      ]];

(**
 * INPUT: x, y: variables
 *        p: minimal polynomial of y over k(x), of the form y^n - rat(x)
 *
 * OUTPUT: an integral basis for k[x,y] over k[x]. An error is generated 
 *        if y is not a simple radical extension over k(x).
 **)
IntegralBasisRadicals[x_, y_, p_] :=
  Module[ {n, w, A, D, F, H, p0, factors, a, i, j},

    p0 = Numerator[Canonic[p]];
    n = Exponent[p0, y];

    D = Coefficient[p0, y, n];
    A = -Canonic[p0 - D*y^n];

    If[ !FreeQ[A, y], 
      A = y^Range[0, n - 1]; 
      Print["Warning: Assuming integral basis " <> ToString[InputForm[A]]];
      Return[A];
    ];

    factors = FactorSquareFreeList[A*D^(n - 1)];
    F = Times@@(factors /. {a_, k_Integer} -> a^Quotient[k, n]);
    H = Times@@(factors /. {a_, k_Integer} -> a^Mod[k, n]);

    factors = FactorSquareFreeList[H];

    Table[(y*D/F)^(i-1)/Times@@(factors /. {a_, j_Integer} -> a^Floor[(i-1)j/n]), {i, 1, n}](*;*)
  ];

(**
 * INPUT: x, y: variables
 *        p: polynomial in x and y
 * 
 * OUTPUT: discriminant of k(x,y) over k(x).
 **)
Discriminant[x_, y_, p_, w_] :=
  Module[ {sp, e, d},

    sp[e_] := sp[e]=Sum[(#[[1,i]]/#[[2]])&[CanonicAlg[e*w[[i]], x, y, p, w]], {i, 1, Length[w]}];
    Return[Det[Map[sp, Outer[Times, w, w], {2}]]]
  ];

(**
 * Solves the linear equation system A.u==b in k[x]/<p>.
 * It is assumed that the system has a unique solution. If this is not the case,
 * an exception is thrown. 
 **)
LinearSolveModPoly[A_, b_, x_, p_] :=
  Module[ {det, At, i},

    (* cramer's rule *)
    det = PolynomialRemainder[Det[A], p, x];
    det = First[ExtendedEuclidean[det, p, 1, x]]; (* == 1/det *)
    At = Transpose[A];

    Table[
      PolynomialRemainder[det*Det[ReplacePart[At, b, i]], p, x]    
    , {i, 1, Length[b]}](*;*)
  ];

(**
 * performs Hermite reduction for the given integrand f 
 * 
 * INPUT: a, d: polynomials in x
 *        x: integration variable, 
 *        y: algebraic over x, 
 *        p: minimal polynomial of y over k(x), monic in y,
 *        w: integral basis for k(x,y) over k[x].
 *
 * The function f:=(a[1]w[1]+..+a[n]w[n])/d must not have a pole or a 
 * branch point at infinity.
 *
 * OUTPUT: g, h in k(x,y) such that 
 * int(f) = g + int(h) and h has a square free denominator.
 **)
HermiteReduce[a_List, d_, x_, y_, p_, w_] :=
  Module[ {n, b, i, j, lcm, factors, U, V, E, D, A, B, M, k, sys, sol},

    n = Length[w];
    E = PolynomialLCM@@(Last /@ (CanonicAlg[#, x, y, p, w]&) /@ (Diff[#, x, y, p]&) /@ w);
    M = Table[First[CanonicAlg[E*Diff[w[[i]], x, y, p], x, y, p, w]], {i, 1, n}];

    G = 0; (* accumulated algebraic part of the integral *)

    lcm = UnivariatePolynomialLCM[d, E, x]; (* ensure E | D *)
    A = PolynomialQuotient[lcm, d, x]*a; D = lcm;
    factors = FactorSquareFreeList[D]; 
    k = Max[Last /@ factors] - 1; (* maximum exponent *)

    While[ k > 0,

      V = Times@@Cases[factors, {i_, k+1} -> i];
      U = Times@@Cases[factors, ({i_, j_} /; j <= k) -> i^j];
      T = PolynomialQuotient[U*V, E, x];

      (* now D = U*V^(k+1) with gcd(V, V') = 1, gcd(U, V) = 1, E | U*V, E*T == U*V *)
      sys = T*M - DiagonalMatrix[Table[k*U*Diff[V, x, y, p], {n}]];
      B = LinearSolveModPoly[sys, A, x, V];

      G += (B.w)/V^k;
      {A, D} = CanonicAlg[(A.w)/(U*V^(k+1))-Diff[(B.w)/V^k,x,y,p], x,y,p,w];

      (* prepare for the next iteration *)
      lcm = UnivariatePolynomialLCM[D, E, x];
      A = PolynomialQuotient[lcm, D, x]*A;     
      D = lcm;
      factors = FactorSquareFreeList[D];
      k = Max[Last /@ factors] - 1;
    ];

    Return[ {G, (A.w)/D} ];
  ];

(***** logarithmic part via groebner bases *****)

LogarithmicPart[num_, den_, u_List, v_List] := Module[ {G, t, factors, f, r, i, F, p},
  G = SingularGroebner[Join[{den, num - t D[den, First[v]]}, u], v, {t}, MonomialOrder->"dp"];
  factors = DeleteCases[Rest[FactorList[First[G]]], {t, _}];
  F[p_, g_] := (t/#2*Log[#1])&@@PrincipalPower[g, Append[u, p], v, {t}];
  Plus@@Apply[r[f[#1], f[#2 F[#1, G]]]&, factors, {1}] /. {t -> #, f -> Function, r -> RootSum}];

PrincipalDivisor[gb_, u_, v__] := Module[{G},
  G = SingularGroebner[Join[gb, u], v, MonomialOrder->"dp"];
  First[Append[Select[G, (SingularGroebner[Append[u, #], v, MonomialOrder->"dp"]===G)&, 1], 1]]];

PrincipalPower[gb_, u_, v__List, bound_Integer:10] := Module[ {id = gb, p, n = 1},
  While[n < bound && (p = PrincipalDivisor[id, u, v])===1, n++; Print[n];
    id = SingularGroebner[Join[SingularTimes[id, gb, v], u], v, MonomialOrder->"dp"]]; Return[{p, n}]];
(**
 * Integrate algebraic functions
 * 
 * INPUT: f in k(x,y)
 *        x,y: field generators
 *        p: minimal polynomial of y over k(x).
 *
 * OUTPUT: {g, h} such that int(f, x) = g + int(h, x) 
 *         and h==0 iff int(f,x) is elementary
 **)
IntegrateAlgebraic::usage = "IntegrateAlgebraic[f,x,y,p] produces expressions g and h such that
f and Diff[g,x,y,p] + h are identical when considered as elements of k(x)[y]/<p>. Example: IntegrateAlgebraic[x^2+y,x,y,y^2-(x^2+1)] gives {x^3/3 + (x*y)/2 - Log[-x + y]/4 + Log[x + y]/4, 0}";
IntegrateAlgebraic[f_, x_, y_, p_] :=
  Module[ {a, b, c, ord, z, w, w0, d, n, f0, p0, x0, y0, den, dis, g, h, am, rad, vars, id, r, ff},

    d = Exponent[p, y];

    (* substitute poles and branch points at infinity away, if there are some. *)
    ord = Orders[(f /. x -> 1/x)*(-1/x^2), 0, x, y, p /. x -> 1/x];
    w = y^Range[0, d - 1];
    {a, den} = CanonicAlg[f, x, y, p, w]; 

    If[ Min[ord] < 0 || !FreeQ[ord, Rational],

      dis = Times@@(First/@FactorSquareFreeList[Discriminant[x, y, p, w]]);
      dis = Times@@NumDen[dis*den, x];
      For[n = 0, Canonic[dis /. x -> n]===0, n++];
      Print["substituting x=", n, " to infinity..."];
      {f0, p0} = ({f, p} /. x -> n + 1/x)*{-1/x^2, 1};
      p0 = ClearContent[Numerator[Canonic[p0]], y];
      am = Coefficient[p0, y, d];
      {f0, p0} = {f0, p0} /. y -> y/am;
      p0 = ClearContent[Numerator[Canonic[p0]], y];
      rad = y/IntegralBasisRadicals[x, y, p0][[2]];
      {f0, p0} = {f0, p0} /. y -> rad*y;
      p0 = ClearContent[p0, y];
      {g, h} = IntegrateAlgebraic[f0, x, y, p0];
      {g, h} = {g, h*(-x^2)} /. y -> y*am/rad /. x -> 1/(x-n);
      g = Canonic[g /. Log[_] -> 0] + g - (g /. Log[_] -> 0) /. Log[z_] :> Log[Canonic[z]]
             /. (Log[a_ b_.] /; FreeQ[a, x]&&FreeQ[a, y]) -> Log[b]
             /. a_?NumberQ*b_Plus :> Distribute[a*b]
             //.  {c_. Log[a_ b_] -> c Log[a] + c Log[b], Log[a_^b_Integer] -> b Log[a]} 
             /. {RootSum -> r, Function -> ff} /. {r -> RootSum, ff -> Function};
      Return[{g, Canonic[h]}];
    ];

    (* compute an integral basis for k[x,y] over k[x] which is normal at infinity *)
    w = IntegralBasisRadicals[x, y, p];

    (* hermite reduction *)
    {a, den} = CanonicAlg[f, x, y, p, w];
    {g, h} = HermiteReduce[a, den, x, y, p, w];

    (* early termination *)
    h = Canonic[h];
    If[ h === 0 || Min[Orders[(h /. x ->1/x)*(-1/x^2), 0, x, y, p /. x -> 1/x]] < 0,
      Print["early termination."];
      Return[ {g, h} ];
    ];    

    id = {p}; vars = {x, y}; w0 = w;
    Do[
      If[ Denominator[w0[[i]]] === 1,
        y0[i] = w0[[i]]
      ,
        id = Join[id, {Numerator[Together[y0[i] - w[[i]]]], 
                       ClearContent[MinPoly[w[[i]], x, y, p], y] /. y -> y0[i]}];
        AppendTo[vars, y0[i]];
        Do[w0[[j]] = w0[[j]]*y0[i]/w[[i]], {j, i, Length[w0]}];
      ]
    , {i, 1, Length[w]}];

    h = CanonicAlg[h, x, y, p, w];
    f0 = {First[h].w0, Last[h]};

    (* logarithmic part *)
    g += LogarithmicPart[First[f0], Last[f0], id, vars] /. y0[i_] :> w[[i]];
    Return[{g, 
        Canonic[(#1.y^Range[0,d-1])/#2&@@CanonicAlg[f - Diff[g, x, y, p], x, y, p, y^Range[0, d-1]]]}
             /. (Log[a_ b_.] /; FreeQ[a, x]&&FreeQ[a, y]) -> Log[b]
             /. a_?NumberQ*b_Plus :> Distribute[a*b]
             //.  {c_. Log[a_ b_] -> c Log[a] + c Log[b], Log[a_^b_Integer] -> b Log[a]} 
        ];
  ];