Weird C program behaviour

Question

I have the following C program:

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

int main() {
    const int opt_count = 2;

    int oc = 30;

    int c = 900;

    printf("%d %f\n", c, pow(oc, opt_count));
    assert(c == (int)(pow(oc, opt_count)));
}

I'm running MinGW on Windows 8.1. Gcc version 4.9.3. I compile my program with:

gcc program.c -o program.exe

When I run it I get this output:

$ program
900 900.000000
Assertion failed: c == (int)(pow(oc, opt_count)), file program.c, line 16

This application has requested the Runtime to terminate it in an unusual way.
Please contact the application's support team for more information.

What is going on? I expect the assertion to pass because 900 == 30^2.

Thanks!

Edit

I'm not using any fractions or decimals. I'm only using integers.

Possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) — too honest for this site, Sep 09 '16 at 18:00
@Olaf I'm aware that computers can't represent 0.1 in binary floating point. But they certainly can represent 3, 30, 900, etc all exactly. — herugar5, Sep 09 '16 at 18:09
@herugar5 And what about the intermediate calculations? Also *how* do you know `900` can be represented exactly? — Eugene Sh., Sep 09 '16 at 18:12
@herugar5: But `log(30)` is neither integer nor exactly representable, thus `pow(30,2)=899.9999999...` and truncating to integer results in `899`. — Lutz Lehmann, Sep 09 '16 at 18:14
@EugeneSh. `Also how do you know 900 can be represented exactly?` Assuming IEEE 754: [Which is the first integer that an IEEE 754 float is incapable of representing exactly?](http://stackoverflow.com/q/3793838/1606345) and yes, the problem comes with intermediate calculations — David Ranieri, Sep 09 '16 at 18:20
@herugar5: You are wrong - in a general sense. It is corect, it can be represented in C floating points, becuase the language mandates this. However, a smaller format might not be able to. In general, it is a bad idea to rely on floating point to be able to precisely represent a value. Less if operations are involved like `pow`. — too honest for this site, Sep 09 '16 at 18:40
@AlterMann: C does not mandate using IEEE754 floating point with all its implications. — too honest for this site, Sep 09 '16 at 18:40
@Olaf I named the compiler I was using. It does in fact use IEEE754. Can you name any C floating-point implementations that can't exactly represent small integer values? — herugar5, Sep 09 '16 at 18:44
@herugar5: gcc uses the format for the targt platform. I'Äm not sure if it does support other platforms, though, or relies on the host and target to havce the same float implementation (I'm not in the mood to verify that right now). Must C targets are not x86 or even ARM. Especially for smaller MCUs, there are floating point formats which cannot represent that value. As I wrote, they are not strictily conforming. Anyway, my comment was a general one, as other readers likely don't hav exactly this value. Let that apart, YOu missed the actual point in my comment! — too honest for this site, Sep 09 '16 at 18:49

Lutz Lehmann · Accepted Answer · 2016-11-10T15:34:15.873

1

This happens when the implementation of pow is via

pow(x,y) = exp(log(x)*y)

Other library implementations first reduce the exponent by integer powers, thus avoiding this small floating point error.

More involved implementations contain steps like

pow(x,y) {
    if(y<0) return 1/pow(x, -y);

    n = (int)round(y);
    y = y-n;
    px = x; powxn = 1;
    while(n>0) {
        if(n%2==1) powxn *= px;
        n /=2; px *= px;
    }
    return powxn * exp(log(x)*y);
}

with the usual divide-n-conquer resp. halving-n-squaring approach for the integer power powxn.

edited Nov 10 '16 at 15:34

answered Sep 09 '16 at 17:59

Lutz Lehmann

25,219
2
22
51

The fuinction has no compliant declaration. As written, it returns an `int` result. And the cast to `int` is wrong. An `int` cannot represent all values a `double` can. FYI: `y` has to be an integer value anyway, see the standard. – too honest for this site Sep 09 '16 at 18:52

score 1 · Answer 2 · answered Sep 09 '16 at 18:48

1

You have a nice answer (and solution) from @LutzL, another solution is comparing the difference with an epsilon, e.g.: 0.00001, in this way you can use the standard function pow included in math.h

#define EPSILON 0.0001
#define EQ(a, b) (fabs(a - b) < EPSILON)

assert(EQ((double)c, pow(oc, opt_count)));

answered Sep 09 '16 at 18:48

David Ranieri

39,972
7
52
94

There's no good way to determine what epsilon value is appropriate, and you could easily mask an actual small difference of values. If you want to raise a value to the 2nd power, just multiply it by itself rather than calling `pow()`. – Keith Thompson Sep 09 '16 at 22:42
@KeithThompson, I agree in this case `30. * 30.` is enough, but the power value can come from an user input – David Ranieri Sep 09 '16 at 23:31

Weird C program behaviour

2 Answers2