pow() cast to integer, unexpected result

Question

I have some problems using an integer cast for the pow() function in the C programming language. The compiler I'm using is the Tiny C Compiler (tcc version 0.9.24) for the Windows platform. When executing the following code, it outputs the unexpected result 100, 99:

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

int main(void)
{
    printf("%d, ", (int) pow(10, 2));
    printf("%d", (int) pow(10, 2));
    return 0;
}

However, at this online compiler the output is as expected: 100, 100. I don't know what is causing this behavior. Any thoughts? Programming error from me, compiler bug?

possible duplicate of [What is a simple example of floating point/rounding error?](http://stackoverflow.com/questions/249467/what-is-a-simple-example-of-floating-point-rounding-error) — , Aug 24 '13 at 10:58
Guys, why is the first value correct? Shouldn't it be `99, 99` if it was the case for the usual imprecise-then-truncate issue? — ppeterka, Aug 24 '13 at 11:02
@ppeterka Conspiracy theory #183742: `printf("%d", some_integer)` is constant-folded at compile time. If the constant-folding algorithm in the compiler is defectious, then the code may very well be changed to `puts("99");`. The `pow()` implementation seems to be honest and correct (in the sense that it pays attention to integer powers), though. But we would **really** need the assembly the compiler generated to prove this. — , Aug 24 '13 at 11:05
@H2CO3 same result (`100`, `99`) with `tcc` with `printf("%d, %d", (int) pow(10, 2), (int) pow(10, 2));`. This does not answer why the two calls results are treated differently in `tcc`. — ouah, Aug 24 '13 at 11:09
@ouah: and with `printf("%d, %d, ", (int) pow(10, 2), (int) pow(10, 2));`? (added `, ` to second `%d` too) — ppeterka, Aug 24 '13 at 11:16
I'm guessing it is as @H2CO3 said: some sort of folding issue. Unfortunately, why it happens is unknown since the best you can do is disassemble the resulting object file because the Tiny C Compiler doesn't output any assembler code, which presumably is a part of the reason why it is so fast (writing opcode octets with values is faster than writing text instructions). — , Aug 24 '13 at 11:40
where do you get your math lib? I do `tcc test.c -lm` and get `100 100` — Antti Haapala -- Слава Україні, Aug 24 '13 at 12:08
or is this implicit? My tcc 0.9.25 on linux does not even link in pow by itself. — Antti Haapala -- Слава Україні, Aug 24 '13 at 12:09
@H2CO3: No, it's not a programming error on his side. It's arguably a programming error in his standard library. — tmyklebu, Aug 24 '13 at 21:52
@tmyklebu Apparently, yes. Rare exception. I should have judged more carefully... — , Aug 24 '13 at 21:56
I don't have to link the math library on the Windows platform as far as I know. Everything else (i.e. all other mathematical operations) works as expected. — Jori, Aug 25 '13 at 09:59

masoud · Answer 1 · 2013-08-24T23:33:01.380

6

Some investigation in assembly code. (OllyDbg)

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

int main(void)
{
    int x1 = (int) pow(10, 2);
    int x2 = (int) pow(10, 2);
    printf("%d %d", x1, x2);
    return 0;
}

The related assembly section:

FLD QWORD PTR DS:[402000]   // Loads 2.0 onto stack
SUB ESP,8
FSTP QWORD PTR SS:[ESP]
FLD QWORD PTR DS:[402008]   // Loads 10.0 onto stack
SUB ESP,8
FSTP QWORD PTR SS:[ESP]
CALL <JMP.&msvcrt.pow>      // Calls pow API
                            // Returned value 100.00000000000000000
...


FLDCW WORD PTR DS:[402042]  //   OH! LOOK AT HERE
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
...

FLD QWORD PTR DS:[402010]   // Loads 2.0 onto stack
SUB ESP,8
FSTP QWORD PTR SS:[ESP]
FLD QWORD PTR DS:[402018]   // Loads 10.0 onto stack
SUB ESP,8
FSTP QWORD PTR SS:[ESP]
CALL <JMP.&msvcrt.pow>      // Calls pow API again
                            // Returned value 99.999999999999999990

The generated code for two calls is the same, but the outputs are different. I don't know why tcc put FLDCW there. But the main reason of two different values are that line.

Before that line the round Mantissa Precision Control Bits is 53bit (10), but after execution of that line (it loads FPU register control) it will set to 64bits (11). On the other hand Rounding Control is nearest so 99.999999999999999990 is the result. Read more...

enter image description here

Solution:

After using (int) to cast a float to an int, you should expect this numeric error, because this casting truncates the values between [0, 1) to zero.

Assume the 10² is 99.9999999999. After that cast, the result is 99.

Try to round the result before casting it to integer, for example:

printf("%d", (int) (floor(pow(10, 2) + 0.5)) );

edited Aug 24 '13 at 23:33

answered Aug 24 '13 at 18:57

masoud

55,379
16
141
208

1

How are you inspecting the values 100.000... and 99.999...90? Printing floats is non-trivial, and that could possibly be stateful/buggy too. – Quadrescence Aug 24 '13 at 20:23
1

@Quadrescence: There is a panel in OllyDbg (Registers FPU) which inspects the values. _(I added an image)_ – masoud Aug 24 '13 at 21:05
Yes, but when I round, how can I be sure that it isn't `floor`'d to 99.9999 and is still truncated to 99 by the integer conversion? As we just saw that 100 is represented by 99.9999. – Jori Aug 25 '13 at 10:28
It's the reason we add `0.5` before `floor` as written in the above code. – masoud Aug 25 '13 at 12:03
It is not necessary to round exactly to the nearest integer to fix the problem here, fortunately, but adding `0.5` does not provide rounding to the nearest. http://blog.frama-c.com/index.php?post/2013/05/02/nearbyintf1 – Pascal Cuoq Aug 25 '13 at 12:48
@MM. Sorry, I think my question was unclear. What I really wanted to ask is, can we, as programmers, be sure that when we `floor` a floating point value and then cast it to an `int`, that it always gives the `int` we wanted. I can't see why this holds, as we just saw that the nearest floating point representation of `100.0` is `99.9999`, and thus is truncated to `99`. `floor` could only possibly work when the implementation has this requirement in mind, and thus doesn't round to the nearest fp representation, but to the nearest fp representation that is larger than the int value. Am I correct? – Jori Aug 25 '13 at 14:26
Also, could you provide a link with some basic explanation of how to use OllyDbg, I'm interested how to obtain you results. Thanks in advance! – Jori Aug 25 '13 at 14:27
@Jori: Just download [OllyDbg](http://www.ollydbg.de/download.htm) open the .exe file and use F8 and F7 to trace. ;-) – masoud Aug 25 '13 at 14:33
@MM. Please look at my previous post. I think you missed it ;) – Jori Aug 26 '13 at 08:41
@Jori: It is generally not possible to ensure that `floor(x)` gives the `int` programmers want, because “what programmers want” is not a definition of a function that can be computed and because different programs want different things. If you want a function that computes what you “want”, you must define what you want, preferably mathematically, but at least in some clear and specific form. – Eric Postpischil Aug 26 '13 at 13:06
2

@Jori: It is clear that in `floor(pow(10, 2))`, you would want the result to be 100. However, `floor` is not passed `pow(10, 2)`. It receives only a number, and it cannot determine intent from this number. When `floor` receives 99.9999…, it must return 99; there is no other correct choice. The actual problem here is that `pow` is failing; the exact mathematical result of `pow(10, 2)` is representable; it is 100; but `pow` is returning a different value. It is `pow` that is broken (or the C implementation generally), not `floor`. – Eric Postpischil Aug 26 '13 at 13:07
Yes I understand that, but what I meant is that if we called `floor(100.5)` and casted it to an int is doesn't give `99`, as `floor` returns a fp and the nearest fp representation of `100` is `99.99999`, but `100`. – Jori Aug 26 '13 at 15:36
And does the same hold for `round` (e.g. `(int) round(a)` is always `int(a)` or `int(a+1)`)? – Jori Aug 26 '13 at 15:41

score 2 · Accepted Answer · answered Aug 28 '13 at 21:14

You found a bug in tcc. Thanks for that. The patch has just been commited to the repository. It will be included in the next release, but that might take a while. You can of course pull the source and build it yourself. The patch is here

http://repo.or.cz/w/tinycc.git/commitdiff/73faaea227a53e365dd75f1dba7a5071c7b5e541

score 0 · Answer 3 · edited Aug 24 '13 at 20:38

It seems that the rounding method may change, thus requiring an ASM instruction finit to reset the FPU. In FreeBASIC on Windows, I get 99.9999 even in the first try, and so I think for you after the first try it would be a consistent 99.9999. (But I call this undefined behavior indeed, more than a bug in the C runtime's pow().)

So my advice is not do the conversion with round down. To avoid such issues, use, for example:

int x1 = (int)(pow(10, 2)+.5);

pow() cast to integer, unexpected result

3 Answers3

Linked

Related