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