hgrsolve.red 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. module hypergeomRsolve;
  2. fluid '(!*tracefps);
  3. algebraic procedure hypergeomRsolve (r,k,a0);
  4. % solves the recurrence equation
  5. %
  6. % a(k+1) = r(k) * a(k), a(0) = a0
  7. begin scalar Re,NNN,DDD,c,p,q,ak,sols,II;
  8. P := {}; Q := {};
  9. C := 1;
  10. Re := R * (k + 1);
  11. NNN := old_factorize num Re; DDD := old_factorize den re;
  12. foreach nn in NNN do
  13. if freeof (nn,k) then c := c * nn
  14. else if deg(nn,k) =1 then
  15. << C:= c*coeffn(nn,k,1);
  16. P:= append (p,list(coeffn(nn,k,0)/coeffn(nn,k,1)))>>
  17. else if deg(nn,k) =2 then
  18. << c := c * lcof(nn,k);
  19. sols := solve(nn,k);
  20. for each s in sols do
  21. << for i:=1:first multiplicities!* do
  22. P:= (- rhs s) . p;
  23. multiplicities!* := rest multiplicities!*;
  24. >>
  25. >>
  26. else rederr(" hypergeomRsolve failed");
  27. foreach dd in DDD do
  28. if freeof (DD,k) then c := c / dd
  29. else if deg(DD,k) =1 then
  30. << C:= C/coeffn(dd,k,1);
  31. Q:= append (Q,list(coeffn(dd,k,0)/coeffn(dd,k,1)))>>
  32. else if deg(DD,k) =2 then
  33. << c := c / lcof(dd,k);
  34. sols := solve(dd,k);
  35. for each s in sols do
  36. << for i:=1:first multiplicities!* do
  37. Q:= (- rhs s) . Q;
  38. multiplicities!* := rest multiplicities!*;
  39. >>;
  40. >>
  41. else rederr(" hypergeomRsolve failed");
  42. RSOLVE := infinite;
  43. for each s in P do if fixp s and s < 0 then RSOLVE := finite;
  44. if symbolic !*traceFPS then write "RSOLVE = ",RSOLVE;
  45. P := for each s in P product pochhammer(s,k);
  46. Q := for each s in Q product pochhammer(s,k);
  47. ak := a0 * (C^k) * P/(q * factorial k);
  48. % Do additional simplification here??
  49. return ak;
  50. end;
  51. % A special ruleset for powerseries; nn has a special meaning here and
  52. % should be treated as integer
  53. algebraic <<
  54. hgspec_pochhammer :=
  55. { Pochhammer(~kk,~nn) => 1 when nn=0,
  56. Pochhammer(~kk,nn) => 0 when kk = 0,
  57. Pochhammer(~kk,nn) => (-1)^nn * factorial(-kk)/factorial(-kk-nn)
  58. when fixp(kk) and kk <=0,
  59. Pochhammer(~kk,nn) => factorial(kk+nn-1)/factorial(kk-1) when fixp kk,
  60. Pochhammer(~kk,~w*nn) =>
  61. factorial(kk+w*nn-1)/factorial(kk-1) when fixp kk}
  62. >>;
  63. endmodule;
  64. end;