560 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SF_ERF_INV_HPP8#define BOOST_MATH_SF_ERF_INV_HPP9 10#ifdef _MSC_VER11#pragma once12#pragma warning(push)13#pragma warning(disable:4127) // Conditional expression is constant14#pragma warning(disable:4702) // Unreachable code: optimization warning15#endif16 17#include <boost/math/tools/config.hpp>18 19#ifndef BOOST_MATH_HAS_NVRTC20 21#include <type_traits>22 23namespace boost{ namespace math{24 25namespace detail{26//27// The inverse erf and erfc functions share a common implementation,28// this version is for 80-bit long double's and smaller:29//30template <class T, class Policy>31BOOST_MATH_GPU_ENABLED T erf_inv_imp(const T& p, const T& q, const Policy&, const std::integral_constant<int, 64>&)32{33 BOOST_MATH_STD_USING // for ADL of std names.34 35 T result = 0;36 37 if(p <= 0.5)38 {39 //40 // Evaluate inverse erf using the rational approximation:41 //42 // x = p(p+10)(Y+R(p))43 //44 // Where Y is a constant, and R(p) is optimised for a low45 // absolute error compared to |Y|.46 //47 // double: Max error found: 2.001849e-1848 // long double: Max error found: 1.017064e-2049 // Maximum Deviation Found (actual error term at infinite precision) 8.030e-2150 //51 // LCOV_EXCL_START52 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.0891314744949340820313f;53 BOOST_MATH_STATIC const T P[] = {54 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000508781949658280665617),55 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00836874819741736770379),56 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0334806625409744615033),57 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0126926147662974029034),58 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0365637971411762664006),59 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0219878681111168899165),60 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00822687874676915743155),61 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00538772965071242932965)62 };63 BOOST_MATH_STATIC const T Q[] = {64 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),65 BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),66 BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),67 BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),68 BOOST_MATH_BIG_CONSTANT(T, 64, 0.662328840472002992063),69 BOOST_MATH_BIG_CONSTANT(T, 64, -0.71228902341542847553),70 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0527396382340099713954),71 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0795283687341571680018),72 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00233393759374190016776),73 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000886216390456424707504)74 };75 // LCOV_EXCL_STOP76 T g = p * (p + 10);77 T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);78 result = g * Y + g * r;79 }80 else if(q >= 0.25)81 {82 //83 // Rational approximation for 0.5 > q >= 0.2584 //85 // x = sqrt(-2*log(q)) / (Y + R(q))86 //87 // Where Y is a constant, and R(q) is optimised for a low88 // absolute error compared to Y.89 //90 // double : Max error found: 7.403372e-1791 // long double : Max error found: 6.084616e-2092 // Maximum Deviation Found (error term) 4.811e-2093 //94 // LCOV_EXCL_START95 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 2.249481201171875f;96 BOOST_MATH_STATIC const T P[] = {97 BOOST_MATH_BIG_CONSTANT(T, 64, -0.202433508355938759655),98 BOOST_MATH_BIG_CONSTANT(T, 64, 0.105264680699391713268),99 BOOST_MATH_BIG_CONSTANT(T, 64, 8.37050328343119927838),100 BOOST_MATH_BIG_CONSTANT(T, 64, 17.6447298408374015486),101 BOOST_MATH_BIG_CONSTANT(T, 64, -18.8510648058714251895),102 BOOST_MATH_BIG_CONSTANT(T, 64, -44.6382324441786960818),103 BOOST_MATH_BIG_CONSTANT(T, 64, 17.445385985570866523),104 BOOST_MATH_BIG_CONSTANT(T, 64, 21.1294655448340526258),105 BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)106 };107 BOOST_MATH_STATIC const T Q[] = {108 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),109 BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),110 BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),111 BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),112 BOOST_MATH_BIG_CONSTANT(T, 64, -20.1432634680485188801),113 BOOST_MATH_BIG_CONSTANT(T, 64, 48.5609213108739935468),114 BOOST_MATH_BIG_CONSTANT(T, 64, 10.8268667355460159008),115 BOOST_MATH_BIG_CONSTANT(T, 64, -22.6436933413139721736),116 BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724)117 };118 // LCOV_EXCL_STOP119 T g = sqrt(-2 * log(q));120 T xs = q - 0.25f;121 T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);122 result = g / (Y + r);123 }124 else125 {126 //127 // For q < 0.25 we have a series of rational approximations all128 // of the general form:129 //130 // let: x = sqrt(-log(q))131 //132 // Then the result is given by:133 //134 // x(Y+R(x-B))135 //136 // where Y is a constant, B is the lowest value of x for which137 // the approximation is valid, and R(x-B) is optimised for a low138 // absolute error compared to Y.139 //140 // Note that almost all code will really go through the first141 // or maybe second approximation. After than we're dealing with very142 // small input values indeed: 80 and 128 bit long double's go all the143 // way down to ~ 1e-5000 so the "tail" is rather long...144 //145 T x = sqrt(-log(q));146 if(x < 3)147 {148 // LCOV_EXCL_START149 // Max error found: 1.089051e-20150 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.807220458984375f;151 BOOST_MATH_STATIC const T P[] = {152 BOOST_MATH_BIG_CONSTANT(T, 64, -0.131102781679951906451),153 BOOST_MATH_BIG_CONSTANT(T, 64, -0.163794047193317060787),154 BOOST_MATH_BIG_CONSTANT(T, 64, 0.117030156341995252019),155 BOOST_MATH_BIG_CONSTANT(T, 64, 0.387079738972604337464),156 BOOST_MATH_BIG_CONSTANT(T, 64, 0.337785538912035898924),157 BOOST_MATH_BIG_CONSTANT(T, 64, 0.142869534408157156766),158 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0290157910005329060432),159 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00214558995388805277169),160 BOOST_MATH_BIG_CONSTANT(T, 64, -0.679465575181126350155e-6),161 BOOST_MATH_BIG_CONSTANT(T, 64, 0.285225331782217055858e-7),162 BOOST_MATH_BIG_CONSTANT(T, 64, -0.681149956853776992068e-9)163 };164 BOOST_MATH_STATIC const T Q[] = {165 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),166 BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),167 BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),168 BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),169 BOOST_MATH_BIG_CONSTANT(T, 64, 2.59301921623620271374),170 BOOST_MATH_BIG_CONSTANT(T, 64, 0.848854343457902036425),171 BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612),172 BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121)173 };174 // LCOV_EXCL_STOP175 T xs = x - 1.125f;176 T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);177 result = Y * x + R * x;178 }179 else if(x < 6)180 {181 // LCOV_EXCL_START182 // Max error found: 8.389174e-21183 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.93995571136474609375f;184 BOOST_MATH_STATIC const T P[] = {185 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0350353787183177984712),186 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00222426529213447927281),187 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0185573306514231072324),188 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00950804701325919603619),189 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00187123492819559223345),190 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000157544617424960554631),191 BOOST_MATH_BIG_CONSTANT(T, 64, 0.460469890584317994083e-5),192 BOOST_MATH_BIG_CONSTANT(T, 64, -0.230404776911882601748e-9),193 BOOST_MATH_BIG_CONSTANT(T, 64, 0.266339227425782031962e-11)194 };195 BOOST_MATH_STATIC const T Q[] = {196 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),197 BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),198 BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),199 BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),200 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0341589143670947727934),201 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00263861676657015992959),202 BOOST_MATH_BIG_CONSTANT(T, 64, 0.764675292302794483503e-4)203 };204 // LCOV_EXCL_STOP205 T xs = x - 3;206 T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);207 result = Y * x + R * x;208 }209 else if(x < 18)210 {211 // LCOV_EXCL_START212 // Max error found: 1.481312e-19213 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.98362827301025390625f;214 BOOST_MATH_STATIC const T P[] = {215 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0167431005076633737133),216 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00112951438745580278863),217 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00105628862152492910091),218 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000209386317487588078668),219 BOOST_MATH_BIG_CONSTANT(T, 64, 0.149624783758342370182e-4),220 BOOST_MATH_BIG_CONSTANT(T, 64, 0.449696789927706453732e-6),221 BOOST_MATH_BIG_CONSTANT(T, 64, 0.462596163522878599135e-8),222 BOOST_MATH_BIG_CONSTANT(T, 64, -0.281128735628831791805e-13),223 BOOST_MATH_BIG_CONSTANT(T, 64, 0.99055709973310326855e-16)224 };225 BOOST_MATH_STATIC const T Q[] = {226 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),227 BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),228 BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),229 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),230 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000964011807005165528527),231 BOOST_MATH_BIG_CONSTANT(T, 64, 0.275335474764726041141e-4),232 BOOST_MATH_BIG_CONSTANT(T, 64, 0.282243172016108031869e-6)233 };234 // LCOV_EXCL_STOP235 T xs = x - 6;236 T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);237 result = Y * x + R * x;238 }239 else if(x < 44)240 {241 // LCOV_EXCL_START242 // Max error found: 5.697761e-20243 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.99714565277099609375f;244 BOOST_MATH_STATIC const T P[] = {245 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0024978212791898131227),246 BOOST_MATH_BIG_CONSTANT(T, 64, -0.779190719229053954292e-5),247 BOOST_MATH_BIG_CONSTANT(T, 64, 0.254723037413027451751e-4),248 BOOST_MATH_BIG_CONSTANT(T, 64, 0.162397777342510920873e-5),249 BOOST_MATH_BIG_CONSTANT(T, 64, 0.396341011304801168516e-7),250 BOOST_MATH_BIG_CONSTANT(T, 64, 0.411632831190944208473e-9),251 BOOST_MATH_BIG_CONSTANT(T, 64, 0.145596286718675035587e-11),252 BOOST_MATH_BIG_CONSTANT(T, 64, -0.116765012397184275695e-17)253 };254 BOOST_MATH_STATIC const T Q[] = {255 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),256 BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),257 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),258 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),259 BOOST_MATH_BIG_CONSTANT(T, 64, 0.145007359818232637924e-4),260 BOOST_MATH_BIG_CONSTANT(T, 64, 0.144437756628144157666e-6),261 BOOST_MATH_BIG_CONSTANT(T, 64, 0.509761276599778486139e-9)262 };263 // LCOV_EXCL_STOP264 T xs = x - 18;265 T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);266 result = Y * x + R * x;267 }268 else269 {270 // LCOV_EXCL_START271 // Max error found: 1.279746e-20272 BOOST_MATH_STATIC_LOCAL_VARIABLE const float Y = 0.99941349029541015625f;273 BOOST_MATH_STATIC const T P[] = {274 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000539042911019078575891),275 BOOST_MATH_BIG_CONSTANT(T, 64, -0.28398759004727721098e-6),276 BOOST_MATH_BIG_CONSTANT(T, 64, 0.899465114892291446442e-6),277 BOOST_MATH_BIG_CONSTANT(T, 64, 0.229345859265920864296e-7),278 BOOST_MATH_BIG_CONSTANT(T, 64, 0.225561444863500149219e-9),279 BOOST_MATH_BIG_CONSTANT(T, 64, 0.947846627503022684216e-12),280 BOOST_MATH_BIG_CONSTANT(T, 64, 0.135880130108924861008e-14),281 BOOST_MATH_BIG_CONSTANT(T, 64, -0.348890393399948882918e-21)282 };283 BOOST_MATH_STATIC const T Q[] = {284 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),285 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),286 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),287 BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),288 BOOST_MATH_BIG_CONSTANT(T, 64, 0.399968812193862100054e-6),289 BOOST_MATH_BIG_CONSTANT(T, 64, 0.161809290887904476097e-8),290 BOOST_MATH_BIG_CONSTANT(T, 64, 0.231558608310259605225e-11)291 };292 // LCOV_EXCL_STOP293 T xs = x - 44;294 T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);295 result = Y * x + R * x;296 }297 }298 return result;299}300 301template <class T, class Policy>302struct erf_roots303{304 boost::math::tuple<T,T,T> operator()(const T& guess)305 {306 BOOST_MATH_STD_USING307 T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));308 T derivative2 = -2 * guess * derivative;309 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);310 }311 erf_roots(T z, int s) : target(z), sign(s) {}312private:313 T target;314 int sign;315};316 317template <class T, class Policy>318T erf_inv_imp(const T& p, const T& q, const Policy& pol, const std::integral_constant<int, 0>&)319{320 //321 // Generic version, get a guess that's accurate to 64-bits (10^-19)322 //323 using tag_type = std::integral_constant<int, 64>;324 T guess = erf_inv_imp(p, q, pol, tag_type());325 T result;326 //327 // If T has more bit's than 64 in it's mantissa then we need to iterate,328 // otherwise we can just return the result:329 //330 if(policies::digits<T, Policy>() > 64)331 {332 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();333 if(p <= 0.5)334 {335 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);336 }337 else338 {339 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);340 }341 policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);342 }343 else344 {345 result = guess;346 }347 return result;348}349 350} // namespace detail351 352template <class T, class Policy>353BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)354{355 typedef typename tools::promote_args<T>::type result_type;356 357 //358 // Begin by testing for domain errors, and other special cases:359 //360 constexpr auto function = "boost::math::erfc_inv<%1%>(%1%, %1%)";361 if((z < 0) || (z > 2))362 return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);363 if(z == 0)364 return policies::raise_overflow_error<result_type>(function, nullptr, pol);365 if(z == 2)366 return -policies::raise_overflow_error<result_type>(function, nullptr, pol);367 //368 // Normalise the input, so it's in the range [0,1], we will369 // negate the result if z is outside that range. This is a simple370 // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)371 //372 result_type p, q, s;373 if(z > 1)374 {375 q = 2 - z;376 p = 1 - q;377 s = -1;378 }379 else380 {381 p = 1 - z;382 q = z;383 s = 1;384 }385 //386 // A bit of meta-programming to figure out which implementation387 // to use, based on the number of bits in the mantissa of T:388 //389 typedef typename policies::precision<result_type, Policy>::type precision_type;390 typedef std::integral_constant<int,391 precision_type::value <= 0 ? 0 :392 precision_type::value <= 64 ? 64 : 0393 > tag_type;394 //395 // Likewise use internal promotion, so we evaluate at a higher396 // precision internally if it's appropriate:397 //398 typedef typename policies::evaluation<result_type, Policy>::type eval_type;399 typedef typename policies::normalise<400 Policy,401 policies::promote_float<false>,402 policies::promote_double<false>,403 policies::discrete_quantile<>,404 policies::assert_undefined<> >::type forwarding_policy;405 406 //407 // And get the result, negating where required:408 //409 return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(410 detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), tag_type()), function);411}412 413template <class T, class Policy>414BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)415{416 typedef typename tools::promote_args<T>::type result_type;417 418 //419 // Begin by testing for domain errors, and other special cases:420 //421 constexpr auto function = "boost::math::erf_inv<%1%>(%1%, %1%)";422 if((z < -1) || (z > 1))423 return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);424 if(z == 1)425 return policies::raise_overflow_error<result_type>(function, nullptr, pol);426 if(z == -1)427 return -policies::raise_overflow_error<result_type>(function, nullptr, pol);428 if(z == 0)429 return 0;430 //431 // Normalise the input, so it's in the range [0,1], we will432 // negate the result if z is outside that range. This is a simple433 // application of the erf reflection formula: erf(-z) = -erf(z)434 //435 result_type p, q, s;436 if(z < 0)437 {438 p = -z;439 q = 1 - p;440 s = -1;441 }442 else443 {444 p = z;445 q = 1 - z;446 s = 1;447 }448 //449 // A bit of meta-programming to figure out which implementation450 // to use, based on the number of bits in the mantissa of T:451 //452 typedef typename policies::precision<result_type, Policy>::type precision_type;453 typedef std::integral_constant<int,454 precision_type::value <= 0 ? 0 :455 precision_type::value <= 64 ? 64 : 0456 > tag_type;457 //458 // Likewise use internal promotion, so we evaluate at a higher459 // precision internally if it's appropriate:460 //461 typedef typename policies::evaluation<result_type, Policy>::type eval_type;462 typedef typename policies::normalise<463 Policy,464 policies::promote_float<false>,465 policies::promote_double<false>,466 policies::discrete_quantile<>,467 policies::assert_undefined<> >::type forwarding_policy;468 //469 // Likewise use internal promotion, so we evaluate at a higher470 // precision internally if it's appropriate:471 //472 typedef typename policies::evaluation<result_type, Policy>::type eval_type;473 474 //475 // And get the result, negating where required:476 //477 return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(478 detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), tag_type()), function);479}480 481template <class T>482BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type erfc_inv(T z)483{484 return erfc_inv(z, policies::policy<>());485}486 487template <class T>488BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type erf_inv(T z)489{490 return erf_inv(z, policies::policy<>());491}492 493} // namespace math494} // namespace boost495 496#else // Special handling for NVRTC497 498namespace boost {499namespace math {500 501template <typename T>502BOOST_MATH_GPU_ENABLED auto erf_inv(T x)503{504 return ::erfinv(x);505}506 507template <>508BOOST_MATH_GPU_ENABLED auto erf_inv(float x)509{510 return ::erfinvf(x);511}512 513template <typename T, typename Policy>514BOOST_MATH_GPU_ENABLED auto erf_inv(T x, const Policy&)515{516 return ::erfinv(x);517}518 519template <typename Policy>520BOOST_MATH_GPU_ENABLED auto erf_inv(float x, const Policy&)521{522 return ::erfinvf(x);523}524 525template <typename T>526BOOST_MATH_GPU_ENABLED auto erfc_inv(T x)527{528 return ::erfcinv(x);529}530 531template <>532BOOST_MATH_GPU_ENABLED auto erfc_inv(float x)533{534 return ::erfcinvf(x);535}536 537template <typename T, typename Policy>538BOOST_MATH_GPU_ENABLED auto erfc_inv(T x, const Policy&)539{540 return ::erfcinv(x);541}542 543template <typename Policy>544BOOST_MATH_GPU_ENABLED auto erfc_inv(float x, const Policy&)545{546 return ::erfcinvf(x);547}548 549} // namespace math550} // namespace boost551 552#endif // BOOST_MATH_HAS_NVRTV553 554#ifdef _MSC_VER555#pragma warning(pop)556#endif557 558#endif // BOOST_MATH_SF_ERF_INV_HPP559 560