123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149 |
- //-----------------------------------------------------------------------------
- //
- // Author : philippe.billet@noos.fr
- //
- // Gamma function gamma(x)
- //
- //-----------------------------------------------------------------------------
- #include "stdafx.h"
- #include "defs.h"
- void gamma(void);
- static void gammaf(void);
- static void gamma_of_sum(void);
- void
- eval_gamma(void)
- {
- push(cadr(p1));
- eval();
- gamma();
- }
- void
- gamma(void)
- {
- save();
- gammaf();
- restore();
- }
- static void
- gammaf(void)
- {
- // double d;
- p1 = pop();
- if (isrational(p1) && MEQUAL(p1->u.q.a, 1) && MEQUAL(p1->u.q.b, 2)) {
- push_symbol(PI);;
- push_rational(1,2);
- power();
- return;
- }
- if (isrational(p1) && MEQUAL(p1->u.q.a, 3) && MEQUAL(p1->u.q.b, 2)) {
- push_symbol(PI);;
- push_rational(1,2);
- power();
- push_rational(1,2);
- multiply();
- return;
- }
-
- // if (p1->k == DOUBLE) {
- // d = exp(lgamma(p1->u.d));
- // push_double(d);
- // return;
- // }
- if (isnegativeterm(p1)) {
- push_symbol(PI);
- push_integer(-1);
- multiply();
- push_symbol(PI);
- push(p1);
- multiply();
- sine();
- push(p1);
- multiply();
- push(p1);
- negate();
- gamma();
- multiply();
- divide();
- return;
- }
-
- if (car(p1) == symbol(ADD)) {
- gamma_of_sum();
- return;
- }
-
-
- push_symbol(GAMMA);
- push(p1);
- list(2);
- return;
- }
- static void
- gamma_of_sum(void)
- {
- p3 = cdr(p1);
- if (isrational(car(p3)) && MEQUAL(car(p3)->u.q.a, 1) && MEQUAL(car(p3)->u.q.b, 1)) {
- push(cadr(p3));
- push(cadr(p3));
- gamma();
- multiply();
- }
- else {
- if (isrational(car(p3)) && MEQUAL(car(p3)->u.q.a, -1) && MEQUAL(car(p3)->u.q.b, 1)) {
- push(cadr(p3));
- gamma();
- push(cadr(p3));
- push_integer(-1);
- add();
- divide();
- }
- else {
- push_symbol(GAMMA);
- push(p1);
- list(2);
- return;
- }
- }
- }
- #if SELFTEST
- static char *s[] = {
- "Gamma(a)",
- "Gamma(a)",
- // "float(gamma(10))",
- // "362880",
- "Gamma(x+1)",
- "x*Gamma(x)",
-
- "Gamma(1/2)",
- "pi^(1/2)",
-
- "Gamma(x-1)-Gamma(x)/(-1+x)",
- "0",
- "Gamma(-x)",
- "-pi/(x*Gamma(x)*sin(pi*x))",
-
- };
- void
- test_gamma(void)
- {
- test(__FILE__, s, sizeof s / sizeof (char *));
- }
- #endif
|