2

I recently wrote some code similar to this:

// Calculate a ^ b
unsigned int result = 1;
for (unsigned int i = 0; i < b; i++) {
    result *= a;
}

I got the comment that I should have used pow from math.h, and I was all ready to point out the floating-point precision issues because pow returns a double. I even found an existing Stack Overflow question that illustrates the problem. Except when I tried to run the code from that other question, it "worked" just fine (i.e. no off-by-1 issues).

Here's the code I've been testing:

#include <math.h>
#include <stdio.h>

int main()
{
    unsigned int a = 10;
    unsigned int result;
    for (unsigned int b = 1; b <= 9; b++) {
        result = (unsigned int)pow(a, b);
        printf("%u\n", result);
    }
    return 0;
}

And here's how I'm compiling and running it (on Ubuntu 18.04.3 with GCC version 7.4.0):

$ gcc -O3 -std=c99 -Wall -Wextra -Werror -Wpedantic -pedantic-errors -o pow_test pow_test.c -lm
$ ./pow_test
10
100
1000
10000
100000
1000000
10000000
100000000
1000000000

So why does (unsigned int)pow(a, b) work? I'm assuming the compiler is doing some magic to prevent the normal floating-point issues? What is it doing and why? Seems like kind of a bad idea to encourage the ignorance of these issues if not all compilers do this. And if all modern compilers do do this, is this not a problem you have to worry about as much anymore?

I've also seen people suggest something like (unsigned int)(pow(a, b) + 0.1) to avoid off-by-1 issues. What are the advantages/disadvantages of doing that vs. my solution with the loop?

jnrbsn
  • 2,498
  • 1
  • 18
  • 25
  • 1
    I think you're just getting lucky that all the errors are in the positve direction rather than negative. So the results are like `100000.00001` rather than `999999.99999`. But it could have gone the other way. – Barmar Jan 13 '20 at 21:39
  • 2
    Whoever told you to use `pow` for integer calculations was wrong. Yet your power function can be improved. – Eugene Sh. Jan 13 '20 at 21:43
  • @Barmar I don't think that's what's happening. I tested it with many other values, and if I print out like 100 decimal places of the `double` using `printf`, it's all zeros. – jnrbsn Jan 13 '20 at 21:50
  • 1
    An implementation is free to treat integer arguments as a special case. But it's not required to do so. You tested with one implementation. Only 99 more to go before you'll have some decent statistics on the percentage of modern implementations that produce exact results for integer arguments. – user3386109 Jan 13 '20 at 21:53
  • @user3386109 Does the standard actually say that? – jnrbsn Jan 13 '20 at 21:58
  • @jnrbsn No, of course not. That would be too easy. Instead the standard has some informative comments like *"The relevant C arithmetic types meet the requirements of LIA−1 types if an implementation adds notification of exceptional arithmetic operations and meets the 1 unit in the last place (ULP) accuracy requirement (LIA−1 subclause 5.2.8)."* To get more information about floating point math, you need the IEEE-754 or IEEE-854 spec, or the IEC-60559 specification. LIA-1 refers to ISO/IEC 10967−1. – user3386109 Jan 13 '20 at 22:22
  • jnrbsn: Quality of implementation: The `pow()` you use is better than others. – chux - Reinstate Monica Jan 14 '20 at 00:46
  • Note that compiler may analyse the code and see `pow()` is always used as `pow(10, some_int)` and then call a different function/code to do the exponentiation. A more challenging test of `pow(a,b)` could pass in `volatile double a,b` to prevent such optimization. – chux - Reinstate Monica Jan 14 '20 at 00:52

2 Answers2

4

pow has double arguments and return value (cppreference.com). double has 53 bit significand precision (Wikipedia). On the other hand maximum value in the example is 1000000000, which has 30 bits, so perfectly fits without loosing precision into the 53 bits.

Update: played with Godbolt and it turned out that compiler could even optimize out all the computations, as in your example a and all the b values are known at compile time, for example clang v9.0. with -O3 option just prints constants (and gcc and msvc do a call to pow)

Renat
  • 7,718
  • 2
  • 20
  • 34
  • 1
    So why is `pow(10, i) == 99` when `i == 2` in the linked question? – Barmar Jan 13 '20 at 21:56
  • I thought the size didn't matter. I thought it was more about the fact that decimal floating-point numbers can't be represented accurately using a binary format. I can find several other SO questions where people had this exact problem with small integers. – jnrbsn Jan 13 '20 at 21:57
  • 1
    @jnrbsn Some decimal *fractions* can't be represented exactly, but whole numbers can be. – Barmar Jan 13 '20 at 21:57
  • @Barmar: The `pow` in the linked question returns a value other than 100 for `pow(10, 2)` because it is not a good implementation of `pow`. It is the same as if `3. * 4.` returned a value other than 12. We would not say that you get 11.9999…5 because of floating-point errors; we would say somebody did a bad job. – Eric Postpischil Jan 13 '20 at 22:20
  • 1
    @EricPostpischil The difference is that there are hardware instructions for floating point multiplication, so no one would even think of using logarithms to implement multiplication. There's no exponentiation instruction in most machines, so the implementor has to decide on an algorithm, and using logarithms for everything is an easy option. – Barmar Jan 13 '20 at 22:26
  • @Barmar. These are irrelevant. Floating-point arithmetic is sometimes implemented on integer-only hardware. If it gets the wrong answer, it is wrong. Similarly, there is no barrier to getting the correct answer for exactly representable results for exponentiation, and using logarithms is no excuse. A sub-ULP `pow` can be implemented with logarithms. If an implementation gets the wrong answer, it is bad. – Eric Postpischil Jan 13 '20 at 23:31
3

This is implementation-dependent.

The implementation that was being used in the linked question apparently uses transcendental functions to implement pow(), even when the arguments are precise integers. As a result, it gets round-off errors, so the result is sometimes slightly lower than the actual integer result.

The implementation you're using apparently checks whether the arguments have zero fractions. In that case, it uses an algorithm that produces accurate integer results (it can be done by shifting the mantissa and some additional multiplications). As long as the result fits within the precision of double, you get a float that will convert back to the precise integer.

Barmar
  • 741,623
  • 53
  • 500
  • 612
  • When exponentiation has an exact result, the cause of round-off errors is not the fact that transcendental functions were used but the fact that the implementation was not good. Essentially every `pow` implementation must use transcendental functions. But, for example, the error in the macOS `pow` function is less than one ULP. When a mathematical result is representable, the only representable value that has error less than an ULP is the mathematical result, so there is no error. – Eric Postpischil Jan 13 '20 at 22:19
  • 2
    @EricPostpischil I don't think it's necessary for it always use transcendental functions, that's my point. It can detect that both parameters have no fraction, and switch to a simple loop. – Barmar Jan 13 '20 at 22:21