]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/tools/igamma_temme_large_coef.cpp
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)
6 #include <boost/math/special_functions/log1p.hpp>
7 #include <boost/math/special_functions/erf.hpp>
8 #include <boost/math/constants/constants.hpp>
15 using namespace boost::math
;
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) ).
23 // See "The Asymptotic Expansion of the Incomplete Gamma Functions"
25 // Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
26 // Coeffient calculation is described from Eq 3.8 (p762) onwards.
32 mp_t
alpha(unsigned k
)
34 static map
<unsigned, mp_t
> data
;
40 map
<unsigned, mp_t
>::const_iterator pos
= data
.find(k
);
44 // OK try and calculate the value:
46 mp_t result
= alpha(k
-1);
47 for(unsigned j
= 2; j
<= k
-1; ++j
)
49 result
-= j
* alpha(j
) * alpha(k
-j
+1);
56 mp_t
gamma(unsigned k
)
58 static map
<unsigned, mp_t
> data
;
60 map
<unsigned, mp_t
>::const_iterator pos
= data
.find(k
);
64 mp_t result
= (k
&1) ? -1 : 1;
66 for(unsigned i
= 1; i
<= (2 * k
+ 1); i
+= 2)
68 result
*= alpha(2 * k
+ 1);
73 mp_t
Coeff(unsigned n
, unsigned k
)
75 map
<unsigned, map
<unsigned, mp_t
> > data
;
77 data
[0][0] = mp_t(-1) / 3;
79 map
<unsigned, map
<unsigned, mp_t
> >::const_iterator p1
= data
.find(n
);
82 map
<unsigned, mp_t
>::const_iterator p2
= p1
->second
.find(k
);
83 if(p2
!= p1
->second
.end())
90 // If we don't have the value, calculate it:
95 mp_t result
= (n
+2) * alpha(n
+2);
100 mp_t result
= gamma(k
) * Coeff(n
, 0) + (n
+2) * Coeff(n
+2, k
-1);
105 void calculate_terms(double sigma
, double a
, unsigned bits
)
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
;
126 for(unsigned n
= 0; n
< 10000; ++n
)
128 double term
= tools::real_cast
<double>(Coeff(n
, 0) * pow(z
, (double)n
));
129 if(fabs(term
) < target
)
135 cout
<< "Max n required: " << max_n
<< endl
;
138 for(unsigned k
= 1; k
< 10000; ++k
)
140 double term
= tools::real_cast
<double>(Coeff(0, k
) * pow(a
, -((double)k
)));
141 if(fabs(term
) < target
)
147 cout
<< "Max k required: " << max_k
<< endl
<< endl
;
150 cout
<< "Print code [0|1]? ";
153 int prec
= 2 + (static_cast<double>(bits
) * 3010LL)/10000;
154 std::cout
<< std::scientific
<< std::setprecision(40);
158 cout
<< " T workspace[" << max_k
+1 << "];\n\n";
159 for(unsigned k
= 0; k
<= max_k
; ++k
)
162 " static const T C" << k
<< "[] = {\n";
163 for(unsigned n
= 0; n
< 10000; ++n
)
165 double term
= tools::real_cast
<double>(Coeff(n
, k
) * pow(a
, -((double)k
)) * pow(z
, (double)n
));
166 if(fabs(term
) < target
)
170 cout
<< " " << Coeff(n
, k
) << "L,\n";
174 " workspace[" << k
<< "] = tools::evaluate_polynomial(C" << k
<< ", z);\n\n";
176 cout
<< " T result = tools::evaluate_polynomial(workspace, 1/a);\n\n";
187 cout
<< "Enter max value for sigma (sigma = |1 - x/a|): ";
190 cout
<< "Enter min value for a: ";
193 cout
<< "Enter number of bits precision required: ";
196 calculate_terms(sigma
, a
, precision
);
198 cout
<< "Try again[0|1]: ";