doc/examples/lanczos.cpp

   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[0] + coeff[1]/(z+1) + coeff[2]/(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  */
  50
  51 #include <ginac/ginac.h>
  52 #include <iostream>
  53
  54 using namespace std;
  55 using namespace GiNaC;
  56
  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 }
  68
  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][1] = 1;
  78         C[2][2] = 1;
  79         for (int i=3; i<size; i+=2)
  80                 C[i][1] = -C[i-2][1];
  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 }
  87
  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 }
  98
  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 }
 121
 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 }
 175
 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[0] + coeff[1]/(z+1) + coeff[2]/(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 }
 197
 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[0] = 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 }
 217
 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[0];
 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 }
 231
 232 void usage(char *progname)
 233 {       cout << "Usage: " << progname << " -n order -D digits" << endl;
 234         exit(0);
 235 }
 236
 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[0]);
 246    }
 247    if(optind!=argc)
 248       usage(argv[0]);
 249 }
 250
 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 }
 259
 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 }