share/contrib/celine.mac

   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
   4 Celine's method is described in Sections 4.1--4.4 of the book "A=B", by Marko Petkovsek, Herbert S. Wilf, and Doron Zeilberger.
   5 This book is available at http://www.math.rutgers.edu/~zeilberg/AeqB.pdf
   6
   7 Let 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
  11 where 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)
  12 is a rational function of n and k, celine returns the empty set. When f involves parameters (variables other than n or k), celine
  13 might make assumptions about these parameters. Using 'put' with a key of 'proviso,' Maxima saves these assumptions on the input
  14 label.
  15
  16 To use this function, first load the package integer_sequence, opsubst, and to_poly_solve.
  17
  18 Examples:
  19
  20   (%i1) celine(n!,n,k,1,0);
  21   (%o1) {fff(n+1,k)-n*fff(n,k)-fff(n,k)}
  22
  23 Check:
  24
  25   (%i2) ratsimp(minfactorial(first(%))),fff(n,k) := n!;
  26   (%o2) 0
  27
  28 An 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
  36 Check:
  37
  38   (%i5) first(%), fff(n,k) := ''(e)$
  39   (%i6) makefact(makegamma(%))$
  40
  41   (%i7) minfactorial(factor(minfactorial(factor(%))));
  42   (%o7) 0
  43
  44 The 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
  52 map('load, ["integer_sequence", "opsubst", "to_poly_solve"]);
  53
  54 celine(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)$