(* Guess. by Manuel Kauers *)

<< LinearSystemSolver.m

Guess`Private`Version = "0.4 2007-10-16"; 

BeginPackage["Guess`", {"DiscreteMath`", "Algebra`SymmetricPolynomials`", 
                        "NumberTheory`NumberTheoryFunctions`", "NumberTheory`AlgebraicNumberFields`",
                        "Algebra`Horner`", "LinearSystemSolver`"}];

SetAttributes[Algebra`Horner, Listable];

`OperatorGCD::usage = "OperatorGCD[rec1, rec2, f[n]] computes the greatest common right divisor of " <>
   "the two given univariate recurrence operators";
`Padd::usage = "Padd[data, u] transforms something like {{1,2},{1,2,3,4,5}} into something like " <>
   "{{1,2,u,u,u},{1,2,3,4,5}}";
`GuessMultRE::usage = "GuessMultRE[data, strset, vars, deg] computes a vector space basis of the space " <>
   "of all recurrences satisfied by the array data.\n\n" <>
   "data... a d-dimensional rectangular array of numbers or rational functions\n" <>
   "strset... a list of expressions of the form f[lc,lc,...] with f being a symbol and " <>
   "the lc being integer-linear combinations of the variables (see below)\n" <>
   "vars... a list of variables, e.g., {n,m,k}, or a list of blocks, e.g., {{n,m},{k}}, indicating " <>
   "that the polynomial coefficients in the recurrence equations should be symmetric wrt. all the " <>
   "variables belonging to the same block. {n,m,k} is equivalent to {{n},{m},{k}}.\n" <>
   "deg... total degree bound for the polynomial coefficients in the recurrence equations, or " <>
   "a list of individual degree bounds for each block.\n\n" <>
   "Options:\n" <>
   "Modulus -> prime number used for solution shape prediction. (0 to skip this step.)\n" <>
   "AdditionalEquations -> number of extra equations to be used, Infinity to take all data.\n" <>
   "StartPoint -> point corresponding to data[[1,1,1...]], default {0,0,0,...}\n" <>
   "Extension -> minimal polynomial of a single parameter, or {} if there is not algebraic extension.\n" <>
   "Except -> list of terms, e.g., n^3 m f[n+2, m+1], that must not occur in the recurrences.\n" <>
   "MustHaveOneOf -> list of terms, e.g., n^3 m f[n+2,m+1], of which at least one must occur in " <>
   "the recurrences, default All\n" <>
   "Factor -> True or False depending on whether coefficients of resulting recurrences should be factored.";

{AdditionalEquations, StartPoint, Except, MustHaveOneOf};

Begin["`Private`"]

(**
 * INPUT: two P-finite recurrence operators
 * OUTPUT: their greatest common right divisor
 **)
OperatorGCD[rec1_, rec2_, f_[n_], opts___] :=
  Module[ {a, b, d, i},
  
    a = rec1 /. n -> n - Min[Cases[{rec1}, f[n + i_.] -> i, Infinity]];
    If[ OperatorOrd[a, f[n]] === 0, Return[ f[n] ]];
    a = Together[#, opts]& /@ Collect[a, f[_]];
    d = PolynomialLCM[#, opts]&@@(Denominator /@ a);
    a = Expand[Cancel[#*d, opts], opts]& /@ a;

    b = rec2 /. n -> n - Min[Cases[{rec2}, f[n + i_.] -> i, Infinity]];
    If[ OperatorOrd[b, f[n]] === 0, Return[ f[n] ]];
    b = Together[#, opts]& /@ Collect[b, f[_]];
    d = PolynomialLCM[#, opts]&@@(Denominator /@ b);
    b = Expand[Cancel[#*d, opts], opts]& /@ b;

    Return[Block[{$RecursionLimit=10^5}, OperatorEuclid[a, b, f[n], opts]]];
  ];
OperatorOrd[a_, f_[n_]] := Max[Cases[{a}, f[n+i_.] -> i, Infinity]];
OperatorEuclid[a_, b_, f_[n_], opts___] := (* brutal euclid *)
  Module[ {da, db, lta, ltb, lcm, c, dc, g, i},

    If[ b === 0, Return[ a ] ];

    da = OperatorOrd[a, f[n]]; db = OperatorOrd[b, f[n]]; 
    If[ da < db, Return[OperatorEuclid[b, a, f[n], opts]] ];
    Print[ {da, db} ];

    lta = Coefficient[a, f[n + da]]; ltb = Coefficient[b, f[n + db]] /. n -> n + da - db;
    lcm = PolynomialLCM[lta, ltb, opts];
    lta = Cancel[lcm/lta, opts]; ltb = Cancel[lcm/ltb, opts];
    c = Sum[Expand[lta*Coefficient[a, f[n+i]] 
                 - ltb*(Coefficient[b, f[n+i-da+db]] /. n -> n + da - db), opts]f[n+i], {i, 0, da}];

    dc = OperatorOrd[c, f[n]];

    Which[
      dc === -Infinity, Return[b],
      dc === 0, Return[f[n]]
    ];

    g = (List@@c) /. f[_] -> 1;
    g = Fold[PolynomialGCD[#1, #2, opts]&, First[g], Rest[g]];
    c = Sum[Cancel[Coefficient[c, f[n+i]]/g, opts]*f[n+i], {i, 0, dc}];
    Return[OperatorEuclid[c, b, f[n], opts]];
  ];

(**
 * Padd a nested list to a full array
 *
 * INPUT: something like {{1,2},{1,2,3,4,5}}
 * OUTPUT: something like {{1,2,u,u,u},{1,2,3,4,5}}
 **)
Padd[data_, u_] := 
  Module[ {dpt, n},
    dpt = Depth[data /. (a_ /; FreeQ[a, List]) -> 1] - 2;
    n = Max[Map[Length, data, {dpt}]];    
    Map[Join[#, Table[u, {n - Length[#]}]]&, data, {dpt}](*;*)
  ];

std2ff[list_] := 
  Module[ {i, j}, 
    Transpose[Table[StirlingS2[i, j], {i, 0, Length[list] - 1}, {j, 0, Length[list] - 1}]] . list
  ];

(* ************************************ multivariate guesser **************************************** *)

(**
 * Multivariate guessing function
 *
 * GIVEN:
 * -> a d-dimensional array data
 * -> a list of variables 
 *    variables may be grouped for indicating symmetry in the data.
 *    then only coefficients with these symmetries will be searched for.
 * -> a structure set sset of the form {f[lc,lc,lc,lc], ...}
 *    with each lc being an integer linear combination of variables 
 *    (f being a fixed symbol)
 * -> a total degree bound deg
 *
 * FIND:
 * -> a linear combination of the structure set elements with polynomial
 *    coefficients that vanishes for the given data
 * 
 **)
Options[GuessMultRE] = {Modulus -> PreviousPrime[2^28], (* 0 to skip hom. precomputation *) 
  Return -> "sol", (* undocumented option for returning the determinant of the system "det" 
                      or the system "sys" itself instead of the solution space "sol", 
                      or "null" to just make a prediction and then return Null *)
  AdditionalEquations -> 100, StartPoint -> Automatic, Extension -> {}, Except -> {},
  MustHaveOneOf -> All, Factor -> True, Sort -> (OrderedQ[{#1,#2}]&)};
GuessMultRE[data_List, sset_, vars_List, deg_, opts:(_Rule|_RuleDelayed)...] :=
  Module[ {i, j, k, mod, dim, hsset, myvars, f, args, terms, sys, sol, p, dom, ret, 
           forced, termorder, 
           actualTerms, e, spols, params, addEqs, blocks, null, ext, except, points},

    {p, addEqs, null, ext, except, ret, fac, forced, termorder} = 
         {Modulus, AdditionalEquations, StartPoint, Extension, Except, Return, Factor, 
          MustHaveOneOf, Sort}     /. {opts} /. Options[GuessMultRE];
    If[ p!=0 && !PrimeQ[p], Throw["Modulus is not a prime!"]];
    If[ p==0 && forced =!= All,
       Throw["modular computation cannot be skipped when there are forced terms."]
    ];
    If[ Length[Complement[First /@ {opts}, First /@ Options[GuessMultRE]]] > 0,
       Throw["unknown options: " <> 
       ToString[Complement[First /@ {opts}, First /@ Options[GuessMultRE]]]];
    ];

    hsset = Complement[sset, {1}]; (* homogeneous terms of the structure set *)

    (* 1. check consistency of input *)
    If[ !MatchQ[Union[Length /@ hsset], {_Integer}],
      Throw["illegal structure set specification"]
    ];
    dim = First[Union[Length /@ hsset]];
    If[ null === Automatic, null = Table[0, {dim}] ];
    myvars = Union[Flatten[Apply[List, Map[Variables, hsset, {2}], 1]]];
    nvars = Length[myvars];

    If[ Sort[Flatten[vars]] =!= Sort[myvars], Throw["unexpected variable declaration"] ];
    myvars = Flatten[vars];
    blocks = If[ MatchQ[Head /@ vars, {List...}], vars, {#}& /@ vars ];

    If[ ListQ[deg] && Length[blocks] =!= Length[deg], 
      Throw["incompatible degree specification"]
    ];

    (* a. check that data is dim-dimensional *)
    If[ Length[Dimensions[data]] =!= dim,
      Throw["data is not a " <> ToString[dim] <> "-dimensional array"];
    ];
    If[ Length[null] =!= Length[myvars],
      Throw["dimensions of data and startpoint don't fit."];
    ];

    (* b. check that the elements of sset have the right form *)
    args = Flatten[Apply[List, hsset, 1]];
    If[ !MatchQ[hsset, {f_[__]..}] ||
        !MatchQ[args /. Thread[myvars -> 0], {_Integer...}] ||
        !(And@@(IntegerQ/@Flatten[Outer[Coefficient, args, myvars]])),

      Throw["illegal structure set specification"];
    ];
    f = hsset[[1, 0]]; (* function symbol in the structure set *)

    (* 2. ansatz *)
    (* a. terms for the polynomial coefficients *)
    terms = 
      If[ IntegerQ[deg],
        (* no symmetries specified *)
        (*
        Flatten[Table@@Join[
           {Times@@(vars^Array[i,Length[vars]])}, 
           Table[{i[j], 0, deg - Sum[i[k], {k, 1, j - 1}]}, {j, 1, nvars}]]],
        *)
        (* use symmetries in ansatz *)
        Flatten[Table@@Join[
           {f@@Thread[SymPolys[blocks, Array[i, Length[blocks]]]]},
           Table[{i[j], 0, deg - Sum[i[k], {k, 1, j - 1}]}, {j, 1, Length[blocks]}]]
           /. f[i__] :> Outer[Times, i]
        ]
      ,
        Flatten[Table@@Join[
           {f@@Thread[SymPolys[blocks, Array[i, Length[blocks]]]]},
           Table[{i[j], 0, deg[[j]]}, {j, 1, Length[blocks]}]]
           /. f[i__] :> Outer[Times, i]
        ]
      ];
       
    (* b. include the structure set *)
    terms = Flatten[Outer[Times, terms, sset]];
    terms = Complement[terms, except]; (* kick out unwanted terms *)
    If[ forced === All, forced = terms ];
    forced = Intersection[terms, forced]; 
    terms = Join[Sort[Complement[terms, forced], termorder], 
                 Sort[forced, termorder]]; (* put forced terms at the end *)
    Print[Length[terms], " terms"]; 

    (* now every term[[i]] has the form n^a k^b ... f[c n + d k + .., ..] *)

    (* 3. determine the domain of points for which all arguments in the
          structure set are within the range supplied in the data array *)
    dom = Table[0 <= hsset[[i,j]] < Length[Nest[First, data, j - 1]], {i, 1, Length[hsset]}, {j, 1, dim}];
    If[ vars=!=myvars,
      dom = Join[dom, Apply[LessEqual, vars, {1}]];
    ];
    dom = CylindricalDecomposition[Flatten[dom], myvars]; 
    If[ dom === False, Throw[ "insufficient input data." ]];
    dom = GenerateTuples[dom, myvars];
    Which[ 
      Count[data, 0, Infinity]/Times@@Dimensions[data] < 1/100,
        Null,
      Length[dom] < 125000,
        e[point__] := data[[Sequence@@({point}+1)]];
        dom = Select[dom, !MatchQ[Apply[e, hsset /. Thread[vars->#], {1}], {0..}]&],
      True,
        Print["Warning: There are many zeros in the given data."]
    ];
    If[ Length[dom] < Length[terms], Throw[ "insufficient input data." ]];
    
    (* 4. construct modular system *)
    params = Union@@(Variables /@ Flatten[data]);
    mod = Thread[params -> Table[Random[Integer, {0, 20}], {i, 1, Length[params]}]];
    ClearAll[e];
    e[point_][t_] := Module[ {j},
      PolynomialMod[t /. Thread[myvars -> point] 
         /. f[j__] :> data[[Sequence@@({j}+1)]]
         /. mod, p] // Return;
    ];

    k = Length[terms] + If[Head[addEqs]===List, First[addEqs], addEqs];
    points = If[ k > Length[dom],
                 Range[Length[dom]], 
                 First[Subsets[Range[Length[dom]], {k}, {Random[Integer, {1, Binomial[Length[dom], k]}]}]]
                ];
    Clear[k];

    If[ ext === {} && ret =!= "det" && p =!= 0, 

      sys = DeleteCases[Table[e[dom[[points[[i]]]]] /@ terms, {i, 1, Length[points]}], {0..}];

      (* 5. solve modular system *)
      Print["modular system: ", Length[sys], " eqns, ", Length[terms], " vars"];
      If[ Length[sys] < Length[terms], Throw["too many zeros."] ];

      sol = NullSpace[sys, Modulus -> p];
      sol = DeleteCases[sol, Prepend[Table[0, {Length[forced]}], ___]];
      If[ Length[sol] == 0, Print["no solution."]; Return[ {} ] ];
    
      Print[Length[sol], " solutions predicted."];
      If[ ret == "null", Return[ Null ] ];

      (* 6. construct refined system *)
      actualTerms = (Min[#, 1]& /@ Plus@@Map[If[#===0, 0, 1]&, sol, {2}]) * terms;
      (* actualTerms = Map[If[#===0,0,1]&, sol[[3]], {1}] * terms; (* SPECIAL *)  *)
      actualTerms = DeleteCases[actualTerms, 0];
    ,
      actualTerms = terms;
    ];

    Which[ 
      Count[data, 0, Infinity]/Times@@Dimensions[data] < 1/100,
        Null,
      Length[dom] < 125000,
        e[point__] := data[[Sequence@@({point}+1)]];
        hsset = Union[Cases[actualTerms, f[__], Infinity]];
        dom = Select[dom, !MatchQ[Apply[e, hsset /. Thread[vars->#], {1}], {0..}]&],
      True, 
        Print["Warning: There are many zeros in the given data."]
    ];
    k = Length[actualTerms] + If[Head[addEqs]===List, Last[addEqs], addEqs];
    points = If[ k > Length[dom],
                 Range[Length[dom]], 
                 First[Subsets[Range[Length[dom]], {k}, {Random[Integer, {1, Binomial[Length[dom], k]}]}]]
                ];
    Clear[k];
    ClearAll[e];
    e[point_][t_] := Module[ {j},
      t /. Thread[myvars -> point] /. f[j__] :> data[[Sequence@@({j}+1)]] // Return;
    ];
    sys = DeleteCases[Table[e[dom[[points[[i]]]]] /@ actualTerms, {i, 1, Length[points]}], {0..}];
    If[ Length[sys] < Length[actualTerms], Throw["too many zeros."] ];

    Which[ 
      ret == "det", 
        While[ Length[sys] > Length[First[sys]], sys = Delete[sys, Random[Integer, Length[sys]]] ];
        Return[Det[sys]],
      ret == "sys",
        Return[sys];
    ];

    (* 7. solve refined system *)
    Print["refined system: ", Length[sys], " eqns, ", Length[actualTerms], " vars"];

    sol = Which[ 
            ext =!= {}, LinearSystemSolver`LinSolve[sys, First[params], ext],
            Length[params] == 0, LinearSystemSolver`LinSolveQ[sys],
            Length[params] == 1, LinearSystemSolver`LinSolveUniv[sys, 0, First[params]],
            NameQ["RISCComb`LinSolve`M"], RISCComb`LinSolve`EANullSpace[sys],
            True, NullSpace[sys]
          ];

    sol = DeleteCases[sol, Prepend[Table[0, {Length[forced]}], ___]];
    Print[Length[sol], " solutions."]; 
    If[ Length[sol] == 0, Return[{}]];

    (* 8. done. *)    
    Return[(If[fac, Factor, Identity]/@ Collect[#, f[__]])& 
            /@ (sol . actualTerms /. Thread[myvars->myvars-null])
            /. f[j__] :> (f[j] /. Thread[myvars->myvars+null])
    ];
  ];

GenerateTuples[dom_, vars_] :=
  GenerateTuples[dom, vars] = (* caching *)
  Module[ {t, int, rest, i, f, a, b, c},

    If[ Length[vars] == 0, Return[{{}}] ];

    If[ Head[dom] === Or, Return[Union @@ (GenerateTuples[#, vars]& /@ dom)]];

    If[ Length[vars] == 1,
      If[ !MemberQ[{Less,LessEqual,Inequality,Equal}, Head[dom]],
        Throw["illegal domain description encountered: ", dom]
      ];
      int = dom;
      rest = True;
    ,
      If[ Head[dom] === And && !MemberQ[{Less,LessEqual,Inequality,Equal}, Head[First[dom]]],
        Throw["illegal domain description encountered: ", dom]
      ];
      int = First[dom];
      rest = Rest[dom];
    ];

    If[ Head[int] === Equal,
      int = Last[int] <= First[int] <= Last[int];
    ];

    If[ MemberQ[{Less, LessEqual}, Head[int]],
      int = Inequality[First[int], Head[int], int[[2]], Head[int], Last[int]];
    ];

    t = Range[Ceiling[First[int]], Floor[Last[int]]];
    If[IntegerQ[First[int]] && int[[2]] === Less, t = Rest[t]];
    If[IntegerQ[Last[int]] && int[[4]] === Less, t = Most[t]];

    Return[
      Union@@Table[
        Prepend[#, t[[i]]]& /@ GenerateTuples[rest /. First[vars] -> t[[i]], Rest[vars]]
      , {i, 1, Length[t]}]
    ];
  ];

(**
 * A basis for the vector space of all homogeneous symmetric polynomials of total degree d
 **)
SymPolys[v_List, d_Integer] := 
  Module[ {e, i, j, k, p},

    If[ d == 0, Return[{1}] ] ;

    e = Table[SymmetricPolynomial[v, i], {i, 1, Length[v]}];
    p = DiscreteMath`IntegerPartitions`IntegerPartitions[d, Min[d, Length[v]]];

    Product[e[[#[[i, 1]]]]^#[[i, 2]], {i, 1, Length[#]}]& /@ p
  ];

End[]

EndPackage[]

If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.796887, 0.789075, 0.871107]]]&,
   Print
  ]["Guess  Package by Manuel Kauers " <> 
    "\[LongDash] \[Copyright] RISC Linz \[LongDash] V " <> Guess`Private`Version];