std::erf versus erf_impl (Abramowitz & Stegun), C++ (gcc)

std::erf versus erf_impl (Abramowitz & Stegun)

#include <chrono>
#include <cstdio>
#include <cmath>
#include <limits>

template <typename T>
inline T erf_impl(T v)
{
   const T t = T(1) / (T(1) + T(0.5) * std::abs(v));

static const T c[] = {
                          T( 1.26551223), T(1.00002368),
                          T( 0.37409196), T(0.09678418),
                          T(-0.18628806), T(0.27886807),
                          T(-1.13520398), T(1.48851587),
                          T(-0.82215223), T(0.17087277)
                        };

T result = T(1) - t * std::exp((-v * v) -
                          c[0] + t * (c[1] + t *
                         (c[2] + t * (c[3] + t *
                         (c[4] + t * (c[5] + t *
                         (c[6] + t * (c[7] + t *
                         (c[8] + t * (c[9]))))))))));

return (v >= T(0)) ? result : -result;
}

template <typename T> const char* type_str();
template<> const char* type_str<double>(){ static const char* s = "double"; return s; };
template<> const char* type_str<float> (){ static const char* s = "float "; return s; };

template <typename T>
void test_precision()
{
   T upper_bound = T(+3);
   T max_error   = T( 0);
   T total_error = T( 0);
   T i           = T(-3);
   T delta       = T(0.000001);

std::size_t cnt = 0;

while (i <= upper_bound)
   {
      T v0 = erf_impl(i);
      T v1 = std::erf(i);
      T curr_error = std::abs(v0 - v1);

if (curr_error > max_error)
        max_error = curr_error;

total_error += curr_error;
      ++cnt;

i += delta;
   }

printf("precision - type: %s max error: %20.17f\tavg error: %20.17f\teps:%20.17f\n",
          type_str<T>(),
          max_error,
          total_error / cnt,
          std::numeric_limits<T>::epsilon());
}

template <typename T>
void test_time()
{
   {
      T upper_bound = T(+3);
      T v_total     = T( 0);
      T i           = T(-3);
      T delta       = T(0.000001);

const auto start = std::chrono::high_resolution_clock::now();

while (i <= upper_bound)
      {
         v_total += erf_impl(i);
         i += delta;
      }

const auto end = std::chrono::high_resolution_clock::now();
      const auto total = std::chrono::duration_cast<std::chrono::duration<double> >(end - start).count();

printf("time - type: %s erf_impl Time: %7.5fsec\tv:%15.10f\n",type_str<T>(), total,v_total);
   }

{
      T upper_bound = T(+3);
      T v_total     = T( 0);
      T i           = T(-3);
      T delta       = T(0.000001);

const auto start = std::chrono::high_resolution_clock::now();

while (i <= upper_bound)
      {
         v_total += std::erf(i);
         i += delta;
      }

const auto end = std::chrono::high_resolution_clock::now();
      const auto total = std::chrono::duration_cast<std::chrono::duration<double> >(end - start).count();

printf("time - type: %s std::erf Time: %7.5fsec\tv:%15.10f\n",type_str<T>(), total,v_total);
   }
}

int main()
{
   test_precision<double>();
   test_precision<float> ();
   test_time<double>     ();
   test_time<float>      ();
   return 0;
}


run \| edit \| history \| help	0