share/contrib/devine.mac

   1 /* devine.usg
   2  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
   3  *                                                                       *
   4  * Copyright (C) 2002 Martin Rubey <Martin.Rubey@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@pap.univie.ac.at>.                  *
  21  * The translation to Maple was done by Jean-Francois Beraud             *
  22  * <Jean-Francois.Beraud@sic.sp2mi.univ-poitiers.fr> and Bruno Gauthier  *
  23  * <Bruno.Gauthier@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 A package to guess closed form for a sequence of numbers.
  31
  32 CALLING SEQUENCE:
  33
  34 guess(l, optional_args);
  35
  36 SYNOPSIS:
  37 - This  package  provides  functions  to find a closed form for a sequence,
  38   of  numbers   within  the  hierarchy  of   expressions   of   the   form,
  39   <rational function>, <product of rational function>, <product of product,
  40   of rational function>, etc.
  41
  42 EXAMPLES:
  43 guess([1,2,3]);
  44                                 [i0]
  45
  46 guess([1,4,9,16]);
  47
  48                                    2
  49                                 [i0 ]
  50
  51 guess([1,2,6,24,120]);
  52
  53                            i0 - 1
  54                            /===\
  55                             ! !
  56                           [ ! !   (i1 + 1)]
  57                             ! !
  58                            i1 = 1
  59
  60 guess(makelist(product(product(gamma(i)*i^2,i,1,j),j,1,k),k,1,8));
  61
  62                       i0 - 1   i1 - 1    i2 - 1
  63                       /===\    /===\     /===\          2
  64                        ! !      ! !       ! !   (i3 + 3)
  65                      [ ! !   4  ! !   18  ! !   ---------]
  66                        ! !      ! !       ! !    i3 + 2
  67                       i1 = 1   i2 = 1    i3 = 1
  68
  69 guess([1,2,7,42,429,7436,218348,10850216]);
  70
  71                     i0 - 1   i1 - 1
  72                     /===\    /===\
  73                      ! !      ! !   3 (3 i2 + 2) (3 i2 + 4)
  74                    [ ! !   2  ! !   -----------------------]
  75                      ! !      ! !   4 (2 i2 + 1) (2 i2 + 3)
  76                     i1 = 1   i2 = 1
  77
  78
  79
  80 guess(makelist(k^3+k^2,k,1,7));
  81
  82
  83 Dependent equations eliminated:  (6)
  84                        i0 - 1
  85                        /===\
  86          2              ! !                       5040
  87       [i0  (i0 + 1), 2  ! !   (- --------------------------------------),
  88                         ! !        4        3         2
  89                        i1 = 1    i1  - 24 i1  + 245 i1  - 1422 i1 + 360
  90
  91                                                       i0 - 1
  92                                                       /===\
  93                                                        ! !   (i1 + 1) (i1 + 2)
  94                                                     2  ! !   -----------------]
  95                                                        ! !            2
  96                                                       i1 = 1        i1
  97
  98 Note that the last example produces three solutions. The first and the last are
  99 equivalent, but the second is different! In this case,
 100
 101 guess(makelist(k^3+k^2,k,1,7),1);
 102
 103 or
 104
 105 guess(makelist(k^3+k^2,k,1,7),"one");
 106
 107                           2
 108 find only the solution i0  (i0 + 1), which is a rational function, and is also
 109 the first function guess finds.
 110
 111 PARAMETERS:
 112   l               - a list of numbers,
 113   level           - an integer (optional),
 114   "one"           - the string "one" (optional),
 115   "nogamma"       - the string "nogamma" (optional),
 116
 117 SYNOPSIS:,
 118 - guess(l) tries to find a closed form for a sequence within the hierarchy,
 119   of expressions  of  the  form  <rational function>, <product of rational,
 120   function>, <product of product of rational function>, etc.
 121
 122 - guess(l,level) does the same thing as guess(l) but it searches only
 123   within the first 'level' levels
 124
 125 - guess(l,"one") does the same thing as guess(l) but it returns the first
 126   solution it finds.
 127
 128 - guess(l,"nogamma") does the same thing as guess(l) but it returns
 129   expressions without gamma functions. In fact, there is not much difference
 130   just at the moment, because Maxima doesn't simplify products yet...
 131 */
 132
 133
 134 /* devine.mac -*- mode: Maxima; -*-
 135  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 136  *                                                                       *
 137  * Copyright (C) 2002 Martin Rubey <Martin.Rubey@LaBRI.fr>               *
 138  *                                                                       *
 139  * This file is part of GNU Maxima.                                      *
 140  *                                                                       *
 141  * This program is free software; you can redistribute it and/or         *
 142  * modify it under the terms of the GNU General Public License as        *
 143  * published by the Free Software Foundation; either version 2 of        *
 144  * the License, or (at your option) any later version.                   *
 145  *                                                                       *
 146  * This program is distributed in the hope that it will be               *
 147  * useful, but WITHOUT ANY WARRANTY; without even the implied            *
 148  * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR               *
 149  * PURPOSE. See the GNU General Public License for more details.         *
 150  *                                                                       *
 151  * History:                                                              *
 152  * This is a translation of the Mathematica package Rate.m               *
 153  * by Christian Krattenthaler <Kratt@pap.univie.ac.at>.                  *
 154  * The translation to Maple was done by Jean-Francois Beraud             *
 155  * <Jean-Francois.Beraud@sic.sp2mi.univ-poitiers.fr> and Bruno Gauthier  *
 156  * <Bruno.Gauthier@univ-mlv.fr>                                          *
 157  *                                                                       *
 158  * All features of this package are due to C. Krattenthaler              *
 159  * The help text is due to Jean-Francois Beraud and Bruno Gauthier       *
 160  *                                                                       *
 161  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 162  */
 163
 164 /*
 165  * Rational Interpolation
 166  * Gives the rational interpolant to the set of points defined by xlist and
 167  * ylist, where m and k are respectively the degrees of the numerator and
 168  * denominator, and xlist is a list of m+k+1 abscissas of the interpolation
 169  * points, x is the variable the result is supposed to be a function of.
 170  */
 171 rationalinterpolation(xlist, ylist, x, m, k) :=
 172 block([tempvec : makelist(1, i, 1, m+k+1), /* contains the new row of mat */
 173        rowlist,                               /* first set of rows of mat */
 174        mat,            /* matrix that describes the interpolation problem */
 175        varlist : makelist('x[i], i, 1, m+k+2)],
 176   mode_declare([tempvec,rowlist,varlist,mat],list,[m,k],fixnum),
 177
 178   if max(m, k) > 0
 179   then rowlist : cons(tempvec,
 180                       makelist(tempvec : tempvec * xlist, i, 1, max(m, k)))
 181   else rowlist : [tempvec],
 182
 183   mat : transpose(apply(matrix,
 184                         append(rest(rowlist, -(max(m, k) - m) ),
 185                                -1 * makelist(rowlist[i] * ylist,
 186                                              i, 1, k + 1)))),
 187   mat : ev(mat . varlist, scalarmatrixp : false),
 188
 189 /* not sure whether it is save to modify xlist... */
 190   xlist : linsolve(makelist(mat[i, 1], i, 1, (m+k)+1), varlist),
 191   if length(xlist) = 0
 192 /* something went wrong */
 193   then 'null
 194 /* use the solution to define the interpolating rational function */
 195   else factor(subst(xlist, sum('x[i+1]*x^i, i, 0, m)
 196                            /sum('x[(i+m)+2]*x^i, i, 0, k))));
 197
 198
 199 /* Intermediate function */
 200 guesscons(l, t) :=
 201 block([lsize : length(l), res : [], x, ri],
 202   mode_declare(lsize, fixnum, res, list, ri, any),
 203
 204   for i : 0 thru lsize-2 do
 205      (ri : rationalinterpolation(makelist(k, k, 1, lsize-1), rest(l,-1),
 206                                  x, (lsize-2)-i, i),
 207       if ri # 'null
 208       then if (subst(x=lsize, denom(ri)) # 0)
 209               and
 210               (subst(x=lsize, ri)-last(l) = 0)
 211               and
 212               not member(subst(x=t, ri), res)
 213            then res : cons(subst(x=t, ri), res)),
 214   res);
 215
 216 /*
 217  * Main function of the package
 218  * it tries to find a closed form  for a sequence
 219  * within the hierarchy of expressions of the
 220  * form <rational function>, <product of rational functions>,
 221  * <product of product of rational functions>, etc. It may
 222  * give several answers
 223  */
 224 guess(l, [optargs]) :=
 225 block([lsize : length(l),
 226        tempres, maxlevel,
 227        maxlevellist : sublist(optargs, numberp),
 228        res : [],
 229        onep : member("one", optargs),
 230        unevp : member("nogamma", optargs), g],
 231       mode_declare([lsize, maxlevel], fixnum,
 232                    [tempres, maxlevellist, res], list,
 233                    [onep, unevp], boolean, g, any),
 234
 235       optargs : delete("nogamma", delete("one", optargs, 1), 1),
 236       if length(maxlevellist) > 1 or length(optargs)-length(maxlevellist) > 0
 237       then error("Wrong number of optional arguments: ", optargs)
 238       else maxlevel : mode_identity(fixnum, apply(min, cons(lsize-1, maxlevellist)) - 1),
 239
 240       g : make_array('any, maxlevel + 1),
 241
 242       for k : 0 thru maxlevel do
 243          (g[k] : l,
 244           l : makelist(l[i+1]/l[i], i, 1, (lsize-k)-1),
 245
 246           tempres : guesscons(g[k], concat('i, k)),
 247           if tempres # []
 248           then (if k > 0
 249                 then for i : 1 thru k do
 250                          if unevp
 251                          then tempres :
 252                                   block ([j : concat('i, (k-i)+1)],
 253                                         map(lambda([expr],
 254                                                    g[k-i][1] *
 255                                                    apply (nounify (product), [expr, j, 1, concat('i, k-i)-1])),
 256                                             tempres))
 257                          else tempres :
 258                                   block ([j : concat('i, (k-i)+1)],
 259                                         map(lambda([expr],
 260                                                    g[k-i][1] *
 261                                                    apply (verbify (product), [expr, j, 1, concat('i, k-i)-1])),
 262                                             tempres)),
 263                 res : append(res, tempres),
 264                 if onep then return(res))),
 265       res);