Faster Logistic Function

Question

I've noticed that a fairly significant number of cycles in one of my programs are being consumed by the logistic function: $$f(x)=\frac{1}{1+e^{-x}}$$ Is there a good approximation I can use to reduce the cost of this function?

Answer 1

Yes! There are nice approximations of the logistic.

Plot of Approximating Functions

As shown below, several functions approximate the logistic (shown as blue dots). This graph is available interactively here.

Benchmark Results of Approximating Functions

To explore this better, I've benchmarked all of the function using code at the bottom of this post. The results of that are as follows:

                                    name      rms_error        maxdiff        time_us        speedup        samples
                      logistic_with_tanh     5.9496e-02     1.5014e-01         0.0393         0.5076      200000001
                      logistic_with_atan     3.9051e-02     9.6934e-02         0.0321         0.6211      200000001
                       logistic_with_erf     6.5068e-02     1.6581e-01         0.0299         0.6676      200000001
            logistic_fexp_quintic_approx     1.2921e-07     5.9050e-07         0.0246         0.8118      200000001
        logistic_product_approx_float128     8.7032e-04     1.7217e-03         0.0209         0.9523      200000001
           logistic_with_exp_no_overflow     4.7660e-17     1.6653e-16         0.0198         1.0084      200000001
              logistic_product_approx128     8.7032e-04     1.7211e-03         0.0164         1.2187      200000001
         log_w_approx_exp_no_overflow128     8.7193e-04     1.7211e-03         0.0158         1.2640      200000001
                      logistic_with_sqrt     8.3414e-02     1.1086e-01         0.0146         1.3662      200000001
          log_w_approx_exp_no_overflow16     6.9726e-03     1.4074e-02         0.0141         1.4114      200000001
  log_w_approx_exp_no_overflow16_clamped     6.9726e-03     1.4074e-02         0.0141         1.4153      200000001
             logistic_schraudolph_approx     1.5661e-03     8.9906e-03         0.0138         1.4497      200000001
                       logistic_with_abs     6.0968e-02     8.2289e-02         0.0134         1.4936      200000001
                           logistic_orig     0.0000e+00     0.0000e+00         0.0199         ------      200000001

Discussion of Approximating Functions

Let's talk about the fastest few approximations.

logistic_with_abs

This is the fastest, but least accurate function and is given by $$f(x)=\frac{1}{2}\left(1+\frac{x}{1+|x|}\right)$$ It is 1.5x faster than the exact logistic with an RMS error of $6\cdot10^{-2}$ in the range $[-10,10]$.

The Schraudolph Approximation

Is drawn from this paper and relies on dark magic involving the IEEE754 definitions of floating-point numbers. From the paper:

After multiplication, the fractional part of y will spill over into the highest-order bits of the mantissa $m$. This spillover is not only harmless, but in fact is highly desirable—under the IEEE-754 format, it amounts to a linear interpolation between neighboring integer exponents. The technique therefore exponentiates real-valued arguments as well as a lookup table with $2^{11}$ entries and linear interpolation.

This function has a relatively low RMS of $1.6\cdot10^{-3}$ and a 1.4x speed-up versus the exact logistic.

log_w_approx_exp_no_overflow16

This function relies on the approximation $$f(x)=\frac{1}{1+\left(1-\frac{x}{n}\right)^{n}}$$ where we've used $n=16$. The value of $n=16$ is that it boils down to very simple assembly code giving a 1.4x speed-up with a decent RMS of $7.0\cdot10^{-3}$. Increasing $n$ would give a better approximation at the cost of worse performance.

Passing the code

double exp_product_approx16(double x){
  x = 1 + x / 16;
  const auto a = x * x;
  const auto b = a * a;
  const auto c = b * b;
  const auto d = c * c;
  return d;
}

inline double log_w_approx_exp_no_overflow16(double x){
  return 1/(1+exp_product_approx16(-x));
}

through Godbolt we get

.LCPI0_0:
        .quad   0x3fb0000000000000              # double 0.0625
.LCPI0_1:
        .quad   0x3ff0000000000000              # double 1
exp_product_approx16(double):              # @exp_product_approx16(double)
        mulsd   xmm0, qword ptr [rip + .LCPI0_0]
        addsd   xmm0, qword ptr [rip + .LCPI0_1]
        mulsd   xmm0, xmm0
        mulsd   xmm0, xmm0
        mulsd   xmm0, xmm0
        mulsd   xmm0, xmm0
        ret

where we see that each doubling of $n$ introduces only a single additional mulsd instruction.

log_w_approx_exp_no_overflow16_clamped

This function appears to have nearly identical characteristics to log_w_approx_exp_no_overflow16, but is defined as $$f(x)=\begin{cases} 1 & x\ge n \\ \frac{1}{1+\left(1-\frac{x}{n}\right)^{n}} & \textrm{otherwise} \end{cases}$$ if you squint at the function, you find that it increases from $(-\infty,n)$ and then decreases from $(n,\infty)$: This is fine if our inputs are all in the range $(-\infty,n]$, but if the inputs are larger than this it becomes problematic. We can handle this either by increasing $n$ (remember that doubling it only adds a single assembly instruction) or by introducing an unlikely if. We choose to use the if in case the input is adversarial and do so with very minor performance costs:

constexpr double exp_product_approx16(double x){
  x = 1 + x / 16;
  const auto a = x * x;
  const auto b = a * a;
  const auto c = b * b;
  const auto d = c * c;
  return d;
}

inline double log_w_approx_exp_no_overflow16_clamped(double x){
  if(x >= 16) [[unlikely]] {
      return 1;
  }
  return 1/(1+exp_product_approx16(-x));
}

Note that the way this code is written is important! If we were to write $x^{16}$ as:

x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x

the compiler would generate 16 multiplication instructions! The -ffast-math flag avoids this, but at the cost of potentially reducing the accuracy of math throughout your program.

Similarly, using

std::pow(x, 16)

will generate a call to the std::pow function, which will be slower than the doubling method we've used above.

Recommendation

I recommend you benchmark the approximations on your own system and choose a function with the quality/performance trade-off that works best for you. Personally, I've found log_w_approx_exp_no_overflow16_clamped to be sufficient for my needs and prefer it to the black magic of the Schraudolph approximation. This is especially so since the accuracy of the Schraudolph approximation is fixed while adjusting $n$ allows me to easily tune my accuracy/performance trade-off (see, for example, log_w_approx_exp_no_overflow128).

Benchmarking Code

// Compile with: clang++ -O3 test.cpp
// Functions plotted here: https://www.desmos.com/calculator/nkblxiypxh
#include <chrono>
#include <cmath>
#include <endian.h>
#include <iomanip>
#include <iostream>
#include <string>

constexpr double STEP_SIZE = 0.0000001;
// constexpr double STEP_SIZE = 0.001;
constexpr int NAME_LEN = 40;

constexpr double logistic_orig(double x) {
  return x > 0 ? (1.0 / (1.0 + std::exp(-x))) : (1.0 - 1.0 / (1.0 + std::exp(x)));
}

// Based on "A Fast, Compact Approximation of the Exponential Function"
// By Nicol N. Schraudolph
// https://nic.schraudolph.org/pubs/Schraudolph99.pdf
constexpr inline double exp_schraudolph_approx(double x){
  constexpr auto c2_to_the_20th = 1 << 20;
  constexpr auto ln2 = 0.6931471805599453;
  constexpr auto a = c2_to_the_20th / ln2;
  constexpr auto b = 1023 * c2_to_the_20th;
  constexpr auto c = 60801;

  union {
    double d;
// NOTE: This works for a little-endian architecture
#if __BYTE_ORDER == __LITTLE_ENDIAN
    struct {
      int j, i;
    } n;
#else
    struct {
      int i, j;
    } n;
#endif
  } eco = {};

  // Black magic happens here. Read the paper.
  eco.n.i = a * x + (b - c);
  return eco.d;
}

constexpr double logistic_schraudolph_approx(double x) {
  return x > 0 ? (1.0 / (1.0 + exp_schraudolph_approx(-x)))
               : (1.0 - 1.0 / (1.0 + exp_schraudolph_approx(x)));
}

constexpr double exp_product_approx128(double x){
  x = 1 + x / 128;
  const auto a = x * x;
  const auto b = a * a;
  const auto c = b * b;
  const auto d = c * c;
  const auto e = d * d;
  const auto f = e * e;
  const auto g = f * f;
  return g;
}

constexpr double exp_product_approx16(double x){
  x = 1 + x / 16;
  const auto a = x * x;
  const auto b = a * a;
  const auto c = b * b;
  const auto d = c * c;
  return d;
}

constexpr double logistic_product_approx128(double x) {
  return x > 0 ? (1.0 / (1.0 + exp_product_approx128(-x)))
               : (1.0 - 1.0 / (1.0 + exp_product_approx128(x)));
}

constexpr float exp_product_approx_float128(float x){
  x = 1 + x / 128;
  const auto a = x * x;
  const auto b = a * a;
  const auto c = b * b;
  const auto d = c * c;
  const auto e = d * d;
  const auto f = e * e;
  const auto g = f * f;
  return g;
}

constexpr float logistic_product_approx_float128(float x) {
  return x > 0 ? (1.0 / (1.0 + exp_product_approx_float128(-x)))
               : (1.0 - 1.0 / (1.0 + exp_product_approx_float128(x)));
}

inline double fexp_quintic(double x){
  constexpr int64_t mantissa = static_cast<int64_t>(1)<<52;
  constexpr int64_t bias = 1023;
  constexpr int64_t ishift = mantissa*bias;
  constexpr double ln2 = 0.6931471805599453;
  constexpr double s1 = -1.90188191959304e-3;
  constexpr double s2 = -9.01146535969578e-3;
  constexpr double s3 = -5.57129652016652e-2;
  constexpr double s4 = -2.40226506959101e-1;
  constexpr double s5 =  3.06852819440055e-1;
  const double y = x/ln2;
  const double yf = y-std::floor(y);
  const double y2 = y-((((s1*yf+s2)*yf+s3)*yf+s4)*yf+s5)*yf;
  const int64_t i8 = mantissa*y2+ishift;
  return *reinterpret_cast<const double*>(&i8);
}

inline double logistic_fexp_quintic_approx(double x) {
  return x > 0 ? (1.0 / (1.0 + fexp_quintic(-x)))
               : (1.0 - 1.0 / (1.0 + fexp_quintic(x)));
}

inline double logistic_with_abs(double x){
  return 0.5*(1+x/(1+std::abs(x)));
}

inline double logistic_with_tanh(double x){
  return 0.5*(1+std::tanh(x));
}

inline double logistic_with_erf(double x){
  return 0.5*(1+std::erf(std::sqrt(M_PI)*x/2));
}

inline double logistic_with_sqrt(double x){
  const auto temp = 0.5/std::sqrt(1+x*x);
  return (x<0)?temp:1-temp;
}

inline double logistic_with_atan(double x){
  return 0.5*(1+std::atan(M_PI*x/2)*2/M_PI);
}

inline double logistic_with_exp_no_overflow(double x){
  return 1/(1+std::exp(-x));
}

inline double log_w_approx_exp_no_overflow128(double x){
  return 1/(1+exp_product_approx128(-x));
}

inline double log_w_approx_exp_no_overflow16(double x){
  return 1/(1+exp_product_approx16(-x));
}

inline double log_w_approx_exp_no_overflow16_clamped(double x){
  // If you squint at exp_product_approx16 you realize it reaches
  // y=1 at x=16 and must decrease thereafter, so we use that
  // as the clamp value here.
  if(x >= 16) [[unlikely]] {
      return 1;
  }
  return 1/(1+exp_product_approx16(-x));
}

template<typename Func>
double time_it(Func func, const std::string& func_name, const double original_time){
  const auto start = std::chrono::high_resolution_clock::now();
  double diff = 0;
  int count = 0;
  double maxdiff = -std::numeric_limits<double>::infinity();
  for(double x=-10;x<10;x+=STEP_SIZE){
      const auto origval = logistic_orig(x);
      const auto newval = func(x);
      diff += std::pow(origval - newval, 2);
      maxdiff = std::max(maxdiff, std::abs(origval - newval));
      count++;
  }
  const auto end = std::chrono::high_resolution_clock::now();
  const auto time_per_sample = std::chrono::duration_cast<std::chrono::microseconds>(end-start).count()/static_cast<double>(count);
  std::cout<<std::setw(NAME_LEN)<<func_name
           <<std::scientific<<std::setprecision(4)<<std::setw(15)<<std::sqrt(diff/count)
           <<std::scientific<<std::setprecision(4)<<std::setw(15)<<maxdiff
           <<std::fixed<<std::setw(15)<<time_per_sample
           <<std::fixed<<std::setw(15)<<original_time/time_per_sample
           <<std::setw(15)<<count
           <<std::endl;
  return time_per_sample;
}

int main(){
  std::cout<<std::setw(NAME_LEN)<<"name"
           <<std::setw(15)<<"rms_error"
           <<std::setw(15)<<"maxdiff"
           <<std::setw(15)<<"time_us"
           <<std::setw(15)<<"speedup"
           <<std::setw(15)<<"samples"
           <<std::endl;
  double original_speed = 1;
  original_speed = time_it(logistic_orig, "logistic_orig", original_speed);
  time_it(logistic_schraudolph_approx, "logistic_schraudolph_approx", original_speed);
  time_it(logistic_product_approx128, "logistic_product_approx128", original_speed);
  time_it(logistic_fexp_quintic_approx, "logistic_fexp_quintic_approx", original_speed);
  time_it(logistic_product_approx_float128, "logistic_product_approx_float128", original_speed);
  time_it(logistic_with_abs, "logistic_with_abs", original_speed);
  time_it(logistic_with_tanh, "logistic_with_tanh", original_speed);
  time_it(logistic_with_erf, "logistic_with_erf", original_speed);
  time_it(logistic_with_sqrt, "logistic_with_sqrt", original_speed);
  time_it(logistic_with_atan, "logistic_with_atan", original_speed);
  time_it(logistic_with_exp_no_overflow, "logistic_with_exp_no_overflow", original_speed);
  time_it(log_w_approx_exp_no_overflow128, "log_w_approx_exp_no_overflow128", original_speed);
  time_it(log_w_approx_exp_no_overflow16, "log_w_approx_exp_no_overflow16", original_speed);
  time_it(log_w_approx_exp_no_overflow16_clamped, "log_w_approx_exp_no_overflow16_clamped", original_speed);

  return 0;
}

Answer 2

If only low-accuracy approximations are needed, it is highly advisable to perform all computation in single precision, for example IEEE-754 binary32 format, usually mapped to the float type in C and C++. If the code is vectorized, fixed-width SIMD machines offer twice as many concurrently active SIMD-lanes when computing with single precision versus to double precision. Divisions usually have lower latency in single precision compare to double precision.

The function $\frac{1}{e^{-x}+1}$ is well behaved numerically, and for some platforms it can suffice to code it in the naive fashion and apply the compiler's "fast math" flag. For example, in the CUDA programming environment the computation, when performed in single precision and with -use_fast_math specified will map straight to hardware instructions implementing the reciprocal, $\frac{1}{x}$, and binary exponentiation, $2^{x}$, both of which are accurate to within a couple of ulps.

With the exception of some DSPs and processors compliant with MIL-STD-1750A, all major processor architectures (in particular x86, ARM, and Power) offer support for IEEE-754. For these I would recommend using a straightforward implementation based on exponentiation that uses a variant of the Schraudolph algorithm. The function logistic_fexp_quintic_approx() in asker's self-answer is an example of such an implementation. The algorithm is easily adjusted to various levels of accuracy, is very competitive in the number of instructions required with other approaches, is trivially amenable to the use of fused multiply-add operations, and is easily vectorized on fixed-width SIMD machines.

N. N. Schraudolph. "A fast, compact approximation of the exponential function." Neural Computation, 11(4), May 1999, pp. 853-862 (online).

Schraudolph uses a simple approach matched well to binary floating-point formats. First, $e^{x}$ is turned into $2^{y}$, where $y := x \log_2 e$. $\lfloor y\rfloor$ is the exponent of the result, one needs to add the exponent bias of the floating-point format to it before inserting it at the appropriate position. $f := y - \lfloor y\rfloor$ is a fraction in $[0, 1)$ and the simplest possible approximation to $2^{f}$ on that interval is $1+f$. However, we know that except at the endpoints of the interval this provides overestimate, so more accurate results can be achieved by subtracting a small constant offset. The exact value of this offset can be found by simple binary search. We also note that we do not have to explicitly separate out the fraction $f$ (for a more detailed discussion, see Schraudolph's original paper).

While Schraudolph's code performs the math needed for the manipulation of a floating-point number at bit-level with integer arithmetic, this is not necessary, as observed, for example, by Paul Mineiro. Performing the computation with floating-point arithmetic, then storing the result to an integer may be advantageous, for example when the processor has support for fused multiply-add. The drawback is that this manner of computation forces some number of least significant bits (at least five) to zero. Therefore Mineiro's variant of Schraudolph's algorithm tends to be a tad faster and a tad less accurate, as can also be glimpsed from the data in the table below.

A flexible approach is easily achieved by making the computation of $f$ explicit, and then subtracting a polynomial correction $x\mathrm{P}(x) \approx 2^{x}-(1+x)$ instead of a simple constant from the baseline approximation $1+f$. Minimax polynomial approximations for this can be generated with various common software, for example Maple, Mathematica, or the open-source Sollya tool. For a high-performance implementation, polynomials of degree up to $4$ may be useful, depending on hardware platform and use case.

Below I am showing exemplary ISO-C99 code that demonstrates the recommended approach in action. All variants of $\exp(x)$ used to compute the logistic function are limited to $[-87.33654, 88.72283]$ to avoid having to handle underflow or overflow in the IEEE-754 binary32 format. Single-precision reciprocals are computed by most modern processors with throughput high enough that it rarely makes sense to roll one's own implementation, but I am including code for that just in case (USE_OWN_RCP), with a simple starting approximation followed by Newton-Raphson or Halley steps. The following table lists the relative error for the logistic function on $[-10, 10]$ with the $\exp$ variants made available in the code.

Relative error of logistic_function() on [-10,10]

                     builtin reciprc  own reciprocal   own reciprocal   own reciprocal
                                      high accuracy    medium accuracy  low accuracy  
expf implementation  max err RMS err  max err RMS err  max err RMS err  max err RMS err
---------------------------------------------------------------------------------------
Mineiro              3.44e-2 9.26e-3  3.44e-2 9.26e-3  3.42e-2 9.14e-3  3.68e-2 7.07e-3
flex Mineiro, deg=0  3.44e-2 9.26e-3  3.44e-2 9.26e-3  3.43e-2 9.14e-3  3.68e-2 7.07e-3
flex Mineiro, deg=1  3.05e-3 2.89e-4  3.05e-3 2.89e-4  3.04e-3 3.17e-4  5.18e-3 2.53e-3
flex Mineiro, deg=2  2.67e-4 2.02e-5  2.62e-4 2.11e-5  3.84e-4 1.29e-4  2.62e-3 2.50e-3
flex Mineiro, deg=3  1.44e-5 8.31e-7  1.81e-5 6.39e-6  1.39e-4 1.26e-4  2.56e-3 2.50e-3
flex Mineiro, deg=4  8.50e-6 4.93e-7  1.45e-5 6.35e-6  1.36e-4 1.26e-4  2.56e-3 2.50e-3

Schraudolph          3.44e-2 9.26e-3  3.44e-2 9.26e-3  3.42e-2 9.14e-3  3.68e-2 7.07e-3
flex Schraudolph, 0  3.44e-2 9.26e-3  3.44e-2 9.26e-3  3.42e-2 9.14e-3  3.68e-2 7.07e-3
flex Schraudolph, 1  3.04e-3 2.89e-4  3.05e-3 2.89e-4  3.04e-3 3.17e-4  5.17e-3 2.53e-3
flex Schraudolph, 2  2.61e-4 2.02e-5  2.57e-4 2.11e-5  3.79e-4 1.29e-4  2.61e-3 2.50e-3
flex Schraudolph, 3  7.93e-6 6.70e-7  1.29e-5 6.37e-6  1.34e-4 1.26e-4  2.56e-3 2.50e-3
flex Schraudolph, 4  6.72e-7 3.15e-8  7.11e-6 6.34e-6  1.29e-4 1.26e-4  2.55e-3 2.50e-3

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

#define TEST_LOGISTIC_FUNCTION (1)
#define fast_expf              flex_schraudolph_expf
#define DEGREE                 (2) // flex_mineiro_expf, flex_schraudolph_expf

#define USE_OWN_RCP            (0) // use own reciprocal implementation
#define RCP_ACCURACY           (0) // accuracy: 0 = low, 1 = medium, 2 = high

uint32_t float_as_uint32 (float a);
float uint32_as_float (uint32_t a);

// max rel err = 3.56027983e-2  RMS rel err = 1.82537677e-2
float schraudolph_expf (float x)
{
    const int FP32_MANT_BITS = 23;
    const int FP32_EXPO_BIAS = 127;
    const float L2E = 1.442695041f;
    const int CORR = -298685;
    float a = L2E * (1 << FP32_MANT_BITS);
    int b = FP32_EXPO_BIAS * (1 << FP32_MANT_BITS) + CORR;
    int i = (int)(a * x) + b;
    return uint32_as_float (i);
}

// max rel err = 3.56061513e-2  RMS rel err = 1.82539111e-2
float mineiro_expf (float x)
{
    const int FP32_MANT_BITS = 23;
    const float FP32_EXPO_BIAS = 127;
    const float L2E = 1.442695041f;
    const float CORR = -298685;
    float c0 = L2E * (1 << FP32_MANT_BITS);
    float c1 = FP32_EXPO_BIAS * (1 << FP32_MANT_BITS) + CORR;
    int i = (int)(c0 * x + c1);
    return uint32_as_float (i);
}

float flex_schraudolph_expf (float x)
{
    const int FP32_MANT_BITS = 23;
    const int FP32_EXPO_BIAS = 127;
    const float L2E = 1.442695041f;
#if DEGREE == 0
    // max rel err = 3.56027983e-2  RMS rel err = 1.82537677e-2
    const int CORR = -298685;
    float ax = x * L2E * (1 << FP32_MANT_BITS);
    int b = FP32_EXPO_BIAS * (1 << FP32_MANT_BITS) + CORR;
#else // DEGREE > 0
    float y = x * L2E;
    float f = y - floorf (y);
    float ax = y * (1 << FP32_MANT_BITS);
#if DEGREE == 1
    // max rel err = 3.06149082e-3  RMS rel err = 5.87195281e-4
    float fcorr = (3.36369932e-1f * f - 3.36369932e-1f) * f;
#elif DEGREE == 2
    // max rel err = 2.64595565e-4  RMS rel err = 4.06953315e-5
    float fcorr = ((7.79931966e-2f * f + 2.27580582e-1f) * f - 
                   3.05573779e-1f) * f;
#elif DEGREE == 3
    // max rel err = 1.13272671e-5  RMS rel err = 1.36108954e-6
    float fcorr = (((1.35335097e-2f * f + 5.22205270e-2f) * f +
                    2.41143385e-1f) * f - 3.06897422e-1f) * f;
#elif DEGREE == 4
    // max rel err = 4.13295586e-6  RMS rel err = 1.80227819e-7
    float fcorr = ((((1.87805750e-3f * f + 9.01247416e-3f) * f + 
                     5.57802263e-2f) * f + 2.40180765e-1f) * f -
                     3.06851522e-1f) * f;
#endif // DEGREE
    int corr = fcorr * (1 << FP32_MANT_BITS);
    int b = FP32_EXPO_BIAS * (1 << FP32_MANT_BITS) + corr;
#endif // DEGREE > 0
    int i = ((int)ax) + b;
    float r = uint32_as_float (i);
    return r;
}

float flex_mineiro_expf (float x)
{
    const int FP32_MANT_BITS = 23;
    const float FP32_EXPO_BIAS = 127;
    const float L2E = 1.442695041f;
    float y = x * L2E;
#if DEGREE == 0
    // max rel err = 3.56061513e-2  RMS rel err = 1.82539012e-2
    const float CORR = -298685;
    float yc = CORR / (1 << FP32_MANT_BITS) + y;
#else // DEGREE > 0 
    float f = y - floorf (y);
#if DEGREE == 1
    // max rel err = 3.06639459e-3  RMS rel err = 5.87191814e-4
    float yc = (3.36369932e-1f * f - 3.36369932e-1f) * f + y;
#elif DEGREE == 2
    // max rel err = 2.72702222e-4  RMS rel err = 4.07118944e-5
    float yc = ((7.79931966e-2f * f + 2.27580582e-1f) * f - 
                3.05573779e-1f) * f + y;
#elif DEGREE == 3
    // max rel err = 1.89591946e-5  RMS rel err = 1.63955571e-6
    float yc = (((1.35335097e-2f * f + 5.22205270e-2f) * f +
                 2.41143385e-1f) * f - 3.06897422e-1f) * f + y;
#elif DEGREE == 4
    // max rel err = 1.51669284e-5  RMS rel err = 9.32432780e-7
    float yc = ((((1.87805750e-3f * f + 9.01247416e-3f) * f + 
                  5.57802263e-2f) * f + 2.40180765e-1f) * f -
                  3.06851522e-1f) * f + y;
#endif // DEGREE 
#endif // DEGREE > 0
    int i = yc * (1 << FP32_MANT_BITS) + FP32_EXPO_BIAS * (1 << FP32_MANT_BITS);
    return uint32_as_float (i);
}

#if USE_OWN_RCP
float fast_rcp (float a)
{
    float x = uint32_as_float (0x7ef311c2 - float_as_uint32 (a));
    float e = 1.0f - a * x;
#if RCP_ACCURACY == 0
    float r = e * x + x;
#elif RCP_ACCURACY == 1
    float r = (e * e + e) * x + x;
#elif RCP_ACCURACY == 2
    x = e * x + x;
    e = 1.0f - a * x;
    float r = e * x + x;
#endif // RCP_ACCURACY
    return r;
}
#endif // USE_OWN_RCP

float logistic_function (float x)
{
#if USE_OWN_RCP
    return fast_rcp (fast_expf (-x) + 1.0f);
#else // USE_OWN_RCP
    return 1.0f / (fast_expf (-x) + 1.0f);
#endif // USE_OWN_RCP
}

double ref_logistic_function (double x)
{
    return 1.0 / (exp (-x) + 1.0);
}

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

#define xstr(a) str(a)
#define str(a) #a

int main (void)
{
#if TEST_LOGISTIC_FUNCTION
    const float MIN_ARG = -10.0f;
    const float MAX_ARG =  10.0f;
#else // TEST_LOGISTIC_FUNCTION
    const float MIN_ARG = -87.33654f;
    const float MAX_ARG =  88.72283f;
#endif // TEST_LOGISTIC_FUNCTION
    float x, res, errloc = INFINITY;
    double ref, err, rmserr, maxerr = 0, sumerrsq = 0;
    int64_t count = 0;

    printf ("Testing %s using %s %s %c\n", TEST_LOGISTIC_FUNCTION ? 
            "logistic function" : "exp",
            xstr (fast_expf), 
            ((fast_expf == flex_mineiro_expf) ||
             (fast_expf == flex_schraudolph_expf)) ? "degree" : "", 
            ((fast_expf == flex_mineiro_expf) ||
             (fast_expf == flex_schraudolph_expf)) ? ('0' + DEGREE) : ' ');
#if TEST_LOGISTIC_FUNCTION && USE_OWN_RCP
    printf ("own reciprocal (%s accuracy)\n", (RCP_ACCURACY > 0) ? 
            ((RCP_ACCURACY > 1) ? "high" : "medium") : "low");
#endif // TEST_LOGISTIC_FUNCTION && USE_OWN_RCP

    for (x = 0.0f; x <= MAX_ARG; x = nextafterf (x, INFINITY)) {
#if TEST_LOGISTIC_FUNCTION
        res = logistic_function (x);
        ref = ref_logistic_function ((double)x);
#else // TEST_LOGISTIC_FUNCTION
        res = fast_expf (x);
        ref = exp ((double)x);
#endif // TEST_LOGISTIC_FUNCTION
        err = fabs (((double)res - ref) / ref);
        if (err > maxerr) { maxerr = err; errloc = x; }
        sumerrsq += err * err;
        count++;
    }
    for (x = -0.0f; x >= MIN_ARG; x = nextafterf (x, -INFINITY)) {
#if TEST_LOGISTIC_FUNCTION
        res = logistic_function (x);
        ref = ref_logistic_function ((double)x);
#else // TEST_LOGISTIC_FUNCTION
        res = fast_expf (x);
        ref = exp ((double)x);
#endif // TEST_LOGISTIC_FUNCTION
        err = fabs (((double)res - ref) / ref);
        if (err > maxerr) { maxerr = err; errloc = x; }
        sumerrsq += err * err;
        count++;
    }
    rmserr = sqrt ((sumerrsq / count));
    printf ("max rel err = %15.8e (@%15.8e)  RMS rel err = %15.8e\n", maxerr, 
            errloc, rmserr);
    return EXIT_SUCCESS;
}