Re: Man page issues: logb, significand, cbrt, log2, log10, exp10

[Morten Welinder <mwelinder@gmail.com>]

Sorry to be a bit of a pain.

Some testing says that the average error from exp10 is 300-500 times
bigger than the average error from pow(10,.). This is consistent
across a large range of arguments.

The maximum error from pow(10,.) in the samples is 1ulp (relative to a
reference rounded value, so the true error is likely less). The
maximum error from exp10 in the samples is 2ulp.

This is not surprising given that exp10 starts out by introducing a
rounding error due to the multiplication before the call to exp.

I didn't bother looking, but it is almost certainly true that there
are arguments for which exp10 is better than pow(10,.).  However, the
numbers below imply that such arguments are very rare compared to the
other way around.


-------------------------- average ----- max
Binade -7      exp10: 0.3378           1
Binade -7  pow(10,.): 0.0007           1
Binade -6      exp10: 0.3429           2
Binade -6  pow(10,.): 0.0007           1
Binade -5      exp10: 0.3532           2
Binade -5  pow(10,.): 0.0008           1
Binade -4      exp10: 0.3774           2
Binade -4  pow(10,.): 0.0008           1
Binade -3      exp10: 0.4402           2
Binade -3  pow(10,.): 0.0010           1
Binade -2      exp10: 0.4118           2
Binade -2  pow(10,.): 0.0009           1
Binade -1      exp10: 0.4228           2
Binade -1  pow(10,.): 0.0009           1
Binade  0      exp10: 0.4204           2
Binade  0  pow(10,.): 0.0009           1
Binade  1      exp10: 0.4221           2
Binade  1  pow(10,.): 0.0009           1
Binade  2      exp10: 0.4204           2
Binade  2  pow(10,.): 0.0009           1
Binade  3      exp10: 0.4222           2
Binade  3  pow(10,.): 0.0009           1
Binade  4      exp10: 0.4209           2
Binade  4  pow(10,.): 0.0009           1
Binade  5      exp10: 0.4200           2
Binade  5  pow(10,.): 0.0009           1
Binade  6      exp10: 0.4210           2
Binade  6  pow(10,.): 0.0009           1
Binade  7      exp10: 0.4210           2
Binade  7  pow(10,.): 0.0009           1

Notes:
1. Only positive arguments tested.
2. powl (long double version of pow) is used as a reference.  I
double-checked with a double-double version of pow (good to about
100ish bits) that this does not matter.




#define _GNU_SOURCE   1
#include <stdio.h>
#include <math.h>
#include <stdint.h>


static uint64_t
murmur64 (uint64_t h)
{
  h ^= h >> 33;
  h *= 0xff51afd7ed558ccdll;
  h ^= h >> 33;
  h *= 0xc4ceb9fe1a85ec53ll;
  h ^= h >> 33;
  return h;
}

static double
exp10ref (double x)
{
  volatile double y = (double)(powl (10.0l, x));
  return y;
}

static double
exp10viapow (double x)
{
  return pow (10, x);
}


static void
test_binade (int b, double (*f) (double), const char *funcname)
{
  uint64_t h = 0x0123456701234567ll;

  double ulps = 0;
  double mulp = 0;
  int N = 1000000;

  for (int i = 0; i < N; i++) {
    h = murmur64 (h);

    double x = ldexp ((h & 0xfffffffffffffll) | 0x10000000000000ll, b - 52);

    double y = f (x);
    double yref = exp10ref (x);
    double dy = fabs (y - yref);

    double ulp = dy / (nextafter (yref, INFINITY) - yref);
    ulps += ulp;
    if (ulp > mulp) mulp = ulp;
  }

  printf ("Binade %2d %10s: %6.4f  %10.0f\n",
      b, funcname, ulps / N, mulp);
}


int
main ()
{
  for (int b = -7; b <= 7; b++) {
    test_binade (b, exp10, "exp10");
    test_binade (b, exp10viapow, "pow(10,.)");
  }

  return 0;
}