share/gamma_simp/psi_simp.mac

   1
   2 /* Author: Barton Willis, copyright 2022
   3
   4 Maxima code for simplifying expressions that involve sums of psi functions.
   5
   6 See LICENSE.md for licensing information.*/
   7
   8 /* We need the function gatherargs--it's in "opsubst." */
   9 load("opsubst");
  10 load("gamma_simp");
  11 load("isequal.lisp");
  12
  13 /* This function is inflag sensitive. */
  14 mtimesp(e) := not mapatom(e) and part(e,0) = "*";
  15 mplusp(e) := not mapatom(e) and part(e,0) = "+";
  16
  17 mysubst(new,old,e) := block([inflag : true, q, prederror : false,listconstvars : false],
  18
  19   if mapatom(e) then (
  20         print(is(equal(e,old))),
  21         if equal(e,old)=true then new else e)
  22
  23   elseif isequal(e,old) then new
  24
  25   elseif mtimesp(e) then (
  26         q : apply('divide, append([e, old], listofvars(old))),
  27         if equal(second(q),0) then (
  28                 new*mysubst(new,old,first(q)))
  29          else e)
  30
  31   elseif mplusp(e) then (
  32         q : apply('divide, append([e, old], listofvars(old))),
  33         if lfreeof(listofvars(old),second(q)) then (
  34                 first(q)*new + second(q))
  35         else e)
  36   else map(lambda([s], mysubst(new, old, s)), e));
  37
  38 make_psi_sum_id(g) := buildq([g,x : gensym(), k : gensym(),n : gensym()],
  39     lambda([x,n], -g(x) + sum(1/(x+k),k,0,n-1) = g(x+n)));
  40
  41 psi_simp(e) := block([a,g : gensym(),n,fn,id],
  42           e : opsubst(g, psi[0],e),
  43           a : map('first, gatherargs(e,g)),
  44           for ak in a do (
  45             for bk in a do (
  46                 n : ratsimp(ak-bk),
  47                         /* see http://dlmf.nist.gov/5.5.E2 */
  48                 if featurep(n,'integer) and n > 0 then (
  49                 fn : make_psi_sum_id(g),
  50                 id : fn(bk,n),
  51                     e : mysubst(lhs(id),rhs(id), e)),
  52
  53                         /* see http://dlmf.nist.gov/5.5.E8 */
  54                         if featurep(2*n,'integer) and n > 0 then (
  55                                 e : mysubst(2*log(2) + 2*g(ak), g(ak)+g(bk),e)),
  56                         n : ratsimp(ak+bk),
  57
  58                         /* see http://dlmf.nist.gov/5.5.E4 */
  59                         if featurep(n,'integer) and n = 1 and gamma_simp_complexity(ak) < gamma_simp_complexity(bk) then (
  60                                 e : mysubst(g(ak)+%pi/tan(%pi*ak), g(bk),e)))),
  61           opsubst(psi[0],g, e));