| // Copyright John Maddock 2006. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #include <boost/math/bindings/rr.hpp> |
| #include <boost/math/special_functions/log1p.hpp> |
| #include <boost/math/special_functions/erf.hpp> |
| #include <boost/math/constants/constants.hpp> |
| #include <map> |
| #include <iostream> |
| #include <iomanip> |
| |
| using namespace std; |
| using namespace NTL; |
| using namespace boost::math; |
| |
| // |
| // This program calculates the coefficients of the polynomials |
| // used for the regularized incomplete gamma functions gamma_p |
| // and gamma_q when parameter a is large, and sigma is small |
| // (where sigma = fabs(1 - x/a) ). |
| // |
| // See "The Asymptotic Expansion of the Incomplete Gamma Functions" |
| // N. M. Temme. |
| // Siam J. Math Anal. Vol 10 No 4, July 1979, p757. |
| // Coeffient calculation is described from Eq 3.8 (p762) onwards. |
| // |
| |
| // |
| // Alpha: |
| // |
| RR alpha(unsigned k) |
| { |
| static map<unsigned, RR> data; |
| if(data.empty()) |
| { |
| data[1] = 1; |
| } |
| |
| map<unsigned, RR>::const_iterator pos = data.find(k); |
| if(pos != data.end()) |
| return (*pos).second; |
| // |
| // OK try and calculate the value: |
| // |
| RR result = alpha(k-1); |
| for(unsigned j = 2; j <= k-1; ++j) |
| { |
| result -= j * alpha(j) * alpha(k-j+1); |
| } |
| result /= (k+1); |
| data[k] = result; |
| return result; |
| } |
| |
| boost::math::ntl::RR gamma(unsigned k) |
| { |
| static map<unsigned, boost::math::ntl::RR> data; |
| |
| map<unsigned, boost::math::ntl::RR>::const_iterator pos = data.find(k); |
| if(pos != data.end()) |
| return (*pos).second; |
| |
| boost::math::ntl::RR result = (k&1) ? -1 : 1; |
| |
| for(unsigned i = 1; i <= (2 * k + 1); i += 2) |
| result *= i; |
| result *= alpha(2 * k + 1); |
| data[k] = result; |
| return result; |
| } |
| |
| boost::math::ntl::RR Coeff(unsigned n, unsigned k) |
| { |
| map<unsigned, map<unsigned, boost::math::ntl::RR> > data; |
| if(data.empty()) |
| data[0][0] = boost::math::ntl::RR(-1) / 3; |
| |
| map<unsigned, map<unsigned, boost::math::ntl::RR> >::const_iterator p1 = data.find(n); |
| if(p1 != data.end()) |
| { |
| map<unsigned, boost::math::ntl::RR>::const_iterator p2 = p1->second.find(k); |
| if(p2 != p1->second.end()) |
| { |
| return p2->second; |
| } |
| } |
| |
| // |
| // If we don't have the value, calculate it: |
| // |
| if(k == 0) |
| { |
| // special case: |
| boost::math::ntl::RR result = (n+2) * alpha(n+2); |
| data[n][k] = result; |
| return result; |
| } |
| // general case: |
| boost::math::ntl::RR result = gamma(k) * Coeff(n, 0) + (n+2) * Coeff(n+2, k-1); |
| data[n][k] = result; |
| return result; |
| } |
| |
| void calculate_terms(double sigma, double a, unsigned bits) |
| { |
| cout << endl << endl; |
| cout << "Sigma: " << sigma << endl; |
| cout << "A: " << a << endl; |
| double lambda = 1 - sigma; |
| cout << "Lambda: " << lambda << endl; |
| double y = a * (-sigma - log1p(-sigma)); |
| cout << "Y: " << y << endl; |
| double z = -sqrt(2 * (-sigma - log1p(-sigma))); |
| cout << "Z: " << z << endl; |
| double dom = erfc(sqrt(y)) / 2; |
| cout << "Erfc term: " << dom << endl; |
| double lead = exp(-y) / sqrt(2 * constants::pi<double>() * a); |
| cout << "Remainder factor: " << lead << endl; |
| double eps = ldexp(1.0, 1 - static_cast<int>(bits)); |
| double target = dom * eps / lead; |
| cout << "Target smallest term: " << target << endl; |
| |
| unsigned max_n = 0; |
| |
| for(unsigned n = 0; n < 10000; ++n) |
| { |
| double term = tools::real_cast<double>(Coeff(n, 0) * pow(z, (double)n)); |
| if(fabs(term) < target) |
| { |
| max_n = n-1; |
| break; |
| } |
| } |
| cout << "Max n required: " << max_n << endl; |
| |
| unsigned max_k; |
| for(unsigned k = 1; k < 10000; ++k) |
| { |
| double term = tools::real_cast<double>(Coeff(0, k) * pow(a, -((double)k))); |
| if(fabs(term) < target) |
| { |
| max_k = k-1; |
| break; |
| } |
| } |
| cout << "Max k required: " << max_k << endl << endl; |
| |
| bool code = false; |
| cout << "Print code [0|1]? "; |
| cin >> code; |
| |
| int prec = 2 + (static_cast<double>(bits) * 3010LL)/10000; |
| RR::SetOutputPrecision(prec); |
| |
| if(code) |
| { |
| cout << " T workspace[" << max_k+1 << "];\n\n"; |
| for(unsigned k = 0; k <= max_k; ++k) |
| { |
| cout << |
| " static const T C" << k << "[] = {\n"; |
| for(unsigned n = 0; n < 10000; ++n) |
| { |
| double term = tools::real_cast<double>(Coeff(n, k) * pow(a, -((double)k)) * pow(z, (double)n)); |
| if(fabs(term) < target) |
| { |
| break; |
| } |
| cout << " " << Coeff(n, k) << "L,\n"; |
| } |
| cout << |
| " };\n" |
| " workspace[" << k << "] = tools::evaluate_polynomial(C" << k << ", z);\n\n"; |
| } |
| cout << " T result = tools::evaluate_polynomial(workspace, 1/a);\n\n"; |
| } |
| } |
| |
| |
| int main() |
| { |
| RR::SetOutputPrecision(50); |
| RR::SetPrecision(1000); |
| |
| bool cont; |
| do{ |
| cont = false; |
| double sigma; |
| cout << "Enter max value for sigma (sigma = |1 - x/a|): "; |
| cin >> sigma; |
| double a; |
| cout << "Enter min value for a: "; |
| cin >> a; |
| unsigned precision; |
| cout << "Enter number of bits precision required: "; |
| cin >> precision; |
| |
| calculate_terms(sigma, a, precision); |
| |
| cout << "Try again[0|1]: "; |
| cin >> cont; |
| |
| }while(cont); |
| |
| |
| return 0; |
| } |
| |