share/contrib/stirling.mac

   1 /* Copyright 2006 by Barton Willis
   2
   3   This is free software; you can redistribute it and/or
   4   modify it under the terms of the GNU General Public License,
   5   http://www.gnu.org/copyleft/gpl.html.
   6
   7  This software has NO WARRANTY, not even the implied warranty of
   8  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9 */
  10
  11 /* See A&S 6.1.37 and 6.1.40
  12
  13 stirling(e) replaces all gamma functions in the expression 'e' with the zeroth
  14 order Stirling expansion.
  15
  16 stirling(e, n) replaces all gamma functions in the expression 'e' with the O(1/x^(2n))
  17 Stirling expansion.
  18
  19 stirling(e, n, pred) replaces all expressions of the form gamma(q) in 'e' with
  20 the O(1/x^(2n)) Stirling expansion provided pred(q) evaluates to true. When
  21 pred(q) is false, the Stirling expansion is not applied to a term of the form
  22 gamma(q). If you want to use the optional argument 'pred,' you must specify
  23 the value of 'n.'
  24
  25 Examples:
  26
  27 (%i1) stirling(gamma(%alpha+x)/gamma(x),1);
  28 (%o1) %e^(-%alpha)*x^(1/2-x)*(x+%alpha)^(x+%alpha-1/2)
  29
  30 (%i2) taylor(%,x,inf,1);
  31 (%o2) x^%alpha+(x^%alpha*%alpha^2-x^%alpha*%alpha)/(2*x)+...
  32
  33 (%i3) map('factor,%);
  34 (%o3) x^%alpha+((%alpha-1)*%alpha*x^(%alpha-1))/2
  35
  36 The function 'stirling' knows the difference between the variable 'gamma' and
  37 the function gamma:
  38
  39 (%i4) stirling(gamma(gamma));
  40 (%o4) sqrt(2)*sqrt(%pi)*gamma^(gamma-1/2)*%e^(-gamma)
  41
  42 To conditionally apply a Stirling expansion, give 'stirling' a third
  43 argument; for example, to apply the expansion only when the argument to
  44 the gamma function is a symbol, you may use:
  45
  46 (%i5) stirling(gamma(a) / gamma(a+b),1, symbolp);
  47 (%o5) (sqrt(2)*sqrt(%pi)*a^(a-1/2)*%e^(-a))/gamma(b+a)
  48
  49 To apply the Stirling expansion only when the argument to gamma does
  50 not depend on the symbol 'a,' you may use:
  51
  52 (%i6) stirling(gamma(a) / gamma(a+b), 5, lambda([s], freeof('a,s)));
  53 (%o6) gamma(a)/gamma(b+a)
  54
  55 To use this code, the file 'opsubst.lisp' must be in your search path.
  56
  57
  58 */
  59
  60 load("opsubst");
  61
  62 stirling(z, [a]) := block([n, s, acc : 0, %m, pred, zerobern : true],
  63
  64        n : if a = [] then 1 else first(a),
  65
  66        pred : if (a = [] or length(a) < 2) then buildq([],lambda([s], true)) else
  67           if length(a) = 2 then second(a) else
  68            error("The function 'stirling' requires one or two arguments"),
  69
  70        if integerp(n) and n > -1 then (
  71           n : n - 1,
  72           for %m from 1 thru n do (
  73                acc : acc + bern(2*%m)/(2*%m *(2*%m-1)*s^(2*%m-1))),
  74           acc : exp(-s)* s^(s-1/2) * sqrt(2*%pi) * exp(acc),
  75           opsubst(buildq([acc,pred], lambda([s], if pred(s) then acc else funmake(gamma, [s]))), 'gamma, z))
  76        else error("The second argument to 'stirling' must be a nonnegative integer"));
  77