(* written by Manuel Kauers and Burkhard Zimmermann *)

(* Research Institute for Symbolic Computation (RISC) *)
(* J. Kepler University Linz, Austria *)

(* http://www.risc.uni-linz.ac.at/research/combinat/risc/software/ *)
(* email: manuel@kauers.de *)

(* ************************************************************************* *)
(*                              Dependencies.m                               *)
(*         Computing Algebraic Dependencies of C-finite Sequences            *)
(* ************************************************************************* *)

(* INSTALLATION:                                                             *)
(*                                                                           *)
(*  Put this file into some directory where Mathematica is able to           *)
(*  find it.                                                                 *)
(*                                                                           *)
(*                                                                           *)
(* BUGS:                                                                     *)
(*                                                                           *)
(*  Please report bugs concerning this package Manuel Kauers                 *)
(*  (manuel@kauers.de)                                                       *)
(*  This package was developped for Mathematica 5. We do not expect          *)
(*  that it works with earlier versions.                                     *)
(*                                                                           *)
(*                                                                           *)
(* 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.                             *)

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

Dependencies`Private`Version = "0.15 (2006-02-22)";

(** TODO ****
 *
 * - FindRelations: Rational functions by homomorphic images and intersection of results. 
 *
 * MMA-Bug: GroebnerBasis with levandovskyy order gives the annihilator {1}
 * upon input Dependencies[{Pell[k^2+n],Pell[n],Pell[k^2-1],(-1)^k}, x]
 * while Faugere doesn't. MMA's DRL also finds a basis =!= {1}, but Lex and
 * block DRL via Matrix do not terminate. 
 * 
 ************)

BeginPackage["RISCComb`"];

 (* declare symbols of special sequences and options *)
 Unprotect[`Pell, `PellLucas, `Lucas, `Perrin, `Where, `Variable];

 Pell::usage = "Pell[n] gives the nth Pell number.";
 PellLucas::usage = "PellLucas[n] gives the nth Pell-Lucas number.";
 Lucas::usage = "Lucas[n] gives the nth Lucas number.";
 Perrin::usage = "Perrin[n] gives the nth Perrin number.";

 Where::usage = "Option to Dependencies"; 
 Variable::usage = "Option to Dependencies";

 Protect[Pell, PellLucas, Lucas, Perrin, Variable]; 

EndPackage[];

BeginPackage["Dependencies`", {"NumberTheory`AlgebraicNumberFields`", "NumberTheory`ContinuedFractions`", "RISCComb`"}];

Unprotect[`Dependencies, Where];

Dependencies::usage = "Dependencies[{f1,f2,...},x, " <>
  "Where -> {definitions}] computes the ideal of algebraic " <>
  "dependencies of f1,f2,...";
Express::usage = "Express[f, {f1,f2,...}, x, " <>
  "Where -> {definitions}] computes representations of f in " <>
  "terms of f1,f2,...";

{`FindRelations}; 

Begin["`Private`"];

DebugPrint[a___] := Null; (* Print[a]; *)

(** define special sequences **)

Unprotect[Pell, PellLucas, Lucas, Perrin];

Pell[0] = 0;
Pell[1] = 1;
Pell[n_Integer] /; n > 0 := (Pell[n] = 2*Pell[n-1] + Pell[n-2]);
Pell[n_Integer] /; n < 0 := (Pell[n] = -2*Pell[n+1] + Pell[n+2]);

PellLucas[0] = 2;
PellLucas[1] = 2;
PellLucas[n_Integer] /; n > 0 := (PellLucas[n] = 2*PellLucas[n-1] + PellLucas[n-2]);
PellLucas[n_Integer] /; n < 0 := (PellLucas[n] = -2*PellLucas[n+1] + PellLucas[n+2]);

Lucas[0] = 2;
Lucas[1] = 1;
Lucas[n_Integer] /; n > 0 := (Lucas[n] = Lucas[n-1] + Lucas[n-2]);
Lucas[n_Integer] /; n < 0 := (Lucas[n] = -Lucas[n+1] + Lucas[n+2]);

Perrin[0] = 3;
Perrin[1] = 0;
Perrin[2] = 2;
Perrin[3] = 3;
Perrin[4] = 2;
Perrin[5] = 5;
Perrin[n_Integer] /; n < 0 := (Perrin[n] = -Perrin[n+1] + Perrin[n+3]);
Perrin[n_Integer] /; n > 0 && EvenQ[n] := Module[{m}, m=n/2; Perrin[n] = 
  (5*Perrin[m]^2 - 12*Perrin[m]*Perrin[1 + m] - 2*Perrin[1 + m]^2 + 
   8*Perrin[m]*Perrin[2 + m] + 18*Perrin[1 + m]*Perrin[2 + m] - 6*Perrin[2 + m]^2)/23];
Perrin[n_Integer] /; n > 0 && OddQ[n] := Module[{m}, m=(n-1)/2; Perrin[n] = 
  (4*Perrin[m]^2 + 18*Perrin[m]*Perrin[1 + m] + 3*Perrin[1 + m]^2 - 
   12*Perrin[m]*Perrin[2 + m] - 4*Perrin[1 + m]*Perrin[2 + m] + 9*Perrin[2 + m]^2)/23];

Protect[Pell, PellLucas, Lucas, Perrin];

(******************************************************)
(* Part 0: Expression manipulation tools              *)
(******************************************************)

(**
 * Given a system of linear recurrences, extract the dependent variable in which it is given.
 * If the dependent variable cannot be determined (e.g. if the system is empty), then a new 
 * symbol will be returned. If more than one dependent variable is found, an exception is thrown.
 *
 * It is assumed that the system is given in the form as required by SolveRecurrence; see
 * the spec of that function for details. 
 **)
ExtractDependentVariable[system_List] :=
  Module[ {n, i},

  If[ system === {}, Return[n] ];

  n = Union[Cases[system, _[n_+Optional[i_Integer]] -> n, Infinity]];
  n = DeleteCases[n, _Integer];

  Which[ 
    Length[n] === 1, Return[First[n]],
    Length[n] < 1, Return[Unique[]],
    Length[n] > 1, Throw["Definitions must be given in one dependent variable only."]
  ];
];

(**
 * Extracts subexpressions of a given expression. 
 * The expression is parsed with respect to the following grammar.
 *
 *   E := E op E | P | C | Ex | r | E^I
 *   op := + | * | - | / 
 *   Ex := r ^ P
 *   C := F[P] 
 *   P := I | S | P^N | P * P | P + P
 *   r := any rational number
 *   F := any expression with Head[F] not Plus, Times, Power
 *   I := any integer
 *   N := any natural number
 *   S := any expression with Head[S] not Plus, Times, Power, Integer
 *
 * (Any Mathematica expression belongs to this language.)
 * The method extracts the values of a specified nonterminal symbol.
 * The list of all subexpressions corresponding to the given symbol will be returned.
 *
 * The symbols "S", "F", "V", "C", "Ex" may be queried. They should
 * be given as strings. 
 **)
ExtractNodes[expr_, node_String] := ParseE[expr, node];

ParseE[expr_, node_] := Module[ {r},

   If[ NumberQ[N[expr]], Return[{}] ];

   If[ MemberQ[{Plus, Times}, Head[expr]], (* E := E + E, E * E *)
      Return[ Union@@(ParseE[#, node]& /@ List@@expr) ];
   ];

   If[ MatchQ[ expr, Power[_, _Integer] ], (* E := E^I *)
      Return[ ParseE[First[expr], node] ]
   ];

   If[ MatchQ[ expr, _[_] ], (* E := C *)

      r = ParseP[First[expr], node];
      If[ node === "C", Return[{expr}], Return[r]];
   ];

   If[ MatchQ[ expr, Power[a_, _] /; NumberQ[N[a]] ], (* E := Ex *)

      r = ParseP[Last[expr], node];
      If[ node === "Ex", Return[{expr}], Return[r]];
   ];

   If[ Element[expr, Rationals], Return[{}] ]; (* E := r *)

   Return[ParseP[expr, node]] (* E := P *)
];

ParseP[expr_, node_] := Module[ {},

   If[ IntegerQ[expr], Return[{}] ]; (* P := I *)

   If[ MatchQ[expr, _^_Integer] && Last[expr] > 0, (* P := P^N *)
       Return[ParseP[First[expr], node]]
   ];

   If[ MemberQ[{Plus, Times}, Head[expr]], (* P := P + P, P * P, *)
       Return[ Union@@(ParseP[#, node]& /@ List@@expr) ];
   ];

   If[ MemberQ[{Rational, Power, Root}, Head[expr]],
     Throw["Arguments and exponents have to be polynomials over the integers."]
   ];

   Return[ If[ node === "S", {expr}, {} ] ]; (* P := S *)
];

(**
 * Given a a list of rational functions of c-finite sequences of multivariate integer polynomial
 * arguments, and a system of recurrence equations defining c-finite functions,
 * find a pair {myseqs, mydefs} of equivalent data in which all "built-in" special functions
 * are made explicit, and their definitions are added.
 **)
Cleanup[seqs_List, defs_List] :=
   Module[ {myseqs, mydefs, undefined, oldundefined, newdefs, defined, rec, ord, A, 
            i, j, n, p, a, b, x, c, f, r, fun, chebU, chebT, exp, fib, luc, pell, pellLuc, per},

   myseqs = seqs; mydefs = defs;

   mydefs = Flatten[mydefs //. Equal[i__, j_, n_?NumberQ] -> {Equal[i, n], Equal[j, n]}];

   n = ExtractDependentVariable[defs];

   (* 0. make uncoupled implicit equations explicit *)
   Do[

     rec = mydefs[[i]];

     If[ !MatchQ[rec, _[n+Optional[_Integer]] == _] && !MatchQ[rec, _[_Integer]==_],

        f = Union[Cases[rec, f_[n+_.] -> f, Infinity]];

        If[ Length[f] =!= 1, Throw["Illegal recurrence: ", rec ]];

        f = First[f];
        r = Max@@Cases[rec, f[n+i_.] -> i, Infinity];

        mydefs[[i]] = f[n+r] == (f[n+r] /. First[Solve[rec, f[n+r]]]);
     ];
   , {i, 1, Length[mydefs]}];

   (* 1. make chebyshev polynomials explicit and add their definitions *)

   (* T[n+2, x] == 2*x*T[n+1, x] - T[n, x],  T[0] == 1, T[1] == x       *)
   (* U[n+2, x] == 2*x*U[n+1, x] - U[n, x],  U[0] == 1, U[1] == 2*x     *)

   myseqs = myseqs /. {ChebyshevT[p_, x_] :> chebT[x][p],
                       ChebyshevU[p_, x_] :> chebU[x][p]};

   mydefs = mydefs /. {ChebyshevT[n+Optional[i_Integer], x_] -> chebT[x][n+i],
                       ChebyshevU[n+Optional[i_Integer], x_] -> chebU[x][n+i]};

   c = Union[Cases[{myseqs,mydefs}, chebT[x_][_] -> chebT[x], Infinity], 
             Cases[{myseqs,mydefs}, chebU[x_][_] -> chebU[x], Infinity]] // Flatten // Union;

   Do[
   
     mydefs = Union[mydefs,
       {c[[i]][n + 2] == 2*c[[i,1]] * c[[i]][n+1] - c[[i]][n],
        c[[i]][0] == 1, c[[i]][1] == If[ Head[c[[i]]] === chebT, c[[i, 1]], 2*c[[i, 1]] ]}
     ];

   , {i, 1, Length[c]}];

   (* 2. add definitions of Fibonacci, Lucas, Pell, PellLucas, if they appear *)

   {myseqs, mydefs} = {myseqs, mydefs} /. {Fibonacci -> fib, Lucas -> luc, Pell -> pell, PellLucas -> pellLuc, Perrin -> per};

   mydefs = Union[mydefs, 
       {fib[n+2]==fib[n+1]+fib[n], fib[0]==Fibonacci[0], fib[1]==Fibonacci[1],
        luc[n+2]==luc[n+1]+luc[n], luc[0]==Lucas[0], luc[1]==Lucas[1],
        pell[n+2]==2*pell[n+1]+pell[n], pell[0]==Pell[0], pell[1]==Pell[1],
        pellLuc[n+2]==2*pellLuc[n+1]+pellLuc[n], pellLuc[0]==PellLucas[0], pellLuc[1]==PellLucas[1],
        per[n+3]==per[n+1]+per[n], per[0]==Perrin[0], per[1]==Perrin[1], per[2]==Perrin[2]
       }];

   (* 3. make exponentials explicit and add their definitions *)

   Ex = Union @@ (ExtractNodes[#, "Ex"]& /@ myseqs);
   
   {myseqs, mydefs} = {myseqs, mydefs} /. (Ex /. (a_^b_ -> (a^b-> exp[a][b])))
                                       /. (a_^b_ /; !NumberQ[b] -> exp[a][b]);
   mydefs = mydefs /. (exp[a_][b_.*n + c_] -> exp[a][b*n]*a^c);

   c = Union[Flatten[Cases[{myseqs,mydefs}, exp[a_][_] -> exp[a], Infinity]]];

   Do[

     mydefs = Union[mydefs, {c[[i]][n+1] == c[[i,1]]*c[[i]][n], c[[i]][0] == 1}];

   , {i, 1, Length[c]}];

   (* 4. reduce high shifts of built-in functions on rhs *)

   Do[

     If[ MatchQ[mydefs[[i]], _[n+_.]==_], (* otherwise, it's an initial value spec. *)
 
        {f, r} = mydefs[[i]] /. a_[n+Optional[j_Integer]] == _ -> {a, j}; (* order and function symbol *)

        fun = Complement[Union[Cases[mydefs[[i]], fun_[n+_.] -> fun, Infinity]], {f}];

        Do[

          rec = Cases[mydefs, fun[[j]][n+_.]==_];

          If[ Length[rec] =!= 1, 
            Throw["Function " <> ToString[fun[[j]]] <> " multiply defined, or not defined at all."] 
          ];

          If[ Complement[Cases[rec, a_[n+_.] -> a], {fun[[j]]}] =!= {},

            ord = -Infinity;
            A := Throw["System illegally coupled."]; (* will be thrown only if the rule below fires *)
          ,
            b = Min[Cases[rec, fun[[j]][n+a_.] -> a, Infinity]];
            ord = Max[Cases[rec, fun[[j]][n+a_.] -> a, Infinity]] - b;
            rec = First[rec];

            A = Rest[IdentityMatrix[ord]];
            AppendTo[A, Table[Coefficient[Last[rec], fun[[j]][n+a]], {a, b, b + ord - 1}]];
          ];

          mydefs[[i]] = mydefs[[i]] /. fun[[j]][n + k_.] /; k >= r
                       :> First[MatrixPower[A, k - r + ord] . Table[ fun[[j]][n + r - ord + a], {a, 0, ord - 1}]];

        , {j, 1, Length[fun]}]; 
     ];

   , {i, 1, Length[mydefs]}];

   (* 5. select the definitions which are actually needed (greedy). *)

   undefined = Head /@ (Union @@ (ExtractNodes[#, "C"]& /@ myseqs));
   newdefs = {};
   defined = {};

   While[ undefined =!= {},

      oldundefined = undefined;
      undefined = {};

      Do[

         f = oldundefined[[i]];

         (* find an element of mydefs that has f[n+r] as its right hand side. *)
         rec = Cases[mydefs, f[n+Optional[j_Integer]] == _];

         Which[ 
           Length[rec] === 0,
             Throw["Undefined function " <> ToString[f] <> " encountered."],
           Length[rec] > 1,
             Throw["Function " <> ToString[f] <> " has multiple definitions."]
         ];

         rec = First[rec];

         newdefs = Union[newdefs, {rec}, Cases[mydefs, f[_Integer]==_]];
         
         undefined = Union[undefined, Head /@ ExtractNodes[Last[rec], "C"]];

         AppendTo[defined, f];

      , {i, 1, Length[oldundefined]}];

      undefined = Complement[undefined, defined];

   ];
   
   (* 5. return all that stuff. *)

   Return[{myseqs, Union[newdefs]}];

];

(**
 * This function converts a given multivariate polynomial 
 * from the standard basis into the falling factorial basis.
 *
 * The result is given as a list of pairs {c, {e1,e2,...}}
 * such that the sum of all c*ff[v1,e1]*ff[v2,e2]*... is equal
 * to poly.
 **)
FallingFactorialRepresentation[poly_, vars_] :=
   Module[ {p, i, j, d, t, v},

   p = Expand[poly];

   Do[

     p = p /. vars[[i]]^Optional[d_Integer] :> Sum[StirlingS2[d, j]*v[i]^j, {j, 0, d}];

   , {i, 1, Length[vars]}];

   p = Expand[p];

   If[ Head[p] === Plus, p = List@@p, p = {p} ];

   Do[
   
     t = p[[i]]; (* current monomial *)
     p[[i]] = {t /. v[_] -> 1, Table[Exponent[t, v[j]], {j, 1, Length[vars]}]};

   , {i, 1, Length[p]}];

   Return[p];
];

FallingFactorial[n_, d_Integer] := Pochhammer[n-d+1, d];

(******************************************************)
(* Part I: solver for systems of c-finite recurrences *)
(******************************************************)

(**
 * Reimplementation of Mathematica's MatrixPower command.
 * This new version does not perform simplifications and returns the result in 
 * a way that is more suitable for further processing.
 *
 * e[base][n] will be used to denote exponentials in the output.
 *
 * It is guaranteed that the result will be a matrix whose entries
 * are linear combinations of terms e[base][n] with polynomial coefficients.
 **)
MyMatrixPower[M_?MatrixQ, n_, e_] :=
   Module[ {A, J, S, d, i, j, exp},

   {A, J} = JordanDecomposition[M]; (* A . J . Inverse[A] == M *)

   (* JordanDecomposition introduces Algebraic objects that are disturbing later
      when A^(-1) is computed. Remove them *)
   A = Map[RootReduce, A, {2}]; 
   d = Length[A];

   (* idea: M^d == A . J^n . Inverse[A] *)
   (* and J^d can be computed for each Jordan block separately. *)

   Do[ exp[i] = J[[i,i]]; J[[i,i]] = e[exp[i]];  , {i, 1, d}];

   Do[ J[[i,i+1]] = J[[i,i+1]] * n * e[exp[i]]/exp[i] , {i, 1, d - 1} ];

   Do[ 
     Do[
       (* compute the value J[[j, j+i]] *)
       If[ J[[j+1,j+i]] =!= 0 && J[[j,j+i-1]] =!= 0,
         J[[j,j+i]] = FallingFactorial[n, i] * e[exp[j]]/exp[j]^i;
       ];
     , {j, 1, d - i}]
   , {i, 2, d - 1}];

   A . J . Inverse[A]
];

(**
 * GIVEN:  A system of equations of the form
 *
 *   _[n + Optional[_Integer]] == rhs            (recurrence)
 *
 * where the right hand side is a linear combination of terms
 * f[n + i] such that f[n + r] for some r > i appears as the left hand side
 * of some equation. 
 *
 * ...or of the form
 * 
 *   _[_Integer] == rhs                          (initial value)
 * 
 * where rhs is free of n and of all f.
 *
 * Each expression f has to appear precisely once in a left hand side of a 
 * recurrence. The expression n must be the same in all equations. 
 *
 * FIND:  Descriptors for the sequences defined by these recurrences.
 * A sequence descriptor is of the form
 *
 *   {f, {r1, r2, ...}, {{c10, c11, ...}, {c20, c21, ...}, ...}}
 *
 * such that f[n] == (c10+c11*n+...)*r1^n + (c20+c21*n+...)*r2^n + ...
 * one solution for each appearing function symbol. 
 **)
SolveRecurrence[{}] = {};
SolveRecurrence[system:HoldPattern[{(_==_)..}]] :=
   Module[ {n, recs, ansatz, functions, init, ord, a, b, c, d, e, f, i, j, k, l, M, v, v0, r,
            sys, sol, phi, p, mult, done},

   (* Remark: constants in the homogenous part are allowed, polynomials not (yet) *)

   DebugPrint["SR 1"];

   n = ExtractDependentVariable[system];

   recs = Cases[system, _[n+Optional[_Integer]] == _];
   init = Cases[system, _[_Integer] == _];

   DebugPrint["SR 2"];

   (* system contains nothing but recurrences and initial value definitions *)
   If[ Complement[system, recs, init] =!= {},

      Throw[StringJoin["Illegal sequence definition(s) encountered: ",
                       ToString[Complement[system, recs, init]]]];
   ];

   If[ !FreeQ[init, n],

      Throw[StringJoin["Specifications of initial values must be free of ",
                       ToString[n], "."]];
   ];

   DebugPrint["SR 2"];

   functions = Union[Cases[recs, f_[_] == _ -> f]];

   (* throw away initial values for functions for which there are no recurrences given *)
   init = Select[init, MemberQ[functions, Part[#, 1, 0]]&];

   (* shift the whole system such that the smallest shift is 0 *)
   ord = Min@@Cases[recs, _[n+Optional[i_Integer]] -> i, Infinity];
   recs = recs /. f_[n+Optional[i_Integer]] -> f[n+i-ord];

   DebugPrint["SR 3"];

   (* ensure that no right hand side contains a term f[n+i] with i > ord[f] *)

   Do[

      f = functions[[i]];
      If[ Max@@Cases[recs, f[n+j_.]==_ -> j] <= Max@@Cases[recs, _ == _.+_.*f[n+j_.] -> j],
      
         Throw[ToString[f] <> " is used on a right hand side with a too high shift."];
      ];

   , {i, 1, Length[functions]}];

   DebugPrint["SR 4"];

   (* next move all recurrences as far as possible "to the right" *)
   
   done = False;

   While[ !done, (* FIXME: this causes problems for Dependencies[{f[n],Fibonacci[n],Fibonacci[n+1]}, x, Where -> {f[n+2]==f[n]+f[n+1]+Fibonacci[n+2],f[0]==0,f[1]==1}] *)

      done = True;

      Do[

         f = recs[[i, 1, 0]];
         r = recs[[i, 1]] /. f[n+j_.] -> j;

         ord = Max@@Cases[Drop[recs, {i}], f[n+j_.] -> j, Infinity];
         k = Min@@Cases[{recs[[i,2]]}, _[n+j_.] -> j, Infinity];

         If[ ord >= 0, 

            (* f appears on the rhs of some other recurrence *)
            k = Max[r - ord, k]; (* >= 0 *)
         ];

         If[ k > 0 && k =!= Infinity,

            recs[[i]] = recs[[i]] /. f_[n+j_.] -> f[n+j-k];
            done = False;
         ];

      , {i, 1, Length[recs]}];

   ];

   DebugPrint["SR 5"];

   Clear[ord, k];

   (* collect the orders of the recurrences *)

   Do[

      f = functions[[i]];
      ord[f] = Max@@Cases[recs, f[n+j_.]==_ -> j];

   , {i, 1, Length[functions]}];

   (* construct the shift matrix *)

   v = Flatten[Table[functions[[i]][n+j], {i, 1, Length[functions]}, {j, 0, ord[functions[[i]]] - 1}]];
   AppendTo[v, 1];

   DebugPrint["SR 6"];

   M = Table[

     f = Head[v[[i]]];

     If[ v[[i]] === f[n + ord[f] - 1],

        (* convert the recurrence for f into a row of the matrix *)
        rhs = First[Cases[recs, f[n+ord[f]] == j_ -> j]];

        Append[Table[Coefficient[rhs, v[[j]]], {j, 1, Length[v] - 1}], rhs /. _[n+_.] -> 0]
     ,
        (* simple shift *)
        rhs = Table[0, {Length[v]}];
        rhs[[i+1]] = 1;
        
        rhs
     ]
   , {i, 1, Length[v] - 1}
   ];   
   AppendTo[M, Append[Table[0,{Length[v]-1}], 1]];

   If[ !FreeQ[M, n],
      Throw["Sequences are not defined by homogenous linear recurrences."]
   ];

   DebugPrint["SR 7"];

   (* find the involved exponentials *)
   M = MyMatrixPower[M, n, e] /. e[0] -> 0;

   (* now v == (M /. e[x_] -> x^n) . (v /. n -> 0) *)
   v0 = v /. n -> 0 /. (init /. Equal -> Rule);

   (* Compute further initial values from the initial values specified outside the range *)
   sys = Select[init, !MemberQ[(v/.n->0), First[#]]&];

   Do[

      f = sys[[i,1,0]]; k = sys[[i,1,1]];
      j = Position[v, f[n], 1][[1,1]];

      sys[[i]] = C[sys[[i,2]], (M . v0)[[j]]] /. e[a_] -> a^n /. n -> k;

   , {i, 1, Length[sys]}];

   sol = Intersection[Cases[v0, _[_Integer]], Cases[sys, _[_Integer], Infinity]];

   If[ Length[sol] > Length[sys], sol = Take[sol, Length[sys]]];

   sol = Solve[sys /. C -> Equal, sol];

   If[ sol === {},
      Throw["Initial value specification is inconsistent."];
   ];

   DebugPrint["SR 8"];

   v0 = v0 /. First[sol];

   (* create the descriptors for the solution *)

   DebugPrint["SR 9"];

   Table[

      f = functions[[i]];
      c = (M . v0)[[Position[v, f[n], 1][[1,1]] ]]; (* closed form of f *)

      phi = Union[Cases[{c}, e[j_] -> j, Infinity]]; 

      {f, phi, CoefficientList[#, n]& /@ (Coefficient[c, #]& /@ (e /@ phi))}

   , {i, 1, Length[functions]}]

];

(**
 * converts a descriptor as delivered by the solver into a closed form expression.
 * example: {f, {a, b}, {{1,2},{3,4}}} ==> (1+2*n)*a^n + (3+4*n)*b^n
 **)
Descriptor2Expression[{f_Symbol, phi_List, coeffs_List}, n_] :=   
  Module[ {i, j},

  Sum[phi[[i]]^n * Sum[coeffs[[i,j]] * If[j==1, 1, n^(j-1)], {j, 1, Length[coeffs[[i]]]}],
     {i, 1, Length[phi]}]  

];

(*******************************************************)
(* Part II: lattice computation routines               *)
(*******************************************************)

(* FindRelations and ZModuleIntersect are defined in the second part of this file *)

(**
 * Given a basis of a lattice L and an integer d,
 * compute a basis of the lattice L' containing all
 * points v such that v*d belongs to L.
 **)
LatticeDivide[lattice_?MatrixQ, d_Integer] :=
  Module[ {n},

  n = Length[First[lattice]];
  ZModuleIntersect[lattice, d*IdentityMatrix[n]] / d

];
LatticeDivide[{}, _] = {};

(**
 * Converts a lattice to an ideal.
 * The variable corresponding to the i-th component is x[i].
 **)
Lattice2Ideal[{}, _] = {};

(*
Lattice2Ideal[b_, x_] := Module[ {base},
   base = Transpose[b];
   Numerator /@ Together /@ (
    Product[x[j]^base[[j]], {j, 1, Length[base]}] - 1
   )
];
*)

Lattice2Ideal[b_, x_] := Module[ {e, y},

  base = Table[ 
    Product[
      e = b[[i,j]]; If[e > 0, x, y][j]^Abs[e]
      , {j, 1, Length[b[[i]]]}] - 1,
    {i, 1, Length[b]}
  ];
  base = Join[base, Table[x[j]*y[j] - 1, {j, 1, Length[First[b]]}]];
  base = EliminationIdeal[base, Table[y[i], {i, 1, Length[First[b]]}], Table[x[i], {i, 1, Length[First[b]]}]];

  Return[base];
];

(*******************************************************)
(* Part III: User interface                            *)
(*******************************************************)

(**
 * Given multivariate sequences f and f1,f2,... and a variable x (opt),
 * compute a representation of f in terms of f1, f2, ... if possible.
 **)
Options[Express] := {Where -> {}, Variables -> Automatic};
Express[f_, terms_List, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {x, i},

  Express[f, terms, x, opts] /. x[i_] :> terms[[i]]
];
Express[f_, terms_List, x_, opts:((_Rule|_RuleDelayed)...)] :=
  Module[ {id, i},

  id = Dependencies[Prepend[terms, f], x, opts, Cleanup -> False] /. x[i_] -> x[i-1];
  id = EliminationIdeal[id, {x[0]}, Table[x[i], {i, 1, Length[terms]}], Keep -> True, Cleanup -> False];
  id = Select[id, !FreeQ[#, x[0]]&];
  id = Flatten[(x[0] /. Solve[#==0, x[0]])& /@ id];

  Return[id];  
];

(**
 * Given multivariate sequences f1,f2,..., and a variable x,
 * compute an ideal basis of kernel(Phi), where 
 * Phi maps x[1] to f1, x[2] to f2, etc.
 **)
Options[Dependencies] = {Where -> {}, Variables -> Automatic};
Dependencies[seqs_List, opts:((_Rule|_RuleDelayed)...)] :=
   Module[ {x, i},

   Dependencies[seqs, x, opts] /. x[i_] :> seqs[[i]]
];
Dependencies[seqs_List, x_, opts:((_Rule|_RuleDelayed)...)] := 
   Module[ {myseqs, C, V, id, p, exp, a, b, c, d, e, f, i, j, n, L, defs, used, vars, basis, u},

   (* I. Preprocessing *)
   DebugPrint["I."];

   (* 1. bring all the data into a standardized form *)
   DebugPrint["I.1"];

   defs = Where /. {opts} /. Options[Dependencies];

   {myseqs, defs} = Cleanup[seqs, defs];

   defs = SolveRecurrence[defs];

   (* 2. extract the exponentials that actually occur *)
   DebugPrint["I.2"];

   phi = Union@@Table[defs[[i, 2]], {i, 1, Length[defs]}];

   (* 4. replace the corresponding entries in defs by pointers to phi *)
   DebugPrint["I.4"];
   
   Do[
      defs[[i, 2, j]] = Position[phi, defs[[i, 2, j]]][[1,1]];
   , {i, 1, Length[defs]}, {j, 1, Length[defs[[i, 2]]]}];

   (* 5. compute the first lattice now and the others on demand *)
   DebugPrint["I.5"];

   L[0] = {};
   L[1] = FindRelations[phi];
   L[d_Integer] := (L[d] = LatticeDivide[L[1], d!]);   

   (* II. Build up the ideal basis *)
   DebugPrint["II."];

   basis = {};

   (* we use c[i][x,y,z,...] as variable for 
      phi[[i]] ^ (ff[V[[1]], x] * ff[V[[2]], y] * ff[V[[3]], z] * ...)

      and c[i] as variable for C[[i]]
      and the elements of V for themselves
      and c[-i] for the denominator of myseqs[[i]] (if nontrivial)
    *)

   (* 1. collect the closed forms of the involved c-finite functions *)
   DebugPrint["II.1"];

   C = Union@@(ExtractNodes[#, "C"]& /@ myseqs);

   V = Variables /. {opts} /. Options[Dependencies];
   If[ Head[V] === List,
     V = Intersection[V, Union@@(ExtractNodes[#, "S"]& /@ myseqs)],
     V = Union@@(ExtractNodes[#, "S"]& /@ myseqs)
   ];

   Do[ 

      p = C[[i, 1]]; (* argument polynomial *)

      ffp = FallingFactorialRepresentation[p, V]; 

      (* get the position of the closed form *)
      d = First[Cases[defs, {C[[i, 0]], _, _}]];
      AppendTo[basis, -c[i] + def2poly[c, d, p, ffp]];

   , {i, 1, Length[C]}];

   basis = basis /. c[i_][0..] :> phi[[i]]; (* a^(n^0)= a^1 = a *)   

   (* 2. add the algebraic dependencies of the exponentials *)
   DebugPrint["II.2"];

   (* which exponent vectors do actually occur in the argument polynomials? *)
   e = Union[Flatten[Cases[basis, c[_][i___] -> a[i], Infinity]]] /. a -> List;

   Do[

     b = Union@@Table[L[e[[i,j]]], {j, 1, Length[e[[i]]]}]; 
     basis = Join[basis, Lattice2Ideal[b, u] /. u[j_] -> c[j]@@e[[i]]];

   , {i, 1, Length[e]}]; 

   (* 3. add the generators of the given expression *)
   DebugPrint["II.3"];

   basis = Union[basis, 
     Table[x[i] - (myseqs[[i]] /. Table[C[[j]] -> c[j], {j, 1, Length[C]}]),
        {i, 1, Length[myseqs]}]
   ];

   basis = Numerator /@ Together /@ basis; (* clear denominators *)

   (* 4. add Rabinovitz variables for nontrivial denominators in myseqs *)
   DebugPrint["II.4"];

   Do[

     d = Denominator[Together[myseqs[[i]]]];

     If[ !NumberQ[d], AppendTo[basis, c[-i]*d - 1 /. Table[C[[j]] -> c[j], {j, 1, Length[C]}]] ];

   , {i, 1, Length[myseqs]}];   

   (* 5. Eliminate all dummy variables *)
   DebugPrint["II.5"];

   basis = EliminationIdeal[basis, 
         Union[Flatten[Cases[basis, c[_][___], Infinity]],
               Flatten[Cases[basis, c[_], Infinity]],
               V],
         Table[x[i], {i, 1, Length[myseqs]}],
         opts
      ];

   (* ...done... *)

   Return[RootReduce /@ basis];
];

polypower[c_, fpp_] := Module[ {i}, 

  Product[ c@@fpp[[i, 2]]^fpp[[i, 1]], {i, 1, Length[fpp]} ]
];

def2poly[c_, d_, p_, fpp_] := Module[ {i,j,k},

  Sum[
    Sum[ d[[3, i, j]]*p^(j-1), {j, 1, Length[d[[3, i]]]}]
      * polypower[c[d[[2, i]]], fpp], 
  {i, 1, Length[d[[2]]]}]
];

(**
 * Compute an Elimination ideal.
 *
 * {p1,p2,...}... a list of multivariate polynomials in x1,x2,...,y1,y2,...
 * {x1,x2,...}... variables to be eliminated
 * {y1,y2,...}... variables to be kept.
 **)
Options[EliminationIdeal] = {`Keep -> False, `Cleanup -> True};
EliminationIdeal[basis_List, {}, stay_List, opts___] := basis;
EliminationIdeal[basis_List, away_List, stay_List, opts___] :=
   Module[ {M, G, mystay, alg, t, i, a, minpoly, keep, cleanup},

   {keep, cleanup} = {Keep, Cleanup} /. {opts} /. Options[EliminationIdeal];

   (* use a fast elimination ordering instead of Lex *)

   (* matrix orders do not work over function fields *)
   mystay = Join[stay, Complement[Variables[basis], away, stay]];

   (* we make algebraic numbers explicit in hope for better performance *)
   alg = Union[Cases[basis, (Algebraic[_,_,_]|Root[_,__]|Complex[_,_]|Power[_,_Rational]|Sqrt[_]), Infinity]];

   alg = Apply[Rule, Transpose[{alg, ToCommonField[alg]}], {1}];
   G = basis /. alg /. Algebraic[_, a_List, _] :> Sum[a[[i]]*t^(i-1),{i,1,Length[a]}];

   If[ Length[alg] > 0, 
     AppendTo[G, alg[[1,2,1]][t]];
     AppendTo[mystay, t];
   ];

   Which[
     NameQ["Singular`SingularVersion"],

      (* use Singular if it is available *)
      If[ keep, 
        G = Singular`SingularStd[G, away, mystay, MonomialOrder -> "dp"],
        G = Singular`SingularEliminate[G, away, mystay]
      ];

    , NameQ["Gb`GbVersion"],

      (* use Faugere if it is available *)
      G = Gb`GBasis[G, away, mystay, 
             MonomialOrder -> DegreeReverseLexicographic,
             Eliminate -> !keep];
    , True,
      (* otherwise, try Mathematica's builtin *)
      If[ keep,
        G = GroebnerBasis[G, Join[away, mystay]],
        G = GroebnerBasis[G, mystay, away]
      ];
      (*
      M = -IdentityMatrix[Length[away] + Length[mystay]];

      M[[1]] = Join[Table[1, {Length[away]}], Table[0, {Length[mystay]}]];
      M[[2]] = 1 - M[[1]];

      If[ Length[away] > 1, M[[2,2]] = -1 ];

      G = GroebnerBasis[G, Join[away, mystay], MonomialOrder -> M];
      G = Select[G, (Intersection[Variables[#], away] == {}) &];
      *)
   ];

   (* undo algebraic numbers *)
   If[ Length[alg] > 0, 
      G = G /. t -> Algebraic[alg[[1,2,1]],{0,1},alg[[1,2,3]]] // RootReduce
   ]; 

   If[ cleanup,
     G = GroebnerBasis[G, stay, MonomialOrder -> DegreeReverseLexicographic]; (* cleanup *)
   ];

   Return[G];
];

Protect[Dependencies];

End[];

(* ###################################################################################### *)
(* ###################################################################################### *)
(* ## Code below provides implementation of Ge's algorithm ############################## *)
(* ###################################################################################### *)
(* ###################################################################################### *)

(* This code should be kept independent from the rest so that it can be pasted into an independent 
   package, if desired *)

(* needs: NumberTheory`AlgebraicNumberFields`, NumberTheory`ContinuedFractions` *)

Unprotect[FindRelations];

FindRelations::usage = 
  StringJoin["FindRelations[{a[1],a[2],...}] computes a basis of the lattice ",
             "of all integer vectors {v[1],v[2]...} such that a[1]^v[1]*a[2]^v[2]*...=1. ",
             "The numbers a[i] have to be algebraic over the rationals, given as Root or ",
             "Algebraic objects. If they are not given as Algebraic objects of the same ",
             "field, then the splitting field has to be computed and this computation may ",
             "dominate all other steps in runtime."];

Begin["`Private`"];

(**
 * ZModuleIntersect[base1, base2];
 *
 * Given two Z-modules, compute a basis for their intersection. 
 **)
ZModuleIntersect[base1_List, base2_List] := 
  Module[ {d, b, a, sol},

  (* like in linear algebra, but with ZNullSpace instead of NullSpace *)
  If[ base1 == {} || base2 == {}, Return[ {} ]];

  sol = ZNullSpace[Join[base1, base2]];

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

  LatticeReduce[Transpose[Transpose[base1] . Transpose[Take[sol, Length[sol], Length[base1]]]]]
];

(**
 * Computes the integer kernel of an integer matrix.
 **)
ZNullSpace[mat_] := 
  Module[ {U, A},    

  {U, A} = Developer`HermiteNormalForm[mat]; 

  (* kernel is generated by the rows of U that correspond to zero rows in A *)
  (* cf. H. Cohen, prop. 2.4.9. *)

  Extract[U, Position[A, {0..}]] 
];

(**
 * Given a list of algebraic numbers, a[1], ..., a[k], this method returns a number B
 * such that the lattice of all vectors {v[1],...,v[k]} with 
 *
 *  a[1]^v[1]*...*a[k]^v[k] == 1
 *
 * has a basis whose basis vectors have norm at most B.
 **)
MasserBound[roots_] := 
  Module[ {k, D, h, h0, p, i, eta, x},

  (* k = number of roots *)
  k = Length[roots]; (* number of roots *)

  (* D = degree of the field extension *)
  p = Union[Cases[roots, Algebraic[i_, __] -> i]];

  Which[ 
    Length[p] == 0 && FreeQ[roots, Sqrt] && FreeQ[roots, Root],
      D = 1,
    Length[p] == 1, 
      D = Exponent[First[p][x], x],
    True,           
      Throw["Input should be specified as Algebraic objects of one single algebraic extension of Q"]
  ];

  (* h = maximum height of the a[i]. The height of an algebraic number is the sum of the degree and 
     the binary length of all coefficients in the defining equation over Q *)
  h = 0;
  For[i = 1, i <= Length[roots], i++,

     p = MinimalPolynomial[ roots[[i]], x];
     h0 = Ceiling[Exponent[p, x] + Plus@@Log[2, Abs[DeleteCases[CoefficientList[p, x], 0, 1]]]];

     If[ h0 > h, h = h0 ];
  ];

  (* Masser's bound *)
  Ceiling[ D^2*(4*h*k*D*(Log[2+D]/Log[Log[2+D]])^3)^(k-1) + 1 ]

];

(**
 * FindRelations[roots]
 *
 * roots... a list of algebraic numbers, given as Algebraic or as Root objects
 * 
 * returns a list of integer vectors {b[1],...,b[d]} generating the lattice of 
 * all vectors v with with 
 *
 *   Product[roots[i]^v[i], {i}] == 1 .
 * 
 * If the roots are given as members of the same algebraic extension of Q, then 
 * the procedure has polynomial runtime. Otherwise, a common algebraic extension
 * field of the given roots is computed, and this may take long.
 **)
FindRelations[roots_List] := 
   Module[ {i, z, a, bound, reLog, imLog, m, m1, m2}, 

   (* first treat zeros in the root list. Assume 0^0 = 1 *)

   If[ MatchQ[ roots, {0...} ], Return[ {} ] ];

   If[ !FreeQ[ roots, 0, 1], 

       ex = Position[roots, 0, 1]; (* positions of the zeroes *)
       B = FindRelations[ DeleteCases[roots, 0] ]; (* basis for the nonzero numbers *)
       Do[ B = Insert[#, 0, ex[[i, 1]] ]& /@ B, {i, 1, Length[ex]}]; (* insert new dimensions *)

       Return[ B ];
   ];

   DebugPrint["FR 1"];

   a = ToCommonField[roots]; (* does nothing if the roots all belong to the same field *)

   DebugPrint["FR 2"];

   reLog = Re[Log[a]];
   imLog = Append[Im[Log[a]], 2*Pi];

   DebugPrint["FR 3"];

   (* replace implizit zeros by explicit ones. *)
   For[i = 1, i <= Length[a], i++,

      z = a[[i]];

      (* Abs[z] == 1 *)
      If[ RootOfUnityQ[z], reLog[[i]] = 0 ];

      (* z is real and z >= 0 *)
      If[ Element[z, Reals] && Element[Sqrt[z], Reals], 
        imLog[[i]] = 0, 
        Null,
         (* try harder if no decision was found *)
         If[ Element[RootReduce[z], Reals] && Element[Sqrt[RootReduce[z]], Reals], imLog[[i]] = 0 ];
      ];
   ];

   DebugPrint["FR 4"];

   (* compute a bound for the coefficients of the generators *)
   bound = MasserBound[a]; 

   DebugPrint["FR 5"];

   m = IdentityMatrix[Length[a]];

   (* successively refine the approximation until only valid generators are returned *)
   For[ level = Ceiling[N[Log[2, bound]] + 1], !(And@@(CheckRelation[a,#]&/@m)), level++, 

      DebugPrint["FR 6"];

      m1 = GetBasis[reLog, level, bound];
      m2 = Drop[#, -1]& /@ GetBasis[imLog, level, bound];
      m = ZModuleIntersect[m1, m2];
   ];

   DebugPrint["FR 7"];

   m
];

(** 
 * Computes the integer relations with coefficients bounded by bound of the real numbers in l based on
 * continued fraction approximations of at least the specified level. Arbitrary expressions are allowed
 * as elements of l, as long as N is able to map them to a float. However, it is assumed that an 
 * expression that is different from 0 corresponds to a number that is different from zero. 
 *
 * The returned basis is guaranteed to contain all integer relations for l with coefficients bounded by
 * bound, but it may contain additional, wrong, "relations", if the approximation level is chosen too 
 * small. For sufficiently large level, only true relations will be returned.
 **)
GetBasis[l_List, level_Integer, bound_Integer] :=
   Module[ {B, i, j, c1, c2, d, e, ex, n, lev},

    n = Length[l];

    (* first treat zero elements in l as special cases *)

    If[ MatchQ[l, {0..}], (* l has zeroes only *)

       Return[IdentityMatrix[n]] 
    ]; 

    If[ !FreeQ[l, 0, 1], (* l contains at least one zero *)

       ex = Position[l, 0, 1]; (* positions of the zeroes *)
       B = GetBasis[ DeleteCases[l, 0], level, bound]; (* basis for the nonzero numbers *)
       Do[ B = Insert[#, 0, ex[[i, 1]] ]& /@ B, {i, 1, Length[ex]}]; (* insert new dimensions *)
       B = Join[B, IdentityMatrix[n][[Flatten[ex]]] ]; (* add unit vectors at the zero positions *)

       Return[ B ];
    ];

    (* now we assume that l is a list of nonzero real numbers *)

    c1 = Last[Convergents[#, level]]& /@ l;
    c2 = Last[Convergents[#, level + 1]]& /@ l;

    e = Table[ 1/Denominator[c1[[i]]*c2[[i]]] , {i, 1, n}];
    (* cfrac theorem says: | log[i] - c1[i] | <= e[i] *)

    (* refine the approximation such that all errors are smaller than the smallest
       element of l in absolute value *)

    lev = level + 1;
    While[ (* exists i, j st e[i] >= c1[j], replace c1[j] by better approximation. *)
      Select[c1, (Max@@e >= Abs[#])&] =!= {} && lev < level + 5,

      ex = Flatten[Position[e, Max@@e]]; (* indices with greatest errors *)

      lev++;
      For[ i=1, i <= Length[ex], i++,
         j = ex[[i]];
         c1[[j]] = c2[[j]];
         c2[[j]] = Last[Convergents[ l[[j]], lev]];
         e[[j]] = If[ c1[[j]] == l[[j]], 0, 1/Denominator[c1[[j]] * c2[[j]]] ];
      ];
    ];

    (* now: max e[i] < min |c1[i]| *)

    (* this bound guarantees that generators whose norm is at most bound will
       appear in the LLL-reduced basis *)
    d = Ceiling[N[2^((Length[c1]-1)/2)*bound/(Min@@(Abs/@c1) - Max@@(Abs/@e))]];

    B = Transpose[Append[IdentityMatrix[n], c1*d]];
    B = LatticeReduce[B];

    (* Vectors whose right hand side is greater than the errors permit can be 
       discarded; they cannot correspond to integer relations. *)
    B = Select[B, (Abs[Last[#]] <= d*Sum[Abs[#[[j]]*e[[j]]], {j, 1, n}])&];

    (* the following old version was sufficient for some time but once became
       incomplete for some strange reason:
    B = Select[B, (Abs[Last[#]] <= d*Abs[Sum[#[[j]]*e[[j]], {j, 1, n}]])&];
    *)

    (* all remaining vectors are returned as candidates *)
    Drop[#, -1]& /@ B
];

(**
 * given a list of algebraic numbers {a[1],a[2],...}, all belonging to the same algebraic 
 * extension of Q, and an integer vector {e[1],e[2],...} this function returns
 * True if a[1]^e[1] * a[2]^e[2] * ... = 1, and False otherwise.
 **)
CheckRelation[a_List, e_List] := (Times@@(a^e) == 1);

Protect[FindRelations];

End[];

(* RISC header *)
If[$Notebooks,
   CellPrint[Cell[#, "Print", FontColor -> RGBColor[0, 0, 0], 
                              CellFrame -> 0.5, 
                              Background -> RGBColor[0.796887, 0.789075, 0.871107]]]&,
   Print
  ]["Dependencies.m by Manuel Kauers and Burkhard Zimmermann " <> 
    "\[LongDash] \[Copyright] RISC Linz \[LongDash] V " <> Dependencies`Private`Version];

EndPackage[];
