share/contrib/ggf.mac

   1 /*
   2 ggf.mac v1.0 for Maxima (tested with Maxima 5.9.1).
   3
   4 Compute the generating function (if it is a fraction of two
   5 polynomials) of a sequence, its first terms being given.
   6
   7 Usage:
   8   ggf(Terms);
   9     Terms - list of first terms of the sequence (integer or rational)
  10
  11 Examples:
  12   ggf(makelist(fib(n),n,0,40));
  13   ggf(makelist(2*fib(n+1)-fib(n),n,0,40));
  14
  15 Flags:
  16   GGFINFINITY
  17     default: 3 - when computing the continued fraction of the
  18     generating function, a partial quotient having a degree
  19     (strictly) greater than GGFINFINITY will be discarded and
  20     the current convergent will be considered as the exact value
  21     of the generating function; most often the degree of all
  22     partial quotients will be 0 or 1; if you use a greater value,
  23     then you should give enough terms in order to make the
  24     computation accurate enough.
  25   GGFCFMAX
  26     default: 24 - when computing the continued fraction of the
  27     generating function, if no good result has been found (see
  28     the GGFINFINITY flag) after having computed GGFCFMAX partial
  29     quotients, the generating function will be considered as
  30     not being a fraction of two polynomials and the function will
  31     exit. Put freely a greater value for more complicated
  32     generating functions.
  33
  34 Results:
  35   The solution is returned as a fraction of two polynomials.
  36   If no solution has been found, it returns with DONE.
  37
  38 History:
  39   2005-09 Thomas Baruchel - version 1.0
  40
  41 --
  42
  43 Copyright (C) 2005  Thomas Baruchel
  44
  45 This library is free software; you can redistribute it and/or
  46 modify it under the terms of the GNU Lesser General Public
  47 License as published by the Free Software Foundation; either
  48 version 2.1 of the License, or (at your option) any later version.
  49
  50 This library is distributed in the hope that it will be useful,
  51 but WITHOUT ANY WARRANTY; without even the implied warranty of
  52 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  53 Lesser General Public License for more details.
  54
  55 You should have received a copy of the GNU Lesser General Public
  56 License along with this library; if not, write to the Free Software
  57 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  58 */
  59
  60 GGFINFINITY: 3;
  61 GGFCFMAX: 24;
  62
  63 ggf(l) := block([p0:0, q0:1, p1:1, q1:0, p, q:1, i], local(l),
  64   p : sum(l[i]*'x^(i-1),i,1,length(l)),
  65   for i:1 thru GGFCFMAX do if block([j:lopow(p,'x),k:lopow(q,'x),a,p2,q2],
  66     if abs(j-k) > GGFINFINITY then return(true)
  67       else a : 'x^(j-k) * coeff(p,'x,j)/coeff(q,'x,k),
  68     p2 : a*p1 + p0, q2 : a*q1 + q0,
  69     p0 : p1, q0 : q1, p1 : p2, q1 : q2,
  70     p : rat(p - a*q),
  71     if p = 0 then return(true),
  72     p2 : p, p : q, q : p2, false)
  73     then return(ratsimp(p1/q1)));