* ginac/registrar.h: dtor is inlined now.
[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.exadd(_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
114
115 //////////
116 // true Gamma function
117 //////////
118
119 static ex tgamma_evalf(const ex & x)
120 {
121         BEGIN_TYPECHECK
122                 TYPECHECK(x,numeric)
123         END_TYPECHECK(tgamma(x))
124         
125         return tgamma(ex_to_numeric(x));
126 }
127
128
129 /** Evaluation of tgamma(x), the true Gamma function.  Knows about integer
130  *  arguments, half-integer arguments and that's it. Somebody ought to provide
131  *  some good numerical evaluation some day...
132  *
133  *  @exception pole_error("tgamma_eval(): simple pole",0) */
134 static ex tgamma_eval(const ex & x)
135 {
136         if (x.info(info_flags::numeric)) {
137                 // trap integer arguments:
138                 if (x.info(info_flags::integer)) {
139                         // tgamma(n) -> (n-1)! for postitive n
140                         if (x.info(info_flags::posint)) {
141                                 return factorial(ex_to_numeric(x).sub(_num1()));
142                         } else {
143                                 throw (pole_error("tgamma_eval(): simple pole",1));
144                         }
145                 }
146                 // trap half integer arguments:
147                 if ((x*2).info(info_flags::integer)) {
148                         // trap positive x==(n+1/2)
149                         // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
150                         if ((x*_ex2()).info(info_flags::posint)) {
151                                 numeric n = ex_to_numeric(x).sub(_num1_2());
152                                 numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1()));
153                                 coefficient = coefficient.div(pow(_num2(),n));
154                                 return coefficient * pow(Pi,_ex1_2());
155                         } else {
156                                 // trap negative x==(-n+1/2)
157                                 // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
158                                 numeric n = abs(ex_to_numeric(x).sub(_num1_2()));
159                                 numeric coefficient = pow(_num_2(), n);
160                                 coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));;
161                                 return coefficient*power(Pi,_ex1_2());
162                         }
163                 }
164                 //  tgamma_evalf should be called here once it becomes available
165         }
166         
167         return tgamma(x).hold();
168 }
169
170
171 static ex tgamma_deriv(const ex & x, unsigned deriv_param)
172 {
173         GINAC_ASSERT(deriv_param==0);
174         
175         // d/dx  tgamma(x) -> psi(x)*tgamma(x)
176         return psi(x)*tgamma(x);
177 }
178
179
180 static ex tgamma_series(const ex & arg,
181                         const relational & rel,
182                         int order,
183                         unsigned options)
184 {
185         // method:
186         // Taylor series where there is no pole falls back to psi function
187         // evaluation.
188         // On a pole at -m use the recurrence relation
189         //   tgamma(x) == tgamma(x+1) / x
190         // from which follows
191         //   series(tgamma(x),x==-m,order) ==
192         //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1);
193         const ex arg_pt = arg.subs(rel);
194         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
195                 throw do_taylor();  // caught by function::series()
196         // if we got here we have to care for a simple pole at -m:
197         numeric m = -ex_to_numeric(arg_pt);
198         ex ser_denom = _ex1();
199         for (numeric p; p<=m; ++p)
200                 ser_denom *= arg+p;
201         return (tgamma(arg+m+_ex1())/ser_denom).series(rel, order+1, options);
202 }
203
204
205 REGISTER_FUNCTION(tgamma, eval_func(tgamma_eval).
206                           evalf_func(tgamma_evalf).
207                           derivative_func(tgamma_deriv).
208                           series_func(tgamma_series));
209
210
211 //////////
212 // beta-function
213 //////////
214
215 static ex beta_evalf(const ex & x, const ex & y)
216 {
217         BEGIN_TYPECHECK
218                 TYPECHECK(x,numeric)
219                 TYPECHECK(y,numeric)
220         END_TYPECHECK(beta(x,y))
221         
222         return tgamma(ex_to_numeric(x))*tgamma(ex_to_numeric(y))/tgamma(ex_to_numeric(x+y));
223 }
224
225
226 static ex beta_eval(const ex & x, const ex & y)
227 {
228         if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) {
229                 // treat all problematic x and y that may not be passed into tgamma,
230                 // because they would throw there although beta(x,y) is well-defined
231                 // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y)
232                 numeric nx(ex_to_numeric(x));
233                 numeric ny(ex_to_numeric(y));
234                 if (nx.is_real() && nx.is_integer() &&
235                         ny.is_real() && ny.is_integer()) {
236                         if (nx.is_negative()) {
237                                 if (nx<=-ny)
238                                         return pow(_num_1(), ny)*beta(1-x-y, y);
239                                 else
240                                         throw (pole_error("beta_eval(): simple pole",1));
241                         }
242                         if (ny.is_negative()) {
243                                 if (ny<=-nx)
244                                         return pow(_num_1(), nx)*beta(1-y-x, x);
245                                 else
246                                         throw (pole_error("beta_eval(): simple pole",1));
247                         }
248                         return tgamma(x)*tgamma(y)/tgamma(x+y);
249                 }
250                 // no problem in numerator, but denominator has pole:
251                 if ((nx+ny).is_real() &&
252                         (nx+ny).is_integer() &&
253                         !(nx+ny).is_positive())
254                          return _ex0();
255                 // everything is ok:
256                 return tgamma(x)*tgamma(y)/tgamma(x+y);
257         }
258         
259         return beta(x,y).hold();
260 }
261
262
263 static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param)
264 {
265         GINAC_ASSERT(deriv_param<2);
266         ex retval;
267         
268         // d/dx beta(x,y) -> (psi(x)-psi(x+y)) * beta(x,y)
269         if (deriv_param==0)
270                 retval = (psi(x)-psi(x+y))*beta(x,y);
271         // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y)
272         if (deriv_param==1)
273                 retval = (psi(y)-psi(x+y))*beta(x,y);
274         return retval;
275 }
276
277
278 static ex beta_series(const ex & arg1,
279                       const ex & arg2,
280                       const relational & rel,
281                       int order,
282                       unsigned options)
283 {
284         // method:
285         // Taylor series where there is no pole of one of the tgamma functions
286         // falls back to beta function evaluation.  Otherwise, fall back to
287         // tgamma series directly.
288         const ex arg1_pt = arg1.subs(rel);
289         const ex arg2_pt = arg2.subs(rel);
290         GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
291         const symbol *s = static_cast<symbol *>(rel.lhs().bp);
292         ex arg1_ser, arg2_ser, arg1arg2_ser;
293         if ((!arg1_pt.info(info_flags::integer) || arg1_pt.info(info_flags::positive)) &&
294             (!arg2_pt.info(info_flags::integer) || arg2_pt.info(info_flags::positive)))
295                 throw do_taylor();  // caught by function::series()
296         // trap the case where arg1 is on a pole:
297         if (arg1.info(info_flags::integer) && !arg1.info(info_flags::positive))
298                 arg1_ser = tgamma(arg1+*s).series(rel, order, options);
299         else
300                 arg1_ser = tgamma(arg1).series(rel,order);
301         // trap the case where arg2 is on a pole:
302         if (arg2.info(info_flags::integer) && !arg2.info(info_flags::positive))
303                 arg2_ser = tgamma(arg2+*s).series(rel, order, options);
304         else
305                 arg2_ser = tgamma(arg2).series(rel,order);
306         // trap the case where arg1+arg2 is on a pole:
307         if ((arg1+arg2).info(info_flags::integer) && !(arg1+arg2).info(info_flags::positive))
308                 arg1arg2_ser = tgamma(arg2+arg1+*s).series(rel, order, options);
309         else
310                 arg1arg2_ser = tgamma(arg2+arg1).series(rel,order);
311         // compose the result (expanding all the terms):
312         return (arg1_ser*arg2_ser/arg1arg2_ser).series(rel, order, options).expand();
313 }
314
315
316 REGISTER_FUNCTION(beta, eval_func(beta_eval).
317                         evalf_func(beta_evalf).
318                         derivative_func(beta_deriv).
319                         series_func(beta_series));
320
321
322 //////////
323 // Psi-function (aka digamma-function)
324 //////////
325
326 static ex psi1_evalf(const ex & x)
327 {
328         BEGIN_TYPECHECK
329                 TYPECHECK(x,numeric)
330         END_TYPECHECK(psi(x))
331         
332         return psi(ex_to_numeric(x));
333 }
334
335 /** Evaluation of digamma-function psi(x).
336  *  Somebody ought to provide some good numerical evaluation some day... */
337 static ex psi1_eval(const ex & x)
338 {
339         if (x.info(info_flags::numeric)) {
340                 numeric nx = ex_to_numeric(x);
341                 if (nx.is_integer()) {
342                         // integer case 
343                         if (nx.is_positive()) {
344                                 // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - Euler
345                                 numeric rat(0);
346                                 for (numeric i(nx+_num_1()); i.is_positive(); --i)
347                                         rat += i.inverse();
348                                 return rat-Euler;
349                         } else {
350                                 // for non-positive integers there is a pole:
351                                 throw (pole_error("psi_eval(): simple pole",1));
352                         }
353                 }
354                 if ((_num2()*nx).is_integer()) {
355                         // half integer case
356                         if (nx.is_positive()) {
357                                 // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - Euler - 2log(2)
358                                 numeric rat(0);
359                                 for (numeric i((nx+_num_1())*_num2()); i.is_positive(); i-=_num2())
360                                                                           rat += _num2()*i.inverse();
361                                                                           return rat-Euler-_ex2()*log(_ex2());
362                         } else {
363                                 // use the recurrence relation
364                                 //   psi(-m-1/2) == psi(-m-1/2+1) - 1 / (-m-1/2)
365                                 // to relate psi(-m-1/2) to psi(1/2):
366                                 //   psi(-m-1/2) == psi(1/2) + r
367                                 // where r == ((-1/2)^(-1) + ... + (-m-1/2)^(-1))
368                                 numeric recur(0);
369                                 for (numeric p(nx); p<0; ++p)
370                                         recur -= pow(p, _num_1());
371                                 return recur+psi(_ex1_2());
372                         }
373                 }
374                 //  psi1_evalf should be called here once it becomes available
375         }
376         
377         return psi(x).hold();
378 }
379
380 static ex psi1_deriv(const ex & x, unsigned deriv_param)
381 {
382         GINAC_ASSERT(deriv_param==0);
383         
384         // d/dx psi(x) -> psi(1,x)
385         return psi(_ex1(), x);
386 }
387
388 static ex psi1_series(const ex & arg,
389                       const relational & rel,
390                       int order,
391                       unsigned options)
392 {
393         // method:
394         // Taylor series where there is no pole falls back to polygamma function
395         // evaluation.
396         // On a pole at -m use the recurrence relation
397         //   psi(x) == psi(x+1) - 1/z
398         // from which follows
399         //   series(psi(x),x==-m,order) ==
400         //   series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x==-m,order);
401         const ex arg_pt = arg.subs(rel);
402         if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
403                 throw do_taylor();  // caught by function::series()