Handling large number calculation with pow() in C

Question

I tried to calculate pow(.1, -120) in C and got:

999999999999993386194947375938605300558731199397053728064304453022541669059515488061520680536169238451435883569728192512.000000

instead of 10^120. Is there a reason why? Thanks in advance for your reply!

I don't understand the downvotes tbh. +1. I have a couple theories but I don't know the answer. — David, Aug 03 '16 at 18:49
`pow` is not very accurate. For that matter, nor is much floating point arithmetic, please see [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken). — Weather Vane, Aug 03 '16 at 18:55
@WeatherVane: "not very accurate"??? `pow` is probably off by no more than a couple of ULPs. But of course, you start with an approximation to .1 and raising that to a large power increases the error factor a bit. All the same, the end result is within one in a (European) billion (i.e. 10^12) — rici, Aug 03 '16 at 18:58
@rici: I am not so certain that it is only a few ULPs. I guess it depends on the implementation. Although this value looks as if it is only off by a few ULPs. — Rudy Velthuis, Aug 03 '16 at 19:25
@RudyVelthuis: Certainly it depends on the implementation, and an implementation which sacrifices accuracy for speed is conceivable; the C standard does not prohibit it. In practice, though, there is a tendency to go for accuracy over speed. And, of course, bugs are possible. Historically, there have been issues with glibc related to certain computations with very large exponents, but they were definitely treated as bugs and, afaik, eventually fixed. — rici, Aug 03 '16 at 19:45
@Rudy ... this value is off by a lot more than a few ULPs if you compare it with 10^12. but it is within a few if you compare it with the inverse of the 120th power of the value actually being used to approximate 0.1. The difference is important because it gives a clue about how to deal with the problem. For example, replacing the `pow` call in this case by a sequence of multiplications is likely to *increase* the error, not decrease it, even though one's intuition might be that multiplication is "more exact". — rici, Aug 03 '16 at 19:47
Glibc table of math library accuracy: https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html#Errors-in-Math-Functions — rici, Aug 03 '16 at 19:57
If you take the hex representation, it is clear that the difference is indeed more than a few ULPs: `result=0x1.8c8dac6a033fcp+398 10^120=0x1.8c8dac6a0342ap+398`. It is off by ~11 bits, AFAICT. — Rudy Velthuis, Aug 03 '16 at 20:29
@Rudy Velthuis By my reckoning in [answer](http://stackoverflow.com/a/38754586/2410359), the result of `pow(0.1, -120)` is within 0.5 ULP. — chux - Reinstate Monica, Aug 03 '16 at 22:37
fact is that the difference between 10^120 and pow(.1, -120) is far more than a few ulp, as the hex representations show. That does not mean that pow() is inaccurate, only that 0.1 is the problem. — Rudy Velthuis, Aug 03 '16 at 22:51

dbush · Answer 1 · 2016-08-03T19:03:28.040

This is due to the inexact nature of floating point. There are not enough significant digits in a double to express that number exactly since numbers are stored internally in binary.

For example, if you run this:

printf("result=%.10e\n", pow(.1, -120));
printf("result=%.20e\n", pow(.1, -120));
printf("result=%f\n", pow(.1, -120));
printf("10^120=%.10e\n", 1e120);
printf("10^120=%.20e\n", 1e120);
printf("10^120=%f\n", 1e120);

You'll get this:

result=1.0000000000e+120
result=9.99999999999993386195e+119
result=999999999999993386194947375938605300558731199397053728064304453022541669059515488061520680536169238451435883569728192512.000000
10^120=1.0000000000e+120
10^120=9.99999999999999980003e+119
10^120=999999999999999980003468347394201181668805192897008518188648311830772414627428725464789434929992439754776075181077037056.000000

The reason pow(.1, -120) is slightly different than 1e120 is that .1 cannot be represented exactly either, so the error magnifies with each iteration of the power.

If on the other hand you printed the result of pow(2, 350), you would get an exact answer.

I think it would help to also print the approximation of .1 which is used in that calculation. It's not so much that pow(x,120) is inaccurate (although of course it is), but that the `x` was inaccurate to start with. — rici, Aug 03 '16 at 18:59

score 1 · Answer 2 · edited May 23 '17 at 12:14

0.1 is typically not exactly representable as a double. Typical double (binary64) can exactly represent about 2^64 different numbers. 0.1 is not one of them. The closest double is 0.100000000000000005551... @rici

So the question is what to expect of pow(0.100000000000000005551..., -120)?

Using Binomial theorem, this would be

pow(exact_one_tenth, -120) 
    - 120*pow(exact_one_tenth, -119)*(0.100000000000000005551... - exact_one_tenth)
    + other smaller terms

1e120 - 120*1e119*(0.000000000000000005551...) + other smaller terms

9.99999999999993_3386...e119

// OP result
9.99999999999993_3861...e119

Looks like a very accurate answer (matches to 16 significant digts.) The next smaller double from OP's answer is below. Since the correct answer is bounded by OP's result and the next representable number and closer to the OP's answer, I assert the calculation is a good as it gets, within 0.5 ULP.

9.99999999999993_2429...e+119

Handling large number calculation with pow() in C

2 Answers2