exp10 different to pow(10)

Question

Firstly, I realise that most base 10 numbers cannot be precisely expressed in base 2, and so my question isn't really about the deficiencies of floating point arithmetic.

I am trying to write a function that will attempt to correct a double tainted by cumulative rounding error by checking the last 6 meaningful digits are within some tolerance and changing it to the next representable above some supposed exact value (only for display purposes - unless it is an integer or power of two).

A component of my function that surprises me is the output of exp10 though; As far as I'm aware, so long as the spacing between two doubles is less than 2 then integer values stored as doubles should be exact - and though 10^14 is pushing it, this should be an exact integer (since 10^14 =~ 2^46.507 < 2^53). However this is not what my testing shows.

An excerpt of my debugging efforts (nothing stands out as obvious) and output is as follows:

double test = 0.000699;
double tmp = fabs(test);
double exp = 10.0 - floor(log10(tmp));
double powTen = exp10(10.0 - floor(log10(tmp)));
double powTen2 = exp10(exp);
double powTen3 = exp10((int)exp);
double powTen4 = exp10(exp);
double powTen5 = pow(10, exp);

printf("exp: %.16lf\n", exp);
printf("powTen: %.16lf\n", powTen);
printf("powTen2: %.16lf\n", powTen2);
printf("powTen3: %.16lf\n", powTen3);
printf("powTen4: %.16lf\n", powTen4);

//these two are exact
printf("10^14: %.16lf\n", exp10(14));
printf("powTen5: %.16lf\n", powTen5);
printf("exp == 14.0: %d\n", exp == 14.0);

output:

exp: 14.0000000000000000
powTen: 100000000000000.1250000000000000
powTen2: 100000000000000.1250000000000000
powTen3: 100000000000000.1250000000000000
powTen4: 100000000000000.1250000000000000
10^14: 100000000000000.0000000000000000
powTen5: 100000000000000.0000000000000000
exp == 14.0: 1

pow is getting the answer exact, as is exp10 with a hardcoded int. For all other cases I am adding in 1/8 (the spacing between 10^14 and 10^14 + next representable is 1/64). The documentation says that exp10 should be equivalent to pow. Can anyone see something I'm missing?

Edit - with O3, O2, O1 optimisation I am getting the expected outputs - unless the data cannot be known until runtime. at this point exp10 is still misbehaving.

Answer 1

It is probable that your exp10 implementation is misbehaving. Note that the results it's returning are sometimes off by an ulp (that's 0.125 relative to your 10^14.)

This is a rather heinous error; you've got a case where the correct answer is representable as a double yet exp10 isn't doing so.

I'd echo Ben Voigt's comment that the compiler may sometimes evaluate things itself instead of passing them off to the math library. It's probably doing a better job since it probably links to arbitrary-precision math library. You might experiment with the -fno-builtin option to see whether it changes anything.

Unfortunately, I don't think crlibm has implemented exp10. Otherwise I'd recommend you just use that and stop worrying.

EDIT: The copy of eglibc source I have seems to implement exp10 thus:

double
__ieee754_exp10 (double arg)
{
  /* This is a very stupid and inprecise implementation.  It'll get
     replaced sometime (soon?).  */
  return __ieee754_exp (M_LN10 * arg);
}

Don't expect this to work well.