histogram/detail/erf_inv.hpp

0001 // Copyright 2022 Hans Dembinski, Jay Gohil
0002 //
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_DETAIL_ERF_INF_HPP
0008 #define BOOST_HISTOGRAM_DETAIL_ERF_INF_HPP
0009
0010 #include <cmath>
0011
0012 namespace boost {
0013 namespace histogram {
0014 namespace detail {
0015
0016 // Simple implementation of erf_inv so that we do not depend on boost::math.
0017 // If you happen to discover this, prefer the boost::math implementation,
0018 // it is more accurate for x very close to -1 or 1 and faster.
0019 // The only virtue of this implementation is its simplicity.
0020 template <int Iterate = 3>
0021 double erf_inv(double x) noexcept {
0022   // Strategy: solve f(y) = x - erf(y) = 0 for given x with Newton's method.
0023   // f'(y) = -erf'(y) = -2/sqrt(pi) e^(-y^2)
0024   // Has quadratic convergence. Since erf_inv<0> is accurate to 1e-3,
0025   // we should have machine precision after three iterations.
0026   const double x0 = erf_inv<Iterate - 1>(x); // recursion
0027   const double fx0 = x - std::erf(x0);
0028   const double pi = std::acos(-1);
0029   double fpx0 = -2.0 / std::sqrt(pi) * std::exp(-x0 * x0);
0030   return x0 - fx0 / fpx0; // = x1
0031 }
0032
0033 template <>
0034 inline double erf_inv<0>(double x) noexcept {
0035   // Specialization to get initial estimate.
0036   // This formula is accurate to about 1e-3.
0037   // Based on https://stackoverflow.com/questions/27229371/inverse-error-function-in-c
0038   const double a = std::log((1 - x) * (1 + x));
0039   const double b = std::fma(0.5, a, 4.120666747961526);
0040   const double c = 6.47272819164 * a;
0041   return std::copysign(std::sqrt(-b + std::sqrt(std::fma(b, b, -c))), x);
0042 }
0043
0044 } // namespace detail
0045 } // namespace histogram
0046 } // namespace boost
0047
0048 #endif