1module hypergeomrsolve;
2
3% Redistribution and use in source and binary forms, with or without
4% modification, are permitted provided that the following conditions are met:
5%
6%    * Redistributions of source code must retain the relevant copyright
7%      notice, this list of conditions and the following disclaimer.
8%    * Redistributions in binary form must reproduce the above copyright
9%      notice, this list of conditions and the following disclaimer in the
10%      documentation and/or other materials provided with the distribution.
11%
12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
16% CONTRIBUTORS
17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
23% POSSIBILITY OF SUCH DAMAGE.
24%
25
26
27fluid '(!*tracefps);
28
29algebraic procedure hypergeomrsolve (r,k,a0);
30
31% solves the recurrence equation
32%
33%  a(k+1) = r(k) * a(k), a(0) = a0
34
35  begin scalar re,nnn,ddd,c,p,q,ak,sols,ii;
36
37   p := {}; q := {};
38   c := 1;
39
40   re := r * (k + 1);
41
42   nnn := old_factorize num re; ddd := old_factorize den re;
43
44   foreach nn in nnn do
45   if freeof (nn,k) then c := c * nn
46        else if deg(nn,k) =1 then
47        << c:= c*coeffn(nn,k,1);
48            p:= append (p,list(coeffn(nn,k,0)/coeffn(nn,k,1)))>>
49        else if deg(nn,k) =2 then
50        << c := c *  lcof(nn,k);
51                sols := solve(nn,k);
52                for each s in sols do
53                   << for i:=1:first multiplicities!* do
54                        p:= (- rhs s) . p;
55                      multiplicities!* := rest multiplicities!*;
56                  >>
57        >>
58        else rederr(" hypergeomRsolve failed");
59
60   foreach dd in ddd do
61   if freeof (dd,k) then c := c / dd
62        else if deg(dd,k) =1 then
63        << c:= c/coeffn(dd,k,1);
64            q:= append (q,list(coeffn(dd,k,0)/coeffn(dd,k,1)))>>
65        else if deg(dd,k) =2 then
66             << c := c /  lcof(dd,k);
67                sols := solve(dd,k);
68                for each s in sols do
69                   << for i:=1:first multiplicities!* do
70                        q:= (- rhs s) . q;
71                      multiplicities!* := rest multiplicities!*;
72                  >>;
73             >>
74        else rederr(" hypergeomRsolve failed");
75
76   rsolve := infinite;
77   for each s in p do if fixp s and s < 0 then rsolve := finite;
78   if symbolic !*tracefps then  write "RSOLVE  = ",rsolve;
79
80   p := for each s in p product Pochhammer(s,k);
81   q := for each s in q product Pochhammer(s,k);
82
83   ak := a0 * (c^k) * p/(q * factorial k);
84
85% Do additional simplification here??
86
87   return ak;
88end;
89
90% A special ruleset for powerseries; nn has a special meaning here and
91% should be treated as integer
92
93algebraic <<
94
95hgspec_pochhammer :=
96
97{ Pochhammer(~kk,~nn) => 1 when nn=0,
98  Pochhammer(~kk,nn) => 0 when kk = 0,
99  Pochhammer(~kk,nn) => (-1)^nn * factorial(-kk)/factorial(-kk-nn)
100                        when fixp(kk) and kk <=0,
101  Pochhammer(~kk,nn) => factorial(kk+nn-1)/factorial(kk-1) when fixp kk,
102  Pochhammer(~kk,~w*nn) =>
103     factorial(kk+w*nn-1)/factorial(kk-1) when fixp kk}
104>>;
105
106endmodule;
107end;
108
109
110