// (C) 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) #define L22 //#include "../tools/ntl_rr_lanczos.hpp" //#include "../tools/ntl_rr_digamma.hpp" #include <boost/math/bindings/rr.hpp> #include <boost/math/tools/polynomial.hpp> #include <boost/math/special_functions.hpp> #include <boost/math/special_functions/zeta.hpp> #include <boost/math/special_functions/expint.hpp> #include <cmath> boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant) { static const boost::math::ntl::RR tiny = boost::math::tools::min_value<boost::math::ntl::RR>() * 64; switch(variant) { case 0: { boost::math::ntl::RR x_ = sqrt(x == 0 ? 1e-80 : x); return boost::math::erf(x_) / x_; } case 1: { boost::math::ntl::RR x_ = 1 / x; return boost::math::erfc(x_) * x_ / exp(-x_ * x_); } case 2: { return boost::math::erfc(x) * x / exp(-x * x); } case 3: { boost::math::ntl::RR y(x); if(y == 0) y += tiny; return boost::math::lgamma(y+2) / y - 0.5; } case 4: // // lgamma in the range [2,3], use: // // lgamma(x) = (x-2) * (x + 1) * (c + R(x - 2)) // // Works well at 80-bit long double precision, but doesn't // stretch to 128-bit precision. // if(x == 0) { return boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008") / 3; } return boost::math::lgamma(x+2) / (x * (x+3)); case 5: { // // lgamma in the range [1,2], use: // // lgamma(x) = (x - 1) * (x - 2) * (c + R(x - 1)) // // works well over [1, 1.5] but not near 2 :-( // boost::math::ntl::RR r1 = boost::lexical_cast<boost::math::ntl::RR>("0.57721566490153286060651209008240243104215933593992"); boost::math::ntl::RR r2 = boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008"); if(x == 0) { return r1; } if(x == 1) { return r2; } return boost::math::lgamma(x+1) / (x * (x - 1)); } case 6: { // // lgamma in the range [1.5,2], use: // // lgamma(x) = (2 - x) * (1 - x) * (c + R(2 - x)) // // works well over [1.5, 2] but not near 1 :-( // boost::math::ntl::RR r1 = boost::lexical_cast<boost::math::ntl::RR>("0.57721566490153286060651209008240243104215933593992"); boost::math::ntl::RR r2 = boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008"); if(x == 0) { return r2; } if(x == 1) { return r1; } return boost::math::lgamma(2-x) / (x * (x - 1)); } case 7: { // // erf_inv in range [0, 0.5] // boost::math::ntl::RR y = x; if(y == 0) y = boost::math::tools::epsilon<boost::math::ntl::RR>() / 64; return boost::math::erf_inv(y) / (y * (y+10)); } case 8: { // // erfc_inv in range [0.25, 0.5] // Use an y-offset of 0.25, and range [0, 0.25] // abs error, auto y-offset. // boost::math::ntl::RR y = x; if(y == 0) y = boost::lexical_cast<boost::math::ntl::RR>("1e-5000"); return sqrt(-2 * log(y)) / boost::math::erfc_inv(y); } case 9: { boost::math::ntl::RR x2 = x; if(x2 == 0) x2 = boost::lexical_cast<boost::math::ntl::RR>("1e-5000"); boost::math::ntl::RR y = exp(-x2*x2); // sqrt(-log(x2)) - 5; return boost::math::erfc_inv(y) / x2; } case 10: { // // Digamma over the interval [1,2], set x-offset to 1 // and optimise for absolute error over [0,1]. // int current_precision = boost::math::ntl::RR::precision(); if(current_precision < 1000) boost::math::ntl::RR::SetPrecision(1000); // // This value for the root of digamma is calculated using our // differentiated lanczos approximation. It agrees with Cody // to ~ 25 digits and to Morris to 35 digits. See: // TOMS ALGORITHM 708 (Didonato and Morris). // and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher. // //boost::math::ntl::RR root = boost::lexical_cast<boost::math::ntl::RR>("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871"); // // Actually better to calculate the root on the fly, it appears to be more // accurate: convergence is easier with the 1000-bit value, the approximation // produced agrees with functions.mathworld.com values to 35 digits even quite // near the root. // static boost::math::tools::eps_tolerance<boost::math::ntl::RR> tol(1000); static boost::uintmax_t max_iter = 1000; boost::math::ntl::RR (*pdg)(boost::math::ntl::RR) = &boost::math::digamma; static const boost::math::ntl::RR root = boost::math::tools::bracket_and_solve_root(pdg, boost::math::ntl::RR(1.4), boost::math::ntl::RR(1.5), true, tol, max_iter).first; boost::math::ntl::RR x2 = x; double lim = 1e-65; if(fabs(x2 - root) < lim) { // // This is a problem area: // x2-root suffers cancellation error, so does digamma. // That gets compounded again when Remez calculates the error // function. This cludge seems to stop the worst of the problems: // static const boost::math::ntl::RR a = boost::math::digamma(root - lim) / -lim; static const boost::math::ntl::RR b = boost::math::digamma(root + lim) / lim; boost::math::ntl::RR fract = (x2 - root + lim) / (2*lim); boost::math::ntl::RR r = (1-fract) * a + fract * b; std::cout << "In root area: " << r; return r; } boost::math::ntl::RR result = boost::math::digamma(x2) / (x2 - root); if(current_precision < 1000) boost::math::ntl::RR::SetPrecision(current_precision); return result; } case 11: // expm1: if(x == 0) { static boost::math::ntl::RR lim = 1e-80; static boost::math::ntl::RR a = boost::math::expm1(-lim); static boost::math::ntl::RR b = boost::math::expm1(lim); static boost::math::ntl::RR l = (b-a) / (2 * lim); return l; } return boost::math::expm1(x) / x; case 12: // demo, and test case: return exp(x); case 13: // K(k): { // x = k^2 boost::math::ntl::RR k2 = x; if(k2 > boost::math::ntl::RR(1) - 1e-40) k2 = boost::math::ntl::RR(1) - 1e-40; /*if(k < 1e-40) k = 1e-40;*/ boost::math::ntl::RR p2 = boost::math::constants::pi<boost::math::ntl::RR>() / 2; return (boost::math::ellint_1(sqrt(k2))) / (p2 - boost::math::log1p(-k2)); } case 14: // K(k) { static double P[] = { 1.38629436111989062502E0, 9.65735902811690126535E-2, 3.08851465246711995998E-2, 1.49380448916805252718E-2, 8.79078273952743772254E-3, 6.18901033637687613229E-3, 6.87489687449949877925E-3, 9.85821379021226008714E-3, 7.97404013220415179367E-3, 2.28025724005875567385E-3, 1.37982864606273237150E-4 }; // x = 1 - k^2 boost::math::ntl::RR mp = x; if(mp < 1e-20) mp = 1e-20; if(mp == 1) mp -= 1e-20; boost::math::ntl::RR k = sqrt(1 - mp); return (boost::math::ellint_1(k) - mp/*boost::math::tools::evaluate_polynomial(P, mp)*/) / -log(mp); } case 15: // E(k) { // x = 1-k^2 boost::math::ntl::RR z = 1 - x * log(x); return boost::math::ellint_2(sqrt(1-x)) / z; } case 16: // Bessel I0(x) over [0,16]: { return boost::math::cyl_bessel_i(0, sqrt(x)); } case 17: // Bessel I0(x) over [16,INF] { boost::math::ntl::RR z = 1 / (boost::math::ntl::RR(1)/16 - x); return boost::math::cyl_bessel_i(0, z) * sqrt(z) / exp(z); } case 18: // Zeta over [0, 1] { return boost::math::zeta(1 - x) * x - x; } case 19: // Zeta over [1, n] { return boost::math::zeta(x) - 1 / (x - 1); } case 20: // Zeta over [a, b] : a >> 1 { return log(boost::math::zeta(x) - 1); } case 21: // expint[1] over [0,1]: { boost::math::ntl::RR tiny = boost::lexical_cast<boost::math::ntl::RR>("1e-5000"); boost::math::ntl::RR z = (x <= tiny) ? tiny : x; return boost::math::expint(1, z) - z + log(z); } case 22: // expint[1] over [1,N], // Note that x varies from [0,1]: { boost::math::ntl::RR z = 1 / x; return boost::math::expint(1, z) * exp(z) * z; } case 23: // expin Ei over [0,R] { static const boost::math::ntl::RR root = boost::lexical_cast<boost::math::ntl::RR>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392"); boost::math::ntl::RR z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::min)() : x; return (boost::math::expint(z) - log(z / root)) / (z - root); } case 24: // Expint Ei for large x: { static const boost::math::ntl::RR root = boost::lexical_cast<boost::math::ntl::RR>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392"); boost::math::ntl::RR z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::max)() : 1 / x; return (boost::math::expint(z) - z) * z * exp(-z); //return (boost::math::expint(z) - log(z)) * z * exp(-z); } case 25: // Expint Ei for large x: { return (boost::math::expint(x) - x) * x * exp(-x); } case 26: { // // erf_inv in range [0, 0.5] // boost::math::ntl::RR y = x; if(y == 0) y = boost::math::tools::epsilon<boost::math::ntl::RR>() / 64; y = sqrt(y); return boost::math::erf_inv(y) / (y); } case 28: { // log1p over [-0.5,0.5] boost::math::ntl::RR y = x; if(fabs(y) < 1e-100) y = (y == 0) ? 1e-100 : boost::math::sign(y) * 1e-100; return (boost::math::log1p(y) - y + y * y / 2) / (y); } case 29: { // cbrt over [0.5, 1] return boost::math::cbrt(x); } } return 0; } void show_extra( const boost::math::tools::polynomial<boost::math::ntl::RR>& n, const boost::math::tools::polynomial<boost::math::ntl::RR>& d, const boost::math::ntl::RR& x_offset, const boost::math::ntl::RR& y_offset, int variant) { switch(variant) { default: // do nothing here... ; } }