
(* Guess. by Manuel Kauers *)

(* ------------------------------------------------------------------------- *)

(* LICENCE:                                                                  *)
(*                                                                           *)
(*  This software package is free; you can redistribute it and/or modify it  *)
(*  under the terms of the GNU General Public License as published by the    *)
(*  Free Software Foundation; either version 2 of the License, or (at your   *)
(*  option) any later version. The package is distributed in the hope that   *)
(*  it will be useful, but without any warranty; without even the implied    *)
(*  warranty of merchantability or fitness for a particular purpose. See the *)
(*  GNU General Public License for more details.                             *)

(* ------------------------------------------------------------------------- *)

(* todo: speed up construction of system by caching evaluations of monomials in GuessMultRE *)
(* todo: are larger primes faster? *)

<< LinearSystemSolver.m (* paste before distribution *) 

Guess`Private`Version = "0.58 2020-11-20"
(* for MMA6 from version 0.6 on *)
(* for MMA7 from version 0.32 on *)
(* for MMA8 from version 0.39 on *)
(* for MMA10 from version 0.53 on *)
(* for MMA11 from version 0.54 on *)

BeginPackage["RISCComb`"]

{`Infolevel}

EndPackage[]

BeginPackage["HolonomicFunctions`"]

{`Der}

EndPackage[]

BeginPackage["GeneratingFunctions`"]

{`AdditionalEquations, `SelectPoints, `SelectTerms}

EndPackage[]

Off[General::compat] (* As of MMA8 I know I only use parts of Combinatorica which have not been changed from MMA7 to MMA8 *)

BeginPackage["Guess`", {"GeneratingFunctions`", "RISCComb`", "Combinatorica`", "LinearSystemSolver`", "HolonomicFunctions`"}]

On[General::compat]

`GuessCurveRE::usage = "GuessCurveRE[data, f[n]] computes a list of pairs {order, degree, dim} such that " <>
    "the data satisfies dim linearly independent recurrence equations of the given order and degree and " <>
    "none of the same order but smaller degree or of the same degree but smaller order.";

`DifferentialGCD::usage = "DifferentialGCD[de1, de2, f[x]] computes the greatest common right divisor of the two given univariate differential operators";
`DifferenceGCD::usage = "DifferenceGCD[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}}";
`GuessSymmetricRE::usage = "under construction"
`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[n+int,m+int,...] with f being a symbol, " <>
   "n,m,... being variables and int being integers\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}}. If a variable " <> 
   "is of the form q^n, then the q-shift is applied in direction n.\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" <>
   "AdditionalTerms -> list of terms to be included in the ansatz in addition to those following " <>
   "from structure set and degree specification.\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. " <>
   "If a term is embraced with Cone[ ], e.g., Cone[n^3 m f[n+2, m+1]], all multiplies of it will " <> 
   "be excluded, too.\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" <>
   "Constraints -> logical expression constraining the points from which the sample values should be " <>
   "chosen, e.g., Abs[n-m]<=3\n" <> 
   "Sort -> internal ordering for terms to be used in the ansatz\n" <> 
   "Factor -> True or False depending on whether coefficients of resulting recurrences should be factored.";
`GuessUnivRE::usage = "GuessUnivRE[list, f[n]] a vector space basis of the space of all recurrences " <>
   "satisfied by the given data.\n\n" <>
   "list... a list of rational numbers\n" <>
   "f[n]... results should be phrased in terms of f and n\n\n" <>
   "Options:\n" <>
   "Order -> maximum order of recurrences to be searched for (default: 5)\n" <>
   "Degree -> maximum polynomial degree of recurrences to be searched for (default: 5)\n" <>
   "AdditionalEquations -> number of extra equations to be used, Infinity to take all data.";
`GuessUnivDE::usage = "GuessUnivDE[list, f[x]] a vector space basis of the space of all differential " <>
   "equations satisfied by power series whose first coefficients match the given data.\n\n" <>
   "list... a list of rational numbers\n" <>
   "f[x]... results should be phrased in terms of f and x\n\n" <>
   "Options:\n" <>
   "Order -> maximum order of diff eqs to be searched for (default: 5)\n" <>
   "Degree -> maximum polynomial degree of diff eqs to be searched for (default: 5)\n" <>
   "AdditionalEquations -> number of extra equations to be used, Infinity to take all data.";
`GuessMinRE::usage = "GuessMinRE[list, f[n]] finds the least order recurrence satisfied by the data in the given list.\n\n" <> 
 "Options:\n" <> 
 "Max -> maximum number of elements to be used for guessing (the rest is used as backup)\n" <>
 "Modulus -> prime to be used for precomputations\n" <> 
 "Return -> \"rec\" to return the final recurrence\n" <> 
 "StartPoint -> n0 such that the first element of the list corresponds to f[n0]\n" <>
 "Involevel -> integer governing verbosity of the computation";
`GuessMinDE::usage = "GuessMinDE[list, f[x]] finds the least order differential equation satisfied by a power series whose first coefficients are given in the list.\n\n" <> 
 "Options:\n" <> 
 "Max -> maximum number of elements to be used for guessing (the rest is used as backup)\n" <>
 "Modulus -> prime to be used for precomputations\n" <> 
 "Return -> \"rec\" to return the final recurrence\n" <> 
 "StartPoint -> n0 such that the first element of the list corresponds to f[n0]\n" <>
 "Involevel -> integer governing verbosity of the computation";

`GuessUnivAE::usage = "GuessUnivAE[list, f[x], dy, dx, mod] finds an algebraic equation of degree dy with polynomial coefficients of degree x over either QQ or the finite field with mod elements.";

`GuessMinAE::usage =  "GuessMinAE[list, f[x]] finds the least degree algebraic equation satisfied by a power series whose first coefficients are given in the list.\n\n" <> 
 "Options:\n" <> 
 "Max -> maximum number of elements to be used for guessing (the rest is used as backup)\n" <>
 "Modulus -> prime to be used for precomputations\n" <> 
 "Return -> \"bool\" to just indicate if something exists, \"shape\" to return bidegree of the equation, \"ae\" (default) to return the algebraic equation itself\n" <> 
 "StartPoint -> n0 such that the first element of the list corresponds to f[n0]\n" <>
 "Involevel -> integer governing verbosity of the computation";

`GuessMultDE::usage = "tba"

{StartPoint, Except, MustHaveOneOf, Constraints, Cone, AdditionalTerms};

Begin["`Private`"]

InfoLevel = 1;
DebugPrint[i_Integer, opts___][data__] := 
  Module[ {level},
    level = Infolevel /. {opts} /. Infolevel -> InfoLevel;
    If[i < level, Print[data]]
  ];

(**** Curve Guessing ****)

ClearAll[GuessCurveRE];
GuessCurveRE[data_] :=
   Module[ {rec, rG, dG, r, d, corners, f, n, i, p, q},
      rec = GuessMinRE[data, f[n], Modulus -> 2^31-1]; corners = {{DifferenceOrd[rec, f[n]], Exponent[rec, n]}};
      p = Coefficient[rec, f[n]]; q = Function[t, Evaluate[Coefficient[rec, f[n+corners[[1,1]]]] /. n -> n + t]];
      Print["minimum point: ", First[corners]];
      disp = Max[Select[Range[50], Exponent[PolynomialGCD[p, q[#], Modulus -> 2^31-1], n]>0&]];
      Print["dispersion = ", disp];
      Print["maximal relevant order therefore: ", corners[[1,1]] + disp corners[[1,2]]];
      Catch[Do[
        If[Mod[r, 10]==0, Print["checking order ", r]];
        For[d = corners[[-1, 2]] - 1, (dim = GuessUnivRE[data, f[n], Order -> r, Degree -> d, Return -> "dim"])>0, d--];
        AppendTo[corners, {r, d+1}];
      , {r, corners[[1,1]] + 1, Min[75, corners[[1,1]] + disp corners[[1,2]]]}]];
      Return[corners];
   ];

(**************************************** Operator GCD routines **********************************)

DifferentialOrd[a_, f_[x_]] :=
  Max[0, Cases[{a}, Derivative[i_][f][_] -> i, Infinity]];

DifferenceOrd[a_, f_[n_]] := Max[Cases[{a}, f[n + i_.] -> i, Infinity]]

(**
 * INPUT: two P-finite differential operators
 * OUTPUT: their greatest common right divisor
 **)
DifferentialGCD[de1_, de2_, f_[x_], opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {g, i, m},
    m = Modulus /. {opts} /. Modulus -> 0;
    g = OreGCD[Table[Coefficient[de1, D[f[x],{x,i}]], {i, 0, DifferentialOrd[de1, f[x]]}],
               Table[Coefficient[de2, D[f[x],{x,i}]], {i, 0, DifferentialOrd[de2, f[x]]}],
               x, m, Identity, D[#, x]&];
    Return[g.Table[D[f[x],{x,i}],{i,0,Length[g]-1}]];
  ];

(**
 * INPUT: two P-finite recurrence operators
 * OUTPUT: their greatest common right divisor
 **)
DifferenceGCD[rec1_, rec2_, f_[n_], opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {g, i, m},

    m = Modulus /. {opts} /. Modulus -> 0;
    g = OreGCD[Table[Coefficient[rec1, f[n+i]], {i, 0, DifferenceOrd[rec1, f[n]]}], 
               Table[Coefficient[rec2, f[n+i]], {i, 0, DifferenceOrd[rec2, f[n]]}], 
               n, m, (#/.n->n+1)&, 0&];
    Return[g.(f/@(n+Range[0,Length[g]-1]))];
  ];

OreGCD[rec1_List, rec2_List, n_, p_Integer, sigma_, delta_] := 
  Block[ {$RecursionLimit=10^5}, 
    OperatorEuclid[OperatorNormal[rec1, n, p], OperatorNormal[rec2, n, p], n, p, sigma, delta]
  ];

OperatorNormal[op_, n_, p_Integer] := Module[ {op0, a, ord, den},

  op0 = op /. {a___, 0..} -> {a}; (* trim leading zeros *)
  If[ (ord = Length[op0]) === 0, Return[{}] ]; (* zero operator *)  
  (* 
  op0 = op0 /. {0.., a__} -> {a}; (* trim trailing zeros *)
  op0 = Together[op0 /. n -> n + Length[op0] - ord, Modulus -> p]; (* adjust coeffs *)
  *)
  den = PolynomialLCM@@Append[Denominator /@ op0, Modulus -> p]; (* common denominator *)
  op0 = Together[op0 den, Modulus -> p]; (* clear denominators *)
  con = PolynomialGCD@@Append[op0, Modulus -> p]; (* content *)
  op0 = Expand[Together[op0/con, Modulus -> p], Modulus -> p]; (* clear content *)
  Return[op0];
];

OperatorEuclid[a_, b_, n_, p_Integer, sigma_, delta_] := (* iterative code... *)
  Module[ {A, B, C, da, db, lta, ltb, lcm, g, i},

    If[ Length[b] === 0, Return[ a ] ];
    If[ Length[a] < Length[b], Return[OperatorEuclid[b, a, n, p, sigma, delta]] ];

    A = a; B = b; 

    While[ True, 

      da = Length[A]; db = Length[B];
      lcm = PolynomialLCM[lta = Last[A], ltb = Nest[sigma, Last[B], da-db], Modulus -> p];
      lta = Together[lcm/lta, Modulus -> p]; ltb = Together[lcm/ltb, Modulus -> p];
      C = Nest[Expand[Prepend[sigma /@ #, 0] + Append[delta /@ #, 0], Modulus -> p]&, B, da-db];
      C = Expand[lta*Most[A] - ltb*Most[C], Modulus -> p];
      Clear[lcm, lta, ltb];
      While[Length[C] > 0 && Last[C]==0, C=Most[C]];
      If[Length[C]==0, Break[]];
      A=B; Clear[B];
      B=Together[C/Fold[PolynomialGCD[#1, #2, Modulus -> p]&, First[C], Rest[C]], Modulus -> p];
      Clear[C];
      If[ Length[A] < Length[B], {A,B}={B,A} ];
    ];

    If[ Length[A] < Length[B], A, B]
  ];

(**
 * 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
  ];

(* ************************************ univariate guesser(s) *************************************** *)

Options[GuessUnivDE] = {Order -> 5, Degree -> 5, AdditionalEquations -> 100, Modulus -> 0,
  Return -> "sol", Infolevel -> 0, Catch -> True};
GuessUnivDE[coeffs_, f_[x_], opts___] := 
  Module[ {ord, deg, add, M, terms, sol, p, mod, ret, eqs, subs, filter, t},

    {ord, deg, add, mod, ret, filter} = {Order, Degree, AdditionalEquations, Modulus, Return, Catch} 
               /. {opts} /. Options[GuessUnivDE];
    Which[
      !MatchQ[coeffs,{(_?(Element[#, Rationals]&))...}], 
         Throw["coefficients must be rational numbers"],
      filter && Length[coeffs] < (ord+1)(deg+2)-1 , 
         Throw["insufficient data. Provide more coefficients, or decrease order and/or degree"] 
    ];
    eqs = Min[(ord+1)(deg+1) + add, Length[coeffs] - ord];
    subs = First[Subsets[Range[Length[coeffs] - ord], {eqs},
                   {Random[Integer, {1, Binomial[Length[coeffs] - ord, eqs]}]}]];
    subs = Take[Union[Range[Min[10, Length[coeffs] - ord]], subs], eqs];
    subs = Range[Length[coeffs] - ord];
    terms = Flatten[Table[x^j D[f[x], {x, i}], {i,0,ord}, {j,0,deg}]];

    M[p_] := Part[Transpose[
            Take[#, Length[coeffs] - ord]& /@ 
            Flatten[Table[PadLeft[#, j + Length[#]], {j, 0, deg}]& /@ 
            NestList[PolynomialMod[Rest[#]Range[Length[#]-1], p]&, PolynomialMod[coeffs, p], ord], 1]
    ], subs];

    p = If[ mod === 0, NextPrime[Developer`$MaxMachineInteger, -1], mod ];

    sol = NullSpace[M[p], Modulus -> p];
    If[ ret == "dim", Return[Length[sol]] ];
    If[ Length[sol] == 0, Return[ {} ]];
    If[ mod > 0, Return[sol.terms]];

    DebugPrint[1, opts][Length[sol], " solutions predicted."];
    Block[ {LinearSystemSolver`Private`PrintLevel = -1 + Infolevel /. {opts} /. Infolevel -> 0}, 
    sol = LinSolveQ[M];
    ];
    ClearAll[M];

    Return[If[Length[sol]==0, {}, Collect[sol.terms, _[x]]]];
  ];

Options[GuessUnivRE] = {Order -> 5, Degree -> 5, AdditionalEquations -> 100, StartPoint -> 0, Modulus -> 0, Infolevel -> 0, Return -> "sol", Catch -> True};
GuessUnivRE[coeffs_, f_[n_], opts___] := 
  Module[ {ord, deg, add, M, terms, sol, eqs, subs, p, start, mod, ret, filter},

    {ord, deg, add, start, mod, ret, filter} = {Order, Degree, AdditionalEquations, StartPoint, Modulus, Return, Catch} 
                             /. {opts} /. Options[GuessUnivRE];

    Which[ 
      !MatchQ[coeffs,{(_?(Element[#, Rationals]&))...}], 
         Throw["coefficients must be rational numbers"],
      filter && Length[coeffs] < (ord+1)(deg+2)-1, 
         Throw["insufficient data. Provide more coefficients, or decrease order and/or degree"] 
    ];
    eqs = Min[(ord+1)(deg+1) + add, Length[coeffs] - ord];
    subs = First[Subsets[Range[Length[coeffs] - ord], {eqs},
                   {Random[Integer, {1, Binomial[Length[coeffs] - ord, eqs]}]}]];
    terms = Flatten[Table[(n-start)^j f[n+i], {i,0,ord}, {j,0,deg}]];

    M[p_] := Part[Transpose[
      Take[#, Length[coeffs] - ord]& /@ 
      Flatten[Table[PolynomialMod[# PowerMod[Range[0, Length[#]-1], j, p], p], {j, 0, deg}]& /@ 
      NestList[Rest, PolynomialMod[coeffs, p], ord], 1]
    ], subs];

    p = If[ mod === 0, NextPrime[Developer`$MaxMachineInteger, -1], mod ];

    sol = NullSpace[M[p], Modulus -> p];
    If[ ret == "dim", Return[Length[sol]] ];
    If[ Length[sol] == 0, Return[ {} ]];
    If[ mod > 0, Return[sol.terms]];

    DebugPrint[1, opts][Length[sol], " solutions predicted."];
    Block[ {LinearSystemSolver`Private`PrintLevel = -1 + Infolevel /. {opts} /. Infolevel -> 0}, 
    sol = LinSolveQ[M];
    ];
    ClearAll[M];

    Return[If[Length[sol]==0, {}, Collect[sol.terms, f[_]]]];
  ];

(* a guesser which goes for the gcd; experience: gcd has shorter coefficients, so it needs
 fewer primes *)

T[1] = 0; (* search time *)
T[2] = 0; (* linear solver time *)
T[3] = 0; (* gcd time *)
T[4] = 0; (* CRA and RR time *)
T[5] = 0; (* termination check time *)
T[6] = 0; (* number of terms needed *)
T[7] = 0; (* order of final recurrence *)
T[8] = 0; (* degree of final recurrence *)
ClearAll[Time];
Time[i_Integer, expr_] := Module[ {r}, T[i] += First[Timing[r = expr]]; r ]; 
SetAttributes[Time, HoldAll];

GuessUnivOp[data_, f_[n_], ord1_Integer, deg1_Integer, ord2_Integer, deg2_Integer, p_Integer, guessfun_, gcdfun_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {sol},

    Time[2,
    sol = guessfun[data, f[n], Order -> ord1, Degree -> deg1, Modulus -> p, opts];
    sol = If[ Length[sol] > 1, sol, Join[sol, 
            guessfun[data, f[n], Order -> ord2, Degree -> deg2, Modulus -> p, opts]
          ]];
    ];

    If[ Length[sol] < 2, Throw["Too few solutions. Try different degrees and/or orders."] ];

    sol = Time[3, gcdfun[sol[[1]], sol[[2]], f[n], Modulus -> p]];

    Return[sol];
  ];

RRdefault[n_, q_] := LinearSystemSolver`Private`ReconstructQ0[n, q];
RRsmalldenom[n_, q_] := If[ n==0, 0,
  Module[ {r, s, rr, ss, quo},
    Return[RRdefault[n, q]];
    {r, s, rr, ss} = {q, 0, n, 1};
    While[ Abs[r]^11 > q^10,
      quo = Quotient[r, rr];
      {r, s, rr, ss} = {rr, ss, r - quo rr, s - quo ss}
    ];
    r/s
  ]];
List2Poly[coeffs_List, var_] := coeffs . var^(Range[Length[coeffs]]-1);
Poly2List[p_, var_, deg_] := Module[ {l}, 
   l = CoefficientList[p, var]; 
   If[ Length[l] <= deg + 1, 
       Return[Join[l, Table[0, {deg + 1 - Length[l]}]]],
       Return[Take[l, deg + 1]]
   ]];

(**
 * polynomial division
 **)
QuoRem[p_, q_, x_, m_] :=
  Module[ {p0, q0, quo, c, i, j},

    p0 = CoefficientList[p, x]; q0 = CoefficientList[q, x]; quo = {};

    If[ Length[p0] < Length[q0], Return[ {0, p} ] ];

    While[ Length[p0] >= Length[q0],

      c = PolynomialMod[Last[p0]/Last[q0], m];
      PrependTo[quo, c];

      j = Length[p0] - Length[q0] + 1;
      Do[ p0[[j++]] = p0[[j]] - c q0[[i]], {i, 1, Length[q0]} ]; 
      p0 = Most[p0];
    ];

    Return[ {List2Poly[quo, x], List2Poly[p0, x]} ]
  ];

MyPolynomialQuotient[p_, q_, x_, m_] := First[QuoRem[p, q, x, m]];
MyPolynomialRemainder[p_, q_, x_, m_] := Last[QuoRem[p, q, x, m]];
  
RRpoly[g_, mod_, x_, p_] :=
  Module[ {t, f, r0, r1, t0, t1, n, d, q, c, df, drat, dr1, dt1, unique},

    c = PolynomialMod[#, p]&;
    If[ c[g] === 0, Return[ 0 ] ];

    f = If[ ListQ[mod], Product[x - mod[[i]], {i, 1, Length[mod]}], mod ];
    df = Exponent[f, x];
    {r0, r1, t0, t1} = {f, g, 0, 1}; {n, d} = {r1, t1};
    drat = Exponent[n, x] + Exponent[d, x];
    unique = True;

    While[ r1 =!= 0,

      t = First[Timing[
      dr1 = Exponent[r1, x]; dt1 = Exponent[t1, x];
      If[ dr1 < df && t1 =!= 1 && drat == dr1 + dt1, unique = False];
      If[ dr1 < df && drat > dr1 + dt1,
        {n, d} = {r1, t1}; drat = dr1 + dt1;
        unique = True;
      ]; 
      ]];

      {t, q} = Timing[MyPolynomialQuotient[r0, r1, x, p]];
      {t, {r0, r1, t0, t1}} = Timing[{r1, c[r0 - q r1], t1, c[t0 - q t1]}];
    ];
    If[ !unique, Throw["ambigous reconstruction result"] ];
    
    {t, {n, d}} = Timing[CoefficientList[#, x]& /@ {n, d}];
    {t, {n, d}} = Timing[Map[c, {n, d}/Last[d], {2}]];
    {n, d} = List2Poly[#, x]& /@ {n, d};

    Return[Together[n/d, Modulus -> p]];
  ];

MyDate[] := Sequence@@(If[Head[#]===Integer && #<10, "0" <> ToString[#], #]& /@ 
                       Riffle[Floor[Date[]], {"-", "-", " ", ":", ":"}]);

GuessUnivOp[data_, f_[n_], ord1_Integer, deg1_Integer, ord2_Integer, deg2_Integer, guessfun_, gcdfun_, opfun_, ordfun_, checkfun_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {gcdq, gcd, gcdp, p, q, ord, deg, den, i, j, k, xstart, done, mod, ret, lc, vars, t, op, new, rr, r, d},

    {mod, ret} = {Modulus, Return} /. {opts} /. {Modulus -> 0, Return -> "rec"};
    If[ MatchQ[data, {0...}], Return[f[n]] ];

    vars = Variables[data];
    If[ Length[vars] == 1,
      t = vars[[1]]; rr = op = f[n]; p = 1; k = 3; 
      While[!checkfun[rr /. {t->k}, f[n], 45007, {t->k}],
        new = GuessUnivOp[data /. t -> k, f[n], ord1, deg1, ord2, deg2, guessfun, gcdfun, opfun, ordfun, checkfun[#1,#2,#3,{t->k}]&, opts];
        r = Max[ordfun[op,f[n]], ordfun[new,f[n]]]; d = Max[Exponent[new, t], Exponent[op, t]];
        new = Expand[new/Last[CoefficientList[Coefficient[new,opfun[f[n],r]], n]]];
        op = Expand[op - p*((op - new)/p /. t -> k)];
        p *= t - (k++);
        rr = Catch[Sum[Sum[RRpoly[Coefficient[Coefficient[op, opfun[f[n],i]],n, j], p, t, 0]n^j, {j, 0, d}]opfun[f[n],i], {i, 0, r}]];
      ];
      Return[rr];
    ];

    If[ mod =!= 0,
      gcdp = GuessUnivOp[data,f[n],ord1,deg1,ord2,deg2, mod, guessfun, gcdfun, opts];
      lc = Coefficient[gcdp, f[n + ordfun[gcdp, f[n]]]];
      If[ lc===0, lc = Coefficient[gcdp, D[f[n], {n, ordfun[gcdp, f[n]]}]] ];
      gcdp = Collect[Expand[gcdp/Last[CoefficientList[lc, n]], Modulus -> mod], f[_]];
      Return[If[ ret == "shape", {ordfun[gcdp, f[n]], Max[Exponent[gcdp, n]]}, gcdp] ];
    ];

    gcdq = {1}; q = 1; 
    p = NextPrime[10^150, -1]; 
    p = NextPrime[Developer`$MaxMachineInteger, -1];

    While[True,

       DebugPrint[3, opts]["...", StringTake[ToString[p], -5], "  ", MyDate[]];
       gcdp = GuessUnivOp[data, f[n], ord1, deg1, ord2, deg2, p, guessfun, gcdfun, opts];
       gcdp = Table[Coefficient[gcdp, opfun[f[n], i]], {i, 0, ord = ordfun[gcdp, f[n]]}];
       gcdp = PolynomialMod[gcdp/Coefficient[Last[gcdp], n, Last[deg = Exponent[gcdp, n]]], p];

       If[ ret == "shape", Return[ {ord, Max[deg]} ] ];

       Which[ 
         ord === 0, Throw["Too little data?"], (* give up *)
         q === 1, (* initialize *)
           DebugPrint[1, opts]["order ", ord, " and degree ", Max[deg], " predicted."];
           T[7] = ord; T[8] = Max[deg]; 
           gcdq = gcdp; q = p, 
         Length[gcdq] =!= ord + 1 || Exponent[gcdq, n] =!= deg, (* unlucky prime *)
           DebugPrint[0, opts]["unlucky prime discarded."]; 
           Continue[], 
         True, 
           gcdq = Time[4, Table[Sum[
                ChineseRemainder[{Coefficient[gcdp[[j]], n, i], Coefficient[gcdq[[j]], n, i]}, {p, q}]n^i, 
                {i, 0, deg[[j]]}], {j, 1, ord + 1}]];
           q *= p;
       ];
       Clear[gcdp];

       Time[4,       
       gcd = Sum[Sum[(*RRsmalldenom*)RRdefault[Coefficient[gcdq[[j + 1]], n, i], q]n^i, {i, 0, deg[[j+1]]}]opfun[f[n], j], {j, 0, ord}];
       ];       

       p = NextPrime[p, -1]; 
       Time[5, 
         gcdp = PolynomialMod[gcd /. a_Rational :> PolynomialMod[a, p], p];
         If[ checkfun[gcdp, f[n], p], Break[] ];
       ];

       Clear[gcd];
    ];

    Return[gcd];
  ];

Options[GuessMinRE] = {Max -> Automatic, Modulus -> NextPrime[Developer`$MaxMachineInteger, -1], Return -> "rec", StartPoint -> 0, Path -> Automatic, Infolevel -> 0, Order -> Infinity, Degree -> Infinity};

GuessMinRE[data_List, f_[n_], opts___] :=
  Module[ {max, mod, start, ret, p, path, o, d, minord=0, maxord, mindeg=0, maxdeg},

    {max, mod, start, ret, path, maxord, maxdeg} = {Max, Modulus, StartPoint, Return, Path, Order, Degree} 
                             /. {opts} /. Options[GuessMinRE];
    If[ !IntegerQ[max], max = Length[data] ];

    If[ ListQ[maxord], {minord, maxord} = maxord ];
    If[ ListQ[maxdeg], {mindeg, maxdeg} = maxdeg ];

    dataeval = data /. Thread[# -> Table[1/Prime[10+k],{k,1,Length[#]}]]&[Variables[data]];
    p[o_Integer, d_Integer] := 
        Time[1, Module[ {xo, xd, xdim},
        If[ d<mindeg || o<minord || d>maxdeg || o>maxord, Return[0]];
        DebugPrint[3, opts][{o, d, (o+1)(d+1)}];
        xdim = GuessUnivRE[Take[dataeval, Min[Length[data], (o+1)(d+2)+25+Count[data, 0]]], 
               f[n], Order -> o, Degree -> d, Return -> "dim", Modulus -> mod, AdditionalEquations -> 500];
        p[o, d] = xdim;
        If[ xdim == 0, Do[p[xo, xd] = 0, {xo, 0, o}, {xd, 0, d}] ];
        xdim
        ]];

    If[ path === Automatic, 
        path = Join[Table[{ord2, Floor[(max - 2 ord2 - 2)/(ord2 + 1)]}, {ord2, 1, max}],
                    Table[{Floor[(max - deg2 - 2)/(deg2 + 2)], deg2}, {deg2, 1, max}]];
        path = Select[path, First[#]>0&&Last[#]>=0&];
        path = {Max[Cases[path, {p_, #} -> p]], #}& /@ Union[Last /@ path];
        path = {#, Max[Cases[path, {#, p_} -> p]]}& /@ Union[First /@ path];
        path = Sort[path, #1[[1]] < #2[[1]]&]; (* by order *)
        path = Sort[path, (#1[[1]]+1)(#1[[2]]+1) < (#2[[1]]+1)(#2[[2]]+1)&]; (* by size *)
        If[ ret == "fullpath", Return[path] ];
        path = Table[path[[p]], {p, 1, Length[path], Max[1, Floor[Length[path]/10]]}]; 
        path = Join[Most[NestWhileList[Floor[1.7#]&, {2, 3}, (Times@@(#+2)<Length[data])&]], 
                   {Floor[{#, 2#}&[(Sqrt[2 Length[data]]-2)/2]]}, path];
    ];

    If[ ret == "path", Return[path] ];

    check[o_, fun_, p_, hom_:{}] :=
      Module[{i}, 
           For[ i = 1, i <= Min[5, Length[data] - DifferenceOrd[o, fun] - 1]
           && PolynomialMod[PolynomialMod[o /. n -> start + i - 1, p] 
               /. Head[fun][j_] :> PolynomialMod[data[[j - start + 1]] /. hom, p], p] === 0,
           i++];
           i > Min[5, Length[data] - DifferenceOrd[o, fun] - 1]
         ];    

    GuessMinOp[data, f[n], path, mod, ret, p, 
      GuessUnivRE, DifferenceGCD, (Head[#1][First[#1]+#2])&, DifferenceOrd, check, opts]
  ];


Options[GuessMinDE] = {Max -> Automatic, Modulus -> NextPrime[Developer`$MaxMachineInteger, -1], Return -> "rec", StartPoint -> 0, Infolevel -> 0, Path -> Automatic, Order -> Infinity, Degree -> Infinity};

GuessMinDE[data_List, f_[x_], opts___] :=
  Module[ {max, mod, start, ret, p, path, ord, deg, minord=0, maxord, mindeg=0, maxdeg},

    {max, mod, start, ret, path, maxord, maxdeg} = {Max, Modulus, StartPoint, Return, Path, Order, Degree} 
                             /. {opts} /. Options[GuessMinDE];

    If[ ListQ[maxord], {minord, maxord} = maxord ];
    If[ ListQ[maxdeg], {mindeg, maxdeg} = maxdeg ];

    If[ !IntegerQ[max], max = Length[data] ];

    dataeval = data /. Thread[# -> Table[1/Prime[10+k],{k,1,Length[#]}]]&[Variables[data]];
    p[o_Integer, d_Integer] := 
        Time[1, Module[ {xo, xd, xdim},
        If[ d<mindeg || o<minord || d>maxdeg || o>maxord, Return[0] ];
        DebugPrint[3, opts][{o, d, (o+1)(d+1)}];
        xdim = GuessUnivDE[Take[dataeval, Min[Length[data], (o+1)(d+2)+25+Count[data, 0]]], 
               f[x], Order -> o, Degree -> d, Return -> "dim", Modulus -> mod];
        p[o, d] = xdim;
        If[ xdim == 0, Do[p[xo, xd] = 0, {xo, 0, o}, {xd, 0, d}] ];
        xdim
        ]];

    If[ path === Automatic, 
        path = Join[Table[{ord2, Floor[(max - 2 ord2 - 2)/(ord2 + 1)]}, {ord2, 1, max}],
                    Table[{Floor[(max - deg2 - 2)/(deg2 + 2)], deg2}, {deg2, 1, max}]];
        (* 
        path = Join[Table[{ord2, Floor[(max - ord2 - 1)/(ord2 + 2)]}, {ord2, 1, max}],
                    Table[{Floor[(max - 2 deg2 - 1)/(deg2 + 1)], deg2}, {deg2, 1, max}]];
         *)
        path = Select[path, First[#]>0&&Last[#]>0&];
        path = {Max[Cases[path, {p_, #} -> p]], #}& /@ Union[Last /@ path];
        path = {#, Max[Cases[path, {#, p_} -> p]]}& /@ Union[First /@ path];
        path = Sort[path, (#1[[1]]+1)(#1[[2]]+1) < (#2[[1]]+1)(#2[[2]]+1)&];
        If[ ret == "fullpath", Return[path] ];
        path = Table[path[[p]], {p, 1, Length[path], Max[1, Floor[Length[path]/10]]}];
        path = Join[Most[NestWhileList[Floor[1.7*#]&, {2, 3}, (Times@@(#+2)<Length[data])&]], 
                    {Floor[{#, 2*#}&[(Sqrt[8 Length[data] +9 ]-5)/4]]}, path];
    ];

    If[ ret == "path", Return[path] ];

    check[o_, fun_, p_, hom_:{}] :=
         Module[{i, j, ord, acc, rem, deg}, 
           acc = Min[5, Length[data] - (deg = Exponent[o, x]) - (ord = DifferentialOrd[o, fun]) - 1];
           For[ i = acc-1, i <= acc, i++,
             rem = PolynomialRemainder[o /. f -> (Evaluate[Take[data /. hom, ord+deg+i+1].#^Range[0,ord+deg+i]]&),
                      x^(i+deg), x, Modulus -> p];
             If[rem=!=0, Break[]]
           ];
           i > acc
         ];

    GuessMinOp[data, f[x], path, mod, ret, p, 
      GuessUnivDE, DifferentialGCD, D[#1,{First[#1],#2}]&, DifferentialOrd, check, opts]
  ];

GuessMinOp[data_List, f_[n_], path_List, mod_Integer, ret_String, p_, 
   guessfun_, gcdfun_, opfun_, ordfun_, checkfun_, opts___] :=
  Module[ {mid, mypath, ord1, deg1, ord2, deg2, dim, o, d, dd, max, o0, d0},

    {o0, d0} = {Order, Degree} /. {opts};

    Do[T[i]=0, {i, 1, 8}];
    max = Length[data];

    (* point on the hyperbola closest to the form (x,2x), or {o0,d0}, if given *)
    If[ FreeQ[{o0, d0}, Automatic],
      mid = {o0, d0},
      mid = Apply[{#1, Abs[#2-2*#1]}&, path, {1}];
      mid = First[Extract[path, Position[mid, {_, Min[Last /@ mid]}]]];
    ];
    mypath = Rest[Reverse[NestWhileList[Floor[#/2]&, mid, Min[#]>1&]]];
    mypath = Join[mypath, path];

    d = 1; dim = 0;
    While[ dim == 0 && d <= Length[mypath],
      {ord1, deg1} = mypath[[d++]]; dim = p[ord1, deg1]
    ];
    If[ ret === "bool", 
      Return[ dim > 0 ];
    ];
    If[ dim == 0, (* give up. *)
      Throw["No recurrence found."];
    ];
    While[ True, (* minimize degree while order is kept fixed. *)
      For[d = Max[0, deg1 - Max[1, Floor[Sqrt[dim]/2]]], 
         (dim = p[ord1, d]) == 0 && d < deg1, d = Ceiling[(deg1+d)/2]];
      If[d == deg1, Break[], deg1 = d];
    ];
    While[ True, (* minimize order while degree is kept fixed. *)
      For[d = Max[1, ord1 - Max[1, Floor[Sqrt[dim]/2]]], 
         (dim = p[d, deg1]) == 0 && d < ord1, d = Ceiling[(ord1+d)/2]];
      If[d == ord1, Break[], ord1 = d];
    ];
    If[ dim > 1 || ord1 == 1,
      (* nullspace is at least two dimensional, no further pair needed *)
      ord2 = ord1; deg2 = deg1, 
      (* nullspace is one dimensional, search for another pair *)
      mypath = Riffle[Table[{ord1 - 1, d}, {d, deg1+1, Floor[(max+1)/Max[1,ord1]-2]}],
                     Table[{d, deg1 - 1}, {d, ord1+1, Floor[(max+1)/Max[1,deg1]-2]}]];
      (* if the configuration for the second solution corresponds to a system
         of more than twice the size of the first, assume minimality. *)
      mypath = Select[mypath, (Length[#]==2 && (#[[1]]+1)(#[[2]]+1) <= 2(ord1+1)(deg1+1))&];
      For[mypath = DeleteCases[mypath, {_, deg1}], 
          Length[mypath] > 0 && (dim = p@@First[mypath]) == 0,
          mypath = Rest[mypath]];
      If[ Length[mypath] == 0, (* recurrence unique? *)
        If[ ret == "shape", Return[{ord1,deg1}] ];
        Catch[Return[First[guessfun[data, f[n], Order -> ord1, Degree -> deg1, opts]]]];
        Return[GuessUnivOp[data, f[n], ord1, deg1, ord1, deg1, guessfun, gcdfun, opfun, ordfun, checkfun, opts]]; (* parameters *)
      ];
      {ord2, deg2} = First[mypath]; 
      While[ True, (* minimize degree while order is kept fixed. *)
        For[d = Max[0, deg2 - Max[1, Floor[Sqrt[dim]/2]]], 
           (dim = p[ord2, d]) == 0 && d < deg2, d = Ceiling[(deg2+d)/2]];
        If[d == deg2, Break[], deg2 = d];
      ];
    ];
    If[ p[ord2, deg2] > 1, {ord1, deg1} = {ord2, deg2} ];
    If[ ret == "pregcdshape", Return[ {ord1, deg1, ord2, deg2} ]];
    DebugPrint[1, opts]["calling gcrd guesser with parameters ", {ord1, deg1, ord2, deg2}];
    T[6] = Max[(ord1+2)(deg1+1)-1, (ord2+2)(deg2+1)-1];
    d = GuessUnivOp[data, f[n], ord1, deg1, ord2, deg2, guessfun, gcdfun, opfun, ordfun, checkfun, opts];
    DebugPrint[2, opts]["timings +", Array[T, 5], " = ", Plus@@Array[T, 5], "s"];
    Return[d];
  ];

(* algebraic guesser *)
GuessUnivAE[data_List, f_[x_], degf_Integer, degx_Integer, mod_Integer, opts___] :=
  Module[ {sys, d, sol, i, j},

    If[Length[data] < (degf+1)(degx+1), Throw["not enough data."]];
    sys = NestList[
       Table[Sum[PolynomialMod[#[[k+1]]data[[n-k+1]], mod],{k,0,n}], 
                 {n,0,Length[data]-1}]&,
       Prepend[Table[0,{Length[data] - 1}], 1], degf];
    Do[sys = Join[sys, Prepend[Most[#], 0]& /@ Take[sys, -(degf + 1)]], {d, 1, degx}];
    Block[ {LinearSystemSolver`Private`PrintLevel = -1 + Infolevel /. {opts} /. Infolevel -> 0}, 
    sol = If[mod===0, LinSolveQ[Transpose[sys]], 
                      NullSpace[Transpose[sys], Modulus -> mod]];
    ];
    If[ Length[sol]==0, {}, 
      Flatten[Table[x^i f[x]^j,{i,0,degx}, {j,0,degf}]].#& /@ sol
    ]
  ];

Options[GuessMinAE] = {
 Modulus -> NextPrime[Developer`$MaxMachineInteger, -1], 
 Return -> "ae", Infolevel -> 0,
 Order -> Infinity,
 Degree -> Infinity,
 Path -> Automatic
};
GuessMinAE[data_List, f_[x_], opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {d, dx, dy, dim, p, Dx, Dy, mod, ret, path, mindeg=0, maxdeg, minord=0, maxord},

    {mod, ret, maxdeg, maxord, path} = {Modulus, Return, Degree, Order, Path} /. {opts} /. Options[GuessMinAE];

    If[ ListQ[maxord], {minord, maxord} = maxord ];
    If[ ListQ[maxdeg], {mindeg, maxdeg} = maxdeg ];

    Dx = Dy = -1;

    p[dx_Integer, dy_Integer] := p[dx,dy] = Module[{r, i, j},
       DebugPrint[3, opts][{dx,dy}];
       If[ dx < mindeg || dy < minord || dx > maxdeg || dy > maxord, Return[False] ];
       r = GuessUnivAE[Take[data, Min[Length[data],(dy+1)(dx+1)+20]], 
                       f[x], dy, dx, mod, opts];
       If[ Length[r]==0, Do[p[i,j]=False,{i,0,dx},{j,0,dy}]; Return[False] ];
       r = PolynomialGCD@@Append[r, Modulus -> mod];
       If[ r===1, Do[p[i,j]=False,{i,0,dx},{j,0,dy}]; Return[False] ];
       Dx = Exponent[r, x]; Dy = Exponent[r, f[x]];
       Return[True];
    ];

    For[ d = 5, Length[data]>=(d+1)^2+1 && !p[d, d], d = Floor[d 3/2] ];
    If[ Dx < 0, d = Floor[Sqrt[Length[data]]] - 1; p[d, d] ];
    If[ Dx < 0,
       d = Select[Table[{dx, Floor[Length[data]/(dx+2)-1]}, {dx, 2, Length[data]-2}], Min[#]>1&];
       d = Select[Union[d, Reverse /@ d], FreeQ[#, 0|1]&];
       d = {Max[Cases[d, {p_, #} -> p]], #}& /@ Union[Last /@ d];
       d = {#, Max[Cases[d, {#, p_} -> p]]}& /@ Union[First /@ d];
       d = Reverse[Sort[d]];
       d = Table[d[[p]], {p, 1, Length[d], Max[1, Floor[Length[d]/10]]}];
       While[Length[d] > 0 && !p@@First[d], d = Rest[d]];
    ];

    If[ret == "bool", Return[ Dx >= 0 ]];
    If[ret == "shape", Return[ {Dx, Dy} ]];
    If[ret == "mod" && Dx >= 0, Return[ GuessUnivAE[data, f[x], Dy, Dx, mod, opts] ] ];
    If[Dx < 0, Throw["No algebraic equation found."] ];
    d = GuessUnivAE[data, f[x], Dy, Dx, 0, opts];
    If[Length[d] == 0, Throw["No algebraic equation found."]];
    Return[Collect[First[d], f[x]]];
  ];

(* ***************************************** other tools ******************************************** *)

Global`DetermineStaircase[data_, start_] :=
  Module[ {r, d, p},
    p = 2^31 - 1; {r, dUp} = start; dold = Infinity;
    While[ r < Length[data], 
      While[ Length[data] < (r + 1)(dUp + 2), Print["cleared point: ", {r, dUp}]; dold = dUp-- ];
      Print["current point: ", {r, dUp}];
      If[ GuessUnivRE[data, f[n], Order -> r, Degree -> dUp, Return -> "dim"] == 0, r++; Continue[] ];
      Print["current point: ", {r, 0}];
      If[ GuessUnivRE[data, f[n], Order -> r, Degree -> 0, Return -> "dim"] > 0, 
          Print["CORNER POINT: ", {r, 0}]; Return[] ];
      dLow = 0;
      While[ dUp - dLow > 1,
        d = Floor[(dUp + dLow)/2]; 
        Print["current point: ", {r, d}];
        {dUp, dLow} = If[ GuessUnivRE[data, f[n], Order -> r, Degree -> d, Return -> "dim"] > 0,
          {d, dLow}, {dUp, d}]
      ];
      If[ dUp < dold, Print[ "CORNER POINT: ", {r, dUp} ]; dold = dUp ];
      r++; dUp--
    ]
  ];

(* ************************************ multivariate guessers **************************************** *)

InverseCantor[z_Integer] := Module[{w, t}, w = Floor[(Sqrt[8z+1]-1)/2]; t = w(w+1)/2; {w-z+t,z-t}];
Cantor[{x_Integer,y_Integer}] := 1/2(x+y)(x+y+1)+y;
Cantor[{-Infinity,-Infinity}] = -1;

(**
 * 
 **)
Options[BivariateHermitePade] = {Modulus -> 0, Break -> Infinity};
BivariateHermitePade[data_List, sset_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {dims, P, f, vars, x, y, z, term, n, k, sdata, dn, dk, degs, pivot, p, t, stop},

     dims = Dimensions[data];
     {p, stop} = {Modulus, Break} /. {opts} /. Options[BivariateHermitePade];
     f = Union[Head /@ sset];
     If[ Length[f] =!= 1, Throw["illegal structure set"] ];
     f = First[f];
     vars = Union[sset /. _Integer -> 0 /. f -> List];
     If[ Dimensions[vars] =!= {1, 2}, Throw["illegal structure set"] ];
     {n, k} = List@@First[vars];
     dk = Max[0, Cases[sset, f[_, k+Optional[i_Integer]] -> i, Infinity]];
     dn = Max[0, Cases[sset, f[n+Optional[i_Integer], _] -> i, Infinity]];

     P = IdentityMatrix[Length[sset]];
     degs = Table[0, {Length[sset]}]; (* cant-inv of exp. vecs of leading terms *)

     Do[

       If[ z > stop, Break[ ]];

       {x, y} = InverseCantor[z];

       If[ x >= dims[[1]] - n && y >= dims[[2]] - k, Break[] ];
       If[ x >= dims[[1]] - n || y >= dims[[2]] - k, Continue[] ];
       sdata = sset /. f[n+i_.,k+j_.] :> data[[x+i+1,y+j+1]];
       If[ MatchQ[sdata, {0..}], Continue[] ];
       cdata = PolynomialMod[(P /. {n -> x, k -> y}).sdata, p];
       lambda = Complement[Range[Length[P]], Flatten[Position[cdata, 0]]];
       If[ lambda=={}, Continue[] ];
       pivot = First[lambda];
       Do[If[ degs[[lambda[[i]]]] < degs[[pivot]], pivot = lambda[[i]] ], {i, 2, Length[lambda]}];
       (P[[#]] = Expand[P[[#]]cdata[[pivot]] - P[[pivot]]cdata[[#]], Modulus -> p])&
            /@ DeleteCases[lambda, pivot];
       AppendTo[P, Expand[(n-x)P[[pivot]], Modulus -> p]];
       AppendTo[degs, Cantor[InverseCantor[degs[[pivot]]] + {1, 0}]];
       P[[pivot]] = Expand[(k-y)P[[pivot]], Modulus -> p];
       degs[[pivot]] = Cantor[InverseCantor[degs[[pivot]]] + {0, 1}];

       P = Singular`SingularGroebner[P, {n,k}, Modulus -> p, MonomialOrder -> {"dp", "c"}];
       degs = Max /@ Map[Cantor[Exponent[#, {n,k}]]&,
          (Singular`SingularLead[#, {n,k}, MonomialOrder -> {"dp", "c"}]&) /@ P, {2}];

       Print[degs];

     , {z, 0, Times@@dims}];
     
     Return[P.sset];
  ];

Options[GuessBivRE] = {Infolevel -> 0, First -> All, AdditionalEquations -> 40, Modulus -> 2^31-1, Return -> "sol", AdditionalTerms -> {}};
GuessBivRE[data_, sset_List, f_[n_, k_], dn_Integer, dk_Integer, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {dim, nP, info, imgs, a, sys, sol, p, ret, lc, Sn, Sk, supp, ordN},

    {info, nP, add, p, ret, supp} = {Infolevel, First, AdditionalEquations, Modulus, Return, AdditionalTerms}
         /. {opts} /. Options[GuessBivRE];

    ordN = Max[0, Cases[sset, f[n+i_., _] -> i, Infinity]];

    dim = Dimensions[data]; 
    If[ Length[dim] =!= 2 || First[dim] <= 2 dn || Last[dim] <= (dk+1) Length[sset] + add, 
      Throw["not enough data."] ];

    If[ IntegerQ[nP],
      If[ Length[supp] == 0, 
        terms = Flatten[Outer[Times, NestList[# k&, 1, dk], sset /. n -> nP]];
        DebugPrint[2, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]];
        sys = Table[Flatten[Outer[Times, NestList[PolynomialMod[# kP, p]&, 1, dk], 
             sset /. {n -> nP, k -> kP} /. f[i_,j_] :> data[[i+1,j+1]]]]
        , {kP, 15, 15 + Length[terms] + add}];
        DebugPrint[2, opts]["starting to solve modular system of size ", Dimensions[sys], " ", MyDate[], ""];
        sol = NullSpace[sys, Modulus -> p];
        DebugPrint[2, opts][Length[sol], " solutions predicted ", MyDate[], ""];
        If[ ret == "dim", Return[Length[sol]] ];
        If[ Length[sol] == 0, Return[ {} ] ];
        If[ ret == "mod", Return[Collect[PolynomialMod[#.terms, p], f[__]]& /@ sol] ];
        terms = DeleteCases[(Min[#, 1]& /@ Plus@@Map[If[#===0, 0, 1]&, sol, {2}]) * terms, 0];
        If[ ret == "support", Return[terms /. Thread[(sset /. n -> nP) -> sset] ] ];
      , terms = supp /. n -> nP
      ];
      DebugPrint[2, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]];
      sys = Table[terms /. k -> kP /. f[i_,j_] :> data[[i+1,j+1]], {kP, 15, 15 + Length[terms] + add}];
      DebugPrint[2, opts]["starting to solve refined system of size ", Dimensions[sys], " ", MyDate[], ""];
      sol = LinSolveQ[sys];
      DebugPrint[2, opts][Length[sol], " solutions ", MyDate[], ""];
      sol = Collect[#.terms, f[__]]& /@ sol;
      lc = Last[CoefficientList[Last[sol] /. f[__] -> 1, k]];
      Return[#/lc& /@ sol];      
    ];

    DebugPrint[1, opts]["computing support ", MyDate[], ""];
    supp = GuessBivRE[data, sset, f[n,k], dn, dk, First -> Length[data] - ordN - 1, Return -> "support", opts];
    DebugPrint[1, opts]["done. ", MyDate[], ""];
    If[ Length[supp] == 0, Return[supp] ];
    DebugPrint[1, opts]["computing images..."];
    imgs = Table[
       GuessBivRE[data, sset, f[n,k], dn, dk, First -> nP, AdditionalTerms -> supp, opts]
       /. Thread[(sset /. n -> nP) -> sset] /. f[n+i_.,k+j_.] -> Sn^i Sk^j
    , {nP, Length[data] - ordN - 2 dn - 5, Length[data] - ordN - 1}];
    If[ Length[Union[Length /@ imgs]]=!=1, Throw["unlucky evaluation points: " <> ToString[Length /@ imgs]]];
    points = Range[Length[data] - ordN - 2 dn - 5, Length[data] - ordN - 1];
    DebugPrint[1, opts]["done. ", MyDate[], ""];
    DebugPrint[1, opts]["Lifting..."];
    sol = Table[Lift[imgs[[All, i]], points, {Sn,Sk,k}, n, 0], {i, 1, Length[imgs[[1]]]}];
    DebugPrint[1, opts]["done. ", MyDate[], ""];
    sol = CoefficientList[#, {Sn,Sk}]& /@ sol;
    sol = (Plus@@Flatten[# Table[f[n+i,k+j],{i,1,Length[#]},{j,1,Length[#[[1]]]}]])& /@ sol;
    sol = Collect[Numerator[Together[#]], f[__]]& /@ sol;
    Return[ sol ];
  ];


(**
 * Guessing of symmetric recurrence equations
 * 
 * GIVEN:
 * -> a d-dimensional array data with symmetric content 
 *    (i.e. data[[i,j,k,...]]==data[[any-permuation[i,j,k,...]]])
 * -> a term f[n,m,k...] specifying function symbol and variable names
 * -> a bound on the (individual variable) order
 * -> a bound on the (total) degree
 *
 * FIND:
 * -> all multivariate recurrences of order and degree not exceeding the 
 *    specified bounds and totally symmetric in all variables and satisfied
 *    by the given data.
 **)
Options[GuessSymmetricRE] = {Return -> "sol", AdditionalEquations -> 100, Modulus -> 2^31-1};
GuessSymmetricRE[data_, f_[nn__], order_, degree_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {dims, vars, coeffs, sset, terms, add, i, j, ret, sys, sol},

    dims = Dimensions[data]; vars = {nn}; d = Length[vars];
    If[ !OrderedQ[dims], Throw["please transpose such that Dimensions[data] is descending."] ];
    If[ Length[dims] =!= Length[vars], Throw["dimension missmatch"] ];
    If[ !IntegerQ[order] || !IntegerQ[degree], Throw["order and degree must be specific integers"] ];
    {ret, add, mod} = {Return, AdditionalEquations, Modulus} /. {opts} /. Options[GuessSymmetricRE];

    coeffs = Table@@(Prepend[Table[{j[i],0,degree-Sum[j[k],{k,1,i-1}]},{i,1,d}], Times@@(vars^Array[j, d])]);
    coeffs = Flatten[coeffs];
    sset = Table@@(Prepend[Table[{j[i],0,j[i-1]},{i,1,d}], Times@@(vars^Array[j, d])] /. j[0] -> order); 
    sset = Expand[sset Times@@vars] /. Times@@(vars^Table[Optional[Pattern[
     Evaluate[Symbol["j" <> ToString[j]]], Blank[]]], {j, 1, d}]) -> f@@(vars + Table[Symbol["j" <> ToString[j]], {j, 1, d}] - 1);
    sset = Flatten[sset];
    terms = Flatten[Outer[Times, coeffs, sset]];
    terms = MySymmetrize[#, f[nn]]& /@ terms;
    terms = Union[terms];

    DebugPrint[1, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]]; 

    If[ret == "terms", Return[terms]];

    points = Flatten[Table@@(Prepend[Table[{j[i], 0, dims[[i]] - order - 1}, {i, 1, d}], f@@Array[j, d]])];

    DebugPrint[1, opts][Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]]; 

    If[ Length[terms] + add > Length[points], Throw["not enough data"] ];

    points = RandomKSubset[points /. f -> List, Length[terms] + add];

    DebugPrint[1, opts]["keeping ", Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]]; 

    If[ret == "points" , Return[points]];

    sys = Table[
      PolynomialMod[terms /. Thread[vars -> points[[i]]], mod] /. f[j__] :> PolynomialMod[data[[##]]&@@Sort[{j}+1], mod]
    , {i, 1, Length[points]}]; 

    DebugPrint[1, opts]["Starting to solve modular system ", If[Length[terms] > 1000, MyDate[], ""]]; 

    sol = NullSpace[sys, Modulus -> mod];

    DebugPrint[1, opts][Length[sol], " solutions predicted."];

    If[ ret == "dim", Return[Length[sol]] ];

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

    If[ ret == "mod", Return[Collect[PolynomialMod[#.terms, mod], f[__]]& /@ sol] ];

    terms = DeleteCases[(Min[#, 1]& /@ Plus@@Map[If[#===0, 0, 1]&, sol, {2}]) * terms, 0];

    If[ ret == "support", Return[terms] ];

    sys = Table[
      terms /. Thread[vars -> points[[i]]] /. f[j__] :> (data[[##]]&@@Sort[{j}+1])
    , {i, 1, Length[points]}];

    If[ ret == "sys", Return[sys] ];

    DebugPrint[1, opts]["Starting to solve system ", If[Length[terms] > 1000, MyDate[], ""]]; 

    (* TODO: select appropriate solver *)

    sol = LinSolveQ[sys];

    DebugPrint[1, opts][Length[sol], " solutions "];

    Return[Collect[#.terms, f[__]]& /@ sol];
  ];

(* input: smth like n*f[n+1,m]. output: smth like n*f[n+1,m] + m*f[n,m+1] *)
MySymmetrize[term_, f_[nn__]] := 
  Module[ {vars, d, P, res, i},
    P = Permutations[vars = {nn}]; d = Length[vars];
    Sum[term /. Thread[vars -> P[[i]]], {i, 1, Length[P]}] /. f[i__] :> f@@Sort[{i}]
  ];

(**
 * Multivariate Guessing with order and degree 
 * 
 * GIVEN:
 * -> a d-dimensional array data 
 * -> a term f[n,m,k...] specifying function symbol and variable names
 * -> a bound on the (total) order
 * -> a bound on the (total) degree
 *
 * FIND:
 * -> all multivariate recurrences of order and degree not exceeding the 
 *    specified bounds and satisfied by the given data.
 **)
GuessMultRE[data_, f_[nn__], order:({__Integer}|_Integer), degree_Integer, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {params, dims, vars, coeffs, sset, terms, add, i, j, ret, sys, sol, mTimes, ord, deg, filter, except, evalvar, null, ff}, 

    dims = Dimensions[data]; vars = {nn}; d = Length[vars]; params = Variables[data];
    If[ Length[dims] =!= Length[vars], Throw["dimension missmatch"] ];
    {ret, add, mod, otype, dtype, filter, except, null} = 
     {Return, AdditionalEquations, Modulus, Order, Degree, SelectPoints, SelectTerms, StartPoint} /. {opts} /. Options[GuessMultRE];
    otype = If[otype=="total", 1, 0];
    dtype = If[dtype=="total", If[Length[#]>0,Most[#],#]&, Identity];
    ord = Which[
      ListQ[order] && otype==0, order,
      ListQ[order] && otype==1, Throw["nonsense: individual orders in total order mode."],
      IntegerQ[order], Table[order, {d}]
    ];

    coeffs = Flatten[Fold[Function[{coeffs, i}, NestList[DeleteCases[Map[vars[[i]]dtype[#]&, #, {i-2}], {}]&, 
           coeffs, degree]], NestList[#*vars[[1]]&, 1, degree], Range[2, d]]];

    sset = Flatten[Table@@(Prepend[Table[{j[i],0,ord[[i]]-otype Sum[j[u], {u,1,i-1}]},{i,1,d}], f@@(vars+Array[j, d])])]; 
    sset = Select[sset, except];
    terms = Flatten[Outer[Times, coeffs, sset]];

    DebugPrint[1, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]]; 

    If[ret == "terms", Return[terms]];

    points = Flatten[Table@@(Prepend[Table[{j[i], 0, dims[[i]] - ord[[i]] - 1}, {i, 1, d}], ff@@Array[j, d]])];
    points = Select[Rest /@ Union[Prepend[#, Plus@@#]& /@ points /. ff -> List], filter];

    DebugPrint[1, opts][Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]];

    If[ Length[terms] + add > Length[points], Throw["not enough data"] ];

    points = Take[points, Length[terms] + add];

    DebugPrint[1, opts]["keeping ", Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]]; 

    If[ret == "points" , Return[points]];

    mTimes = PolynomialMod[#1*#2, mod]&;
    mdata = PolynomialMod[data /. Thread[params -> 123 + Range[Length[params]]], mod];
    sys = Union[Table[
      Flatten[Outer[mTimes, 
        Flatten[Fold[Function[{coeffs, i}, NestList[DeleteCases[Map[PolynomialMod[points[[k, i]]dtype[#], mod]&,
          #, {i - 2}], {}]&, coeffs, degree]], PolynomialMod[NestList[#*points[[k, 1]]&, 1, degree], mod], 
          Range[2, d]]], sset /. Thread[{nn} -> points[[k]]] /. f[j__] :> (mdata[[##]]&@@({j}+1))]]
    , {k, 1, Length[points]}]];

    If[ Length[sys] < Length[First[sys]], Throw["Underdetermined system: " <> ToString[Dimensions[sys]]]];

    DebugPrint[1, opts]["Starting to solve modular system of size ", Length[sys], "x", Length[First[sys]], "   ", If[Length[terms] > 1000, MyDate[], ""]]; 

    sol = NullSpace[sys, Modulus -> mod];

    DebugPrint[1, opts][Length[sol], " solutions predicted."];

    If[ ret == "dim", Return[Length[sol]] ];

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

    If[ ret == "mod", Return[Collect[PolynomialMod[#.terms, mod], f[__]]& /@ sol] ];

    terms = DeleteCases[(Min[#, 1]& /@ Plus@@Map[If[#===0, 0, 1]&, sol, {2}]) * terms, 0];
    
    If[ ret == "support", Return[terms] ];

    sys = Table[
      terms /. Thread[vars -> points[[i]]] /. f[j__] :> (data[[##]]&@@({j}+1))
    , {i, 1, Length[points]}];

    If[ ret == "sys", Return[sys] ];

    DebugPrint[1, opts]["Starting to solve system ", If[Length[terms] > 1000, MyDate[], ""]]; 

    (* TODO: select appropriate solver *)
    Which[ 
      Length[params] == 0, sol = LinSolveQ[sys],
      Length[params] == 1, sol = LinSolveUniv[sys, 0, params[[1]]],
      True, Throw["no solver available"];
    ];

    DebugPrint[1, opts][Length[sol], " solutions "];

    Return[(Collect[#.terms, f[__]]& /@ (sol /. Thread[{nn}->{nn}-null])) /. f[j__] :> (f[j] /. Thread[{nn}->{nn}+null]) /. Automatic -> 0]
  ];

(**
 * Multivariate Guessing with order and degree 
 * 
 * GIVEN:
 * -> a d-dimensional array data 
 * -> a term f[n,m,k...] specifying function symbol and variable names
 * -> a bound on the (total) order
 * -> a bound on the (total) degree
 *
 * FIND:
 * -> all multivariate pdes of order and degree not exceeding the 
 *    specified bounds and satisfied by the given data.
 **)
GuessMultDE[data_, f_[xx__], order:({__Integer}|_Integer), degree:({__Integer}|_Integer), opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {params, dims, vars, coeffs, sset, terms, add, i, j, ret, sys, sol, mTimes, ord, deg, filter, except, evalvar, null, ff, poch, ii, jj, kk}, 

    dims = Dimensions[data]; vars = {xx}; d = Length[vars]; params = Variables[data];
    If[ Length[dims] =!= Length[vars], Throw["dimension missmatch"] ];
    {ret, add, mod, otype, dtype, filter, except, null} = 
     {Return, AdditionalEquations, Modulus, Order, Degree, SelectPoints, SelectTerms, StartPoint} /. {opts} /. Options[GuessMultRE];
    otype = If[otype=="total", 1, 0];
    ord = Which[
      ListQ[order] && otype==0, order,
      ListQ[order] && otype==1, Throw["nonsense: individual orders in total order mode."],
      IntegerQ[order], Table[order, {d}]
    ];
    dtype = If[dtype=="total", 1, 0];
    deg = Which[
      ListQ[degree] && dtype==0, degree,
      ListQ[degree] && dtype==1, Throw["nonsense: individual orders in total order mode."],
      IntegerQ[degree], Table[degree, {d}]
    ];

    (*
    coeffs = Flatten[Fold[Function[{coeffs, i}, NestList[DeleteCases[Map[vars[[i]]dtype[#]&, #, {i-2}], {}]&, 
           coeffs, degree]], NestList[#vars[[1]]&, 1, degree], Range[2, d]]];
     *)
    coeffs = Flatten[Table@@(Prepend[Table[{j[i],0,deg[[i]]-dtype Sum[j[u],{u,1,i-1}]},{i,1,d}], Product[vars[[i]]^j[i],{i,1,d}]])];

    sset = Flatten[Table@@(Prepend[Table[{j[i],0,ord[[i]]-otype Sum[j[u], {u,1,i-1}]},{i,1,d}], Product[Der[vars[[i]]]^j[i], {i, 1, d}]])]; 
    sset = Select[sset, except];
    terms = Flatten[Outer[Times, coeffs, sset]];

    DebugPrint[1, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]]; 

    If[ret == "terms", Return[terms]];

    points = Flatten[Table@@(Prepend[Table[{j[i], 0, dims[[i]] - ord[[i]] - 1}, {i, 1, d}], ff@@Array[j, d]])];
    points = Select[Rest /@ Union[Prepend[#, Plus@@#]& /@ points /. ff -> List], filter];

    DebugPrint[1, opts][Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]];

    If[ Length[terms] + add > Length[points], Throw["not enough data"] ];

    points = Take[points, Length[terms] + add];

    DebugPrint[1, opts]["keeping ", Length[points], " points ", If[Length[points] > 1000, MyDate[], ""]]; 

    If[ret == "points" , Return[points]];

    mTimes = PolynomialMod[#1#2, mod]&;
    mdata = PolynomialMod[data /. Thread[params -> 123 + Range[Length[params]]], mod];
    poch[i_Integer, j_Integer] := If[j==0, 1, poch[i, j] = mTimes[i, poch[i - 1, j - 1]]];
    poch[i_List, j_List] := poch[i, j] = Product[poch[i[[k]], j[[k]]], {k, 1, d}];
    sys = DeleteCases[Union[Table[
      kk = points[[k]];
      Table[
        ii = Exponent[terms[[i]], vars];
        jj = Exponent[terms[[i]], Der /@ vars];
        If[ Min[kk + jj - ii] < 0, 0, (mdata[[##]]&@@(kk + jj - ii + 1))*poch[kk + jj - ii, jj]]
      , {i, 1, Length[terms]}]
    , {k, 1, Length[points]}]], {0..}];
    ClearAll[poch];

    If[ Length[sys] < Length[First[sys]], Throw["Underdetermined system: " <> ToString[Dimensions[sys]]]];

    DebugPrint[1, opts]["Starting to solve modular system of size ", Length[sys], "x", Length[First[sys]], "   ", If[Length[terms] > 1000, MyDate[], ""]]; 

    sol = NullSpace[sys, Modulus -> mod];

    DebugPrint[1, opts][Length[sol], " solutions predicted."];

    If[ ret == "dim", Return[Length[sol]] ];

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

    If[ ret == "mod", Return[PolynomialMod[#.terms, mod]& /@ sol] ];

    terms = DeleteCases[(Min[#, 1]& /@ Plus@@Map[If[#===0, 0, 1]&, sol, {2}]) * terms, 0];
    
    If[ ret == "support", Return[terms] ];

    poch[i_Integer, j_Integer] := If[j==0, 1, poch[i, j] = i*poch[i - 1, j - 1]];
    poch[i_List, j_List] := poch[i, j] = Product[poch[i[[k]], j[[k]]], {k, 1, d}];
    sys = DeleteCases[Union[Table[
      kk = points[[k]];
      Table[
        ii = Exponent[terms[[i]], vars];
        jj = Exponent[terms[[i]], Der /@ vars];
        If[ Min[kk + jj - ii] < 0, 0, (data[[##]]&@@(kk + jj - ii + 1))*poch[kk + jj - ii, jj] ]
      , {i, 1, Length[terms]}]
    , {k, 1, Length[points]}]], {0..}];
    ClearAll[poch];

    If[ ret == "sys", Return[sys] ];

    DebugPrint[1, opts]["Starting to solve system ", If[Length[terms] > 1000, MyDate[], ""]]; 

    (* TODO: select appropriate solver *)
    Which[ 
      Length[params] == 0, sol = LinSolveQ[sys],
      Length[params] == 1, sol = LinSolveUniv[sys, 0, params[[1]]],
      True, Throw["no solver available"];
    ];

    DebugPrint[1, opts][Length[sol], " solutions "];

    Return[(#.terms&) /@ sol];
  ];

(**
 * 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 -> NextPrime[Developer`$MaxMachineInteger,-1], 
  (* 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,
                      or "terms" to return the list of terms in the original (unrefined) ansatz,
                      or "support" to return the predicted structure set,
                      or "except" to return the complement of the predicted structre set,
                      or "dim" to return the predicted dimension of the solution space,
                      or "mod" to return the modular result (possible without parameters only). *)
  AdditionalEquations -> 100, StartPoint -> Automatic, Extension -> {}, Except -> {},
  AdditionalTerms -> {}, Infolevel -> 0, 
  Order -> "individual", Degree -> "total", Point -> All, 
  SelectTerms -> (True&), (* for the other version *) SelectPoints -> (True&), (* for the other version *)
  MustHaveOneOf -> All, Factor -> True, Sort -> (OrderedQ[{#1,#2}]&), Constraints -> True};
GuessMultRE[data_List, sset_List, vars_List, deg:(_Integer|_List), opts:(_Rule|_RuleDelayed)...] :=
  Module[ {i, j, k, mod, dim, hsset, myvars, f, args, terms, sys, sol, p, dom, ret, 
           forced, termorder, s, mdata, mterms, termeval, inhomog,
           actualTerms, e, spols, params, addEqs, blocks, null, ext, except, points, constr},

    {p, addEqs, addTerms, null, ext, except, ret, fac, forced, termorder, constr} = 
         {Modulus, AdditionalEquations, AdditionalTerms, StartPoint, Extension, Except, 
          Return, Factor, MustHaveOneOf, Sort, Constraints} /. {opts} /. Options[GuessMultRE];
    If[ p!=0 && !PrimeQ[p], Throw["Modulus is not a prime!"]];
    If[ p==0 && (forced =!= All || MemberQ[{"dim", "null", "mod", "except", "support"}, ret]),
       Throw["Forced terms and some special return values are inconsistent with modulus 0."]
    ];
    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[ !ListQ[vars], Throw["unexpected variable declaration"] ];

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

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

    If[ !MatchQ[deg, {_Integer...}] && !IntegerQ[deg],
      Throw["degree specification must be an integer or a list of integers."];
    ];

    (* 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 *)
    (* a structure set must be given even if all terms are supplied via the AdditionalTerms option.
       in the latter case, one may declare a negative degree. *)
    If[ !MatchQ[hsset, {f_[__]..}] ||
        Union[Length /@ hsset] =!= {Length[myvars]} || 
        !MatchQ[Flatten[(#-myvars)& /@ Apply[List, hsset, {1}]], {_Integer...}],
      Throw["illegal structure set specification"];
    ];
    f = hsset[[1, 0]]; (* function symbol in the structure set *)

    (* c. check that the elements of addTerms have the right form *)
    If[ !And@@(MatchQ[#, a_. f[__] /; PolynomialQ[a, Flatten[vars]]]& /@ addTerms) || 
        !MatchQ[Flatten[(#-myvars)& /@ Apply[List, Flatten[Cases[addTerms, f[a__], Infinity]], {1}]], 
                {_Integer...}],
      Throw["illegal additional terms specified"];
    ];

    (* 2. ansatz *)
    (* a. terms for the polynomial coefficients *)
    terms = 
      If[ IntegerQ[deg],
        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 = Union[Flatten[Outer[Times, terms, sset]], addTerms];
    terms = Complement[terms, except]; (* kick out unwanted terms *)
    
    inhomog = Select[terms, FreeQ[#, f]&];

    (* kick out unwanted cones [todo: this code only works if no symmetries are exploited] *)
    {terms, except} = {terms, except} /. f[a__] :> Times@@((s /@ myvars)^({a} - myvars));
    Do[
      If[ Head[except[[i]]] === Cone, terms = Select[terms, (Denominator[#/except[[i, 1]]]=!=1)&] ]
    , {i, 1, Length[except]}];
    terms = (#*f@@(myvars + Exponent[#, s /@ myvars])& /@ terms) /. s[_] -> 1;

    terms = Join[inhomog, terms];
    
    If[ forced === All, forced = terms ];
    forced = Intersection[terms, forced]; 
    If[ Length[forced] == 0, Throw[ "cannot force terms from empty set."] ];
    terms = Join[Sort[Complement[terms, forced], termorder], 
                 Sort[forced, termorder]]; (* put forced terms at the end *)
    DebugPrint[1, opts][Length[terms], " terms ", If[Length[terms] > 1000, MyDate[], ""]]; 

    If[ ret == "terms", Return[terms] ];

    (* now every term[[i]] has the form n^a k^b ... f[c n + d k + .., ..] *)
    (* where n^a etc may also be q^(n a) etc., or elementary polynomials in those, if symmetries were 
       specified *)

    (* 3. determine the domain of points for which all arguments in the
          structure set are within the range supplied in the data array *)
    hsset = Union[hsset, Cases[addTerms, f[__], Infinity]];
    dom = Table[0 <= hsset[[i,j]] < Length[Nest[First, data, j - 1]], {i, 1, Length[hsset]}, {j, 1, dim}];
    If[ myvars=!=(vars/.a_^b_->b),
      dom = Join[dom, Apply[LessEqual, vars, {1}]] /. a_^b_ -> b;
    ];
    dom = CylindricalDecomposition[Append[Flatten[dom], constr], myvars]; 
    If[ dom === False, Throw[ "insufficient input data." ]];
    dom = GenerateTuples[dom, myvars];
    Which[ 
      Sum[KroneckerDelta[data[[Sequence@@(1+dom[[i]])]]], {i,1,Length[dom]}]/Length[dom] < 1/100,
        Null,
      Length[dom] < 125000,
        DebugPrint[1, opts]["collecting nonzero points..."];
        e[point__] := data[[Sequence@@({point}+1)]];
        dom = Select[dom, !MatchQ[Apply[e, hsset /. Thread[myvars->#], {1}], {0..}]&],
      True,
        DebugPrint[0, opts]["Warning: There are many zeros in the given data."];
        DebugPrint[0, opts]["Consider using the option Constraints."]
    ];
    If[ Length[dom] < Length[terms], Throw[ "insufficient input data." ]];
    
    (* 4. construct modular system *)
    params = Union@@(Variables /@ Flatten[data]);

    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, 

      mod = Thread[params -> Table[Random[Integer, {1, 20}], {i, 1, Length[params]}]];
      mdata = PolynomialMod[data /. mod, p];
      mterms = Reverse[Sort[Union[terms /. f[__] -> 1, vars]]] /. mod;
      mterms1 = If[p === 0, mterms, mterms /. Power[a1_, a2_] :> PowerMod[a1, a2, p]];

      sys = DeleteCases[
        (terms /. mod
               /. Thread[mterms -> PolynomialMod[mterms1 /. Thread[myvars -> dom[[#]]], p]]
               /. Thread[myvars -> dom[[#]]]
               /. f[j__] :> mdata[[Sequence@@({j}+1)]])& /@ points, {0..}];

      (* 5. solve modular system *)
      DebugPrint[1, opts]["modular system: ", Length[sys], " eqns, ", Length[terms], " vars ",
            If[Length[terms] > 1000, MyDate[], ""]];
      If[ Length[sys] < Length[terms], Throw["too many zeros."] ];

      sol = NullSpace[sys, Modulus -> p];
      If[ Length[terms] > 1000, DebugPrint[1, opts]["solved: ", MyDate[]] ];
      sol = DeleteCases[sol, Prepend[Table[0, {Length[forced]}], ___]];
      If[ ret == "dim", Return[ Length[sol] ] ];
      If[ Length[sol] == 0, DebugPrint[1, opts]["no solution."]; Return[ {} ] ];
    
      DebugPrint[1, opts][Length[sol], " solutions predicted. ", If[Length[actualTerms] > 1000, MyDate[], ""]];
      If[ ret == "null", Return[ Null ] ];

      If[ ret == "mod",
        Return[Collect[#, f[__]]& /@ (sol . terms /. Thread[myvars->myvars-null])
                               /. f[j__] :> (f[j] /. Thread[myvars->myvars+null])
        ];
      ];

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

      If[ ret == "except", Return[Complement[terms, actualTerms]]];
      If[ ret == "support", Return[actualTerms] ];
    ,
      actualTerms = terms;
    ];

    Which[ 
      Sum[KroneckerDelta[data[[Sequence@@(1+dom[[i]])]]], {i,1,Length[dom]}]/Length[dom] < 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[myvars->#], {1}], {0..}]&],
      True, 
        DebugPrint[0, opts]["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 *)
    DebugPrint[1, opts]["refined system: ", Length[sys], " eqns, ", Length[actualTerms], " vars ", 
        If[Length[actualTerms] > 1000, MyDate[], ""]];

    Block[ {LinearSystemSolver`Private`PrintLevel = -3 + Infolevel /. {opts} /. Infolevel -> InfoLevel}, 
    sol = Which[ 
            ext =!= {}, LinearSystemSolver`LinSolveAlg[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]}], ___]];
    DebugPrint[1, opts][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[i, v], {i, 1, Length[v]}];
    p = IntegerPartitions[d, Infinity, Range[Min[d, Length[v]]]];

    Product[e[[#[[i]]]], {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];

