special_functions/detail/erf_inv.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SF_ERF_INV_HPP
0007 #define BOOST_MATH_SF_ERF_INV_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #pragma warning(push)
0012 #pragma warning(disable:4127) // Conditional expression is constant
0013 #pragma warning(disable:4702) // Unreachable code: optimization warning
0014 #endif
0015
0016 #include <type_traits>
0017
0018 namespace boost{ namespace math{
0019
0020 namespace detail{
0021 //
0022 // The inverse erf and erfc functions share a common implementation,
0023 // this version is for 80-bit long double's and smaller:
0024 //
0025 template <class T, class Policy>
0026 T erf_inv_imp(const T& p, const T& q, const Policy&, const std::integral_constant<int, 64>*)
0027 {
0028    BOOST_MATH_STD_USING // for ADL of std names.
0029
0030    T result = 0;
0031
0032    if(p <= 0.5)
0033    {
0034       //
0035       // Evaluate inverse erf using the rational approximation:
0036       //
0037       // x = p(p+10)(Y+R(p))
0038       //
0039       // Where Y is a constant, and R(p) is optimised for a low
0040       // absolute error compared to |Y|.
0041       //
0042       // double: Max error found: 2.001849e-18
0043       // long double: Max error found: 1.017064e-20
0044       // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
0045       //
0046       static const float Y = 0.0891314744949340820313f;
0047       static const T P[] = {
0048          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000508781949658280665617),
0049          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00836874819741736770379),
0050          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0334806625409744615033),
0051          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0126926147662974029034),
0052          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0365637971411762664006),
0053          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0219878681111168899165),
0054          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00822687874676915743155),
0055          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00538772965071242932965)
0056       };
0057       static const T Q[] = {
0058          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0059          BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),
0060          BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),
0061          BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),
0062          BOOST_MATH_BIG_CONSTANT(T, 64, 0.662328840472002992063),
0063          BOOST_MATH_BIG_CONSTANT(T, 64, -0.71228902341542847553),
0064          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0527396382340099713954),
0065          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0795283687341571680018),
0066          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00233393759374190016776),
0067          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000886216390456424707504)
0068       };
0069       T g = p * (p + 10);
0070       T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
0071       result = g * Y + g * r;
0072    }
0073    else if(q >= 0.25)
0074    {
0075       //
0076       // Rational approximation for 0.5 > q >= 0.25
0077       //
0078       // x = sqrt(-2*log(q)) / (Y + R(q))
0079       //
0080       // Where Y is a constant, and R(q) is optimised for a low
0081       // absolute error compared to Y.
0082       //
0083       // double : Max error found: 7.403372e-17
0084       // long double : Max error found: 6.084616e-20
0085       // Maximum Deviation Found (error term) 4.811e-20
0086       //
0087       static const float Y = 2.249481201171875f;
0088       static const T P[] = {
0089          BOOST_MATH_BIG_CONSTANT(T, 64, -0.202433508355938759655),
0090          BOOST_MATH_BIG_CONSTANT(T, 64, 0.105264680699391713268),
0091          BOOST_MATH_BIG_CONSTANT(T, 64, 8.37050328343119927838),
0092          BOOST_MATH_BIG_CONSTANT(T, 64, 17.6447298408374015486),
0093          BOOST_MATH_BIG_CONSTANT(T, 64, -18.8510648058714251895),
0094          BOOST_MATH_BIG_CONSTANT(T, 64, -44.6382324441786960818),
0095          BOOST_MATH_BIG_CONSTANT(T, 64, 17.445385985570866523),
0096          BOOST_MATH_BIG_CONSTANT(T, 64, 21.1294655448340526258),
0097          BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)
0098       };
0099       static const T Q[] = {
0100          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0101          BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),
0102          BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),
0103          BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),
0104          BOOST_MATH_BIG_CONSTANT(T, 64, -20.1432634680485188801),
0105          BOOST_MATH_BIG_CONSTANT(T, 64, 48.5609213108739935468),
0106          BOOST_MATH_BIG_CONSTANT(T, 64, 10.8268667355460159008),
0107          BOOST_MATH_BIG_CONSTANT(T, 64, -22.6436933413139721736),
0108          BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724)
0109       };
0110       T g = sqrt(-2 * log(q));
0111       T xs = q - 0.25f;
0112       T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0113       result = g / (Y + r);
0114    }
0115    else
0116    {
0117       //
0118       // For q < 0.25 we have a series of rational approximations all
0119       // of the general form:
0120       //
0121       // let: x = sqrt(-log(q))
0122       //
0123       // Then the result is given by:
0124       //
0125       // x(Y+R(x-B))
0126       //
0127       // where Y is a constant, B is the lowest value of x for which
0128       // the approximation is valid, and R(x-B) is optimised for a low
0129       // absolute error compared to Y.
0130       //
0131       // Note that almost all code will really go through the first
0132       // or maybe second approximation.  After than we're dealing with very
0133       // small input values indeed: 80 and 128 bit long double's go all the
0134       // way down to ~ 1e-5000 so the "tail" is rather long...
0135       //
0136       T x = sqrt(-log(q));
0137       if(x < 3)
0138       {
0139          // Max error found: 1.089051e-20
0140          static const float Y = 0.807220458984375f;
0141          static const T P[] = {
0142             BOOST_MATH_BIG_CONSTANT(T, 64, -0.131102781679951906451),
0143             BOOST_MATH_BIG_CONSTANT(T, 64, -0.163794047193317060787),
0144             BOOST_MATH_BIG_CONSTANT(T, 64, 0.117030156341995252019),
0145             BOOST_MATH_BIG_CONSTANT(T, 64, 0.387079738972604337464),
0146             BOOST_MATH_BIG_CONSTANT(T, 64, 0.337785538912035898924),
0147             BOOST_MATH_BIG_CONSTANT(T, 64, 0.142869534408157156766),
0148             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0290157910005329060432),
0149             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00214558995388805277169),
0150             BOOST_MATH_BIG_CONSTANT(T, 64, -0.679465575181126350155e-6),
0151             BOOST_MATH_BIG_CONSTANT(T, 64, 0.285225331782217055858e-7),
0152             BOOST_MATH_BIG_CONSTANT(T, 64, -0.681149956853776992068e-9)
0153          };
0154          static const T Q[] = {
0155             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0156             BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),
0157             BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),
0158             BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),
0159             BOOST_MATH_BIG_CONSTANT(T, 64, 2.59301921623620271374),
0160             BOOST_MATH_BIG_CONSTANT(T, 64, 0.848854343457902036425),
0161             BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612),
0162             BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121)
0163          };
0164          T xs = x - 1.125f;
0165          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0166          result = Y * x + R * x;
0167       }
0168       else if(x < 6)
0169       {
0170          // Max error found: 8.389174e-21
0171          static const float Y = 0.93995571136474609375f;
0172          static const T P[] = {
0173             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0350353787183177984712),
0174             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00222426529213447927281),
0175             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0185573306514231072324),
0176             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00950804701325919603619),
0177             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00187123492819559223345),
0178             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000157544617424960554631),
0179             BOOST_MATH_BIG_CONSTANT(T, 64, 0.460469890584317994083e-5),
0180             BOOST_MATH_BIG_CONSTANT(T, 64, -0.230404776911882601748e-9),
0181             BOOST_MATH_BIG_CONSTANT(T, 64, 0.266339227425782031962e-11)
0182          };
0183          static const T Q[] = {
0184             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0185             BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),
0186             BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),
0187             BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),
0188             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0341589143670947727934),
0189             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00263861676657015992959),
0190             BOOST_MATH_BIG_CONSTANT(T, 64, 0.764675292302794483503e-4)
0191          };
0192          T xs = x - 3;
0193          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0194          result = Y * x + R * x;
0195       }
0196       else if(x < 18)
0197       {
0198          // Max error found: 1.481312e-19
0199          static const float Y = 0.98362827301025390625f;
0200          static const T P[] = {
0201             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0167431005076633737133),
0202             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00112951438745580278863),
0203             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00105628862152492910091),
0204             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000209386317487588078668),
0205             BOOST_MATH_BIG_CONSTANT(T, 64, 0.149624783758342370182e-4),
0206             BOOST_MATH_BIG_CONSTANT(T, 64, 0.449696789927706453732e-6),
0207             BOOST_MATH_BIG_CONSTANT(T, 64, 0.462596163522878599135e-8),
0208             BOOST_MATH_BIG_CONSTANT(T, 64, -0.281128735628831791805e-13),
0209             BOOST_MATH_BIG_CONSTANT(T, 64, 0.99055709973310326855e-16)
0210          };
0211          static const T Q[] = {
0212             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0213             BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),
0214             BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),
0215             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),
0216             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000964011807005165528527),
0217             BOOST_MATH_BIG_CONSTANT(T, 64, 0.275335474764726041141e-4),
0218             BOOST_MATH_BIG_CONSTANT(T, 64, 0.282243172016108031869e-6)
0219          };
0220          T xs = x - 6;
0221          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0222          result = Y * x + R * x;
0223       }
0224       else if(x < 44)
0225       {
0226          // Max error found: 5.697761e-20
0227          static const float Y = 0.99714565277099609375f;
0228          static const T P[] = {
0229             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0024978212791898131227),
0230             BOOST_MATH_BIG_CONSTANT(T, 64, -0.779190719229053954292e-5),
0231             BOOST_MATH_BIG_CONSTANT(T, 64, 0.254723037413027451751e-4),
0232             BOOST_MATH_BIG_CONSTANT(T, 64, 0.162397777342510920873e-5),
0233             BOOST_MATH_BIG_CONSTANT(T, 64, 0.396341011304801168516e-7),
0234             BOOST_MATH_BIG_CONSTANT(T, 64, 0.411632831190944208473e-9),
0235             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145596286718675035587e-11),
0236             BOOST_MATH_BIG_CONSTANT(T, 64, -0.116765012397184275695e-17)
0237          };
0238          static const T Q[] = {
0239             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0240             BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),
0241             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),
0242             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),
0243             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145007359818232637924e-4),
0244             BOOST_MATH_BIG_CONSTANT(T, 64, 0.144437756628144157666e-6),
0245             BOOST_MATH_BIG_CONSTANT(T, 64, 0.509761276599778486139e-9)
0246          };
0247          T xs = x - 18;
0248          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0249          result = Y * x + R * x;
0250       }
0251       else
0252       {
0253          // Max error found: 1.279746e-20
0254          static const float Y = 0.99941349029541015625f;
0255          static const T P[] = {
0256             BOOST_MATH_BIG_CONSTANT(T, 64, -0.000539042911019078575891),
0257             BOOST_MATH_BIG_CONSTANT(T, 64, -0.28398759004727721098e-6),
0258             BOOST_MATH_BIG_CONSTANT(T, 64, 0.899465114892291446442e-6),
0259             BOOST_MATH_BIG_CONSTANT(T, 64, 0.229345859265920864296e-7),
0260             BOOST_MATH_BIG_CONSTANT(T, 64, 0.225561444863500149219e-9),
0261             BOOST_MATH_BIG_CONSTANT(T, 64, 0.947846627503022684216e-12),
0262             BOOST_MATH_BIG_CONSTANT(T, 64, 0.135880130108924861008e-14),
0263             BOOST_MATH_BIG_CONSTANT(T, 64, -0.348890393399948882918e-21)
0264          };
0265          static const T Q[] = {
0266             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0267             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),
0268             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),
0269             BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),
0270             BOOST_MATH_BIG_CONSTANT(T, 64, 0.399968812193862100054e-6),
0271             BOOST_MATH_BIG_CONSTANT(T, 64, 0.161809290887904476097e-8),
0272             BOOST_MATH_BIG_CONSTANT(T, 64, 0.231558608310259605225e-11)
0273          };
0274          T xs = x - 44;
0275          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0276          result = Y * x + R * x;
0277       }
0278    }
0279    return result;
0280 }
0281
0282 template <class T, class Policy>
0283 struct erf_roots
0284 {
0285    boost::math::tuple<T,T,T> operator()(const T& guess)
0286    {
0287       BOOST_MATH_STD_USING
0288       T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
0289       T derivative2 = -2 * guess * derivative;
0290       return boost::math::make_tuple(((sign > 0) ? static_cast<T>(boost::math::erf(guess, Policy()) - target) : static_cast<T>(boost::math::erfc(guess, Policy())) - target), derivative, derivative2);
0291    }
0292    erf_roots(T z, int s) : target(z), sign(s) {}
0293 private:
0294    T target;
0295    int sign;
0296 };
0297
0298 template <class T, class Policy>
0299 T erf_inv_imp(const T& p, const T& q, const Policy& pol, const std::integral_constant<int, 0>*)
0300 {
0301    //
0302    // Generic version, get a guess that's accurate to 64-bits (10^-19)
0303    //
0304    T guess = erf_inv_imp(p, q, pol, static_cast<std::integral_constant<int, 64> const*>(nullptr));
0305    T result;
0306    //
0307    // If T has more bit's than 64 in it's mantissa then we need to iterate,
0308    // otherwise we can just return the result:
0309    //
0310    if(policies::digits<T, Policy>() > 64)
0311    {
0312       std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0313       if(p <= 0.5)
0314       {
0315          result = tools::halley_iterate(detail::erf_roots<typename std::remove_cv<T>::type, Policy>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
0316       }
0317       else
0318       {
0319          result = tools::halley_iterate(detail::erf_roots<typename std::remove_cv<T>::type, Policy>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (policies::digits<T, Policy>() * 2) / 3, max_iter);
0320       }
0321       policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);
0322    }
0323    else
0324    {
0325       result = guess;
0326    }
0327    return result;
0328 }
0329
0330 template <class T, class Policy>
0331 struct erf_inv_initializer
0332 {
0333    struct init
0334    {
0335       init()
0336       {
0337          do_init();
0338       }
0339       static bool is_value_non_zero(T);
0340       static void do_init()
0341       {
0342          // If std::numeric_limits<T>::digits is zero, we must not call
0343          // our initialization code here as the precision presumably
0344          // varies at runtime, and will not have been set yet.
0345          if(std::numeric_limits<T>::digits)
0346          {
0347             boost::math::erf_inv(static_cast<T>(0.25), Policy());
0348             boost::math::erf_inv(static_cast<T>(0.55), Policy());
0349             boost::math::erf_inv(static_cast<T>(0.95), Policy());
0350             boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
0351             // These following initializations must not be called if
0352             // type T can not hold the relevant values without
0353             // underflow to zero.  We check this at runtime because
0354             // some tools such as valgrind silently change the precision
0355             // of T at runtime, and numeric_limits basically lies!
0356             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130))))
0357                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());
0358
0359             // Some compilers choke on constants that would underflow, even in code that isn't instantiated
0360             // so try and filter these cases out in the preprocessor:
0361 #if LDBL_MAX_10_EXP >= 800
0362             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800))))
0363                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
0364             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900))))
0365                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
0366 #else
0367             if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800))))
0368                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800)), Policy());
0369             if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900))))
0370                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900)), Policy());
0371 #endif
0372          }
0373       }
0374       void force_instantiate()const{}
0375    };
0376    static const init initializer;
0377    static void force_instantiate()
0378    {
0379       initializer.force_instantiate();
0380    }
0381 };
0382
0383 template <class T, class Policy>
0384 const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer;
0385
0386 template <class T, class Policy>
0387 BOOST_NOINLINE bool erf_inv_initializer<T, Policy>::init::is_value_non_zero(T v)
0388 {
0389    // This needs to be non-inline to detect whether v is non zero at runtime
0390    // rather than at compile time, only relevant when running under valgrind
0391    // which changes long double's to double's on the fly.
0392    return v != 0;
0393 }
0394
0395 } // namespace detail
0396
0397 template <class T, class Policy>
0398 typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
0399 {
0400    typedef typename tools::promote_args<T>::type result_type;
0401
0402    //
0403    // Begin by testing for domain errors, and other special cases:
0404    //
0405    static const char* function = "boost::math::erfc_inv<%1%>(%1%, %1%)";
0406    if((z < 0) || (z > 2))
0407       return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
0408    if(z == 0)
0409       return policies::raise_overflow_error<result_type>(function, nullptr, pol);
0410    if(z == 2)
0411       return -policies::raise_overflow_error<result_type>(function, nullptr, pol);
0412    //
0413    // Normalise the input, so it's in the range [0,1], we will
0414    // negate the result if z is outside that range.  This is a simple
0415    // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
0416    //
0417    result_type p, q, s;
0418    if(z > 1)
0419    {
0420       q = 2 - z;
0421       p = 1 - q;
0422       s = -1;
0423    }
0424    else
0425    {
0426       p = 1 - z;
0427       q = z;
0428       s = 1;
0429    }
0430    //
0431    // A bit of meta-programming to figure out which implementation
0432    // to use, based on the number of bits in the mantissa of T:
0433    //
0434    typedef typename policies::precision<result_type, Policy>::type precision_type;
0435    typedef std::integral_constant<int,
0436       precision_type::value <= 0 ? 0 :
0437       precision_type::value <= 64 ? 64 : 0
0438    > tag_type;
0439    //
0440    // Likewise use internal promotion, so we evaluate at a higher
0441    // precision internally if it's appropriate:
0442    //
0443    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0444    typedef typename policies::normalise<
0445       Policy,
0446       policies::promote_float<false>,
0447       policies::promote_double<false>,
0448       policies::discrete_quantile<>,
0449       policies::assert_undefined<> >::type forwarding_policy;
0450
0451    detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
0452
0453    //
0454    // And get the result, negating where required:
0455    //
0456    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
0457       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(nullptr)), function);
0458 }
0459
0460 template <class T, class Policy>
0461 typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
0462 {
0463    typedef typename tools::promote_args<T>::type result_type;
0464
0465    //
0466    // Begin by testing for domain errors, and other special cases:
0467    //
0468    static const char* function = "boost::math::erf_inv<%1%>(%1%, %1%)";
0469    if((z < -1) || (z > 1))
0470       return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
0471    if(z == 1)
0472       return policies::raise_overflow_error<result_type>(function, nullptr, pol);
0473    if(z == -1)
0474       return -policies::raise_overflow_error<result_type>(function, nullptr, pol);
0475    if(z == 0)
0476       return 0;
0477    //
0478    // Normalise the input, so it's in the range [0,1], we will
0479    // negate the result if z is outside that range.  This is a simple
0480    // application of the erf reflection formula: erf(-z) = -erf(z)
0481    //
0482    result_type p, q, s;
0483    if(z < 0)
0484    {
0485       p = -z;
0486       q = 1 - p;
0487       s = -1;
0488    }
0489    else
0490    {
0491       p = z;
0492       q = 1 - z;
0493       s = 1;
0494    }
0495    //
0496    // A bit of meta-programming to figure out which implementation
0497    // to use, based on the number of bits in the mantissa of T:
0498    //
0499    typedef typename policies::precision<result_type, Policy>::type precision_type;
0500    typedef std::integral_constant<int,
0501       precision_type::value <= 0 ? 0 :
0502       precision_type::value <= 64 ? 64 : 0
0503    > tag_type;
0504    //
0505    // Likewise use internal promotion, so we evaluate at a higher
0506    // precision internally if it's appropriate:
0507    //
0508    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0509    typedef typename policies::normalise<
0510       Policy,
0511       policies::promote_float<false>,
0512       policies::promote_double<false>,
0513       policies::discrete_quantile<>,
0514       policies::assert_undefined<> >::type forwarding_policy;
0515    //
0516    // Likewise use internal promotion, so we evaluate at a higher
0517    // precision internally if it's appropriate:
0518    //
0519    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0520
0521    detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
0522    //
0523    // And get the result, negating where required:
0524    //
0525    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
0526       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(nullptr)), function);
0527 }
0528
0529 template <class T>
0530 inline typename tools::promote_args<T>::type erfc_inv(T z)
0531 {
0532    return erfc_inv(z, policies::policy<>());
0533 }
0534
0535 template <class T>
0536 inline typename tools::promote_args<T>::type erf_inv(T z)
0537 {
0538    return erf_inv(z, policies::policy<>());
0539 }
0540
0541 } // namespace math
0542 } // namespace boost
0543
0544 #ifdef _MSC_VER
0545 #pragma warning(pop)
0546 #endif
0547
0548 #endif // BOOST_MATH_SF_ERF_INV_HPP
0549