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 <vector>
  52 #include <cstddef> // for size_t
  53 #include <iostream>
  54 #include <cln/io.h>
  55 #include <cln/number.h>
  56 #include <cln/integer.h>
  57 #include <cln/rational.h>
  58 #include <cln/real.h>
  59 #include <cln/real_io.h>
  60 #include <cln/complex.h>
  61
  62 using namespace std;
  63 using namespace cln;
  64
  65 /*
  66  * Chebyshev polynomial coefficient matrix as far as is required for
  67  * the Lanczos approximation.
  68  */
  69 void calc_chebyshev(vector<vector<cl_I> >& C, const size_t size)
  70 {
  71         C.reserve(size);
  72         for (size_t i=0; i<size; ++i)
  73                 C.push_back(vector<cl_I>(size, 0));
  74         C[1][1] = cl_I(1);
  75         C[2][2] = cl_I(1);
  76         for (size_t i=3; i<size; i+=2)
  77                 C[i][1] = -C[i-2][1];
  78         for (size_t i=3; i<size; ++i)
  79                 C[i][i] = 2*C[i-1][i-1];
  80         for (size_t j=2; j<size; ++j)
  81                 for (size_t i=j+2; i<size; i+=2)
  82                         C[i][j] = 2*C[i-1][j-1] - C[i-2][j];
  83 }
  84
  85 /*
  86  * The coefficients p_n(g) that occur in the Lanczos approximation.
  87  */
  88 const cl_F p(const size_t k, const cl_F& g, const vector<vector<cln::cl_I> >& C)
  89 {
  90         const float_format_t prec = float_format(g);
  91         const cl_F sqrtPi = sqrt(pi(prec));
  92         cl_F result = cl_float(0, prec);
  93
  94         result = result +
  95                 (2*C[2*k+1][1])*recip(sqrtPi)*
  96                 recip(sqrt(2*g+1))*exp(g+cl_RA(1)/2);
  97
  98         for (size_t a=1; a <= k; ++a)
  99                 result = result +
 100                         (2*C[2*k+1][2*a+1]*doublefactorial(2*a-1))*recip(sqrtPi)*
 101                         recip(expt(2*cl_I(a)+2*g+1, a))*recip(sqrt(2*cl_I(a)+2*g+1))*
 102                         exp(cl_I(a)+g+cl_RA(1)/2);
 103         return result;
 104 }
 105
 106 /*
 107  * Calculate coefficients that occurs in the expression for the
 108  * Lanczos approximation of order n.
 109  */
 110 void calc_lanczos_coeffs(vector<cl_F>& lanc, const cln::cl_F& g)
 111 {
 112         const size_t n = lanc.size();
 113         vector<vector<cln::cl_I> > C;
 114         calc_chebyshev(C, 2*n+2);
 115
 116         // \Pi_{i=1}^n (z-i+1)/(z+i) = \Pi_{i=1}^n (1 - (2i-1)/(z+i))
 117         // Such a product can be rewritten as multivariate polynomial \in
 118         // Q[1/(z+1),...1/(z+n)] of degree 1. To store coefficients of this
 119         // polynomial we use vector<cl_I>, so the set of such polynomials is
 120         // stored as
 121         vector<vector<cln::cl_I> > fractions(n);
 122
 123         // xi = 1/(z+i)
 124         fractions[0] = vector<cl_I>(1);
 125         fractions[0][0] = cl_I(1);
 126         fractions[1] = fractions[0];
 127         fractions[1].push_back(cl_I(-1)); // 1 - 1/(z+1)
 128
 129         for (size_t i=2; i<n; ++i)
 130         {
 131                 // next = previous*(1 - (2*i-1)*xi) =
 132                 // previous - (2*i-1)*xi - (2i-1)*(previous-1)*xi
 133                 // Such representation is convenient since previous contains only
 134                 // only x1...x[i-1]
 135
 136                 fractions[i] = fractions[i-1];
 137                 fractions[i].push_back(-cl_I(2*i-1)); // - (2*i-1)*xi
 138                 // Now add -(2i+1)*(previous-1)*xi using formula
 139                 // 1/(z+j)*1/(z+i) = 1/(i-j)*(1/(z+j) - 1/(z+i)), that is
 140                 // xj*xi = 1/(i-j)(xj - xi)
 141                 for (size_t j=1; j < i; j++) {
 142                         cl_I curr = - exquo(cl_I(2*i-1)*fractions[i-1][j], cl_I(i-j));
 143                         fractions[i][j] = fractions[i][j] + curr;
 144                         fractions[i][i] = fractions[i][i] - curr;
 145                         // N.B. i-th polynomial has (i+1) non-zero coefficients
 146                 }
 147         }
 148
 149         vector<cl_F> p_cache(n);
 150         for (size_t i=0; i < n; i++)
 151                 p_cache[i] = p(i, g, C);
 152
 153         lanc[0] = p_cache[0]/2;
 154
 155         float_format_t prec = float_format(g);
 156         // A = p(i, g, C)*fraction[i] = p(i, g, C)*F[i][j]*xj,
 157         for (size_t j=1; j < n; ++j) {
 158                 lanc[j] = cl_float(0, prec);
 159                 lanc[0] = lanc[0] + p_cache[j];
 160                 for (size_t i=j; i < n; i++)
 161                         lanc[j] = lanc[j] + p_cache[i]*fractions[i][j];
 162         }
 163 }
 164
 165 /*
 166  * Calculate Gamma(z) using the Lanczos approximation with parameter g and
 167  * coefficients stored in the exvector coeffs.
 168  */
 169 const cl_N calc_gamma(const cl_N& z, const cl_F& g, const vector<cl_F>& lanc)
 170 {
 171         const cl_F thePi = pi(float_format(g));
 172         // XXX: need to check if z is floating-point and precision of z
 173         // matches one of g
 174         if (realpart(z) < 0.5)
 175                 return (thePi/sin(thePi*z))/calc_gamma(1-z, g, lanc);
 176         cl_N A = lanc[0];
 177         for (size_t i=1; i < lanc.size(); ++i)
 178                 A = A + lanc[i]/(z-1+i);
 179         cl_N result = sqrt(2*thePi)*expt(z+g-cl_RA(1)/2, z-cl_RA(1)/2)*
 180                             exp(-(z+g-cl_RA(1)/2))*A;
 181         return result;
 182 }
 183
 184 void usage(char *progname)
 185 {       cout << "Usage: " << progname << " -n order -D digits" << endl;
 186         exit(0);
 187 }
 188
 189 void read_options(int argc, char**argv, int &order, int& digits)
 190 {  int c;
 191    while((c=getopt(argc,argv,"n:D:"))!=-1)
 192    {  if(c=='n')
 193          order = atoi(optarg);
 194       else if (c=='D')
 195          digits = atoi(optarg);
 196       else
 197          usage(argv[0]);
 198    }
 199    if(optind!=argc)
 200       usage(argv[0]);
 201 }
 202
 203 int main(int argc, char *argv[])
 204 {
 205 /*
 206  * Handle command line options.
 207  */
 208         int order = 10;
 209         int digits_ini = 17;
 210         read_options(argc, argv, order, digits_ini);
 211         float_format_t prec = float_format(digits_ini);
 212         const cl_F thePi = pi(prec);
 213 /*
 214  * Calculate coefficients.
 215  */
 216         const cl_F g_val = cl_float(order, prec);
 217         vector<cl_F> coeffs(order);
 218         calc_lanczos_coeffs(coeffs, g_val);
 219 /*
 220  * Determine the accuracy by comparing Gamma(1/2) to sqrt(Pi).
 221  */
 222         cl_N gamma_half = calc_gamma(cl_float(cl_RA(1)/2, prec), g_val, coeffs);
 223         cl_F digits = (ln(thePi)/2 - ln(abs(gamma_half - sqrt(thePi))))/ln(cl_float(10, prec));
 224         int i_digits = cl_I_to_int(floor1(digits));
 225         if (digits < cl_I(1))
 226                 cout << "Forget it, this is waaaaaaay too inaccurate." << endl;
 227         else
 228                 cout << "Reliable digits: " << i_digits << endl;
 229 /*
 230  * Print the coefficients.
 231  */
 232         for (size_t i=0; i<coeffs.size(); ++i) {
 233                 cout << "coeffs_" << order << "[" << i << "] = \"";
 234                 cout << coeffs[i];
 235                 cout << "\";" << endl;
 236         }
 237         return 0;
 238 }