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