1
$\begingroup$

Wolfram isn't helping me much so I am curious if there are other programs out there.

I don't know what degree it is, but I have a series of numbers and I'd like to determine the linear recurrence function/coefficients so I can find any value I want via matrix exponentiation.

Are there such programs for interpolating, and how would I do it?

  • 0
    @Manos such as what? Solver? Sometimes Solver can't find an answer even when I know one exists (to which it'll find the solution if I bootstrap the numbers first)2012-11-18

3 Answers 3

3

You want the book The Book of Numbers by Conway and Guy. A description of the algorithm is at ALGORITHM but is not typeset.

  • 0
    When I first had seen Conways method (somewhere in the 70ies or 80ies in the german version of the Scientific American, I think) I was much impressed, but did never check it thoroughly and so never really used it. I'd like it, if it could be expressed using the matrix-ansatz in such a simple translation. Perhaps tomoorow I'll take a deeper look into it...2012-11-19
2

It can easily be done using for instance the freely available program Pari/GP. Here is a short excerpt of an example session.
We define the linear recurrence $x_{k+1} = B x_k + A x_{k-1} \quad $
For an example we set $[A,B]=[1,4] \quad $, then produce testdata beginning with $[x_1,x_2]=[1,2] \quad $. This is the vector of example data $[9, 38, 161, 682, 2889, 12238] \quad $ for which we want to find the recurrence-coefficients A and B (and possibly more if we do not know the order of recurrence beforehand).

Given for instance 6 numbers which are consecutive iterations, we could restate the problem as a 3x3-matrix equation: $ \small \begin{array} {} & \cdot & \begin{bmatrix} A \\ B \\ C \end{bmatrix} \\ \begin{bmatrix} 9 & 38 & 161 \\ 38&161&682 \\ 161&682&2889 \end{bmatrix} & = & \begin{bmatrix} 682 \\ 2889 \\ 12238 \end{bmatrix}\end{array}$ If the lhs-matrix has smaller rank than 3, then we need only two coefficients $[A,B] \quad $ instead of three $[A,B,C] \quad $.

The Pari/GP code to solve the problem:

\\pari/gp   \\ leading ? indicates a command to Pari/GP,  \\ leading % and a number indicates the (numbered) answer   \\ we introduce the given vector of numbers ? v= [9, 38, 161, 682, 2889, 12238] %141 = [9, 38, 161, 682, 2889, 12238]  \\ of six consecutive values we can generate a 3x4-matrix ? M=matrix(3,4,r,c,v[r+c-1]) %142 =  [9 38 161 682] [38 161 682 2889] [161 682 2889 12238]  \\ the rank of the matrix gives the order of the linear recurrence ? rk=matrank(M)  %143 = 2  \\ so we'll have two coefficients and need a 2x2 and a 2x1 matrix to solve ? MX=matrix(rk,rk,r,c,M[r,c])  %144 =  [9 38] [38 161]  ? MY=matrix(rk,1,r,c,M[r,rk+1])  %145 =  [161] [682]  \\ the following command solves the problem: ? AB = MX^-1 * MY  %146 =  [1] [4] 

The resulting matrix contains the coefficients A and B.

As a complete function this could be written:

{solvelinrec(v)=local(rs,cs,M,rk,MX,MY,AB);   rs=#v \ 2; cs = rs+1 ;   M=matrix(rs,cs,r,c,v[r+c-1]);   rk=matrank(M);   MX=matrix(rk,rk,r,c,M[r,c]);   MY=vectorv(rk,r,M[r,rk+1]);   AB = matsolve(MX,MY); \\ MX^-1 * MY ;  return(AB); } 

and then

? solvelinrec([9, 38, 161, 682, 2889, 12238]) %173 = [1, 4]~ 

or, receiving a order-three-recurrent sequence:

 ? solvelinrec([4, 5, 13, 11, 46, 11, 184, -94, 841, -1033, 4303, -7594])  %178 = [1, 4, -1]~ 

thus in this case we have $x_k = -1 x_{k-1} + 4 x_{k-2} + 1 x_{k-3}$

0

There are more powerful tools in Pari/GP, for example the LLL algorithm. Here is a small program that eats a vector of integers and spits out the rational generating function or 0 if it can't find it. Note that the coefficients of the denominator polynom determine the coefficients of the linear recurrence:

ggf(v)=local(l,m,p,q,B);l=length(v);B=floor(l/2);if(B<3,return(0));m=matrix(B,B,x,y,v[x-y+B+1]);q=qflll(m,4)[1];if(length(q)==0,return(0));p=sum(k=1,B,x^(k-1)*q[k,1]);q=Pol(Pol(vector(l,n,v[l-n+1]))*p+O(x^(B+1)));if(polcoeff(p,0)<0,q=-q;p=-p);q=q/p;p=Ser(q+O(x^(l+1)));for(m=1,l,if(polcoeff(p,m-1)!=v[m],return(0)));q 

Example:

? ggf([1,3,7,15,31,63,127,255]) 1/(2*x^2 - 3*x + 1) 

So, $a_{n}=3a_{n-1}-2a_{n-2},$ starting $a_0=1, a_2=3$.

Note also that my program needs $2d$ values to determine the coefficients of a recurrence of order $d$. This is just to have a check on the result and could perhaps be optimized if really necessary.

Have fun.