
(* started by MK 2019-12-30 *)

(**
 * Input:
 * -- ideal : a list of bivariate polynomials over QQ generating an ideal of dimension 0.
 * -- x, y : the variables with respect to which the polynomials are given
 *
 * Output:
 * -- a list of generators of the the algebra of separated polynomials 
 **)
Separate0D[ideal_List, x_, y_] :=
  Module[ {xpure, ypure, terms, rems, G, vars, a, mixed},
    If[ PolynomialGCD@@ideal =!= 1, Throw["not zero dimensional"] ];
    {xpure, ypure} = First[GroebnerBasis[ideal, {First[#]}, {Last[#]}]]& /@ {{x,y},{y,x}};
    If[ xpure === 1, Return[{{1,0},{x,0},{0,1},{0,y}}] ];
    terms = Join[{0,#}& /@ Reverse[y^Range[0, Exponent[ypure,y]-1]], {#,0}& /@ Reverse[x^Range[0, Exponent[xpure,x]-1]]];
    vars = Array[a, Length[terms]]; G = GroebnerBasis[ideal, {x,y}];
    mixed = Transpose[terms].#& /@ NullSpace[Outer[Coefficient,
       Flatten[CoefficientList[(Last[PolynomialReduce[First[DeleteCases[{1,-1}*#, 0]], G, {x,y}]]&
          /@ terms).vars, {x,y}]], vars]];
    Return[Join[mixed, {#,0}& /@ (xpure x^Range[0, Exponent[xpure,x]-1]), {0,#}& /@ (ypure y^Range[0, Exponent[ypure,y]-1])]]
  ]

(**
 * Input:
 * -- A0 : a list of generators of the algebra of separated polynomials of a bivariate zero-dimensional ideal, as produced by the function Separate0D
 * -- x, y : the variables with respect to which the polynomials are given
 *
 * Output:
 * -- a linear function which applied to a pair of polynomials produces a vector (a finite list of numbers) that is zero if and only if the input pair belongs to A0 
 **)
MakeReductor[A0_, x_, y_] :=
  Module[ {p, q, B, i, f},
    p = First[Cases[A0, {p_, 0} :> p]]; q = First[Cases[A0, {0, q_} :> q]];
    B = (#/Last[CoefficientList[First[#], x]])& /@ Sort[DeleteCases[A0, {_, 0}|{0, _}], (Exponent[First[#1], x] > Exponent[First[#2], x])&];
    f[{u_, v_}] := Module[{u0, v0, t},
      u0 = PolynomialRemainder[u, p, x]; v0 = PolynomialRemainder[v, q, y];
      Do[{u0, v0} = Expand[{u0, v0} - Coefficient[u0, x, Exponent[First[b], x]]*b] + {t x^Exponent[First[b], x], 0}, {b, B}];
      Return[Join[DeleteCases[Most[CoefficientList[u0 + x^Exponent[p, x], x]], t],
                  Most[CoefficientList[v0 + y^Exponent[q, y], y]]]];
    ];
    Return[f];
  ];

(**
 * Input:
 * -- poly : a univariate polynomial
 * -- x : the variable in which the polynomial is stated
 *
 * Output:
 * -- a positive integer n such that poly = const*Cyclotomic[n, x], or -1 if no such n exists 
 **)
CyclotomicQ[poly_, x_] :=
  Module[ {phin, bound, n},
    phin = Exponent[poly, x];
    If[ Coefficient[poly, x, phin] =!= 1, Return[CyclotomicQ[poly/Coefficient[poly, x, phin], x] ]];
    If[ Expand[poly - (x-1)] === 0, Return[1]];
    If[ Coefficient[poly, x, 0] =!= 1, Return[-1]];
    For[n = 0, n < 3 || n/(Exp[EulerGamma]Log[Log[n]] + 3/Log[Log[n]]) <= phin, n += 1, 
      If[ EulerPhi[n] == phin && Expand[Cyclotomic[n, x] - poly]===0, Return[n] ]];
    Return[-1];
  ];

(**
 * Input:
 * -- f : a bivariate polynomial
 * -- x, y : the two variables in which f is stated
 *
 * Output:
 * -- a list of generators of the algebra of separated polynomials of the ideal generated by f 
 **)
SeparatePrincipal[f_, x_, y_] :=
  Module[ {h, t, n, alpha, index, a, b, i, j, sol, p, q},
    Which[ Head[f] === List && Length[f] == 1, Return[SeparatePrincipal[First[f], x, y]], 
      Expand[f] === 0, Return[ {{1,1}} ],
      Expand[f] === 1, Return[ {{1,0},{x,0},{0,1},{0,y}} ],
      FreeQ[Expand[f], x], Return[ Prepend[{0,y^# f}& /@ Range[0, Exponent[f,y]-1], {1,1}] ],
      FreeQ[Expand[f], y], Return[ Prepend[{x^# f,0}& /@ Range[0, Exponent[f,x]-1], {1,1}] ]
    ];
    p = Exponent[f /. y -> 0, x]; q = Exponent[f /. x -> 0, y]; (* f contains x^p and y^q *)
    If[ !IntegerQ[p] || !IntegerQ[q], Return[ {{1,1}} ] ]; (* f has x or y as factor *)
    n = Exponent[h = Expand[Last[CoefficientList[f /. {x -> x t^q, y -> y t^p}, t]]], x]; (* omega(x^i*y^j) = q*i + p*j *)
    If[ Expand[(h/.y->0)*(h/.x->0)] === 0, Return[ {{1,1}} ] ]; (* nontrivial Newton polygon *)
    h = Expand[h/Coefficient[h, y, q]]; alpha = ToNumberField[Coefficient[h, x, p]^(1/p)]; (* normalize *)
    If[ Length[Complement[Variables[h], {x,y}]] > 0, Return[ {{1,1}}] ]; (* no parameters beyond this point *)
    h = CyclotomicQ[#, x]& /@ (MinimalPolynomial[Root[Function[x, First[#]], 1], x]&) /@
       Select[FactorList[h /. {x -> x/alpha, y -> 1}, Extension -> alpha], !FreeQ[#, x]&];
    If[ MemberQ[h, -1], Return[{{1, 1}}]];
    sol = {Null, Null}; n = LCM@@h + 1; (* now n bounds the order of the roots of unity *)
    While[ Length[sol] > 1,
      n -= 1; vars = Join[Table[a[i], {i,0,n}], Table[b[i],{i,1,n*q/p}]];
      sol = Expand[Join[x^Range[0,n], y^Range[n*q/p]].#]& /@
        NullSpace[Outer[Coefficient, Flatten[CoefficientList[PolynomialRemainder[vars.Join[x^Range[0,n],y^Range[n*q/p]], f, x], {x,y}]], vars]];
    ];
    Return[Prepend[Expand[{#/.y->0,(#/.y->0)-#}]& /@ sol, {1,1}]]
  ];

(**
 * Input:
 * -- ideal : a list of bivariate polynomials
 * -- x, y : the two variables in which the ideal generators are stated
 *
 * Output:
 * -- a list of generators of the algebra of separated polynomials of the input ideal
 **)
Separate[ideal_, x_, y_] :=
  Module[ {id, A0, A1, f, d, a, t, g, G, Delta, S, s, p},
    A1 = SeparatePrincipal[id[1] = PolynomialGCD@@ideal, x, y];
    A0 = Separate0D[id[0] = GroebnerBasis[Together[ideal/id[1]], {x,y}], x, y]; 
    Which[ id[1] === 1, Return[A0], id[0] === {1} || A1 === {{1,1}}, Return[A1],
    FreeQ[A1, y], p = First[GroebnerBasis[ideal, {x}, {y}]]; Return[{x^# p, 0}& /@ Range[0, Exponent[p, x] - 1]],
    FreeQ[A1, x], p = First[GroebnerBasis[ideal, {y}, {x}]]; Return[{0, y^# p}& /@ Range[0, Exponent[p, y] - 1]]
    ];
    f = MakeReductor[A0, x, y]; d = Length[f[{1,1}]]; a = A1[[2]]; G = {{1,1}}; g = 0; Delta = {}; 
    While[ True, 
      Which[
        Length[Delta] == 0, S = Range[d+1],
        Length[Delta] == 1, S = Complement[Range[g*(d+1)], g*Range[d+1]],
        g != 1, S = Select[Range[g*(d+1)], Length[FrobeniusSolve[Delta, #, 1]]===0&],
        True, S = Select[Range[FrobeniusNumber[Delta]], Length[FrobeniusSolve[Delta, #, 1]]===0&];
      ];
      If[ Length[S] > d, S = Take[S, d+1] ];
      If[ Length[S] == 0, Break[] ];
      p = NullSpace[Transpose[Table[f[a^s], {s, S}]]];
      If[ Length[p] == 0, Break[], p = (#*LCM@@(Denominator/@#))&[Last[p]] ];
      AppendTo[G, Expand[Sum[p[[i]] a^S[[i]], {i, 1, Length[S]}]]];
      AppendTo[Delta, Exponent[p.(t^S), t]];
      g = GCD[g, Last[Delta]]; 
    ];
    Return[G];    
  ];

(* the functions below are test functions *)

CheckSeparate0D[ideal_List, x_, y_] :=
  Module[ {A0, G, u, v, a, b, dx},
    A0 = Separate0D[ideal, x, y];
    G = GroebnerBasis[ideal, {x,y}];
    If[ !MatchQ[Last[PolynomialReduce[#, G, {x,y}]]& /@ (A0 /. {u_, v_} :> u - v), {0..}],
      Throw["incorrect!"];
    ];
  ];

CheckSeparatePrincipal[f_, x_, y_] :=
  Module[ {A1, dx, dy},
    A1 = SeparatePrincipal[f, x, y];
    If[ Length[A1] > 1 && !FreeQ[Denominator[Together[(A1[[2, 1]] - A1[[2, 2]])/f]], x|y],
      Throw["incorrect!"];
    ];
    dx = If[ Length[A1] == 1, 10, Exponent[A1[[2,1]], x] - 1];
    dy = If[ Length[A1] == 1, 10, Exponent[A1[[2,2]], y] - 1];
    If[ Length[DeleteCases[First[Solve[Flatten[CoefficientList[Last[PolynomialReduce[Sum[a[i]x^i, {i,0,dx}] + Sum[b[i]y^i, {i,1,dy}], {f}, {x,y}]], {x,y}]]==0]], _ -> 0]] > 0,
      Throw["incomplete!"];
    ];
  ];

CheckSeparate[ideal_, x_, y_] := 
  Module[ {A, dx, dy, u, v},
    A = Separate[ideal, x, y]; G = GroebnerBasis[ideal, {x,y}];
    If[ !MatchQ[Last[PolynomialReduce[#, G, {x,y}]]& /@ (A /. {u_, v_} :> u - v), {0..}],
      Throw["incorrect!"];
    ];
    dx = If[ Length[A] == 1, 10, Exponent[A[[2,1]], x] - 1];
    dy = If[ Length[A] == 1, 10, Exponent[A[[2,2]], y] - 1];
    If[ Length[DeleteCases[First[Solve[Flatten[CoefficientList[Last[PolynomialReduce[Sum[a[i]x^i, {i,0,dx}] + Sum[b[i]y^i, {i,1,dy}], {f}, {x,y}]], {x,y}]]==0]], _ -> 0]] > 0,
      Throw["incomplete!"];
    ];
  ];
