c4f4f1a9aec248c2b1c1aaba4549925188e829de
[ginac.git] / ginac / inifcns_gamma.cpp
1 /** @file inifcns_gamma.cpp
2  *
3  *  Implementation of Gamma-function, Beta-function, Polygamma-functions, and
4  *  some related stuff. */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <vector>
25 #include <stdexcept>
26
27 #include "inifcns.h"
28 #include "constant.h"
29 #include "pseries.h"
30 #include "numeric.h"
31 #include "power.h"
32 #include "relational.h"
33 #include "symbol.h"
34 #include "utils.h"
35
36 namespace GiNaC {
37
38 //////////
39 // Logarithm of Gamma function
40 //////////
41
42 static ex lgamma_evalf(const ex & x)
43 {
44         BEGIN_TYPECHECK
45                 TYPECHECK(x,numeric)
46         END_TYPECHECK(lgamma(x))
47         
48         return lgamma(ex_to_numeric(x));
49 }
50
51
52 /** Evaluation of lgamma(x), the natural logarithm of the Gamma function.
53  *  Knows about integer arguments and that's it.  Somebody ought to provide
54  *  some good numerical evaluation some day...
55  *
56  *  @exception GiNaC::pole_error("lgamma_eval(): logarithmic pole",0) */
57 static ex lgamma_eval(const ex & x)
58 {
59         if (x.info(info_flags::numeric)) {
60                 // trap integer arguments:
61                 if (x.info(info_flags::integer)) {
62                         // lgamma(n) -> log((n-1)!) for postitive n
63                         if (x.info(info_flags::posint))
64                                 return log(factorial(x + _ex_1()));
65                         else
66                                 throw (pole_error("lgamma_eval(): logarithmic pole",0));
67                 }
68                 //  lgamma_evalf should be called here once it becomes available
69         }
70         
71         return lgamma(x).hold();
72 }
73
74
75 static ex lgamma_deriv(const ex & x, unsigned deriv_param)
76 {
77         GINAC_ASSERT(deriv_param==0);
78         
79         // d/dx  lgamma(x) -> psi(x)
80         return psi(x);
81 }
82
83
84 static ex lgamma_series(const ex & arg,
85                         const relational & rel,
86                         int order,
87                         unsigned options)
88 {
89         // method:
90         // Taylor series where there is no pole falls back to psi function
91         // evaluation.
92         // On a pole at -m we could use the recurrence relation
93         //   lgamma(x) == lgamma(x+1)-log(x)
94         // from which follows
95         //   series(lgamma(x),x==-m,order) ==
96         //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order);
97         const ex arg_pt = arg.subs(rel);
98         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
99                 throw do_taylor();  // caught by function::series()
100         // if we got here we have to care for a simple pole of tgamma(-m):
101         numeric m = -ex_to_numeric(arg_pt);
102         ex recur;
103         for (numeric p; p<=m; ++p)
104                 recur += log(arg+p);
105         return (lgamma(arg+m+_ex1())-recur).series(rel, order, options);
106 }
107
108
109 REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval).
110                           evalf_func(lgamma_evalf).
111                           derivative_func(lgamma_deriv).
112                           series_func(lgamma_series).
113                           latex_name("\\log \\Gamma"));
114
115
116 //////////
117 // true Gamma function
118 //////////
119
120 static ex tgamma_evalf(const ex & x)
121 {
122         BEGIN_TYPECHECK
123                 TYPECHECK(x,numeric)
124         END_TYPECHECK(tgamma(x))
125         
126         return tgamma(ex_to_numeric(x));
127 }
128
129
130 /** Evaluation of tgamma(x), the true Gamma function.  Knows about integer
131  *  arguments, half-integer arguments and that's it. Somebody ought to provide
132  *  some good numerical evaluation some day...
133  *
134  *  @exception pole_error("tgamma_eval(): simple pole",0) */
135 static ex tgamma_eval(const ex & x)
136 {
137         if (x.info(info_flags::numeric)) {
138                 // trap integer arguments:
139                 if (x.info(info_flags::integer)) {
140                         // tgamma(n) -> (n-1)! for postitive n
141                         if (x.info(info_flags::posint)) {
142                                 return factorial(ex_to_numeric(x).sub(_num1()));
143                         } else {
144                                 throw (pole_error("tgamma_eval(): simple pole",1));
145                         }
146                 }
147                 // trap half integer arguments:
148                 if ((x*2).info(info_flags::integer)) {
149                         // trap positive x==(n+1/2)
150                         // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
151                         if ((x*_ex2()).info(info_flags::posint)) {
152                                 numeric n = ex_to_numeric(x).sub(_num1_2());
153                                 numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1()));
154                                 coefficient = coefficient.div(pow(_num2(),n));
155                                 return coefficient * pow(Pi,_ex1_2());
156                         } else {
157                                 // trap negative x==(-n+1/2)
158                                 // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
159                                 numeric n = abs(ex_to_numeric(x).sub(_num1_2()));
160                                 numeric coefficient = pow(_num_2(), n);
161                                 coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));;
162                                 return coefficient*power(Pi,_ex1_2());
163                         }
164                 }
165                 //  tgamma_evalf should be called here once it becomes available
166         }
167         
168         return tgamma(x).hold();
169 }
170
171
172 static ex tgamma_deriv(const ex & x, unsigned deriv_param)
173 {
174         GINAC_ASSERT(deriv_param==0);
175         
176         // d/dx  tgamma(x) -> psi(x)*tgamma(x)
177         return psi(x)*tgamma(x);
178 }
179
180
181 static ex tgamma_series(const ex & arg,
182                         const relational & rel,
183                         int order,
184                         unsigned options)
185 {
186         // method:
187         // Taylor series where there is no pole falls back to psi function
188         // evaluation.
189         // On a pole at -m use the recurrence relation
190         //   tgamma(x) == tgamma(x+1) / x
191         // from which follows
192         //   series(tgamma(x),x==-m,order) ==
193         //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1);
194         const ex arg_pt = arg.subs(rel);
195         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
196                 throw do_taylor();  // caught by function::series()
197         // if we got here we have to care for a simple pole at -m:
198         numeric m = -ex_to_numeric(arg_pt);
199         ex ser_denom = _ex1();
200         for (numeric p; p<=m; ++p)
201                 ser_denom *= arg+p;
202         return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1, options);
203 }
204
205
206 REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval).
207                           evalf_func(tgamma_evalf).
208                           derivative_func(tgamma_deriv).
209                           series_func(tgamma_series).
210                           latex_name("\\Gamma"));
211
212
213 //////////
214 // beta-function
215 //////////
216
217 static ex beta_evalf(const ex & x, const ex & y)
218 {
219         BEGIN_TYPECHECK
220                 TYPECHECK(x,numeric)
221                 TYPECHECK(y,numeric)
222         END_TYPECHECK(beta(x,y))
223         
224         return tgamma(ex_to_numeric(x))*tgamma(ex_to_numeric(y))/tgamma(ex_to_numeric(x+y));
225 }
226
227
228 static ex beta_eval(const ex & x, const ex & y)
229 {
230         if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) {
231                 // treat all problematic x and y that may not be passed into tgamma,
232                 // because they would throw there although beta(x,y) is well-defined
233                 // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y)
234                 numeric nx(ex_to_numeric(x));
235                 numeric ny(ex_to_numeric(y));
236                 if (nx.is_real() && nx.is_integer() &&
237                         ny.is_real() && ny.is_integer()) {
238                         if (nx.is_negative()) {
239                                 if (nx<=-ny)
240                                         return pow(_num_1(), ny)*beta(1-x-y, y);
241                                 else
242                                         throw (pole_error("beta_eval(): simple pole",1));
243                         }
244                         if (ny.is_negative()) {
245                                 if (ny<=-nx)
246                                         return pow(_num_1(), nx)*beta(1-y-x, x);
247                                 else
248                                         throw (pole_error("beta_eval(): simple pole",1));
249                         }
250                         return tgamma(x)*tgamma(y)/tgamma(x+y);
251                 }
252                 // no problem in numerator, but denominator has pole:
253                 if ((nx+ny).is_real() &&
254                         (nx+ny).is_integer() &&
255                         !(nx+ny).is_positive())
256                          return _ex0();
257                 // everything is ok:
258                 return tgamma(x)*tgamma(y)/tgamma(x+y);
259         }
260         
261         return beta(x,y).hold();
262 }
263
264
265 static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param)
266 {
267         GINAC_ASSERT(deriv_param<2);
268         ex retval;
269         
270         // d/dx beta(x,y) -> (psi(x)-psi(x+y)) * beta(x,y)
271         if (deriv_param==0)
272                 retval = (psi(x)-psi(x+y))*beta(x,y);
273         // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y)
274         if (deriv_param==1)
275                 retval = (psi(y)-psi(x+y))*beta(x,y);
276         return retval;
277 }
278
279
280 static ex beta_series(const ex & arg1,
281                       const ex & arg2,
282                       const relational & rel,
283                       int order,
284                       unsigned options)
285 {
286         // method:
287         // Taylor series where there is no pole of one of the tgamma functions
288         // falls back to beta function evaluation.  Otherwise, fall back to
289         // tgamma series directly.
290         const ex arg1_pt = arg1.subs(rel);
291         const ex arg2_pt = arg2.subs(rel);
292         GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
293         const symbol *s = static_cast<symbol *>(rel.lhs().bp);
294         ex arg1_ser, arg2_ser, arg1arg2_ser;
295         if ((!arg1_pt.info(info_flags::integer) || arg1_pt.info(info_flags::positive)) &&
296             (!arg2_pt.info(info_flags::integer) || arg2_pt.info(info_flags::positive)))
297                 throw do_taylor();  // caught by function::series()
298         // trap the case where arg1 is on a pole:
299         if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))
300                 arg1_ser = tgamma(arg1+*s).series(rel, order, options);
301         else
302                 arg1_ser = tgamma(arg1).series(rel,order);
303         // trap the case where arg2 is on a pole:
304         if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))
305                 arg2_ser = tgamma(arg2+*s).series(rel, order, options);
306         else
307                 arg2_ser = tgamma(arg2).series(rel,order);
308         // trap the case where arg1+arg2 is on a pole:
309         if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))
310                 arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options);
311         else
312                 arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
313         // compose the result (expanding all the terms):
314         return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel, order, options).expand();
315 }
316
317
318 REGISTER_FUNCTION(beta, eval_func(beta_eval).
319                         evalf_func(beta_evalf).
320                         derivative_func(beta_deriv).
321                         series_func(beta_series).
322                         latex_name("\\mbox{B}"));
323
324
325 //////////
326 // Psi-function (aka digamma-function)
327 //////////
328
329 static ex psi1_evalf(const ex & x)
330 {
331         BEGIN_TYPECHECK
332                 TYPECHECK(x,numeric)
333         END_TYPECHECK(psi(x))
334         
335         return psi(ex_to_numeric(x));
336 }
337
338 /** Evaluation of digamma-function psi(x).
339  *  Somebody ought to provide some good numerical evaluation some day... */
340 static ex psi1_eval(const ex & x)
341 {
342         if (x.info(info_flags::numeric)) {
343                 numeric nx = ex_to_numeric(x);
344                 if (nx.is_integer()) {
345                         // integer case 
346                         if (nx.is_positive()) {
347                                 // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - Euler
348                                 numeric rat(0);
349                                 for (numeric i(nx+_num_1()); i.is_positive(); --i)
350                                         rat += i.inverse();
351                                 return rat-Euler;
352                         } else {
353                                 // for non-positive integers there is a pole:
354                                 throw (pole_error("psi_eval(): simple pole",1));
355                         }
356                 }
357                 if ((_num2()*nx).is_integer()) {
358                         // half integer case
359                         if (nx.is_positive()) {
360                                 // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - Euler - 2log(2)
361                                 numeric rat(0);
362                                 for (numeric i((nx+_num_1())*_num2()); i.is_positive(); i-=_num2())
363                                                                           rat += _num2()*i.inverse();
364                                                                           return rat-Euler-_ex2()*log(_ex2());
365                         } else {
366                                 // use the recurrence relation
367                                 //   psi(-m-1/2) == psi(-m-1/2+1) - 1 / (-m-1/2)
368                                 // to relate psi(-m-1/2) to psi(1/2):
369                                 //   psi(-m-1/2) == psi(1/2) + r
370                                 // where r == ((-1/2)^(-1) + ... + (-m-1/2)^(-1))
371                                 numeric recur(0);
372                                 for (numeric p(nx); p<0; ++p)
373                                         recur -= pow(p, _num_1());
374                                 return recur+psi(_ex1_2());
375                         }
376                 }
377                 //  psi1_evalf should be called here once it becomes available
378         }
379         
380         return psi(x).hold();
381 }
382
383 static ex psi1_deriv(const ex & x, unsigned deriv_param)
384 {
385         GINAC_ASSERT(deriv_param==0);
386         
387         // d/dx psi(x) -> psi(1,x)
388         return psi(_ex1(), x);
389 }
390
391 static ex psi1_series(const ex & arg,
392                       const relational & rel,
393                       int order,
394                       unsigned options)
395 {
396         // method:
397         // Taylor series where there is no pole falls back to polygamma function
398         // evaluation.
399         // On a pole at -m use the recurrence relation
400         //   psi(x) == psi(x+1) - 1/z
401         // from which follows
402         //   series(psi(x),x==-m,order) ==
403         //   series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x==-m,order);
404         const ex arg_pt = arg.subs(rel);
405         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
406                 throw do_taylor();  // caught by function::series()
407         // if we got here we have to care for a simple pole at -m:
408         numeric m = -ex_to_numeric(arg_pt);
409         ex recur;
410         for (numeric p; p<=m; ++p)
411                 recur += power(arg+p,_ex_1());
412         return (psi(arg+m+_ex1())-recur).series(rel, order, options);
413 }
414
415 const unsigned function_index_psi1 =
416         function::register_new(function_options("psi").
417                                eval_func(psi1_eval).
418                                evalf_func(psi1_evalf).
419                                derivative_func(psi1_deriv).
420                                series_func(psi1_series).
421                                overloaded(2));
422
423 //////////
424 // Psi-functions (aka polygamma-functions)  psi(0,x)==psi(x)
425 //////////
426
427 static ex psi2_evalf(const ex & n, const ex & x)
428 {
429         BEGIN_TYPECHECK
430                 TYPECHECK(n,numeric)
431                 TYPECHECK(x,numeric)
432         END_TYPECHECK(psi(n,x))
433         
434         return psi(ex_to_numeric(n), ex_to_numeric(x));
435 }
436
437 /** Evaluation of polygamma-function psi(n,x). 
438  *  Somebody ought to provide some good numerical evaluation some day... */
439 static ex psi2_eval(const ex & n, const ex & x)
440 {
441         // psi(0,x) -> psi(x)
442         if (n.is_zero())
443                 return psi(x);
444         // psi(-1,x) -> log(tgamma(x))
445         if (n.is_equal(_ex_1()))
446                 return log(tgamma(x));
447         if (n.info(info_flags::numeric) && n.info(info_flags::posint) &&
448                 x.info(info_flags::numeric)) {
449                 numeric nn = ex_to_numeric(n);
450                 numeric nx = ex_to_numeric(x);
451                 if (nx.is_integer()) {
452                         // integer case 
453                         if (nx.is_equal(_num1()))
454                                 // use psi(n,1) == (-)^(n+1) * n! * zeta(n+1)
455                                 return pow(_num_1(),nn+_num1())*factorial(nn)*zeta(ex(nn+_num1()));
456                         if (nx.is_positive()) {
457                                 // use the recurrence relation
458                                 //   psi(n,m) == psi(n,m+1) - (-)^n * n! / m^(n+1)
459                                 // to relate psi(n,m) to psi(n,1):
460                                 //   psi(n,m) == psi(n,1) + r
461                                 // where r == (-)^n * n! * (1^(-n-1) + ... + (m-1)^(-n-1))
462                                 numeric recur(0);
463                                 for (numeric p(1); p<nx; ++p)
464                                         recur += pow(p, -nn+_num_1());
465                                 recur *= factorial(nn)*pow(_num_1(), nn);
466                                 return recur+psi(n,_ex1());
467                         } else {
468                                 // for non-positive integers there is a pole:
469                                 throw (pole_error("psi2_eval(): pole",1));
470                         }
471                 }
472                 if ((_num2()*nx).is_integer()) {
473                         // half integer case
474                         if (nx.is_equal(_num1_2()))
475                                 // use psi(n,1/2) == (-)^(n+1) * n! * (2^(n+1)-1) * zeta(n+1)
476                                 return pow(_num_1(),nn+_num1())*factorial(nn)*(pow(_num2(),nn+_num1()) + _num_1())*zeta(ex(nn+_num1()));
477                         if (nx.is_positive()) {
478                                 numeric m = nx - _num1_2();
479                                 // use the multiplication formula
480                                 //   psi(n,2*m) == (psi(n,m) + psi(n,m+1/2)) / 2^(n+1)
481                                 // to revert to positive integer case
482                                 return psi(n,_num2()*m)*pow(_num2(),nn+_num1())-psi(n,m);
483                         } else {
484                                 // use the recurrence relation
485                                 //   psi(n,-m-1/2) == psi(n,-m-1/2+1) - (-)^n * n! / (-m-1/2)^(n+1)
486                                 // to relate psi(n,-m-1/2) to psi(n,1/2):
487                                 //   psi(n,-m-1/2) == psi(n,1/2) + r
488                                 // where r == (-)^(n+1) * n! * ((-1/2)^(-n-1) + ... + (-m-1/2)^(-n-1))
489                                 numeric recur(0);
490                                 for (numeric p(nx); p<0; ++p)
491                                         recur += pow(p, -nn+_num_1());
492                                 recur *= factorial(nn)*pow(_num_1(), nn+_num_1());
493                                 return recur+psi(n,_ex1_2());
494                         }
495                 }
496                 //  psi2_evalf should be called here once it becomes available
497         }
498         
499         return psi(n, x).hold();
500 }    
501
502 static ex psi2_deriv(const ex & n, const ex & x, unsigned deriv_param)
503 {
504         GINAC_ASSERT(deriv_param<2);
505         
506         if (deriv_param==0) {
507                 // d/dn psi(n,x)
508                 throw(std::logic_error("cannot diff psi(n,x) with respect to n"));
509         }
510         // d/dx psi(n,x) -> psi(n+1,x)
511         return psi(n+_ex1(), x);
512 }
513
514 static ex psi2_series(const ex & n,
515                       const ex & arg,
516                       const relational & rel,
517                       int order,
518                       unsigned options)
519 {
520         // method:
521         // Taylor series where there is no pole falls back to polygamma function
522         // evaluation.
523         // On a pole at -m use the recurrence relation
524         //   psi(n,x) == psi(n,x+1) - (-)^n * n! / x^(n+1)
525         // from which follows
526         //   series(psi(x),x==-m,order) == 
527         //   series(psi(x+m+1) - (-1)^n * n! * ((x)^(-n-1) + (x+1)^(-n-1) + ...
528         //                                      ... + (x+m)^(-n-1))),x==-m,order);
529         const ex arg_pt = arg.subs(rel);
530         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
531                 throw do_taylor();  // caught by function::series()
532         // if we got here we have to care for a pole of order n+1 at -m:
533         numeric m = -ex_to_numeric(arg_pt);
534         ex recur;
535         for (numeric p; p<=m; ++p)
536                 recur += power(arg+p,-n+_ex_1());
537         recur *= factorial(n)*power(_ex_1(),n);
538         return (psi(n, arg+m+_ex1())-recur).series(rel, order, options);
539 }
540
541 const unsigned function_index_psi2 =
542         function::register_new(function_options("psi").
543                                eval_func(psi2_eval).
544                                evalf_func(psi2_evalf).
545                                derivative_func(psi2_deriv).
546                                series_func(psi2_series).
547                            latex_name("\\psi").
548                                overloaded(2));
549
550
551 } // namespace GiNaC