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