special_functions/detail/erf_inv.hpp

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