share/gamma_simp/gamma_simp.mac

   1 /* Author: Barton Willis, copyright 2022
   2
   3 Maxima code for simplifying expressions that involve products of
   4 gamma functions.
   5
   6 See LICENSE.md for licensing information.
   7 */
   8
   9  /* The CL translator suggests adding the following: */
  10  eval_when(translate,declare_translated(gamma_arg_increase,gamma_arg_decrease,
  11          gamma_product_simp,homogenize_gamma,gamma_magnitude_simp))$
  12
  13 /* Load opsubst so we can use gatherargs. */
  14  load("opsubst");
  15
  16 /* This function is inflag sensitive. */
  17 mtimesp(e) := not mapatom(e) and part(e,0) = "*";
  18
  19   /* Return a lambda form that generates the Gauss gamma function
  20  product identity; see http://dlmf.nist.gov/5.5.E6 */
  21  make_gamma_product_id() := buildq([x : gensym(), k : gensym(),n : gensym()],
  22     lambda([x,n], prod(gamma(x+k/n),k,0,n-1) = gamma(n*x)/((2*%pi)^((1-n)/2)*n^(n*x-1/2))));
  23
  24 /* Apply ratsimp only to the arguments of gamma. Doing this is an effort to allow
  25    ratsubst to make more matches to gamma expressions without dispatching
  26    ratsimp on the entire expression.*/
  27 ratsimp_gamma_args(e) := block([keepfloat : true],
  28   if mapatom(e) then e
  29   elseif inpart(e,0) = 'gamma then gamma(ratsimp(first(e)))
  30   else map('ratsimp_gamma_args, e));
  31
  32  /* The functions gamma_simp and factorial_simp are the only two user level
  33 functions in this file.
  34
  35 The function ratsubst does not respect the value of keepfloat--this
  36 is its documented behavior. The gamma_simp code is a heavy user of
  37 ratsubst--fixing this bugs that are a result of ratsubst converting
  38 floats to rationals likely isn't worthwhile. The local assignments
  39 keepfloat : true don't prevent all float to rational conversions.
  40
  41 When radsubstflag is true, apparently ratsubst sometimes introduces some
  42 complex exponentials--as far as I know, the answers are OK, just messy.
  43 Instead of limiting a users freedom by locally setting radsubstflag to
  44 false, I'm going to make that setting in the rtest_gamma_simp.*/
  45
  46 gamma_simp(e) := block([prederror : false, opsubst : true,keepfloat : true,
  47    inflag : true],
  48   e : ratsimp_gamma_args(e),
  49   e : homogenize_gamma(e),
  50   e : subst("*" = lambda([[q]], gamma_magnitude_simp(xreduce("*",q))),e),
  51   e : gamma_product_simp(e),
  52   e : facsum_gamma_terms(e),
  53   e : gamma_arg_increase(gamma_arg_decrease(e)));
  54
  55 /* This function simplifies a product of two or more gamma functions. */
  56  gamma_product_simp(e) := block([a,N,id,ee,dn,fn,keepfloat : true],
  57     ee : e,
  58     a : map('first, gatherargs(e,'gamma)),
  59     N : length(a),
  60
  61    /* Do gamma(X) gamma(n-X) ==> pochhammer(1-X,n-1) %pi/sin(%pi X), when
  62        n is a declared integer and X isn't a declared integer. */
  63     for ak in a do (
  64         for bk in a do (
  65         dn : ratsimp(ak+bk),
  66         if featurep(dn,'integer) and dn >= 0 and (not featurep(ak,'integer)) then (
  67             e : ratsubst(pochhammer(1-ak, dn-1) *%pi/(sin(%pi*ak)),
  68                 gamma(ak)*gamma(bk),e)))),
  69
  70     /* Do prod(gamma(x+k/n),k,0,n-1) ==> gamma(n*x)/((2*%pi)^((1-n)/2)*n^(n*x-1/2))
  71     */
  72     fn : make_gamma_product_id(),
  73     for ak in a do (
  74        for n : 2 thru N do (
  75           id : ratsimp(fn(ak,n)),
  76           e : ratsubst(rhs(id), lhs(id), e))),
  77
  78    /* Do gamma(2*x+n)/gamma(x) = pochhammer(2*x,n)*gamma(x+1/2)/(sqrt(%pi)*2^(1-2x))
  79       I think we only want to do this when n is an explict, not declared integer.
  80       Otherwise pochhammer(2*x,n)*gamma(x+1/2) is effectively a product of two
  81       gamma functions, and that's not really a simplification. So we'll test
  82       for an integer using integerp, not featurep.*/
  83     for ak in a do (
  84        for bk in a do (
  85          dn : ratsimp(ak-2*bk),
  86          if integerp(dn) then (
  87             e : ratsubst(pochhammer(2*bk,dn)*gamma(bk+1/2)/(sqrt(%pi)*2^(1-2*bk)),
  88                          gamma(ak)/gamma(bk),e)))),
  89
  90     /* When the expression has changed, recurse; otherwise return e.*/
  91     ee : factor(ee),
  92     if is(factor(e) # ee) then gamma_product_simp(homogenize_gamma(e)) else e)$
  93
  94 /* When n is a declared positive integer and the expression e involves both
  95    gamma(X+n) and gamma(X), do gamma(X+n) ==> pochhammer(X,n) gamma(X).
  96    The expansion of the pochhammer term is controlled by the option variable
  97    pochhammer_max_index (default 100). If the expression e changes after finishing
  98    the double loop, recurse.
  99
 100    Each time a match to gamma(X+n) & gamma(X) is detected, dispatch facsum on the
 101    result. Unfortunately, subsequent processing sometimes expands the factorization
 102    done by facsum.*/
 103
 104 homogenize_gamma(e) := block([a,ee : e, keepfloat : true],
 105    a : map('first, gatherargs(e,'gamma)),
 106    for ak in a do (
 107       for bk in a do (
 108          if featurep(ratsimp(ak-bk),'integer) and ak-bk > 0 then (
 109             e : ratsubst(pochhammer(bk,ak-bk)*gamma(bk), gamma(ak),e),
 110             e : facsum(e,gamma(bk))))),
 111       if is(e # ee) then homogenize_gamma(e) else e)$
 112
 113 /* See http://dlmf.nist.gov/5.4.E3 and http://dlmf.nist.gov/5.4.E4:
 114     gamma(1/2 + %i*y)*gamma(1/2-%i*y) = %pi/cosh(%pi*y)
 115     gamma(%iy)*gamma(-%i*y) = %pi/(y * (sinh(%pi * y))) */
 116
 117 gamma_magnitude_simp(e) :=  block([a, ee : e, y],
 118    a : map('first, gatherargs(e,'gamma)),
 119    for ak in a do (
 120       for bk in a do (
 121          if equal(ak+bk,1) and equal(realpart(ak),1/2) then (
 122             e : ratsubst((%pi/cosh(%pi*imagpart(ak))), gamma(ak)*gamma(bk),e)),
 123
 124          if equal(ak,-bk) and equal(realpart(ak),0) and equal(realpart(bk),0) then (
 125             y : imagpart(ak),
 126             e : ratsubst(%pi/(y * (sinh(%pi * y))), gamma(ak)*gamma(bk),e)))),
 127       if is(e # ee) then gamma_magnitude_simp(e) else e)$
 128
 129 /* Apply facsum to e with respect to all gamma terms. */
 130 facsum_gamma_terms(e) := block([a],
 131   a : map('first, gatherargs(e,'gamma)),
 132   a : map('gamma, a),
 133   if emptyp(a) then e else (
 134      buildq([z : a], lambda([x], facsum(x,splice(z))))(e)));
 135
 136 linearp(e,x) :=
 137   polynomialp(e,[x], lambda([q],freeof(x,q)), lambda([q], is(q<=1))) and
 138   lopow(e,x)=1;
 139
 140 /* Do X*gamma(X) ==> gamma(X+1). Repeat until there are no more such terms.
 141    See http://dlmf.nist.gov/5.5.E1
 142
 143    To match to X*gamma(X), we divide by the numerical coefficient of X.
 144    This allows matchings to, for example, x*gamma(2*x).  The linearity check is
 145    supposed to prevent substitutions such as (x+5)*gamma(x) ==> gamma(x+1)+5*gamma(x).
 146    This result isn't wrong, but it's not what we want, I think. */
 147
 148 gamma_arg_increase(e) := block([a, ee : e, es, g : gensym(),nf,inflag : false],
 149   if mapatom(e) then e
 150   elseif mtimesp(e) then (
 151      a : map('first, gatherargs(e,'gamma)),
 152      for ak in a do (
 153         nf : if numberp(ak) then 1 else numfactor(ak),
 154         es : ratsubst(g/nf, ratsimp(ak/nf)*gamma(ak),e),
 155         if linearp(es,g) then (
 156          e : ratsubst(gamma(ratsimp(ak+1)),g,es),
 157          e : factor(e))),
 158       if is(e # ee) then gamma_arg_increase(e) else e)
 159   else map('gamma_arg_increase, e));
 160
 161
 162 /* A measure of the complexity of an expression.*/
 163 gamma_simp_complexity(e) := block([],
 164  if mapatom(e) then 1 else xreduce("+", map('gamma_simp_complexity, args(e))));
 165
 166 /* Do gamma(X)/(X-1) ==> gamma(X-1). Repeat until there are no more
 167    such terms. */
 168
 169 gamma_arg_decrease(e) := block([a, ee : e,ea],
 170   if mapatom(e) then e
 171   elseif mtimesp(e) then (
 172      a : map('first, gatherargs(e,'gamma)),
 173      for ak in a do (
 174         ea : ratsubst((ak-1)*gamma(ak-1), gamma(ak),e),
 175         ea : factor(ea),
 176         if gamma_simp_complexity(ea) <= gamma_simp_complexity(e) then (
 177          e : ea)),
 178       if is(e # ee) then gamma_arg_decrease(e) else e)
 179   else map('gamma_arg_decrease, e));
 180
 181 /* Convert factorials to gamma functions, dispatch gamma_simp, and finally
 182    convert gamma functions back to factorial form. The gamma functions that
 183    are in the input are protected from participating the the simplification process.
 184    We do this by substituting a gensym for gamma. After converting back to factorials,
 185    this substitution is reverted.*/
 186 factorial_simp(e) := block([g : gensym(),opsubst : true, keepfloat : true],
 187   e : subst('gamma = lambda([x], g(x)),e),
 188   e : subst("!" = lambda([q], gamma(q+1)),e),
 189   e : gamma_simp(e),
 190   e : makefact(e),
 191   subst(g = lambda([x], gamma(x)),e));
 192