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[ ];