1 /*
2  * This program can be used to find the coefficients needed to approximate
3  * the gamma function using the Lanczos approximation.
4  *
5  * Usage: lanczos -n order -D digits
6  *
7  * The order defaults to 10. digits defaults to GiNaCs default.
8  * It is recommended to run the program several times with an increasing
9  * value for digits until numerically stablilty for the values has been
10  * reached. The program will also print the number of digits for which
11  * the approximation is still reliable. This will be lower then the
12  * number of digits given at command line. It is determined by comparing
13  * Gamma(1/2) to sqrt(Pi). Note that the program may crash if the number of
14  * digits is unreasonably small for the given order. Another thing that
15  * can happen if the number of digits is too small is that the program will
16  * print "Forget it, this is waaaaaaay too inaccurate." at the top of the
17  * output.
18  *
19  * The gamma function can be (for real_part(z) > -1/2) calculated using
20  *
21  *    Gamma(z+1) = sqrt(2*Pi)*power(z+g+ex(1)/2, z+ex(1)/2)*exp(-(z+g+ex(1)/2))
22  *                  *A_g(z),
23  * where,
24  *
25  *    A_g(z) = coeff + coeff/(z+1) + coeff/(z+2) + ...
26  *                + coeff[N-1]/(z+N-1).
27  *
28  * The value of g is taken to be equal to the order N.
29  *
30  * More details can be found at Wikipedia:
31  * http://en.wikipedia.org/wiki/Lanczos_approximation.
32  *
33  * (C) 2006 Chris Dams
34  *
35  * This program is free software; you can redistribute it and/or modify
36  * it under the terms of the GNU General Public License as published by
37  * the Free Software Foundation; either version 2 of the License, or
38  * (at your option) any later version.
39  *
40  * This program is distributed in the hope that it will be useful,
41  * but WITHOUT ANY WARRANTY; without even the implied warranty of
42  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
43  * GNU General Public License for more details.
44  *
45  * You should have received a copy of the GNU General Public License
46  * along with this program; if not, write to the Free Software
47  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
48  * MA  02110-1301 USA
49  */
51 #include <ginac/ginac.h>
52 #include <iostream>
54 using namespace std;
55 using namespace GiNaC;
57 /*
58  * Double factorial function.
59  */
60 ex doublefact(int i)
61 {       if (i==0 || i==-1)
62                 return 1;
63         if (i>0)
64                 return i*doublefact(i-2);
65         cout << "function doublefact called with wrong argument" << endl;
66         exit(1);
67 }
69 /*
70  * Chebyshev polynomial coefficient matrix as far as is required for
71  * the Lanczos approximation.
72  */
73 void initialize_C_vec(vector<exvector> &C, int size)
74 {       C.reserve(size);
75         for (int i=0; i<size; ++i)
76                 C.push_back(exvector(size, 0));
77         C = 1;
78         C = 1;
79         for (int i=3; i<size; i+=2)
80                 C[i] = -C[i-2];
81         for (int i=3; i<size; ++i)
82                 C[i][i] = 2*C[i-1][i-1];
83         for (int j=2; j<size; ++j)
84                 for (int i=j+2; i<size; i+=2)
85                         C[i][j] = 2*C[i-1][j-1] - C[i-2][j];
86 }
88 /*
89  * The coefficients p_n(g) that occur in the Lanczos approximation.
90  */
91 ex p(int k, ex g, const vector<exvector> &C)
92 {       ex result = 0;
93         for (int a=0; a<=k; ++a)
94                 result += 2*C[2*k+1][2*a+1]/sqrt(Pi)*doublefact(2*a-1)
95                         *power(2*a+2*g+1, -a-ex(1)/2)*exp(a+g+ex(1)/2);
96         return result;
97 }
99 /*
100  * Check wheter x is of the form 1/(z+n) where z is given by the second
101  * argument and n should be a positive integer that is returned as the
102  * third argument. If x is not of this form, false is returned.
103  */
104 bool is_z_pole(const ex &x, const ex &z, ex &n)
105 {       if (!is_a<power>(x))
106                 return false;
107         if (x.op(1) != -1)
108                 return false;
109         ex denom = x.op(0);
110         if (!is_a<add>(denom))
111                 return false;
112         if (denom.nops() != 2)
113                 return false;
114         if (denom.op(0) != z)
115                 return false;
116         if (!denom.op(1).info(info_flags::posint))
117                 return false;
118         n = denom.op(1);
119         return true;
120 }
122 /*
123  * Simplify the expression x by applying the rules
124  *
125  *    1/(z+n)/(z+m) = 1/(n-m)/(z+m) - 1/(n-m)/(z+n);
126  *    z^m/(z+n) = z^(m-1) - n*z^(m-1)/(z+n)
127  *
128  * as often as possible, where z is given as an argument to the function
129  * and n and m are arbitrary positive numbers.
130  */
131 ex poles_simplify(const ex &x, const ex &z)
132 {       if (is_a<mul>(x))
133         {       for (int i=0; i<x.nops(); ++i)
134                 {       ex arg1;
135                         if (!is_z_pole(x.op(i), z, arg1))
136                                 continue;
137                         for (int j=i+1; j<x.nops(); ++j)
138                         {       ex arg2;
139                                 if (!is_z_pole(x.op(j), z, arg2))
140                                         continue;
141                                 ex result = x/x.op(i)/(arg1-arg2)-x/x.op(j)/(arg1-arg2);
142                                 result = poles_simplify(result, z);
143                                 return result;
144                         }
145                 }
146                 ex result = expand(x);
147                 if (is_a<add>(result))
148                         return poles_simplify(result, z);
149                 for (int i=0; i<x.nops(); ++i)
150                 {       ex arg1;
151                         if (!is_z_pole(x.op(i), z, arg1))
152                                 continue;
153                         for (int j=0; j<x.nops(); ++j)
154                         {       ex expon = 0;
155                                 if (is_a<power>(x.op(j)) && x.op(j).op(0)==z
156                                                 && x.op(j).op(1).info(info_flags::posint))
157                                         expon = x.op(j).op(1);
158                                 if (x.op(j) == z)
159                                         expon = 1;
160                                 if (expon == 0)
161                                         continue;
162                                 ex result = x/x.op(i)/z - arg1*x/z;
163                                 result = poles_simplify(result, z);
164                                 return result;
165                         }
166                 }
167                 return x;
168         }
169         if (is_a<add>(x))
170         {       pointer_to_map_function_1arg<const ex&> mf(poles_simplify, z);
171                 return x.map(mf);
172         }
173         return x;
174 }
176 /*
177  * Calculate the expression A_g(z) that occurs in the expression for the
178  * Lanczos approximation of order n. This is returned as an expression of
179  * the form
180  *
181  *    A_g(z) = coeff + coeff/(z+1) + coeff/(z+2) + ...
182  *                + coeff[N-1]/(z+N-1).
183  */
184 ex A(const ex &g, const ex &z, int n)
185 {       vector<exvector> C;
186         initialize_C_vec(C, 2*n+2);
187         ex result = evalf(p(0, g, C))/2;
188         ex fraction = 1;
189         for (int i=1; i<n; ++i)
190         {       fraction *= (z-i+1)/(z+i);
191                 fraction = poles_simplify(fraction, z);
192                 result += evalf(p(i, g, C))*fraction;
193         }
194         result = poles_simplify(result, z);
195         return result;
196 }
198 /*
199  * The exvector coeffs should contain order elements and these are set to
200  * the coefficients that belong to the Lanczos approximation of the gamma
201  * function at the given order.
202  */
203 void calc_lanczos_coeffs(exvector &coeffs, double g_val, int order)
204 {       symbol g("g"), z("z");
205         ex result = A(g, z, order);
206         result = result.subs(g==g_val);
207         for (int i=0; i<result.nops(); ++i)
208         {       if (is_a<numeric>(result.op(i)))
209                         coeffs = result.op(i);
210                 else
211                 {       ex n;
212                         is_z_pole(result.op(i).op(0), z, n);
213                         coeffs[ex_to<numeric>(n).to_int()] = result.op(i).op(1);
214                 }
215         }
216 }
218 /*
219  * Calculate Gamma(z) using the Lanczos approximation with parameter g and
220  * coefficients stored in the exvector coeffs.
221  */
222 ex calc_gamma(const ex &z, const ex &g, exvector &coeffs)
223 {       if (real_part(evalf(z)) < 0.5)
224                  return evalf(Pi/sin(Pi*z)/calc_gamma(1-z, g, coeffs));
225         ex A = coeffs;
226         for (int i=1; i<coeffs.size(); ++i)
227                 A += evalf(coeffs[i]/(z-1+i));
228         ex result = sqrt(2*Pi)*power(z+g-ex(1)/2, z-ex(1)/2)*exp(-(z+g-ex(1)/2))*A;
229         return evalf(result);
230 }
232 void usage(char *progname)
233 {       cout << "Usage: " << progname << " -n order -D digits" << endl;
234         exit(0);
235 }
237 void read_options(int argc, char**argv, int &order)
238 {  int c;
239    while((c=getopt(argc,argv,"n:D:"))!=-1)
240    {  if(c=='n')
241          order = atoi(optarg);
242       else if (c=='D')
243          Digits = atoi(optarg);
244       else
245          usage(argv);
246    }
247    if(optind!=argc)
248       usage(argv);
249 }
251 /*
252  * Do something stupid to the number x in order to round it to the current
253  * numer of digits. Does somebody know a better way to do this? In that case,
254  * please fix me.
255  */
256 ex round(const ex &x)
257 {       return ex_to<numeric>(x).add(numeric("0.0"));
258 }
260 int main(int argc, char *argv[])
261 {
262 /*
263  * Handle command line options.
264  */
265         int order = 10;
266         read_options(argc, argv, order);
267 /*
268  * Calculate coefficients.
269  */
270         const int g_val = order;
271         exvector coeffs(order);
272         calc_lanczos_coeffs(coeffs, g_val, order);
273 /*
274  * Determine the accuracy by comparing Gamma(1/2) to sqrt(Pi).
275  */
276         ex gamma_half = calc_gamma(ex(1)/2, g_val, coeffs);
277         ex digits = -log(abs((gamma_half - sqrt(Pi)))/sqrt(Pi))/log(10.0);
278         digits = evalf(digits);
279         int i_digits = int(ex_to<numeric>(digits).to_double());
280         if (digits < 1)
281                 cout << "Forget it, this is waaaaaaay too inaccurate." << endl;
282         else
283                 cout << "Reliable digits: " << i_digits << endl;
284 /*
285  * Print the coefficients.
286  */
287         Digits = i_digits + 10; // Don't print too many spurious digits.
288         for (int i=0; i<coeffs.size(); ++i)
289                 cout << "coeffs_" << order << "[" << i << "] = numeric(\""
290                      << round(coeffs[i]) << "\");"<< endl;
291         return 0;
292 }