2 * This program can be used to find the coefficients needed to approximate
3 * the gamma function using the Lanczos approximation.
5 * Usage: lanczos -n order -D digits
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
19 * The gamma function can be (for real_part(z) > -1/2) calculated using
21 * Gamma(z+1) = sqrt(2*Pi)*power(z+g+ex(1)/2, z+ex(1)/2)*exp(-(z+g+ex(1)/2))
25 * A_g(z) = coeff[0] + coeff[1]/(z+1) + coeff[2]/(z+2) + ...
26 * + coeff[N-1]/(z+N-1).
28 * The value of g is taken to be equal to the order N.
30 * More details can be found at Wikipedia:
31 * http://en.wikipedia.org/wiki/Lanczos_approximation.
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.
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.
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,
51 #include <ginac/ginac.h>
55 using namespace GiNaC;
58 * Double factorial function.
64 return i*doublefact(i-2);
65 cout << "function doublefact called with wrong argument" << endl;
70 * Chebyshev polynomial coefficient matrix as far as is required for
71 * the Lanczos approximation.
73 void initialize_C_vec(vector<exvector> &C, int size)
75 for (int i=0; i<size; ++i)
76 C.push_back(exvector(size, 0));
79 for (int i=3; i<size; 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];
89 * The coefficients p_n(g) that occur in the Lanczos approximation.
91 ex p(int k, ex g, const vector<exvector> &C)
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);
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.
104 bool is_z_pole(const ex &x, const ex &z, ex &n)
105 { if (!is_a<power>(x))
110 if (!is_a<add>(denom))
112 if (denom.nops() != 2)
114 if (denom.op(0) != z)
116 if (!denom.op(1).info(info_flags::posint))
123 * Simplify the expression x by applying the rules
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)
128 * as often as possible, where z is given as an argument to the function
129 * and n and m are arbitrary positive numbers.
131 ex poles_simplify(const ex &x, const ex &z)
133 { for (int i=0; i<x.nops(); ++i)
135 if (!is_z_pole(x.op(i), z, arg1))
137 for (int j=i+1; j<x.nops(); ++j)
139 if (!is_z_pole(x.op(j), z, arg2))
141 ex result = x/x.op(i)/(arg1-arg2)-x/x.op(j)/(arg1-arg2);
142 result = poles_simplify(result, z);
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)
151 if (!is_z_pole(x.op(i), z, arg1))
153 for (int j=0; j<x.nops(); ++j)
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);
162 ex result = x/x.op(i)/z - arg1*x/z;
163 result = poles_simplify(result, z);
170 { pointer_to_map_function_1arg<const ex&> mf(poles_simplify, z);
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
181 * A_g(z) = coeff[0] + coeff[1]/(z+1) + coeff[2]/(z+2) + ...
182 * + coeff[N-1]/(z+N-1).
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;
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;
194 result = poles_simplify(result, z);
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.
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[0] = result.op(i);
212 is_z_pole(result.op(i).op(0), z, n);
213 coeffs[ex_to<numeric>(n).to_int()] = result.op(i).op(1);
219 * Calculate Gamma(z) using the Lanczos approximation with parameter g and
220 * coefficients stored in the exvector coeffs.
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));
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);
232 void usage(char *progname)
233 { cout << "Usage: " << progname << " -n order -D digits" << endl;
237 void read_options(int argc, char**argv, int &order)
239 while((c=getopt(argc,argv,"n:D:"))!=-1)
241 order = atoi(optarg);
243 Digits = atoi(optarg);
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,
256 ex round(const ex &x)
257 { return ex_to<numeric>(x).add(numeric("0.0"));
260 int main(int argc, char *argv[])
263 * Handle command line options.
266 read_options(argc, argv, order);
268 * Calculate coefficients.
270 const int g_val = order;
271 exvector coeffs(order);
272 calc_lanczos_coeffs(coeffs, g_val, order);
274 * Determine the accuracy by comparing Gamma(1/2) to sqrt(Pi).
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());
281 cout << "Forget it, this is waaaaaaay too inaccurate." << endl;
283 cout << "Reliable digits: " << i_digits << endl;
285 * Print the coefficients.
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;