
(* Experimental code for multiplication and factorization of C-finite sequences *)
(* C-finite sequences specified by characteristic polynomials, without initial values *)

(* Manuel Kauers, December 2015 *)

(* Main functions:

 CFiniteTimes[poly, poly, variable] -- returns p\otimes q for given p and q. constant fields with 
    parameters might not work well. 

 CFiniteFactor[poly, variable] -- returns a list of pairs {p,q} such that p\otimes q is the given poly.
    case of multiple roots not yet implemented.  

 *)

SquareFreePart[p_] := Times@@(First/@Rest[FactorSquareFreeList[p]])

GuessCFinite[data_List, x_] := (Reverse[Prepend[#, -1]].x^Range[0, Length[#]])&[FindLinearRecurrence[data]]

CFinite2List[p_, x_, init_, n_Integer] :=
  Module[ {rec, L},
    rec = -Most[CoefficientList[p, x]]/Coefficient[p, x, Exponent[p, x]];
    L = init;
    While[ Length[L] < n, AppendTo[L, Together[Take[L, -Length[rec]].rec]] ];
    Return[Take[L, n]]
  ]

MakeMonic[p_, x_] := Expand[(#/Coefficient[#, x, Exponent[p, x]])&[p]]

MakeDoubleMonic[p_, x_] := MakeMonic[p /. x -> Coefficient[p, x, Exponent[p, x] - 1] x, x];  

CFiniteTimes[p1_, p2_, x_] := (* checked *)
  Module[ {p1list, p2list, t, q, d1, d2, q1, q2, qq},
    p1list = Rest[FactorSquareFreeList[p1]];
    p2list = Rest[FactorSquareFreeList[p2]];
    Collect[Numerator[Together[#/x^Exponent[#,x,Min]/(#/.x->0)]]&[PolynomialLCM@@Flatten[Table[
      q1 = p1list[[i, 1]]; d1 = Exponent[q1, x]; q2 = p2list[[j, 1]]; d2 = Exponent[q2, x];
      qq = ToNumberField[Join[CoefficientList[q1, x], CoefficientList[q2, x]]];
      q1 = Take[qq, d1+1].x^Range[0,d1]; q2 = Take[qq, -(d2+1)].x^Range[0,d2]; 
      SquareFreePart[
        GuessCFinite[Expand[CFinite2List[q1, x, Prime /@ (10+Range[d1]), 2*d1*d2 + 1]*
                            CFinite2List[q2, x, Prime /@ (20+Range[d2]), 2*d1*d2 + 1]], x]
      ]^(p1list[[i, 2]] + p2list[[j, 2]] - 1)
      , {i, 1, Length[p1list]}
      , {j, 1, Length[p2list]}      
    ]]], x, Factor]
  ];

CFiniteFactor[pol_, x_] :=
  Module[ {p, v, eq, M, piList, pi, k, i, j},
     p = SquareFreePart[pol] /. x -> #; (* NYI: polynomials with multiple roots *)
     v = Complement[Variables[p], {#}];
     v = Thread[v -> (Prime/@(10+Range[Length[v]]))];
     eq = Abs[N[#1/.v] - N[#2/.v]] < 10^(-6)&;
     roots = Table[Root[Evaluate[p]&, i], {i, 1, Exponent[p, #]}];
     M = Prepend[IdentityMatrix[Length[roots]], Range[Length[roots]]];
     piList = addRow[M, 2, roots, eq];
     Union[Union /@ Table[
       pi = piList[[k]];
       {(x-1)Product[x - roots[[pi[[i,1]]]]/roots[[pi[[1,1]]]], {i, 2, Length[pi]}],
        (x-roots[[pi[[1,1]]]])Product[x - roots[[pi[[1,j]]]], {j, 2, Length[pi[[1]]]}]}
     , {k, 1, Length[piList]}]]
  ];

SetAttributes[addRow, HoldFirst]
addRow[M_, n_, roots_, eq_] :=
  Module[ {out, pi},
     out = {};
     Do[
       M[[n]] = 0*M[[n]]; M[[n,1]] = p;
       Do[
         If[ M[[n - 1, j]] != 0,
           q = Select[Complement[Range[Length[roots]], M[[n]]],
                      eq[roots[[1]]/roots[[p]], roots[[j]]/roots[[#]]] &&
                      eq[roots[[1]]/roots[[j]], roots[[p]]/roots[[#]]]&];
           If[ Length[q] > 0, M[[n,j]] = q[[1]] ];
         ];
       , {j, 2, Length[roots]}];
       pi = Select[Transpose[Take[M, n]], Last[#]!=0&];
       Which[
         Union[Flatten[pi]] == Range[Length[roots]], AppendTo[out, Transpose[pi]],
         Length[DeleteCases[M[[n]], 0]] > n, out = Join[out, addRow[M, n+1, roots, eq]]
       ];
     , {p, M[[n-1,1]] + 1, Length[roots]}];
     Return[out];
  ];

