1 // Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #include <boost/math/special_functions/log1p.hpp>
7 #include <boost/math/special_functions/erf.hpp>
8 #include <boost/math/constants/constants.hpp>
9 #include <map>
10 #include <iostream>
11 #include <iomanip>
12 #include "mp_t.hpp"
13
14 using namespace std;
15 using namespace boost::math;
16
17 //
18 // This program calculates the coefficients of the polynomials
19 // used for the regularized incomplete gamma functions gamma_p
20 // and gamma_q when parameter a is large, and sigma is small
21 // (where sigma = fabs(1 - x/a) ).
22 //
23 // See "The Asymptotic Expansion of the Incomplete Gamma Functions"
24 // N. M. Temme.
25 // Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
26 // Coefficient calculation is described from Eq 3.8 (p762) onwards.
27 //
28
29 //
30 // Alpha:
31 //
alpha(unsigned k)32 mp_t alpha(unsigned k)
33 {
34 static map<unsigned, mp_t> data;
35 if(data.empty())
36 {
37 data[1] = 1;
38 }
39
40 map<unsigned, mp_t>::const_iterator pos = data.find(k);
41 if(pos != data.end())
42 return (*pos).second;
43 //
44 // OK try and calculate the value:
45 //
46 mp_t result = alpha(k-1);
47 for(unsigned j = 2; j <= k-1; ++j)
48 {
49 result -= j * alpha(j) * alpha(k-j+1);
50 }
51 result /= (k+1);
52 data[k] = result;
53 return result;
54 }
55
gamma(unsigned k)56 mp_t gamma(unsigned k)
57 {
58 static map<unsigned, mp_t> data;
59
60 map<unsigned, mp_t>::const_iterator pos = data.find(k);
61 if(pos != data.end())
62 return (*pos).second;
63
64 mp_t result = (k&1) ? -1 : 1;
65
66 for(unsigned i = 1; i <= (2 * k + 1); i += 2)
67 result *= i;
68 result *= alpha(2 * k + 1);
69 data[k] = result;
70 return result;
71 }
72
Coeff(unsigned n,unsigned k)73 mp_t Coeff(unsigned n, unsigned k)
74 {
75 map<unsigned, map<unsigned, mp_t> > data;
76 if(data.empty())
77 data[0][0] = mp_t(-1) / 3;
78
79 map<unsigned, map<unsigned, mp_t> >::const_iterator p1 = data.find(n);
80 if(p1 != data.end())
81 {
82 map<unsigned, mp_t>::const_iterator p2 = p1->second.find(k);
83 if(p2 != p1->second.end())
84 {
85 return p2->second;
86 }
87 }
88
89 //
90 // If we don't have the value, calculate it:
91 //
92 if(k == 0)
93 {
94 // special case:
95 mp_t result = (n+2) * alpha(n+2);
96 data[n][k] = result;
97 return result;
98 }
99 // general case:
100 mp_t result = gamma(k) * Coeff(n, 0) + (n+2) * Coeff(n+2, k-1);
101 data[n][k] = result;
102 return result;
103 }
104
calculate_terms(double sigma,double a,unsigned bits)105 void calculate_terms(double sigma, double a, unsigned bits)
106 {
107 cout << endl << endl;
108 cout << "Sigma: " << sigma << endl;
109 cout << "A: " << a << endl;
110 double lambda = 1 - sigma;
111 cout << "Lambda: " << lambda << endl;
112 double y = a * (-sigma - log1p(-sigma));
113 cout << "Y: " << y << endl;
114 double z = -sqrt(2 * (-sigma - log1p(-sigma)));
115 cout << "Z: " << z << endl;
116 double dom = erfc(sqrt(y)) / 2;
117 cout << "Erfc term: " << dom << endl;
118 double lead = exp(-y) / sqrt(2 * constants::pi<double>() * a);
119 cout << "Remainder factor: " << lead << endl;
120 double eps = ldexp(1.0, 1 - static_cast<int>(bits));
121 double target = dom * eps / lead;
122 cout << "Target smallest term: " << target << endl;
123
124 unsigned max_n = 0;
125
126 for(unsigned n = 0; n < 10000; ++n)
127 {
128 double term = tools::real_cast<double>(Coeff(n, 0) * pow(z, (double)n));
129 if(fabs(term) < target)
130 {
131 max_n = n-1;
132 break;
133 }
134 }
135 cout << "Max n required: " << max_n << endl;
136
137 unsigned max_k;
138 for(unsigned k = 1; k < 10000; ++k)
139 {
140 double term = tools::real_cast<double>(Coeff(0, k) * pow(a, -((double)k)));
141 if(fabs(term) < target)
142 {
143 max_k = k-1;
144 break;
145 }
146 }
147 cout << "Max k required: " << max_k << endl << endl;
148
149 bool code = false;
150 cout << "Print code [0|1]? ";
151 cin >> code;
152
153 int prec = 2 + (static_cast<double>(bits) * 3010LL)/10000;
154 std::cout << std::scientific << std::setprecision(40);
155
156 if(code)
157 {
158 cout << " T workspace[" << max_k+1 << "];\n\n";
159 for(unsigned k = 0; k <= max_k; ++k)
160 {
161 cout <<
162 " static const T C" << k << "[] = {\n";
163 for(unsigned n = 0; n < 10000; ++n)
164 {
165 double term = tools::real_cast<double>(Coeff(n, k) * pow(a, -((double)k)) * pow(z, (double)n));
166 if(fabs(term) < target)
167 {
168 break;
169 }
170 cout << " " << Coeff(n, k) << "L,\n";
171 }
172 cout <<
173 " };\n"
174 " workspace[" << k << "] = tools::evaluate_polynomial(C" << k << ", z);\n\n";
175 }
176 cout << " T result = tools::evaluate_polynomial(workspace, 1/a);\n\n";
177 }
178 }
179
180
main()181 int main()
182 {
183 bool cont;
184 do{
185 cont = false;
186 double sigma;
187 cout << "Enter max value for sigma (sigma = |1 - x/a|): ";
188 cin >> sigma;
189 double a;
190 cout << "Enter min value for a: ";
191 cin >> a;
192 unsigned precision;
193 cout << "Enter number of bits precision required: ";
194 cin >> precision;
195
196 calculate_terms(sigma, a, precision);
197
198 cout << "Try again[0|1]: ";
199 cin >> cont;
200
201 }while(cont);
202
203
204 return 0;
205 }
206
207