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       // LCOV_EXCL_START
0047       static const float Y = 0.0891314744949340820313f;
0048       static const T P[] = {
0049          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000508781949658280665617),
0050          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00836874819741736770379),
0051          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0334806625409744615033),
0052          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0126926147662974029034),
0053          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0365637971411762664006),
0054          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0219878681111168899165),
0055          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00822687874676915743155),
0056          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00538772965071242932965)
0057       };
0058       static const T Q[] = {
0059          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0060          BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),
0061          BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),
0062          BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),
0063          BOOST_MATH_BIG_CONSTANT(T, 64, 0.662328840472002992063),
0064          BOOST_MATH_BIG_CONSTANT(T, 64, -0.71228902341542847553),
0065          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0527396382340099713954),
0066          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0795283687341571680018),
0067          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00233393759374190016776),
0068          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000886216390456424707504)
0069       };
0070       // LCOV_EXCL_STOP
0071       T g = p * (p + 10);
0072       T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
0073       result = g * Y + g * r;
0074    }
0075    else if(q >= 0.25)
0076    {
0077       //
0078       // Rational approximation for 0.5 > q >= 0.25
0079       //
0080       // x = sqrt(-2*log(q)) / (Y + R(q))
0081       //
0082       // Where Y is a constant, and R(q) is optimised for a low
0083       // absolute error compared to Y.
0084       //
0085       // double : Max error found: 7.403372e-17
0086       // long double : Max error found: 6.084616e-20
0087       // Maximum Deviation Found (error term) 4.811e-20
0088       //
0089       // LCOV_EXCL_START
0090       static const float Y = 2.249481201171875f;
0091       static const T P[] = {
0092          BOOST_MATH_BIG_CONSTANT(T, 64, -0.202433508355938759655),
0093          BOOST_MATH_BIG_CONSTANT(T, 64, 0.105264680699391713268),
0094          BOOST_MATH_BIG_CONSTANT(T, 64, 8.37050328343119927838),
0095          BOOST_MATH_BIG_CONSTANT(T, 64, 17.6447298408374015486),
0096          BOOST_MATH_BIG_CONSTANT(T, 64, -18.8510648058714251895),
0097          BOOST_MATH_BIG_CONSTANT(T, 64, -44.6382324441786960818),
0098          BOOST_MATH_BIG_CONSTANT(T, 64, 17.445385985570866523),
0099          BOOST_MATH_BIG_CONSTANT(T, 64, 21.1294655448340526258),
0100          BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)
0101       };
0102       static const T Q[] = {
0103          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0104          BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),
0105          BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),
0106          BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),
0107          BOOST_MATH_BIG_CONSTANT(T, 64, -20.1432634680485188801),
0108          BOOST_MATH_BIG_CONSTANT(T, 64, 48.5609213108739935468),
0109          BOOST_MATH_BIG_CONSTANT(T, 64, 10.8268667355460159008),
0110          BOOST_MATH_BIG_CONSTANT(T, 64, -22.6436933413139721736),
0111          BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724)
0112       };
0113       // LCOV_EXCL_STOP
0114       T g = sqrt(-2 * log(q));
0115       T xs = q - 0.25f;
0116       T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0117       result = g / (Y + r);
0118    }
0119    else
0120    {
0121       //
0122       // For q < 0.25 we have a series of rational approximations all
0123       // of the general form:
0124       //
0125       // let: x = sqrt(-log(q))
0126       //
0127       // Then the result is given by:
0128       //
0129       // x(Y+R(x-B))
0130       //
0131       // where Y is a constant, B is the lowest value of x for which
0132       // the approximation is valid, and R(x-B) is optimised for a low
0133       // absolute error compared to Y.
0134       //
0135       // Note that almost all code will really go through the first
0136       // or maybe second approximation.  After than we're dealing with very
0137       // small input values indeed: 80 and 128 bit long double's go all the
0138       // way down to ~ 1e-5000 so the "tail" is rather long...
0139       //
0140       T x = sqrt(-log(q));
0141       if(x < 3)
0142       {
0143          // LCOV_EXCL_START
0144          // Max error found: 1.089051e-20
0145          static const float Y = 0.807220458984375f;
0146          static const T P[] = {
0147             BOOST_MATH_BIG_CONSTANT(T, 64, -0.131102781679951906451),
0148             BOOST_MATH_BIG_CONSTANT(T, 64, -0.163794047193317060787),
0149             BOOST_MATH_BIG_CONSTANT(T, 64, 0.117030156341995252019),
0150             BOOST_MATH_BIG_CONSTANT(T, 64, 0.387079738972604337464),
0151             BOOST_MATH_BIG_CONSTANT(T, 64, 0.337785538912035898924),
0152             BOOST_MATH_BIG_CONSTANT(T, 64, 0.142869534408157156766),
0153             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0290157910005329060432),
0154             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00214558995388805277169),
0155             BOOST_MATH_BIG_CONSTANT(T, 64, -0.679465575181126350155e-6),
0156             BOOST_MATH_BIG_CONSTANT(T, 64, 0.285225331782217055858e-7),
0157             BOOST_MATH_BIG_CONSTANT(T, 64, -0.681149956853776992068e-9)
0158          };
0159          static const T Q[] = {
0160             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0161             BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),
0162             BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),
0163             BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),
0164             BOOST_MATH_BIG_CONSTANT(T, 64, 2.59301921623620271374),
0165             BOOST_MATH_BIG_CONSTANT(T, 64, 0.848854343457902036425),
0166             BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612),
0167             BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121)
0168          };
0169          // LCOV_EXCL_STOP
0170          T xs = x - 1.125f;
0171          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0172          result = Y * x + R * x;
0173       }
0174       else if(x < 6)
0175       {
0176          // LCOV_EXCL_START
0177          // Max error found: 8.389174e-21
0178          static const float Y = 0.93995571136474609375f;
0179          static const T P[] = {
0180             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0350353787183177984712),
0181             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00222426529213447927281),
0182             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0185573306514231072324),
0183             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00950804701325919603619),
0184             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00187123492819559223345),
0185             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000157544617424960554631),
0186             BOOST_MATH_BIG_CONSTANT(T, 64, 0.460469890584317994083e-5),
0187             BOOST_MATH_BIG_CONSTANT(T, 64, -0.230404776911882601748e-9),
0188             BOOST_MATH_BIG_CONSTANT(T, 64, 0.266339227425782031962e-11)
0189          };
0190          static const T Q[] = {
0191             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0192             BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),
0193             BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),
0194             BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),
0195             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0341589143670947727934),
0196             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00263861676657015992959),
0197             BOOST_MATH_BIG_CONSTANT(T, 64, 0.764675292302794483503e-4)
0198          };
0199          // LCOV_EXCL_STOP
0200          T xs = x - 3;
0201          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0202          result = Y * x + R * x;
0203       }
0204       else if(x < 18)
0205       {
0206          // LCOV_EXCL_START
0207          // Max error found: 1.481312e-19
0208          static const float Y = 0.98362827301025390625f;
0209          static const T P[] = {
0210             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0167431005076633737133),
0211             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00112951438745580278863),
0212             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00105628862152492910091),
0213             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000209386317487588078668),
0214             BOOST_MATH_BIG_CONSTANT(T, 64, 0.149624783758342370182e-4),
0215             BOOST_MATH_BIG_CONSTANT(T, 64, 0.449696789927706453732e-6),
0216             BOOST_MATH_BIG_CONSTANT(T, 64, 0.462596163522878599135e-8),
0217             BOOST_MATH_BIG_CONSTANT(T, 64, -0.281128735628831791805e-13),
0218             BOOST_MATH_BIG_CONSTANT(T, 64, 0.99055709973310326855e-16)
0219          };
0220          static const T Q[] = {
0221             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0222             BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),
0223             BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),
0224             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),
0225             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000964011807005165528527),
0226             BOOST_MATH_BIG_CONSTANT(T, 64, 0.275335474764726041141e-4),
0227             BOOST_MATH_BIG_CONSTANT(T, 64, 0.282243172016108031869e-6)
0228          };
0229          // LCOV_EXCL_STOP
0230          T xs = x - 6;
0231          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0232          result = Y * x + R * x;
0233       }
0234       else if(x < 44)
0235       {
0236          // LCOV_EXCL_START
0237          // Max error found: 5.697761e-20
0238          static const float Y = 0.99714565277099609375f;
0239          static const T P[] = {
0240             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0024978212791898131227),
0241             BOOST_MATH_BIG_CONSTANT(T, 64, -0.779190719229053954292e-5),
0242             BOOST_MATH_BIG_CONSTANT(T, 64, 0.254723037413027451751e-4),
0243             BOOST_MATH_BIG_CONSTANT(T, 64, 0.162397777342510920873e-5),
0244             BOOST_MATH_BIG_CONSTANT(T, 64, 0.396341011304801168516e-7),
0245             BOOST_MATH_BIG_CONSTANT(T, 64, 0.411632831190944208473e-9),
0246             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145596286718675035587e-11),
0247             BOOST_MATH_BIG_CONSTANT(T, 64, -0.116765012397184275695e-17)
0248          };
0249          static const T Q[] = {
0250             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0251             BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),
0252             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),
0253             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),
0254             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145007359818232637924e-4),
0255             BOOST_MATH_BIG_CONSTANT(T, 64, 0.144437756628144157666e-6),
0256             BOOST_MATH_BIG_CONSTANT(T, 64, 0.509761276599778486139e-9)
0257          };
0258          // LCOV_EXCL_STOP
0259          T xs = x - 18;
0260          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0261          result = Y * x + R * x;
0262       }
0263       else
0264       {
0265          // LCOV_EXCL_START
0266          // Max error found: 1.279746e-20
0267          static const float Y = 0.99941349029541015625f;
0268          static const T P[] = {
0269             BOOST_MATH_BIG_CONSTANT(T, 64, -0.000539042911019078575891),
0270             BOOST_MATH_BIG_CONSTANT(T, 64, -0.28398759004727721098e-6),
0271             BOOST_MATH_BIG_CONSTANT(T, 64, 0.899465114892291446442e-6),
0272             BOOST_MATH_BIG_CONSTANT(T, 64, 0.229345859265920864296e-7),
0273             BOOST_MATH_BIG_CONSTANT(T, 64, 0.225561444863500149219e-9),
0274             BOOST_MATH_BIG_CONSTANT(T, 64, 0.947846627503022684216e-12),
0275             BOOST_MATH_BIG_CONSTANT(T, 64, 0.135880130108924861008e-14),
0276             BOOST_MATH_BIG_CONSTANT(T, 64, -0.348890393399948882918e-21)
0277          };
0278          static const T Q[] = {
0279             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0280             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),
0281             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),
0282             BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),
0283             BOOST_MATH_BIG_CONSTANT(T, 64, 0.399968812193862100054e-6),
0284             BOOST_MATH_BIG_CONSTANT(T, 64, 0.161809290887904476097e-8),
0285             BOOST_MATH_BIG_CONSTANT(T, 64, 0.231558608310259605225e-11)
0286          };
0287          // LCOV_EXCL_STOP
0288          T xs = x - 44;
0289          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
0290          result = Y * x + R * x;
0291       }
0292    }
0293    return result;
0294 }
0295
0296 template <class T, class Policy>
0297 struct erf_roots
0298 {
0299    boost::math::tuple<T,T,T> operator()(const T& guess)
0300    {
0301       BOOST_MATH_STD_USING
0302       T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
0303       T derivative2 = -2 * guess * derivative;
0304       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);
0305    }
0306    erf_roots(T z, int s) : target(z), sign(s) {}
0307 private:
0308    T target;
0309    int sign;
0310 };
0311
0312 template <class T, class Policy>
0313 T erf_inv_imp(const T& p, const T& q, const Policy& pol, const std::integral_constant<int, 0>*)
0314 {
0315    //
0316    // Generic version, get a guess that's accurate to 64-bits (10^-19)
0317    //
0318    T guess = erf_inv_imp(p, q, pol, static_cast<std::integral_constant<int, 64> const*>(nullptr));
0319    T result;
0320    //
0321    // If T has more bit's than 64 in it's mantissa then we need to iterate,
0322    // otherwise we can just return the result:
0323    //
0324    if(policies::digits<T, Policy>() > 64)
0325    {
0326       std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0327       if(p <= 0.5)
0328       {
0329          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);
0330       }
0331       else
0332       {
0333          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);
0334       }
0335       policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);
0336    }
0337    else
0338    {
0339       result = guess;
0340    }
0341    return result;
0342 }
0343
0344 } // namespace detail
0345
0346 template <class T, class Policy>
0347 typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
0348 {
0349    typedef typename tools::promote_args<T>::type result_type;
0350
0351    //
0352    // Begin by testing for domain errors, and other special cases:
0353    //
0354    static const char* function = "boost::math::erfc_inv<%1%>(%1%, %1%)";
0355    if((z < 0) || (z > 2))
0356       return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
0357    if(z == 0)
0358       return policies::raise_overflow_error<result_type>(function, nullptr, pol);
0359    if(z == 2)
0360       return -policies::raise_overflow_error<result_type>(function, nullptr, pol);
0361    //
0362    // Normalise the input, so it's in the range [0,1], we will
0363    // negate the result if z is outside that range.  This is a simple
0364    // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
0365    //
0366    result_type p, q, s;
0367    if(z > 1)
0368    {
0369       q = 2 - z;
0370       p = 1 - q;
0371       s = -1;
0372    }
0373    else
0374    {
0375       p = 1 - z;
0376       q = z;
0377       s = 1;
0378    }
0379    //
0380    // A bit of meta-programming to figure out which implementation
0381    // to use, based on the number of bits in the mantissa of T:
0382    //
0383    typedef typename policies::precision<result_type, Policy>::type precision_type;
0384    typedef std::integral_constant<int,
0385       precision_type::value <= 0 ? 0 :
0386       precision_type::value <= 64 ? 64 : 0
0387    > tag_type;
0388    //
0389    // Likewise use internal promotion, so we evaluate at a higher
0390    // precision internally if it's appropriate:
0391    //
0392    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0393    typedef typename policies::normalise<
0394       Policy,
0395       policies::promote_float<false>,
0396       policies::promote_double<false>,
0397       policies::discrete_quantile<>,
0398       policies::assert_undefined<> >::type forwarding_policy;
0399
0400    //
0401    // And get the result, negating where required:
0402    //
0403    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
0404       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(nullptr)), function);
0405 }
0406
0407 template <class T, class Policy>
0408 typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
0409 {
0410    typedef typename tools::promote_args<T>::type result_type;
0411
0412    //
0413    // Begin by testing for domain errors, and other special cases:
0414    //
0415    static const char* function = "boost::math::erf_inv<%1%>(%1%, %1%)";
0416    if((z < -1) || (z > 1))
0417       return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
0418    if(z == 1)
0419       return policies::raise_overflow_error<result_type>(function, nullptr, pol);
0420    if(z == -1)
0421       return -policies::raise_overflow_error<result_type>(function, nullptr, pol);
0422    if(z == 0)
0423       return 0;
0424    //
0425    // Normalise the input, so it's in the range [0,1], we will
0426    // negate the result if z is outside that range.  This is a simple
0427    // application of the erf reflection formula: erf(-z) = -erf(z)
0428    //
0429    result_type p, q, s;
0430    if(z < 0)
0431    {
0432       p = -z;
0433       q = 1 - p;
0434       s = -1;
0435    }
0436    else
0437    {
0438       p = z;
0439       q = 1 - z;
0440       s = 1;
0441    }
0442    //
0443    // A bit of meta-programming to figure out which implementation
0444    // to use, based on the number of bits in the mantissa of T:
0445    //
0446    typedef typename policies::precision<result_type, Policy>::type precision_type;
0447    typedef std::integral_constant<int,
0448       precision_type::value <= 0 ? 0 :
0449       precision_type::value <= 64 ? 64 : 0
0450    > tag_type;
0451    //
0452    // Likewise use internal promotion, so we evaluate at a higher
0453    // precision internally if it's appropriate:
0454    //
0455    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0456    typedef typename policies::normalise<
0457       Policy,
0458       policies::promote_float<false>,
0459       policies::promote_double<false>,
0460       policies::discrete_quantile<>,
0461       policies::assert_undefined<> >::type forwarding_policy;
0462    //
0463    // Likewise use internal promotion, so we evaluate at a higher
0464    // precision internally if it's appropriate:
0465    //
0466    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
0467
0468    //
0469    // And get the result, negating where required:
0470    //
0471    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
0472       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(nullptr)), function);
0473 }
0474
0475 template <class T>
0476 inline typename tools::promote_args<T>::type erfc_inv(T z)
0477 {
0478    return erfc_inv(z, policies::policy<>());
0479 }
0480
0481 template <class T>
0482 inline typename tools::promote_args<T>::type erf_inv(T z)
0483 {
0484    return erf_inv(z, policies::policy<>());
0485 }
0486
0487 } // namespace math
0488 } // namespace boost
0489
0490 #ifdef _MSC_VER
0491 #pragma warning(pop)
0492 #endif
0493
0494 #endif // BOOST_MATH_SF_ERF_INV_HPP
0495