BeginPackage["RCF`"];

PNF::usage = "PNF[R, x] computes the strict polynomial normal form (a, b, c) of R as a rational function of x.";

RNF::usage = "RNF[R, x] computes a strict rational normal form (F, V) of R as a rational function of x.";

RCF1::usage = "RCF1[R, x] computes the first rational canonical form (F, V) of R as a rational function of x, minimizing deg den V.";

RCF2::usage = "RCF2[R, x] computes the second rational canonical form (F, V) of R as a rational function of x, minimizing deg num V.";

RCF3::usage = "RCF3[R, x] computes the third rational canonical form (F, V) of R as a rational function of x, minimizing first deg num V + deg den V, then deg den V.";

RCF4::usage = "RCF4[R, x] computes the fourth rational canonical form (F, V) of R as a rational function of x, minimizing first deg num V + deg den V, then deg num V.";

SHF::usage = "SHF[R, x] computes the shift-homogeneous factorization of R as a rational function of x over the rational number field.";


Begin["`Private`"];

Off[General::spell1];

MyResultant[p_ q_, r_, x_] := MyResultant[p, r, x]MyResultant[q, r, x];
MyResultant[r_, p_ q_, x_] := (-1)^(Exponent[r,
 x](Exponent[p, x] + Exponent[q, x]))MyResultant[p q, r, x];
MyResultant[p_, q_, x_] := Resultant[p, q, x];

PNF[R_, x_Symbol] := Module[{a =
 Factor[Numerator[R]], b = Factor[Denominator[R]], c, h, rez, nicle, s},
      rez = MyResultant[a, b /. x -> x + h, x];
      nicle = h /. Select[Solve[rez == 0, h], ListQ];
      If[nicle != {}, nicle = Union[Select[nicle, IntegerQ[#] && # >= 0&]]];
      For[c = 1, nicle != {}, nicle = Rest[nicle], h = First[nicle];
        s = PolynomialGCD[a, b /. x -> x + h];
        a = Cancel[a/s];
        b = Cancel[b/(s /. x -> x - h)];
        c = c Product[s /. x -> x - j, {j, h}]];
      Return[{a, b, Factor[c]}]];

RNF[R_, x_Symbol] := RCF1[R, x];

RCF1[R_, x_Symbol] := Module[{a, b, c, r, s, d},
      {a, b, c} = PNF[R, x];
      {s, r, d} = PNF[b/a, x];
      Return[{r/s, Cancel[c/d]}]];

RCF2[R_, x_Symbol] := Module[{F, V},
      {F, V} = RCF1[1/R, x];
      Return[{1/F, 1/V}]];

lCoeff[p_, x_] := Coefficient[p, x, Exponent[p, x]];

group[{}, _] = {};
group[{a_, b___}, x_] := Module[{d = Exponent[
      First[a], x], 
        ac, agroup = Select[{a, b}, Exponent[First[#], x] == Exponent[First[
    a], x] &], h},
      ac = Coefficient[First[a], x, d - 1];
      h = (Coefficient[#[[1]], x, d - 1] - ac)/d& /@ agroup;
      agroup = Select[Thread[{agroup, h}], 
    IntegerQ[#[[2]]] && Expand[#[[1, 1]] - (First[a] /. x -> x + #[[2]])] ==
         0&];
      Prepend[
          group[Complement[{b}, First /@ 
            agroup], x], {First[a], Thread[{Last /@ agroup, 
                Last /@ (First /@ agroup)}]}]];

SHF[r_, x_] := Module[{flist = FactorList[r], lc, lclist},
      lc = Select[flist, FreeQ[#, x] &];
      flist = Complement[flist, lc];
      lc = Times @@ (Power @@ # & /@ lc);
      lclist = lCoeff[#, x] & /@ (First /@ flist);
      lc = lc(Times @@ (lclist^(Last /@ flist)));
      flist = Thread[{Expand[(First /@ flist)/lclist], Last /@ flist}];
      Prepend[group[flist, x], lc]
      ];

mshRCF3[{p_, a_, b_, x_}] := Module[{m = Length[
      a], n = Length[b], Nb, Ma, w, y, vars, rules, rangef, F, V},
      Which[m < n, (* Case 1 *)
        Nb = (Plus @@ b) + 1;
        w =
     Table[If[r <= m, Abs[a[[r]] - b[[s]]], 1 - b[[s]]/Nb], {r, n}, {s, n}];
        vars = Array[y, {n, n}];
        rules = Thread[(Join @@ vars) -> LinearProgramming[Join @@ w, 
      Table[If[i <= n, If[(i - 1)n < j <= i n, 1, 0], If[Mod[i - j, n] == 0, 1, 
        0]], {i, 2n}, {j, n^2}], Table[{1, 0}, {2n}]]];
        (* rules = ConstrainedMin[{Plus @@ Plus @@ (w vars),
                    Join[
                    Thread[Plus @@ vars == 1], Thread[Plus @@ Thread[
                      vars] == 1]]}, Join @@ vars][[2]]; *)
        rules = Select[rules, #[[2]] == 1&] /. (y[i_, j_] -> 1) -> (i -> j);
        rangef = Range[m] /. rules;
        F = 1/(Times @@ (p /. 
          x -> x + b[[#]] & /@ Complement[Range[n], rangef]));
        V = RCF1[Product[(p /. x -> x + a[[k]])/(p /. x -> x + 
                  b[[k /. rules]]), {k, m}], x][[2]];
        Return[{F, V}],
        m >= n, (* Case 2 *)
        Ma = (Plus @@ a) + 1;
        w = Table[If[s <= n, Abs[a[[r]] - b[[
        s]]], a[[r]]/Ma], {r, m}, {s, m}];
        vars = Array[y, {m, m}];
        rules = Thread[(Join @@ vars) -> 
        LinearProgramming[
          Join @@ w, Table[If[i <= m, If[(i - 
              1)m < j <= i m, 1, 0], If[Mod[i - j, m] == 0, 
              1, 0]], {i, 2m}, {j, m^2}], Table[{1, 0}, {2m}]]];
        (* rules = 
        ConstrainedMin[Plus @@ Plus @@ (w vars), Join[Thread[Plus @@ vars == 
        1], Thread[Plus @@ Thread[vars] == 1]], Join @@ vars][[2]] *);
        rules = Select[rules, #[[2]] == 1&] /. (y[i_, j_] -> 1) -> (j -> i);
        rangef = Range[n] /. rules;
        F = 
      Times @@ (p /. x -> x + a[[#]] & /@ Complement[Range[m], rangef]);
        V = RCF1[Product[(p /. x -> x + 
a[[k /. rules]])/(p /. x -> x + b[[k]]), {k, n}], x][[2]];
        Return[{F, V}]
        ]
      ];

prepare[{p_, hlist_}] := Module[{pos = 
      Select[hlist, #[[2]] > 0 &], neg = Select[hlist, #[[2]] < 0 &]},
      pos = Join @@ (Table[#[[1]], {#[[2]]}] & /@ pos);
      neg = Join @@ (Table[#[[1]], {-#[[2]]}] & /@ neg);
      Return[{p, pos, neg}]];

RCF3[R_, x_Symbol] := Module[{shf = SHF[R, x], z, classes, FVlist},
      z = First[shf];
      classes = Rest[shf];
      FVlist = mshRCF3[Append[prepare[#], x]] & /@ classes;
      Return[Times @@ Prepend[FVlist, {z, 1}]]
      ];

RCF4[R_, x_Symbol] := Module[{F, V},
      {F, V} = RCF3[1/R, x];
      Return[{1/F, 1/V}]];

Off[General::spell1];



End[ ];

EndPackage[ ];