1/* Maxima implementation of Sister Celine's method
2    Barton Willis wrote this code. It is released under the Creative Commons CC0 license (https://creativecommons.org/about/cc0)
3
4Celine's method is described in Sections 4.1--4.4 of the book "A=B", by Marko Petkovsek, Herbert S. Wilf, and Doron Zeilberger.
5This book is available at http://www.math.rutgers.edu/~zeilberg/AeqB.pdf
6
7Let f = F(n,k). The function celine returns a set of recursion relations for F of the form
8
9    p_0(n) * fff(n,k) + p_1(n) * fff(n+1,k) + ... +  p_p(n) * fff(n+p,k+q),
10
11where p_0 through p_p are polynomials. If Maxima is unable to determine that sum(sum(a(i,j) * F(n+i,k+j),i,0,p),j,0,q) / F(n,k)
12is a rational function of n and k, celine returns the empty set. When f involves parameters (variables other than n or k), celine
13might make assumptions about these parameters. Using 'put' with a key of 'proviso,' Maxima saves these assumptions on the input
14label.
15
16To use this function, first load the package integer_sequence, opsubst, and to_poly_solve.
17
18Examples:
19
20  (%i1) celine(n!,n,k,1,0);
21  (%o1) {fff(n+1,k)-n*fff(n,k)-fff(n,k)}
22
23Check:
24
25  (%i2) ratsimp(minfactorial(first(%))),fff(n,k) := n!;
26  (%o2) 0
27
28An example with parameters:
29
30  (%i3) e : pochhammer(a,k) * pochhammer(-k,n) / (pochhammer(b,k));
31  (%o3) (pochhammer(a,k)*pochhammer(-k,n))/pochhammer(b,k)
32
33  (%i4) recur : celine(e,n,k,2,1);
34  (%o4) {fff(n+2,k+1)-fff(n+2,k)-b*fff(n+1,k+1)+n*(-fff(n+1,k+1)+2*fff(n+1,k)-a*fff(n,k)-fff(n,k))+a*(fff(n+1,k)-fff(n,k))+2*fff(n+1,k)-n^2*fff(n,k)}
35
36Check:
37
38  (%i5) first(%), fff(n,k) := ''(e)$
39  (%i6) makefact(makegamma(%))$
40
41  (%i7) minfactorial(factor(minfactorial(factor(%))));
42  (%o7) 0
43
44The proviso data suggests that setting a = b may result in a lower order recursion
45
46  (%i8) get('%i4,'proviso);
47  (%o8) (-(b-1)*(b-a)*n*(n+a-1)#0) %and ((b-1)*(b-a)*n*(n+a-1)#0) %and (n-b+a#0)
48
49  (%i9) celine(subst(b=a,e),n,k,1,1);
50  (%o9) {fff(n+1,k+1)-fff(n+1,k)+n*fff(n,k)+fff(n,k)} */
51
52map('load, ["integer_sequence", "opsubst", "to_poly_solve"]);
53
54celine(f,n,k,p,q) := block([e, recur, v, sol : set(), fff, ratmx : false, mat,cnd],
55   f : makefact(makegamma(f)),
56   p : n .. (n + p),
57   q : k .. (k + q),
58   v : outermap(lambda([i,j], gensym()), p, q),
59   recur : xreduce("+", outermap('fff, p, q) . v),
60   e : xreduce("+", outermap(lambda([i,j], subst([n=i, k=j], f)), p, q) . v),
61   v : xreduce('append,v),
62   e : minfactorial(factor(e/f)),
63   if polynomialp(ratnum(e),[n,k], lambda([s], freeof(n,k,s))) and polynomialp(ratdenom(e),[n,k],lambda([s], freeof(n,k,s))) then (
64       e : rat(ratnum(e)),
65       e : map(lambda([i], ratcoeff(e,k,i)), 0 .. hipow(e,k)),
66       mat : triangularize(coefmatrix(e,v)),
67       cnd : map(lambda([s], delete(0,s)), args(mat)),
68       cnd : xreduce("%and", map(lambda([s], factor(first(s)) # 0),cnd)),
69       sol :  block([scalarmatrixp : false], algsys(xreduce('append, args(mat . transpose(v))),v)),
70       recur : expand(subst(sol, recur)),
71       recur : setify(map(lambda([s], coeff(recur,s)), %rnum_list)),
72       recur : map(lambda([s], block([ar : gatherargs(s, 'fff)], /* standardize to minimum argument of n & k. */
73                 expand(subst([n = 2*n - lmin(map('first,ar)), k = 2*k - lmin(map('second,ar))],s)))),recur),
74       sol : map(lambda([s], block([ar : gatherargs(s, 'fff)], /* standardize leading coefficient to one.*/
75                   ratsimp(s / coeff(s, funmake('fff,last(sort(ar))))))),recur)),
76  put(first(labels), cnd, 'proviso),
77  sol)$
78