Merge branch 'develop' of https://github.com/boostorg/math into develop

Author: jzmaddock

doc/sf/erf.qbk

@@ -66,14 +66,29 @@ than the one shown will have __zero_error.
 
 [table_erfc]
 
-The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 
-and GCC-7.1/Ubuntu for `long double` and `__float128`.
+The following error plots are based on an exhaustive search of the functions domain, MSVC-16.7.1 at `double` precision, 
+and GCC-10/Mingw64 for `long double` and `__float128`.
 
-[graph erf__double]
+[$../../reporting/accuracy/erf_errors_double.png]
 
-[graph erf__80_bit_long_double]
+[$../../reporting/accuracy/erf_errors_long_double.png]
+
+[$../../reporting/accuracy/erf_errors_float128.png]
+
+In the erfc case, error rates are almost entirely the error in calculating `exp(-x*x)`:
+
+[$../../reporting/accuracy/erfc_errors_double.png]
+
+[$../../reporting/accuracy/erfc_errors_long_double.png]
+
+[$../../reporting/accuracy/erfc_errors_float128.png]
+
+Multiprecision error rates are similar, albeit with a much larger error in calculating the exponent term for erfc:
+
+[$../../reporting/accuracy/erf_errors_cpp_bin_float_50.png]
+
+[$../../reporting/accuracy/erfc_errors_cpp_bin_float_50.png]
 
-[graph erf____float128]
 
 [h4 Testing]
 

include/boost/math/special_functions/erf.hpp

@@ -105,6 +105,37 @@ inline T erf_asymptotic_limit()
    return erf_asymptotic_limit_N(tag_type());
 }
 
+template <class T>
+struct erf_series_near_zero
+{
+   typedef T result_type;
+   T         term;
+   T         zz;
+   int       k;
+   erf_series_near_zero(const T& z) : term(z), zz(-z * z), k(0) {}
+
+   T operator()()
+   {
+      T result = term / (2 * k + 1);
+      term *= zz / ++k;
+      return result;
+   }
+};
+
+template <class T, class Policy>
+T erf_series_near_zero_sum(const T& x, const Policy& pol)
+{
+   //
+   // We need Kahan summation here, otherwise the errors grow fairly quickly.
+   // This method is *much* faster than the alternatives even so.
+   //
+   erf_series_near_zero<T> sum(x);
+   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+   T result = constants::two_div_root_pi<T>() * tools::kahan_sum_series(sum, tools::digits<T>(), max_iter);
+   policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
+   return result;
+}
+
 template <class T, class Policy, class Tag>
 T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
 {
@@ -132,20 +163,12 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
    else
    {
       T x = z * z;
-      if(x < 0.6)
+      if(z < 1.3f)
       {
          // Compute P:
-         result = z * exp(-x);
-         result /= sqrt(boost::math::constants::pi<T>());
-         if(result != 0)
-            result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol);
-      }
-      else if(x < 1.1f)
-      {
-         // Compute Q:
-         invert = !invert;
-         result = tgamma_small_upper_part(T(0.5f), x, pol);
-         result /= sqrt(boost::math::constants::pi<T>());
+         // This is actually good for z p to 2 or so, but the cutoff given seems
+         // to be the best compromise.  Performance wise, this is way quicker than anything else...
+         result = erf_series_near_zero_sum(z, pol);
       }
       else if(x > 1 / tools::epsilon<T>())
       {

reporting/accuracy/Jamfile.v2

@@ -159,3 +159,25 @@ boostbook standalone
         <xsl:param>generate.section.toc.level=10
     ;
 
+lib gmp ;
+lib mpfr ;
+lib quadmath ;
+#
+# Some manual tests that are expensive to run:
+#
+run erf_error_plot.cpp mpfr gmp : : : release <cxxstd>17 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erf_error_plot_double ;
+explicit erf_error_plot_double ;
+run erf_error_plot.cpp mpfr gmp : : : release <cxxstd>17 <define>TEST_TYPE="\"long double\"" [ check-target-builds ../../config//has_mpfr : : <build>no ] : erf_error_plot_long_double ;
+explicit erf_error_plot_long_double ;
+run erf_error_plot.cpp mpfr gmp : : : release <cxxstd>17 <define>TEST_TYPE=cpp_bin_float_50 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erf_error_plot_cpp_bin_float_50 ;
+explicit erf_error_plot_cpp_bin_float_50 ;
+run erf_error_plot.cpp mpfr gmp quadmath : : : release <cxxstd>17 <cxxstd-dialect>gnu <define>TEST_TYPE=float128 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erf_error_plot_float128 ;
+explicit erf_error_plot_cpp_bin_float_50 ;
+run erfc_error_plot.cpp mpfr gmp : : : release <cxxstd>17 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erfc_error_plot_double ;
+explicit erfc_error_plot_double ;
+run erfc_error_plot.cpp mpfr gmp : : : release <cxxstd>17 <define>TEST_TYPE="\"long double\"" [ check-target-builds ../../config//has_mpfr : : <build>no ] : erfc_error_plot_long_double ;
+explicit erfc_error_plot_long_double ;
+run erfc_error_plot.cpp mpfr gmp : : : release <cxxstd>17 <define>TEST_TYPE=cpp_bin_float_50 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erfc_error_plot_cpp_bin_float_50 ;
+explicit erfc_error_plot_cpp_bin_float_50 ;
+run erfc_error_plot.cpp mpfr gmp quadmath : : : release <cxxstd>17 <cxxstd-dialect>gnu <define>TEST_TYPE=float128 [ check-target-builds ../../config//has_mpfr : : <build>no ] : erfc_error_plot_float128 ;
+explicit erfc_error_plot_cpp_bin_float_50 ;

reporting/accuracy/erf_error_plot.cpp

@@ -0,0 +1,75 @@
+
+//  (C) Copyright Nick Thompson 2020.
+//  (C) Copyright John Maddock 2020.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <iostream>
+#include <boost/math/tools/ulps_plot.hpp>
+#include <boost/core/demangle.hpp>
+#include <boost/multiprecision/mpfr.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+#endif
+
+using namespace boost::multiprecision;
+
+#ifndef TEST_TYPE
+#define TEST_TYPE double
+#endif
+
+std::string test_type_name(BOOST_STRINGIZE(TEST_TYPE));
+std::string test_type_filename(BOOST_STRINGIZE(TEST_TYPE));
+
+using boost::math::tools::ulps_plot;
+
+int main() 
+{
+   std::string::size_type n;
+   while ((n = test_type_filename.find_first_not_of("_qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM1234567890")) != std::string::npos)
+   {
+      test_type_filename[n] = '_';
+   }
+
+   using PreciseReal = boost::multiprecision::mpfr_float_100;
+   using CoarseReal = TEST_TYPE;
+
+   typedef boost::math::policies::policy<
+      boost::math::policies::promote_float<false>,
+      boost::math::policies::promote_double<false> >
+      no_promote_policy;
+
+   auto ai_coarse = [](CoarseReal const& x)->CoarseReal {
+      return erf(x);
+   };
+   auto ai_precise = [](PreciseReal const& x)->PreciseReal {
+      return erf(x);
+   };
+
+   std::string filename = "erf_errors_";
+   filename += test_type_filename;
+   filename += ".svg";
+   int samples = 100000;
+   // How many pixels wide do you want your .svg?
+   int width = 700;
+   // Near a root, we have unbounded relative error. So for functions with roots, we define an ULP clip:
+   PreciseReal clip = 40;
+   // Should we perturb the abscissas? i.e., should we compute the high precision function f at x,
+   // and the low precision function at the nearest representable x̂ to x?
+   // Or should we compute both the high precision and low precision function at a low precision representable x̂?
+   bool perturb_abscissas = false;
+   auto plot = ulps_plot<decltype(ai_precise), PreciseReal, CoarseReal>(ai_precise, CoarseReal(-10), CoarseReal(10), samples, perturb_abscissas);
+   // Note the argument chaining:
+   plot.clip(clip).width(width);
+   plot.background_color("white").font_color("black");
+   // Sometimes it's useful to set a title, but in many cases it's more useful to just use a caption.
+   std::string title = "Erf ULP plot at " + test_type_name + " precision";
+   plot.title(title);
+   plot.vertical_lines(6);
+   plot.add_fn(ai_coarse);
+   // You can write the plot to a stream:
+   //std::cout << plot;
+   // Or to a file:
+   plot.write(filename);
+}

reporting/accuracy/erf_errors_cpp_bin_float_50.png

@@ -0,0 +1,75 @@
+
+//  (C) Copyright Nick Thompson 2020.
+//  (C) Copyright John Maddock 2020.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <iostream>
+#include <boost/math/tools/ulps_plot.hpp>
+#include <boost/core/demangle.hpp>
+#include <boost/multiprecision/mpfr.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+#endif
+
+using namespace boost::multiprecision;
+
+#ifndef TEST_TYPE
+#define TEST_TYPE cpp_bin_float_50
+#endif
+
+std::string test_type_name(BOOST_STRINGIZE(TEST_TYPE));
+std::string test_type_filename(BOOST_STRINGIZE(TEST_TYPE));
+
+using boost::math::tools::ulps_plot;
+
+int main() 
+{
+   std::string::size_type n;
+   while ((n = test_type_filename.find_first_not_of("_qwertyuiopasdfghjklzxcvbnmQWERTYUIOPASDFGHJKLZXCVBNM1234567890")) != std::string::npos)
+   {
+      test_type_filename[n] = '_';
+   }
+
+   using PreciseReal = boost::multiprecision::mpfr_float_100;
+   using CoarseReal = TEST_TYPE;
+
+   typedef boost::math::policies::policy<
+      boost::math::policies::promote_float<false>,
+      boost::math::policies::promote_double<false> >
+      no_promote_policy;
+
+   auto ai_coarse = [](CoarseReal const& x)->CoarseReal {
+      return erfc(x);
+   };
+   auto ai_precise = [](PreciseReal const& x)->PreciseReal {
+      return erfc(x);
+   };
+
+   std::string filename = "erfc_errors_";
+   filename += test_type_filename;
+   filename += ".svg";
+   int samples = 100000;
+   // How many pixels wide do you want your .svg?
+   int width = 700;
+   // Near a root, we have unbounded relative error. So for functions with roots, we define an ULP clip:
+   PreciseReal clip = 40;
+   // Should we perturb the abscissas? i.e., should we compute the high precision function f at x,
+   // and the low precision function at the nearest representable x̂ to x?
+   // Or should we compute both the high precision and low precision function at a low precision representable x̂?
+   bool perturb_abscissas = false;
+   auto plot = ulps_plot<decltype(ai_precise), PreciseReal, CoarseReal>(ai_precise, CoarseReal(-10), CoarseReal(30), samples, perturb_abscissas);
+   // Note the argument chaining:
+   plot.clip(clip).width(width);
+   plot.background_color("white").font_color("black");
+   // Sometimes it's useful to set a title, but in many cases it's more useful to just use a caption.
+   std::string title = "Erfc ULP plot at " + test_type_name + " precision";
+   plot.title(title);
+   plot.vertical_lines(6);
+   plot.add_fn(ai_coarse);
+   // You can write the plot to a stream:
+   //std::cout << plot;
+   // Or to a file:
+   plot.write(filename);
+}