| 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