gamma.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Author : philippe.billet@noos.fr
   4 //
   5 //      Gamma function gamma(x)
   6 //
   7 //-----------------------------------------------------------------------------
   8
   9 #include "stdafx.h"
  10 #include "defs.h"
  11 void gamma(void);
  12 static void gammaf(void);
  13 static void gamma_of_sum(void);
  14
  15 void
  16 eval_gamma(void)
  17 {
  18         push(cadr(p1));
  19         eval();
  20         gamma();
  21 }
  22
  23 void
  24 gamma(void)
  25 {
  26         save();
  27         gammaf();
  28         restore();
  29 }
  30
  31 static void
  32 gammaf(void)
  33 {
  34 //      double d;
  35
  36         p1 = pop();
  37
  38         if (isrational(p1) && MEQUAL(p1->u.q.a, 1) && MEQUAL(p1->u.q.b, 2))  {
  39                 push_symbol(PI);;
  40                 push_rational(1,2);
  41                 power();
  42                 return;
  43         }
  44
  45         if (isrational(p1) && MEQUAL(p1->u.q.a, 3) && MEQUAL(p1->u.q.b, 2))  {
  46                 push_symbol(PI);;
  47                 push_rational(1,2);
  48                 power();
  49                 push_rational(1,2);
  50                 multiply();
  51                 return;
  52         }
  53
  54 //      if (p1->k == DOUBLE) {
  55 //              d = exp(lgamma(p1->u.d));
  56 //              push_double(d);
  57 //              return;
  58 //      }
  59
  60         if (isnegativeterm(p1)) {
  61                 push_symbol(PI);
  62                 push_integer(-1);
  63                 multiply();
  64                 push_symbol(PI);
  65                 push(p1);
  66                 multiply();
  67                 sine();
  68                 push(p1);
  69                 multiply();
  70                 push(p1);
  71                 negate();
  72                 gamma();
  73                 multiply();
  74                 divide();
  75                 return;
  76         }
  77
  78         if (car(p1) == symbol(ADD)) {
  79                 gamma_of_sum();
  80                 return;
  81         }
  82
  83
  84         push_symbol(GAMMA);
  85         push(p1);
  86         list(2);
  87         return;
  88 }
  89
  90 static void
  91 gamma_of_sum(void)
  92 {
  93         p3 = cdr(p1);
  94         if (isrational(car(p3)) && MEQUAL(car(p3)->u.q.a, 1) && MEQUAL(car(p3)->u.q.b, 1)) {
  95                 push(cadr(p3));
  96                 push(cadr(p3));
  97                 gamma();
  98                 multiply();
  99         }
 100         else {
 101                 if (isrational(car(p3)) && MEQUAL(car(p3)->u.q.a, -1) && MEQUAL(car(p3)->u.q.b, 1)) {
 102                         push(cadr(p3));
 103                         gamma();
 104                         push(cadr(p3));
 105                         push_integer(-1);
 106                         add();
 107                         divide();
 108                 }
 109                 else    {
 110                         push_symbol(GAMMA);
 111                         push(p1);
 112                         list(2);
 113                         return;
 114                 }
 115         }
 116 }
 117
 118 #if SELFTEST
 119
 120 static char *s[] = {
 121
 122         "Gamma(a)",
 123         "Gamma(a)",
 124
 125 //      "float(gamma(10))",
 126 //      "362880",
 127
 128         "Gamma(x+1)",
 129         "x*Gamma(x)",
 130
 131         "Gamma(1/2)",
 132         "pi^(1/2)",
 133
 134         "Gamma(x-1)-Gamma(x)/(-1+x)",
 135         "0",
 136
 137         "Gamma(-x)",
 138         "-pi/(x*Gamma(x)*sin(pi*x))",
 139
 140 };
 141
 142 void
 143 test_gamma(void)
 144 {
 145         test(__FILE__, s, sizeof s / sizeof (char *));
 146 }
 147
 148 #endif