ceph/src/boost/boost/math/special_functions/detail/erf_inv.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SF_ERF_INV_HPP
   7 #define BOOST_MATH_SF_ERF_INV_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #pragma warning(push)
  12 #pragma warning(disable:4127) // Conditional expression is constant
  13 #pragma warning(disable:4702) // Unreachable code: optimization warning
  14 #endif
  15
  16 #include <type_traits>
  17
  18 namespace boost{ namespace math{
  19
  20 namespace detail{
  21 //
  22 // The inverse erf and erfc functions share a common implementation,
  23 // this version is for 80-bit long double's and smaller:
  24 //
  25 template <class T, class Policy>
  26 T erf_inv_imp(const T& p, const T& q, const Policy&, const std::integral_constant<int, 64>*)
  27 {
  28    BOOST_MATH_STD_USING // for ADL of std names.
  29
  30    T result = 0;
  31
  32    if(p <= 0.5)
  33    {
  34       //
  35       // Evaluate inverse erf using the rational approximation:
  36       //
  37       // x = p(p+10)(Y+R(p))
  38       //
  39       // Where Y is a constant, and R(p) is optimised for a low
  40       // absolute error compared to |Y|.
  41       //
  42       // double: Max error found: 2.001849e-18
  43       // long double: Max error found: 1.017064e-20
  44       // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
  45       //
  46       static const float Y = 0.0891314744949340820313f;
  47       static const T P[] = {
  48          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000508781949658280665617),
  49          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00836874819741736770379),
  50          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0334806625409744615033),
  51          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0126926147662974029034),
  52          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0365637971411762664006),
  53          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0219878681111168899165),
  54          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00822687874676915743155),
  55          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00538772965071242932965)
  56       };
  57       static const T Q[] = {
  58          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
  59          BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),
  60          BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),
  61          BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),
  62          BOOST_MATH_BIG_CONSTANT(T, 64, 0.662328840472002992063),
  63          BOOST_MATH_BIG_CONSTANT(T, 64, -0.71228902341542847553),
  64          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0527396382340099713954),
  65          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0795283687341571680018),
  66          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00233393759374190016776),
  67          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000886216390456424707504)
  68       };
  69       T g = p * (p + 10);
  70       T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
  71       result = g * Y + g * r;
  72    }
  73    else if(q >= 0.25)
  74    {
  75       //
  76       // Rational approximation for 0.5 > q >= 0.25
  77       //
  78       // x = sqrt(-2*log(q)) / (Y + R(q))
  79       //
  80       // Where Y is a constant, and R(q) is optimised for a low
  81       // absolute error compared to Y.
  82       //
  83       // double : Max error found: 7.403372e-17
  84       // long double : Max error found: 6.084616e-20
  85       // Maximum Deviation Found (error term) 4.811e-20
  86       //
  87       static const float Y = 2.249481201171875f;
  88       static const T P[] = {
  89          BOOST_MATH_BIG_CONSTANT(T, 64, -0.202433508355938759655),
  90          BOOST_MATH_BIG_CONSTANT(T, 64, 0.105264680699391713268),
  91          BOOST_MATH_BIG_CONSTANT(T, 64, 8.37050328343119927838),
  92          BOOST_MATH_BIG_CONSTANT(T, 64, 17.6447298408374015486),
  93          BOOST_MATH_BIG_CONSTANT(T, 64, -18.8510648058714251895),
  94          BOOST_MATH_BIG_CONSTANT(T, 64, -44.6382324441786960818),
  95          BOOST_MATH_BIG_CONSTANT(T, 64, 17.445385985570866523),
  96          BOOST_MATH_BIG_CONSTANT(T, 64, 21.1294655448340526258),
  97          BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)
  98       };
  99       static const T Q[] = {
 100          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 101          BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),
 102          BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),
 103          BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),
 104          BOOST_MATH_BIG_CONSTANT(T, 64, -20.1432634680485188801),
 105          BOOST_MATH_BIG_CONSTANT(T, 64, 48.5609213108739935468),
 106          BOOST_MATH_BIG_CONSTANT(T, 64, 10.8268667355460159008),
 107          BOOST_MATH_BIG_CONSTANT(T, 64, -22.6436933413139721736),
 108          BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724)
 109       };
 110       T g = sqrt(-2 * log(q));
 111       T xs = q - 0.25f;
 112       T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 113       result = g / (Y + r);
 114    }
 115    else
 116    {
 117       //
 118       // For q < 0.25 we have a series of rational approximations all
 119       // of the general form:
 120       //
 121       // let: x = sqrt(-log(q))
 122       //
 123       // Then the result is given by:
 124       //
 125       // x(Y+R(x-B))
 126       //
 127       // where Y is a constant, B is the lowest value of x for which
 128       // the approximation is valid, and R(x-B) is optimised for a low
 129       // absolute error compared to Y.
 130       //
 131       // Note that almost all code will really go through the first
 132       // or maybe second approximation.  After than we're dealing with very
 133       // small input values indeed: 80 and 128 bit long double's go all the
 134       // way down to ~ 1e-5000 so the "tail" is rather long...
 135       //
 136       T x = sqrt(-log(q));
 137       if(x < 3)
 138       {
 139          // Max error found: 1.089051e-20
 140          static const float Y = 0.807220458984375f;
 141          static const T P[] = {
 142             BOOST_MATH_BIG_CONSTANT(T, 64, -0.131102781679951906451),
 143             BOOST_MATH_BIG_CONSTANT(T, 64, -0.163794047193317060787),
 144             BOOST_MATH_BIG_CONSTANT(T, 64, 0.117030156341995252019),
 145             BOOST_MATH_BIG_CONSTANT(T, 64, 0.387079738972604337464),
 146             BOOST_MATH_BIG_CONSTANT(T, 64, 0.337785538912035898924),
 147             BOOST_MATH_BIG_CONSTANT(T, 64, 0.142869534408157156766),
 148             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0290157910005329060432),
 149             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00214558995388805277169),
 150             BOOST_MATH_BIG_CONSTANT(T, 64, -0.679465575181126350155e-6),
 151             BOOST_MATH_BIG_CONSTANT(T, 64, 0.285225331782217055858e-7),
 152             BOOST_MATH_BIG_CONSTANT(T, 64, -0.681149956853776992068e-9)
 153          };
 154          static const T Q[] = {
 155             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 156             BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),
 157             BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),
 158             BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),
 159             BOOST_MATH_BIG_CONSTANT(T, 64, 2.59301921623620271374),
 160             BOOST_MATH_BIG_CONSTANT(T, 64, 0.848854343457902036425),
 161             BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612),
 162             BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121)
 163          };
 164          T xs = x - 1.125f;
 165          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 166          result = Y * x + R * x;
 167       }
 168       else if(x < 6)
 169       {
 170          // Max error found: 8.389174e-21
 171          static const float Y = 0.93995571136474609375f;
 172          static const T P[] = {
 173             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0350353787183177984712),
 174             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00222426529213447927281),
 175             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0185573306514231072324),
 176             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00950804701325919603619),
 177             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00187123492819559223345),
 178             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000157544617424960554631),
 179             BOOST_MATH_BIG_CONSTANT(T, 64, 0.460469890584317994083e-5),
 180             BOOST_MATH_BIG_CONSTANT(T, 64, -0.230404776911882601748e-9),
 181             BOOST_MATH_BIG_CONSTANT(T, 64, 0.266339227425782031962e-11)
 182          };
 183          static const T Q[] = {
 184             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 185             BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),
 186             BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),
 187             BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),
 188             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0341589143670947727934),
 189             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00263861676657015992959),
 190             BOOST_MATH_BIG_CONSTANT(T, 64, 0.764675292302794483503e-4)
 191          };
 192          T xs = x - 3;
 193          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 194          result = Y * x + R * x;
 195       }
 196       else if(x < 18)
 197       {
 198          // Max error found: 1.481312e-19
 199          static const float Y = 0.98362827301025390625f;
 200          static const T P[] = {
 201             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0167431005076633737133),
 202             BOOST_MATH_BIG_CONSTANT(T, 64, -0.00112951438745580278863),
 203             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00105628862152492910091),
 204             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000209386317487588078668),
 205             BOOST_MATH_BIG_CONSTANT(T, 64, 0.149624783758342370182e-4),
 206             BOOST_MATH_BIG_CONSTANT(T, 64, 0.449696789927706453732e-6),
 207             BOOST_MATH_BIG_CONSTANT(T, 64, 0.462596163522878599135e-8),
 208             BOOST_MATH_BIG_CONSTANT(T, 64, -0.281128735628831791805e-13),
 209             BOOST_MATH_BIG_CONSTANT(T, 64, 0.99055709973310326855e-16)
 210          };
 211          static const T Q[] = {
 212             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 213             BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),
 214             BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),
 215             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),
 216             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000964011807005165528527),
 217             BOOST_MATH_BIG_CONSTANT(T, 64, 0.275335474764726041141e-4),
 218             BOOST_MATH_BIG_CONSTANT(T, 64, 0.282243172016108031869e-6)
 219          };
 220          T xs = x - 6;
 221          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 222          result = Y * x + R * x;
 223       }
 224       else if(x < 44)
 225       {
 226          // Max error found: 5.697761e-20
 227          static const float Y = 0.99714565277099609375f;
 228          static const T P[] = {
 229             BOOST_MATH_BIG_CONSTANT(T, 64, -0.0024978212791898131227),
 230             BOOST_MATH_BIG_CONSTANT(T, 64, -0.779190719229053954292e-5),
 231             BOOST_MATH_BIG_CONSTANT(T, 64, 0.254723037413027451751e-4),
 232             BOOST_MATH_BIG_CONSTANT(T, 64, 0.162397777342510920873e-5),
 233             BOOST_MATH_BIG_CONSTANT(T, 64, 0.396341011304801168516e-7),
 234             BOOST_MATH_BIG_CONSTANT(T, 64, 0.411632831190944208473e-9),
 235             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145596286718675035587e-11),
 236             BOOST_MATH_BIG_CONSTANT(T, 64, -0.116765012397184275695e-17)
 237          };
 238          static const T Q[] = {
 239             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 240             BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),
 241             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),
 242             BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),
 243             BOOST_MATH_BIG_CONSTANT(T, 64, 0.145007359818232637924e-4),
 244             BOOST_MATH_BIG_CONSTANT(T, 64, 0.144437756628144157666e-6),
 245             BOOST_MATH_BIG_CONSTANT(T, 64, 0.509761276599778486139e-9)
 246          };
 247          T xs = x - 18;
 248          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 249          result = Y * x + R * x;
 250       }
 251       else
 252       {
 253          // Max error found: 1.279746e-20
 254          static const float Y = 0.99941349029541015625f;
 255          static const T P[] = {
 256             BOOST_MATH_BIG_CONSTANT(T, 64, -0.000539042911019078575891),
 257             BOOST_MATH_BIG_CONSTANT(T, 64, -0.28398759004727721098e-6),
 258             BOOST_MATH_BIG_CONSTANT(T, 64, 0.899465114892291446442e-6),
 259             BOOST_MATH_BIG_CONSTANT(T, 64, 0.229345859265920864296e-7),
 260             BOOST_MATH_BIG_CONSTANT(T, 64, 0.225561444863500149219e-9),
 261             BOOST_MATH_BIG_CONSTANT(T, 64, 0.947846627503022684216e-12),
 262             BOOST_MATH_BIG_CONSTANT(T, 64, 0.135880130108924861008e-14),
 263             BOOST_MATH_BIG_CONSTANT(T, 64, -0.348890393399948882918e-21)
 264          };
 265          static const T Q[] = {
 266             BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 267             BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),
 268             BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),
 269             BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),
 270             BOOST_MATH_BIG_CONSTANT(T, 64, 0.399968812193862100054e-6),
 271             BOOST_MATH_BIG_CONSTANT(T, 64, 0.161809290887904476097e-8),
 272             BOOST_MATH_BIG_CONSTANT(T, 64, 0.231558608310259605225e-11)
 273          };
 274          T xs = x - 44;
 275          T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
 276          result = Y * x + R * x;
 277       }
 278    }
 279    return result;
 280 }
 281
 282 template <class T, class Policy>
 283 struct erf_roots
 284 {
 285    boost::math::tuple<T,T,T> operator()(const T& guess)
 286    {
 287       BOOST_MATH_STD_USING
 288       T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
 289       T derivative2 = -2 * guess * derivative;
 290       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);
 291    }
 292    erf_roots(T z, int s) : target(z), sign(s) {}
 293 private:
 294    T target;
 295    int sign;
 296 };
 297
 298 template <class T, class Policy>
 299 T erf_inv_imp(const T& p, const T& q, const Policy& pol, const std::integral_constant<int, 0>*)
 300 {
 301    //
 302    // Generic version, get a guess that's accurate to 64-bits (10^-19)
 303    //
 304    T guess = erf_inv_imp(p, q, pol, static_cast<std::integral_constant<int, 64> const*>(0));
 305    T result;
 306    //
 307    // If T has more bit's than 64 in it's mantissa then we need to iterate,
 308    // otherwise we can just return the result:
 309    //
 310    if(policies::digits<T, Policy>() > 64)
 311    {
 312       std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 313       if(p <= 0.5)
 314       {
 315          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);
 316       }
 317       else
 318       {
 319          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);
 320       }
 321       policies::check_root_iterations<T>("boost::math::erf_inv<%1%>", max_iter, pol);
 322    }
 323    else
 324    {
 325       result = guess;
 326    }
 327    return result;
 328 }
 329
 330 template <class T, class Policy>
 331 struct erf_inv_initializer
 332 {
 333    struct init
 334    {
 335       init()
 336       {
 337          do_init();
 338       }
 339       static bool is_value_non_zero(T);
 340       static void do_init()
 341       {
 342          // If std::numeric_limits<T>::digits is zero, we must not call
 343          // our initialization code here as the precision presumably
 344          // varies at runtime, and will not have been set yet.
 345          if(std::numeric_limits<T>::digits)
 346          {
 347             boost::math::erf_inv(static_cast<T>(0.25), Policy());
 348             boost::math::erf_inv(static_cast<T>(0.55), Policy());
 349             boost::math::erf_inv(static_cast<T>(0.95), Policy());
 350             boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
 351             // These following initializations must not be called if
 352             // type T can not hold the relevant values without
 353             // underflow to zero.  We check this at runtime because
 354             // some tools such as valgrind silently change the precision
 355             // of T at runtime, and numeric_limits basically lies!
 356             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130))))
 357                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());
 358
 359             // Some compilers choke on constants that would underflow, even in code that isn't instantiated
 360             // so try and filter these cases out in the preprocessor:
 361 #if LDBL_MAX_10_EXP >= 800
 362             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800))))
 363                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
 364             if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900))))
 365                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
 366 #else
 367             if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800))))
 368                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800)), Policy());
 369             if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900))))
 370                boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900)), Policy());
 371 #endif
 372          }
 373       }
 374       void force_instantiate()const{}
 375    };
 376    static const init initializer;
 377    static void force_instantiate()
 378    {
 379       initializer.force_instantiate();
 380    }
 381 };
 382
 383 template <class T, class Policy>
 384 const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer;
 385
 386 template <class T, class Policy>
 387 bool erf_inv_initializer<T, Policy>::init::is_value_non_zero(T v)
 388 {
 389    // This needs to be non-inline to detect whether v is non zero at runtime
 390    // rather than at compile time, only relevant when running under valgrind
 391    // which changes long double's to double's on the fly.
 392    return v != 0;
 393 }
 394
 395 } // namespace detail
 396
 397 template <class T, class Policy>
 398 typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
 399 {
 400    typedef typename tools::promote_args<T>::type result_type;
 401
 402    //
 403    // Begin by testing for domain errors, and other special cases:
 404    //
 405    static const char* function = "boost::math::erfc_inv<%1%>(%1%, %1%)";
 406    if((z < 0) || (z > 2))
 407       return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
 408    if(z == 0)
 409       return policies::raise_overflow_error<result_type>(function, 0, pol);
 410    if(z == 2)
 411       return -policies::raise_overflow_error<result_type>(function, 0, pol);
 412    //
 413    // Normalise the input, so it's in the range [0,1], we will
 414    // negate the result if z is outside that range.  This is a simple
 415    // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
 416    //
 417    result_type p, q, s;
 418    if(z > 1)
 419    {
 420       q = 2 - z;
 421       p = 1 - q;
 422       s = -1;
 423    }
 424    else
 425    {
 426       p = 1 - z;
 427       q = z;
 428       s = 1;
 429    }
 430    //
 431    // A bit of meta-programming to figure out which implementation
 432    // to use, based on the number of bits in the mantissa of T:
 433    //
 434    typedef typename policies::precision<result_type, Policy>::type precision_type;
 435    typedef std::integral_constant<int,
 436       precision_type::value <= 0 ? 0 :
 437       precision_type::value <= 64 ? 64 : 0
 438    > tag_type;
 439    //
 440    // Likewise use internal promotion, so we evaluate at a higher
 441    // precision internally if it's appropriate:
 442    //
 443    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
 444    typedef typename policies::normalise<
 445       Policy,
 446       policies::promote_float<false>,
 447       policies::promote_double<false>,
 448       policies::discrete_quantile<>,
 449       policies::assert_undefined<> >::type forwarding_policy;
 450
 451    detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
 452
 453    //
 454    // And get the result, negating where required:
 455    //
 456    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
 457       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(0)), function);
 458 }
 459
 460 template <class T, class Policy>
 461 typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
 462 {
 463    typedef typename tools::promote_args<T>::type result_type;
 464
 465    //
 466    // Begin by testing for domain errors, and other special cases:
 467    //
 468    static const char* function = "boost::math::erf_inv<%1%>(%1%, %1%)";
 469    if((z < -1) || (z > 1))
 470       return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
 471    if(z == 1)
 472       return policies::raise_overflow_error<result_type>(function, 0, pol);
 473    if(z == -1)
 474       return -policies::raise_overflow_error<result_type>(function, 0, pol);
 475    if(z == 0)
 476       return 0;
 477    //
 478    // Normalise the input, so it's in the range [0,1], we will
 479    // negate the result if z is outside that range.  This is a simple
 480    // application of the erf reflection formula: erf(-z) = -erf(z)
 481    //
 482    result_type p, q, s;
 483    if(z < 0)
 484    {
 485       p = -z;
 486       q = 1 - p;
 487       s = -1;
 488    }
 489    else
 490    {
 491       p = z;
 492       q = 1 - z;
 493       s = 1;
 494    }
 495    //
 496    // A bit of meta-programming to figure out which implementation
 497    // to use, based on the number of bits in the mantissa of T:
 498    //
 499    typedef typename policies::precision<result_type, Policy>::type precision_type;
 500    typedef std::integral_constant<int,
 501       precision_type::value <= 0 ? 0 :
 502       precision_type::value <= 64 ? 64 : 0
 503    > tag_type;
 504    //
 505    // Likewise use internal promotion, so we evaluate at a higher
 506    // precision internally if it's appropriate:
 507    //
 508    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
 509    typedef typename policies::normalise<
 510       Policy,
 511       policies::promote_float<false>,
 512       policies::promote_double<false>,
 513       policies::discrete_quantile<>,
 514       policies::assert_undefined<> >::type forwarding_policy;
 515    //
 516    // Likewise use internal promotion, so we evaluate at a higher
 517    // precision internally if it's appropriate:
 518    //
 519    typedef typename policies::evaluation<result_type, Policy>::type eval_type;
 520
 521    detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate();
 522    //
 523    // And get the result, negating where required:
 524    //
 525    return s * policies::checked_narrowing_cast<result_type, forwarding_policy>(
 526       detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), forwarding_policy(), static_cast<tag_type const*>(0)), function);
 527 }
 528
 529 template <class T>
 530 inline typename tools::promote_args<T>::type erfc_inv(T z)
 531 {
 532    return erfc_inv(z, policies::policy<>());
 533 }
 534
 535 template <class T>
 536 inline typename tools::promote_args<T>::type erf_inv(T z)
 537 {
 538    return erf_inv(z, policies::policy<>());
 539 }
 540
 541 } // namespace math
 542 } // namespace boost
 543
 544 #ifdef _MSC_VER
 545 #pragma warning(pop)
 546 #endif
 547
 548 #endif // BOOST_MATH_SF_ERF_INV_HPP
 549