`exp2f` benchmark

test_exp2f.cpp

#include <cstdio>
#include <iostream>
#include <chrono>
#include <cmath>
#include <vector>

using namespace std;

inline float fastexp2 (float p)
{
    // Mineiro's method, from: http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html
    union { float f; uint32_t i; } vp = { p };
    int sign = (vp.i >> 31);
    int w = p;
    float z = p - w + sign;
    union { uint32_t i; float f; } v = { (1 << 23) * (p + 121.2740838f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z) };
    return v.f;
}

static float P[] = {
    2.30933477057345225087E-2f,
    2.02020656693165307700E1f,
    1.51390680115615096133E3f,
};
static float Q[] = {
    /* 1.00000000000000000000E0,*/
    2.33184211722314911771E2f,
    4.36821166879210612817E3f,
};

inline float polevl(float x, float* coef, int N)
{
  float ans;
  int i;
  float *p;

  p = coef;
  ans = *p++;
  i = N;

  do
    ans = ans * x + *p++;
  while (--i);

  return (ans);
};

inline float p1evl(float x, float* coef, int N)
{
  float ans;
  float *p;
  int i;

  p = coef;
  ans = x + *p++;
  i = N - 1;

  do
    ans = ans * x + *p++;
  while (--i);

  return (ans);
};

inline float fastexp2p2 (float x)
{
    // Cephes' method, from: https://github.com/nearform/node-cephes/blob/master/cephes/exp2.c
    float xx, px, n;
    xx = x; /* save x */
    /* separate into integer and fractional parts */
    px = floor(x + 0.5);
    n = px;
    x = x - px;

    /* rational approximation
    * exp2(x) = 1 +  2xP(xx)/(Q(xx) - P(xx))
    * where xx = x**2
    */
    xx = x * x;
    px = x * polevl(xx, P, 2);
    x = px / (p1evl(xx, Q, 2) - px);
    x = 1.0 + ldexp(x, 1);

    /* scale by power of 2 */
    x = ldexp(x, n);
    return (x);
}

inline float fastexp2p3 (float x)
{
    // 3rd order polynomial approximation, from `np.polyfit`
    return 0.05700169f * x * x * x + 0.24858144f * x * x + 0.69282515f * x + 0.9991608f;
}

int main() {
    typedef std::chrono::high_resolution_clock Time;
    typedef std::chrono::milliseconds ms;
    typedef std::chrono::duration<float> fsec;
    const int N = 10000;
    const float low = -0.5f;
    const float upp = 0.5f;

    const float sr = (upp - low) / (float)N;
    float k = low;

    vector<float> vecC(N);

    std::chrono::steady_clock::time_point t0 = Time::now();
    for (int i = 0; i < N; i++) {
        vecC[i] = fastexp2(k);
        k += sr;
    }
    std::chrono::steady_clock::time_point t1 = Time::now();
    fsec fs = t1 - t0;
    std::cout << "Used time: " << fs.count() << " secs for " << N << " calls" << std::endl;

    float avgerr = 0;
    float maxerr = -1000;
    float minerr = 1000;
    k = low;
    for (int i = 0; i < N; i++) {
        float err = abs(vecC[i] - exp2f(k));
        
        avgerr += err;
        maxerr = max(err, maxerr);
        minerr = min(err, minerr);

        k += sr;
    }
    avgerr /= N;
    std::cout << "maxerr: " << maxerr << " minerr: " << minerr << " avgerr: " << avgerr << std::endl;
    return 1;
}