Michael Somos on Wed, 24 Mar 2004 20:18:38 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Solving recurrences.


Marko Riedel wrote :

> I read somewhere that there is a PARI package that can solve
> recurrences of the type
>
> P_2(n) a(n+2) + P_1(n) a(n+1) + P_0(n) a(n) = 0

Yes, that is easily done. I wrote a small function for myself :

/* ------------------------------------------------------------------------- */
/* find recurrence relation for a sequence an=[a1,a2,...] */
{vecrec(n1,n2,an,oo)=
   if(n1*n2>oo,error("myrec: more terms"));
   mm=matrix(n1*n2,oo);
   for(i=0,n1-1,for(j=0,n2-1,
   for(k=0,oo-n2, mm[1+i+n1*j,1+k]=an[k+n2-j]*(k+n2)^i )));
   matker(mm~)~
} /* end vecrec() */
/* ------------------------------------------------------------------------- */
/* print recurrence relation found for sequence ( choice of org=origin ) */
{prrec(n1,n2,v,org=1)=
   for(j=0,n2-1,print("+a(n-"j")*("sum(i=0,n1-1,v[1+i+n1*j]*(n-org+1)^i)")")) }
/* ------------------------------------------------------------------------- */

parisize = 4000000, primelimit = 500000
? an=vector(6,n,(2*n+1)!*(n+1)/n!/2^n)
%1 = [6, 45, 420, 4725, 62370, 945945]
? vecrec(3,2,an,6)
%2 = 
[0 -1/2 0 1/2 3/2 1]

? prrec(3,2,%[1,]*2)
+a(n-0)*(-n)
+a(n-1)*(2*n^2 + 3*n + 1)

It works for me. I know there could be better functions. Shalom, Michael