share/contrib/guess/guess.mac

   1 /* guess.mac -*- mode: Maxima; -*-
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                       *
   4  * Copyright (C) 2002 Martin Rubey <Martin.Rubey at LaBRI.fr>            *
   5  *                                                                       *
   6  * This file is part of GNU Maxima.                                      *
   7  *                                                                       *
   8  * This program is free software; you can redistribute it and/or         *
   9  * modify it under the terms of the GNU General Public License as        *
  10  * published by the Free Software Foundation; either version 2 of        *
  11  * the License, or (at your option) any later version.                   *
  12  *                                                                       *
  13  * This program is distributed in the hope that it will be               *
  14  * useful, but WITHOUT ANY WARRANTY; without even the implied            *
  15  * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR               *
  16  * PURPOSE. See the GNU General Public License for more details.         *
  17  *                                                                       *
  18  * History:                                                              *
  19  * This is a translation of the Mathematica package Rate.m               *
  20  * by Christian Krattenthaler <Kratt at pap.univie.ac.at>.               *
  21  * The translation to Maple was done by Jean-Francois Beraud             *
  22  * <Jean-Francois.Beraud at sic.sp2mi.univ-poitiers.fr> and Bruno Gauthier
  23  * <Bruno.Gauthier at univ-mlv.fr>                                       *
  24  *                                                                       *
  25  * All features of this package are due to C. Krattenthaler              *
  26  * The help text is due to Jean-Francois Beraud and Bruno Gauthier       *
  27  *                                                                       *
  28  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  29  */
  30
  31 /*
  32  * Rational Interpolation
  33  * Gives the rational interpolant to f (a function),
  34  * where m and k are respectively
  35  * the degrees of the numerator and denominator,
  36  * and xlist is a list of m+k+1 abscissas of the
  37  * interpolation points
  38  */
  39
  40 rationalinterpolation(f, x, m, k, xlist) :=
  41 block([fx : map(f, xlist),                 /* values of f at given points */
  42        tempvec : makelist(1, i, 1, m+k+1), /* contains the new row of mat */
  43        rowlist,                               /* first set of rows of mat */
  44        mat,            /* matrix that describes the interpolation problem */
  45        varlist : makelist('x[i], i, 1, m+k+2)],
  46
  47   if max(m, k) > 0
  48   then rowlist : cons(tempvec,
  49                       makelist(tempvec : tempvec * xlist, i, 1, max(m, k)))
  50   else rowlist : [tempvec],
  51
  52   mat : transpose(apply(matrix,
  53                         append(rest(rowlist, -(max(m, k) - m) ),
  54                                -1 * makelist(rowlist[i] * fx, i, 1, k + 1)))),
  55   mat : ev(mat . varlist, scalarmatrixp : false),
  56
  57 /* not sure whether it is save to modify xlist... */
  58   xlist : linsolve(makelist(mat[i, 1], i, 1, m+k+1), varlist),
  59   if length(xlist) = 0
  60 /* something went wrong */
  61   then NULL
  62 /* use the solution to define the interpolating rational function */
  63   else factor(subst(xlist, sum('x[i+1]*x^i, i, 0, m)
  64                            /sum('x[i+m+2]*x^i, i, 0, k))));
  65
  66
  67 /* Intermediate function */
  68 guesscons(l, t) :=
  69 block([lsize : length(l), res : [], x, ri],
  70
  71   for i : 0 thru lsize-2 do
  72      (ri : rationalinterpolation(lambda([x], l[x]), x,
  73                                  (lsize-2)-i, i,
  74                                  makelist(k, k, 1, lsize-1)),
  75       if ri # NULL
  76       then if (subst(x=lsize, denom(ri)) # 0)
  77               and
  78               (subst(x=lsize, ri)-last(l) = 0)
  79               and
  80               not member(subst(x=t, ri), res)
  81            then res : cons(subst(x=t, ri), res)),
  82   res);
  83
  84 /*
  85  * Main function of the package
  86  * it tries to find a closed form  for a sequence
  87  * within the hierarchy of expressions of the
  88  * form <rational function>, <product of rational functions>,
  89  * <product of product of rational functions>, etc. It may
  90  * give several answers
  91  */
  92 guess(l, [optargs]) :=
  93 block([lsize : length(l),
  94        tempres, maxlevel,
  95        maxlevellist : sublist(optargs, numberp),
  96        res : [],
  97        onep : member("one", optargs),
  98        unevp : member("nogamma", optargs), g],
  99
 100       optargs : delete("nogamma", delete("one", optargs, 1), 1),
 101       if length(maxlevellist) > 1 or length(optargs)-length(maxlevellist) > 0
 102       then error("Wrong number of optional arguments: ", optargs)
 103       else maxlevel : apply(min, cons(lsize-1, maxlevellist)) - 1,
 104
 105       g : make_array('flonum, maxlevel + 1),
 106
 107       for k : 0 thru maxlevel do
 108          (g[k] : l,
 109           l : makelist(l[i+1]/l[i], i, 1, lsize-k-1),
 110
 111           tempres : guesscons(g[k], concat('i, k)),
 112           if tempres # []
 113           then (if k > 0
 114                 then for i : 1 thru k do
 115                          if unevp
 116                          then tempres :
 117                                   subst('j = concat('i, (k-i)+1),
 118                                         map(lambda([expr],
 119                                                    g[k-i][1] *
 120                                                    'product(expr, j, 1,
 121                                                             concat('i, k-i)-1)),
 122                                             tempres))
 123                          else tempres :
 124                                   subst('j = concat('i, (k-i)+1),
 125                                         map(lambda([expr],
 126                                                    g[k-i][1] *
 127                                                    product(expr, j, 1,
 128                                                            concat('i, k-i)-1)),
 129                                             tempres)),
 130                 res : append(res, tempres),
 131                 if onep then return())),
 132       res);