e5643d9c36391c0254d03fad8cf417361cc85719
[ginac.git] / check / time_gammaseries.cpp
1 /** @file time_gammaseries.cpp
2  *
3  *  Some timings on series expansion of the Gamma function around a pole. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include "times.h"
24
25 unsigned tgammaseries(unsigned order)
26 {
27         unsigned result = 0;
28         symbol x;
29
30         ex myseries = series(tgamma(x),x==0,order);
31         // compute the last coefficient numerically:
32         ex last_coeff = myseries.coeff(x,order-1).evalf();
33         // compute a bound for that coefficient using a variation of the leading
34         // term in Stirling's formula:
35         ex bound = exp(-.57721566490153286*(order-1))/(order-1);
36         if (abs((last_coeff-pow(-1,order))/bound) > 1) {
37                 clog << "The " << order-1
38                      << "th order coefficient in the power series expansion of tgamma(0) was erroneously found to be "
39                      << last_coeff << ", violating a simple estimate." << endl;
40                 ++result;
41         }
42
43         return result;
44 }
45
46 unsigned time_gammaseries(void)
47 {
48         unsigned result = 0;
49
50         cout << "timing Laurent series expansion of Gamma function" << flush;
51         clog << "-------Laurent series expansion of Gamma function:" << endl;
52
53         vector<unsigned> sizes;
54         vector<double> times;
55         timer omega;
56
57         sizes.push_back(10);
58         sizes.push_back(15);
59         sizes.push_back(20);
60         sizes.push_back(25);
61
62         for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i) {
63                 omega.start();
64                 result += tgammaseries(*i);
65                 times.push_back(omega.read());
66                 cout << '.' << flush;
67         }
68
69         if (!result) {
70                 cout << " passed ";
71                 clog << "(no output)" << endl;
72         } else {
73                 cout << " failed ";
74         }
75         // print the report:
76         cout << endl << "       order: ";
77         for (vector<unsigned>::iterator i=sizes.begin(); i!=sizes.end(); ++i)
78                 cout << '\t' << (*i);
79         cout << endl << "       time/s:";
80         for (vector<double>::iterator i=times.begin(); i!=times.end(); ++i)
81                 cout << '\t' << int(1000*(*i))*0.001;
82         cout << endl;
83         
84         return result;
85 }