CFSum: A Program for Summing C-Finite Sequences

Herbert Wilf
Curtis Greene
December 2004

Introduction

This is a program for evaluating sums S(n) = (∑_i) (Subscript[F, 1])((Subscript[a, 1])n + (Subscript[b, 1])j + (Subscript[c, 1]))⋯(Subscript[F, k])((Subscript[a, k])n + (Subscript[b, k])j + (Subscript[c, k])) in closed form, where the functions (Subscript[F, j]) are C-finite, i.e., each (Subscript[F, j]) satisfies a linear recurrence of degree D(j) < ∞ with constant coefficients.The output is a linear combination of the terms of the form (Subscript[F, 1])(((Subscript[a, 1])+ (Subscript[b, 1]) )n +(Subscript[i, 1]))⋯(Subscript[F, k])(((Subscript[a, k]) + (Subscript[b, k]))n+ i\(TraditionalForm\`\_k\)) and   ψ(n) (Subscript[F, 1])((Subscript[a, 1])n+(Subscript[i, 1]))⋯(Subscript[F, k])((Subscript[a, k])n +(Subscript[i, k])) where ψ(n) is a polynomial, and, for each j, the indices (Subscript[i, j]) range between 0 and D(j)-1.

The method follows that described in our paper,  "Closed Form Summation of C-Finite Sequences", to appear in Trans. American Math. Society.

How the program works

To compute a sum of the form S(n) = (∑_i) (Subscript[F, 1])((Subscript[a, 1])n + (Subscript[b, 1])j + (Subscript[c, 1]))⋯(Subscript[F, k])((Subscript[a, k])n + (Subscript[b, k])j + (Subscript[c, k])) , the only required input is a list of the form

{{a[1], b[1], c[1]}, {a[2], b[2], c[2]},  ... , {a[k], b[k], 
        c[k]}}

and explicit definitions for functions

F[n, 1],  F[n, 2],  ... ,  F[n, k]

which are assumed to be C-finite. The program computes the actual C-finite degree, the multiplicities of the roots, and also a minimal recurrence satisfied by each F[n,i]. The exact format is

CFSum[<abclist>, <newnames>]

It then computes a list of N distinct "target monomials", and solves a system of N linear equations involving S(1), S(2), ..., S(N) to express S(n) as a linear combination of those monomials. By results in our paper, any such linear combination agreeing with S(n) for n = 1,2, ..., N must agree for all values of N.

If the resulting linear combination is not uniquely determined, the expression will contain arbitrary constants denoted xx[1], xx[2], xx[3], ..., etc.

The parameter <newnames> is optional. If omitted, then the names of the functions F[n,i] as they appear on output will be "f[n,i]". Otherwise, <newnames> can be specified as a list of replacement rules, e.g., {F[n_,1]->n,F[n_,2]->Fibonacci[n]}, and these name substitutions will be made on output. Note that if you replace with the name of a user-defined function, as opposed to a built-in function, the output might look unpleasant, with towers of "If" statements, so this should be avoided.

Important: If you define a function directly with a recurrence relation,  be sure to use the Mathematica "f[n_]:=f[n]=" construction so the program will remember previously computed values. This can speed the calculation immensely.

Before running the program, be sure to execute the initalization cells using the Kernel > Evaluation menu.

Example 3.1 (Fibonacci Example)

F[n_, 1] := Fibonacci[n];
F[n_, 2] := Fibonacci[n];
F[n_, 3] := Fibonacci[n];

CFSum[{{0, 1, 0}, {0, 1, 0}, {2, -1, 0}}, {f[n_, _] -> 
          Fib[n]}]

Example 3.2 (Subword Avoidance, A=2)

A = 2;

Clear[G];
G[n_]  :=  G[n]  =  
        If[n == 0, 1, 
            If[n == 1, A, 
              If[n == 2, A^2, A G[n - 1]  -  G[n - 3]]]] // 
          Simplify

F[n_, 1] := G[n];
F[n_, 2] := G[n];

CFSum[{{0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> GG[n]}]

Collect[%, xx[11]]

Example 3.2 (Subword Avoidance, General A)

A =.

Clear[G];
G[n_]  :=  G[n]  =  
        If[n == 0, 1, 
          If[n == 1, A, 
            If[n == 2, A^2, 
              A G[n - 1]  -  G[n - 3] // Expand]]]

F[n_, 1] := G[n];
F[n_, 2] := G[n];

Table[G[n], {n, 0, 10}]

CFSum[{{0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> GG[n]}]

% // Simplify

Example 3.3 (Fibonacci Power Sums)

F[n_, 1] := Fibonacci[n];
F[n_, 2] := Fibonacci[n];
F[n_, 3] := Fibonacci[n];
F[n_, 4] := Fibonacci[n];
F[n_, 5] := Fibonacci[n];
F[n_, 6] := Fibonacci[n];
F[n_, 7] := Fibonacci[n];

CFSum[{{0, 1, 0}}, {f[n_, _] -> f[n]}]

CFSum[{{0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> f[n]}]

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> 
          f[n]}]

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> 
          f[n]}]

Collect[%, xx[18]]

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 
          0}}, {f[n_, _] -> f[n]}]

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 
          0}}, {f[n_, _] -> f[n]}]

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 
          0}, {0, 1, 0}}, {f[n_, _] -> f[n]}]

Example 3.4 (Generic Power Sums)

A =. ; B =.

Clear[H];
H[n_]  :=  H[n]  =  
        If[n == 0, p, 
          If[n == 1, q, 
            A H[n - 1]  +  B H[n - 2] // Expand]]

Table[H[n], {n, 0, 5}]

F[n_, 1] := H[n];
F[n_, 2] := H[n];

CFSum[{{0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> h[n]}]

% // Simplify

Example 9.1 (Mixed Convolution)

Clear[G];
G[n_]  :=  G[n]  =  
        If[n == 0, 1, 
            If[n == 1, 3, 
              If[n == 2, 9, 3 G[n - 1]  -  G[n - 3]]]] // 
          Simplify

F[n_, 1]  :=  n;
F[n_, 2] := Fibonacci[n];
F[n_, 3] := G[n];

CFSum[{{0, 1, 0}, {0, 1, 0}, {1, -1, 0}}, {f[n_, 1] -> n, 
        f[n_, 2] -> f[n], f[n_, 3] -> g[n]}]

Example 9.2 (Independence Property)

Clear[F];
F[n_]  := F[n]  =  
        If[n == 0, p, If[n == 1, 1, F[n - 1]  + F[n - 2]]] // 
          Simplify

Table[F[n], {n, 0, 6}]

F[n_, 1]  := F[n];
F[n_, 2] := F[n];
F[n_, 3] := F[n];

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> 
          f[n]}]

F[n_, 1]  := F[n];
F[n_, 2] := F[n];
F[n_, 3] := F[n];
F[n_, 4] := F[n];

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}}, {f[n_, _] -> 
          f[n]}]

To get it to agree with the expression in our paper:

% /.  xx[18] -> (19/25 - K);
Collect[%, K]

Example 9.3 (Partial summation of series)

F[n_, 1]  := x^n;
F[n_, 2] := Fibonacci[n];

CFSum[{{0, 1, 0}, {0, 1, 0}}, {f[n_, 1] -> x^n, 
        f[n_, 2] -> f[n]}]

F[n_, 1]  := x^n;
F[n_, 2] := Fibonacci[n];
F[n_, 3] := Fibonacci[n];

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}}, {f[n_, 1] -> x^n, 
        f[n_, 2] -> f[n], f[n_, 3] -> f[n]}]

% // Simplify

Example 9.4 (Degree of the polynomial multiplier)

F[n_, 1]  := n^5;

CFSum[{{0, 1, 0}}, {f[n_, 1] -> n^5}]

F[n_, 1]  := n;
F[n_, 2]  := n;
F[n_, 3]  := n;
F[n_, 4]  := n;
F[n_, 5]  := n;

CFSum[{{0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 
          0}}, {f[n_, _] -> n}]

Example 9.5 (Example from the theory of partitions)

Clear[P5]
P5[n_]  :=  P5[n]  = 
        SeriesCoefficient[
          Series[Product[1/(1 - z^i), {i, 1, 5}], {z, 0, 200}], 
          n]

Table[P5[n], {n, 0, 30}]

F[n_, 1]  := P5[n];

CFSum[{{0, 1, 0}}, {f[n_, _] -> p[n]}]

Collect[%, xx[17]]

-1001 - ((3875 n)/12) - ((935 (n^2))/24) - ((25 (n^3))/12) - 
        ((n^4)/24) // Factor

Technical Details (This part can safely be ignored)

\(CFSum[abc_, newnames_: {}] := 
      Module[{k, ff, dst, a, b, c, i, pl, oo, \ outr, eqs, varbs, mxm, pdg, 
          tt, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ kk, ind, pp, 
          psi, sm, zz, ty, ss, nn, h, deg, 
          x}, 

 (*\ Evaluates\ a\ sum, \ 
          from\ j = 0\ to\ n - 1, \ \(of

              \ \ \ F[a[1] n + b[1] j + c[1], 1] 
                  F[a[2] n + b[2] j + c[2], 2]\  ... \)\ F[
                a[k] n + b[k] j + c[k], k]
\ 
              in\ which\ each\ F[m, j]\ is\ C - 
            finite . 
The\ input\ abc\ is\ the\ array\ of\ length\ k\ of\ the\ triples\ {a[j], b[j], 
                  c[j]} . \ 
The\ k\ functions\ F[m, 
                j]\ are\ specified\ externally . \ 
Input\ newnames\ is\ optional . \ If\ omitted, \ 
          then\ the\ names\ of\ the\ functions\ F[n, 
              i]
\ \ \ \ as\ they\ appear\ on\ output\ will\ be\ "f[n,i]" . \ Otherwise, \ 
          newnames\ can\ be
\ \ \ \ specified\ as\ a\ list\ of\ replacement\ rules, \ 
          e . g . , \ {F[n_, 1] -> n, 
                  F[n_, 2] -> 
                    Fibonacci[
                      n]} . 
\ \ \ \ and\ these\ name\ substitutions\ will\ be\ made\ on\ output . \ 
                Note\ that\ if\ you\ replace
\ \ \ \ with\ the\ name\ of\ a\ user - defined\ function, \ 
          as\ opposed\ to\ a\ built - in\ function, \ 
          the
\ \ \ \ \ output\ might\ look\ unpleasant, \ 
          with\ towers\ of\ "If"\ statements, \ 
          so\ this\ should\ be
\ \ \ \ \ avoided . 
This\ routine\ assumes\ that\ the\ routines\ "indices", \ "maxmult", \ "companion", \ "FindDegree", 
\ \ \ \ \ \ and\ "kronecker"\ below\ are\ activated . 
\(Important : \ 
                If\ you\ define\ a\ function\ directly\ with\ a\ recurrence\ relation\), \ 
\ \ \ be\ sure\ to\ use\ the\ Mathematica\ "f[n_]:=f[n]="\ construction\ so\ the\ program
\ \ will\ remember\ previously\ computed\ values . \ 
              This\ can\ speed\ the\ calculation
\ \ \(\(immensely\)\(.\)\)\ 
\ \ *) 

k = Length[abc]; Clear[ff, f, n]; \ ff = 0; 
        Off[Solve::svars]; 

 (*\ 
          Find\ CF - 
            degrees\ of\ input\ sequences\ *) 
deg = 
          Table[0, {k}]; 
For[i = 1, 
          i ≤ k, \(i++\), h[n_] := F[n, i]; 
          deg[\([i]\)] = 
            FindDegree[h]]; 
Print["Degrees are ", 
          deg, "."]; 
dst = Table[0, {k}]; 
        mxm = Table[0, {k}]; 
        outr := {{1}}; 
a = 
          Table[\(abc[\([i]\)]\)[\([1]\)], {i, 1, k}]; 
b = 
          Table[\(abc[\([i]\)]\)[\([2]\)], {i, 1, k}]; 
c = 
          Table[\(abc[\([i]\)]\)[\([3]\)], {i, 1, 
              k}]; 

 (*\ 
          Finds\ polynomials\ associated\ with\ each\ of\ the\ F' 
            s\ *) 
For[i = 1, 
          i ≤ k, \(i++\), 
\ 
          pl = x^\((deg[\([i]\)])\) - 
              Sum[cc[j]\ x^j, {j, 0, deg[\([i]\)] - 1}]; 
\ 
          eqs = Table[
              F[n + deg[\([i]\)], i] == 
                Sum[cc[j]\ F[n + j, i], {j, 0, deg[\([i]\)] - 1}], {n, 0, 
                deg[\([i]\)] - 1}]; 
\ 
          varbs = Table[cc[k], {k, 0, deg[\([i]\)]}]; 
\ 
          pl = \((\(Expand[pl /. Solve[eqs, varbs]]\)[\([1]\)])\) /. 
              Table[cc[j] -> 0, {j, 0, 
                  deg[\([i]\)] - 1}]; 
\  (*\ 
            Set\ up\ array\ of\ numbers\ of\ distinct\ roots\ of\ each\ *) 
\ 
          dst = ReplacePart[dst, 
              Exponent[
                PolynomialQuotient[pl, PolynomialGCD[pl, D[pl, x]], x], x], 
              i]; 
\  (*\ 
            Set\ up\ array\ of\ maximum\ multiplicities\ of\ each\ polynomial\ *) 
\ 
          mxm = ReplacePart[mxm, maxmult[pl, x], 
              i]; 
\ \  (*\ 
            
            Accumulate\ companion\ matrix\ of\ pl\ into\ NoMonom\ matrix\ outer\ product\ \ *) 
\ 
          outr = kronecker[outr, 
              MatrixPower[companion[pl, x], 
                b[\([i]\)]]];
]; 

 (*Test\ whether\ any\ monomials\  = \ 
            1\ *) 
oo = \(Dimensions[outr]\)[\([1]\)]; 
        NoMonom = \((Det[IdentityMatrix[oo] - outr] ≠ 
              0)\); 
If[NoMonom, 
          Print["No monomial product equals 1."], 
          Print["There are monomial products equal to 1."], 
          Print["If ", 
            NoMonom, " then no monomial product equals 1."]]; 

 (*Compute\ max\ degree\ of\ polynomial\ psi*) 
pdg = 
          1 + Sum[If[a[\([i]\)] == 0, mxm[\([i]\)] - 1, 0], {i, 1, 
                k}]; 

 (*\ 
          Assemble\ list\ of\ target\ monomials\ *) 
tt = 
          indices[deg]; 
        ls = {}; 
For[kk = 1, 
          kk ≤ Length[tt], \(kk++\), 
ind = 
            tt[\([kk]\)]; 
ls = 
            Union[ls, \ {Product[
                  If[a[\([i]\)] + b[\([i]\)] == 0, 1, 
                    f[\((a[\([i]\)] + b[\([i]\)])\)\ n + ind[\([i]\)], 
                      i]], {i, 1, k}]}]; 
pp = 
            Product[If[a[\([i]\)] ≠ 0, 
                f[\((a[\([i]\)]\ n + ind[\([i]\)])\), i], 1], {i, 1, 
                k}]; 
ls = 
            Union[ls, {pp}]; 
If[Not[NoMonom], 
            For[h = 1, h ≤ pdg, \(h++\), 
              ls = Union[ls, {n^h\ pp}]]]]; 
jj = 
          Length[ls]; \ 
        Print["There are ", 
          jj, " target monomials."]; 
ff = 
          Sum[xx[j]\ ls[\([j]\)], {j, 1, jj}]; 
For[i = 1, 
          i ≤ k, \(i++\), 
          f[m_, i_] := F[m, i]]; 
sm[nn_] := 
          Sum[Product[
              F[a[\([i]\)]\ nn + b[\([i]\)]\ j + c[\([i]\)], i], {i, 1, 
                k}], {j, 0, nn - 1}]; 
 (*\ 
          cfd = FindDegree[
              sm]; 
Print["The C-finite degree of the sum is ", cfd];\ *) 
zz = 
          Table[xx[j], {j, 1, jj}]; 
ty = 
          Table[sm[nn] == ff /. {n -> nn}, {nn, 0, 
              jj - 1}]; 
ss = \(Solve[ty, 
              zz]\)[\([1]\)]; 
zz = \((zz /. {ss})\)[\([1]\)]; 
 (*\(Print["The evaluation below has been checked for ", 
              jj, " consecutive values of n, and is therefore correct for all n: "];\)*) 
Clear[f]; 
ff = 
          ff /. Table[
              xx[i] -> zz[\([i]\)], {i, 1, 
                jj}]; 
\ \ \ \ \ \ ff = 
          ExpandAll[Simplify[ff]]; 
ff = ff /. newnames; 
        ff = Expand[ff]; 
Return[ff];]\)

indices[vec_] := 
      Flatten[Outer[List, Apply[Sequence, Map[Range[0, # - 1] &, vec]]], 
        Length[vec] - 1]

maxmult[p_, x_] := 
      Module[{pp, t, dg, j, px}, 
  
px = 
          Expand[p]; 
pp = 
          Factor[px] /. {a_^(m_) -> t[m]}; 
        dg = Exponent[px, x]; 
For[j = dg, 
          j ≥ 0, j--, 
          If[(pp /. {t[j] -> 0}) == 0, 
            Return[j]]]; 
Return[1]]

FindDegree[f_] := 
      Module[{mat, i, j, cfd, maxsize}, 
  
maxsize = 1; 
        cfd = 2; 
While[
          cfd ≥ maxsize, 
maxsize = 
            maxsize + 3; 
mat := 
            Table[f[i + j], {i, 0, maxsize}, {j, 0, 
                maxsize}]; 
cfd = 
            maxsize + 1 - Length[NullSpace[mat]]]; 
Return[
          cfd]]

companion[f_, x_] := 
      Module[{n, aa, cc, i, j}, 
  
n = Exponent[f, x]; 
aa = 
          Table[0, {n}, {n}]; 
cc = 
          CoefficientList[f, x]; 
For[i = 2, 
          i ≤ n, i++, 
          aa[[i, i - 1]] = 1]; 
For[j = 1, 
          j ≤ n, j++, 
          aa[[1, j]] = -cc[[n + 1 - j]]/
              cc[[n + 1]]]; 
Return[aa]]

kronecker[M_, N_] := 
      Map[Flatten, 
        Flatten[Transpose[Outer[Times, M, N], {1, 3, 2, 4}], 1]]