What is happening here in pow function?
Asked Answered
C

5

13

I have seen various answer here that depicts Strange behavior of pow function in C.
But I Have something different to ask here.

In the below code I have initialized int x = pow(10,2) and int y = pow(10,n) (int n = 2).

In first case it when I print the result it shows 100 and in the other case it comes out to be 99.

I know that pow returns double and it gets truncated on storing in int, but I want to ask why the output comes to be different.

CODE1

#include<stdio.h>
#include<math.h>
int main()
{
     int n = 2;
     int x;
     int y;
     x = pow(10,2);   //Printing Gives Output 100   
     y = pow(10,n);   //Printing Gives Output 99


     printf("%d %d" , x , y);

}

Output : 100 99

Why is the output coming out to be different. ?

My gcc version is 4.9.2

Update :

Code 2

int main()
    {
         int n = 2;
         int x;
         int y;
         x = pow(10,2);   //Printing Gives Output 100   
         y = pow(10,n);   //Printing Gives Output 99
         double k = pow(10,2);
         double l = pow(10,n);


         printf("%d %d\n" , x , y);
         printf("%f %f\n" , k , l);

    }

Output : 100 99
100.000000 100.000000

Update 2 Assembly Instructions FOR CODE1

Generated Assembly Instructions GCC 4.9.2 using gcc -S -masm=intel :

    .LC1:
    .ascii "%d %d\0"
    .text
    .globl  _main
    .def    _main;  .scl    2;  .type   32; .endef
_main:
    push    ebp
    mov ebp, esp
    and esp, -16
    sub esp, 48
    call    ___main
    mov DWORD PTR [esp+44], 2
    mov DWORD PTR [esp+40], 100      //Concerned Line
    fild    DWORD PTR [esp+44]
    fstp    QWORD PTR [esp+8]
    fld QWORD PTR LC0
    fstp    QWORD PTR [esp]
    call    _pow                    //Concerned Line
    fnstcw  WORD PTR [esp+30]
    movzx   eax, WORD PTR [esp+30]
    mov ah, 12
    mov WORD PTR [esp+28], ax
    fldcw   WORD PTR [esp+28]
    fistp   DWORD PTR [esp+36]
    fldcw   WORD PTR [esp+30]
    mov eax, DWORD PTR [esp+36]
    mov DWORD PTR [esp+8], eax
    mov eax, DWORD PTR [esp+40]
    mov DWORD PTR [esp+4], eax
    mov DWORD PTR [esp], OFFSET FLAT:LC1
    call    _printf
    leave
    ret
    .section .rdata,"dr"
    .align 8
LC0:
    .long   0
    .long   1076101120
    .ident  "GCC: (tdm-1) 4.9.2"
    .def    _pow;   .scl    2;  .type   32; .endef
    .def    _printf;    .scl    2;  .type   32; .endef
Copaiba answered 10/2, 2017 at 16:48 Comment(29)
I have no problem ideone.com/sbaCitArelus
I get 100 both times. What is your build command?Nerve
@EugeneSh. Please read it again , it has different answer n = 2 pow(10,n) and pow(10,2) , in first it is 99 and in second with 100.Copaiba
Cannot reproduce (VS2013 on x64).Puppy
Compiling it with Apple LLVM version 8.0.0 (clang-800.0.42.1), I obtain 100 and 100. With and without #include<math.h>.Chickamauga
Can you print the resulting double values and see any differences? Just as an experiment.Bioscope
BTW, you can look for the source code of pow corresponding to your gcc version, duplicate it and run with the debugger to see what's going on.Bioscope
@Copaiba Just googled version of your compiler and it looks like it can be a bug (if you use optimisation) See my answer for detailsMercurio
It's more likely glibc (or your C standard library) than gcc. gcc only generates the call to pow and links to it; the library itself is the likely cause of rounding errors.Jerboa
The exact same problem had already been found here: #19127309, but with another compiler. Since the compiler is so important in this case, maybe this question shouldn't be considered a duplicate.Sapient
@FabioTurati Actually looks like exact duplicate, but I would rather close that one as a duplicate of this one, as it is higher quality...Bioscope
Simple solution: don't mix floating point functions with integers if you need exact results.Flawy
@Olaf Did you read the question, i am asking why it is happening what happens behind the scene atleast read question first..Copaiba
This question has an answer here: #32353801. The C standard says the precision of math functions like pow() is implementation defined. The glibc library fixed bugs due to rounding mode recently, pow(1.01,1.1) would crash (!) in one case. Discussed here: sourceware.org/bugzilla/show_bug.cgi?id=3976 Fixed here: github.com/lattera/glibc/commit/….Mok
typo... previous comment should read "accuracy" not "precision"Mok
@Mok Sir , What does this mean std::pow<int, int>(int, int)Copaiba
The answer to your question about what does __gnu_cxx::__promote_2<int, int, __gnu_cxx::__promote<int, std::__is_integer<int>::__value>::__type, __gnu_cxx::__promote<int, std::__is_integer<int>::__value>::__type>::__type std::pow<int, int>(int, int) mean is roughly that the implementation calls a template function and converts (promotes) your arguments from their original type to the type required by the std::pow<int, int>(int, int). It isn't immediately obvious why the conversions are deemed necessary.Straightout
@JonathanLeffler That means there is pow function that takes int argument , how can it be ? What is meant by template function ?Copaiba
Amongst other things, it means you're compiling your code with a C++ compiler, which makes the C tag on the question odd. How are you abusing the system so that you're compiling it as C++? You must either be using g++ as the compiler, or possibly using gcc but specifying a file name with an extension like .cpp or .cxx — or perhaps .C (upper-case letter C) as the extension, and gcc is interpreting that as a C++ source file and running appropriate compiler phases for C++ instead of C.Straightout
I have given link, i have compiled using gcc 6.3 in godbolt.orgCopaiba
The IDE says it is C++ source, so your code is being compiled as C++.Straightout
@JonathanLeffler I noticed now , godbolt is doing that , waitCopaiba
@JonathanLeffler I dont have gcc 6.3 , but my friends codeblock has it , maybe that is why his pow function is giving correct answer , why is pow giving wrong answer in pow(10,n) and pow(10,2) it is giving 100. How can i generate assembly online for C, ?Copaiba
You should read the documentation, but the answer is gcc -S.Straightout
@JonathanLeffler I used this answer https://mcmap.net/q/20237/-using-gcc-to-produce-readable-assembly , and now godbolt is treating code as C godbolt.org/g/LCslhxCopaiba
Re. "update 3", it means you used a C++ compiler, which is irrelevant and should be removed since this is a C-only question. You should also remove the C++ parts of "update 2".'Reformer
@Reformer Done , Sir , can you tell me how to generate the complete floating point number of pow(10,n) just so that i could see where it fell shortCopaiba
Same Problem Here #18156383Sentimentalism
Compilers can constant fold away calculations sometimes and this can lead to differences compared to when the calculation is done at run-time, usually using -fno-builtin will disable. Also see log(10.0) can compile but log(0.0) cannot? and Inconsistent strcmp() return value when passing strings as pointers or as literalsBisset
J
9

I know that pow returns double and it gets truncated on storing in int, but I want to ask why the output comes to be different.

You must first, if you haven't already, divest yourself of the idea that floating-point numbers are in any way sensible or predictable. double only approximates real numbers and almost anything you do with a double is likely to be an approximation to the actual result.

That said, as you have realized, pow(10, n) resulted in a value like 99.99999999999997, which is an approximation accurate to 15 significant figures. And then you told it to truncate to the largest integer less than that, so it threw away most of those.

(Aside: there is rarely a good reason to convert a double to an int. Usually you should either format it for display with something like sprintf("%.0f", x), which does rounding correctly, or use the floor function, which can handle floating-point numbers that may be out of the range of an int. If neither of those suit your purpose, like in currency or date calculations, possibly you should not be using floating point numbers at all.)

There are two weird things going on here. First, why is pow(10, n) inaccurate? 10, 2, and 100 are all precisely representable as double. The best answer I can offer is that the C standard library you are using has a bug. (The compiler and the standard library, which I assume are gcc and glibc, are developed on different release schedules and by different teams. If pow is returning inaccurate results, that is probably a bug in glibc, not gcc.)

In the comments on your question, amdn found a glibc bug to do with FP rounding that might be related and another Q&A that goes into more detail about why this happens and how it's not a violation of the C standard. chux's answer also addresses this. (C doesn't require implementation of IEEE 754, but even if it did, pow isn't required to use correct rounding.) I will still call this a glibc bug, because it's an undesirable property.

(It's also conceivable, though unlikely, that your processor's FPU is wrong.)

Second, why is pow(10, n) different from pow(10, 2)? This one is far easier. gcc optimizes away function calls for which the result can be calculated at compile time, so pow(10, 2) is almost certainly being optimized to 100.0. If you look at the generated assembly code, you will find only one call to pow.

The GCC manual, section 6.59 describes which standard library functions may be treated in this way (follow the link for the full list):

The remaining functions are provided for optimization purposes.

With the exception of built-ins that have library equivalents such as the standard C library functions discussed below, or that expand to library calls, GCC built-in functions are always expanded inline and thus do not have corresponding entry points and their address cannot be obtained. Attempting to use them in an expression other than a function call results in a compile-time error.

[...]

The ISO C90 functions abort, abs, acos, asin, atan2, atan, calloc, ceil, cosh, cos, exit, exp, fabs, floor, fmod, fprintf, fputs, frexp, fscanf, isalnum, isalpha, iscntrl, isdigit, isgraph, islower, isprint, ispunct, isspace, isupper, isxdigit, tolower, toupper, labs, ldexp, log10, log, malloc, memchr, memcmp, memcpy, memset, modf, pow, printf, putchar, puts, scanf, sinh, sin, snprintf, sprintf, sqrt, sscanf, strcat, strchr, strcmp, strcpy, strcspn, strlen, strncat, strncmp, strncpy, strpbrk, strrchr, strspn, strstr, tanh, tan, vfprintf, vprintf and vsprintf are all recognized as built-in functions unless -fno-builtin is specified (or -fno-builtin-function is specified for an individual function).

So it would seem you can disable this behavior with -fno-builtin-pow.

Jerboa answered 10/2, 2017 at 17:23 Comment(22)
As you pointed in the comments, the math library is precompiled, so I doubt pow(10,2) can be computed in the compile-time, as the compiler is just not aware of the implementation of pow.Bioscope
@EugeneSh. Try it yourself. Compiling with gcc -S, only one call to pow is in the assembly output. My x86 isn't quite good enough to figure out what happens to the other one, but I assume it gets optimized to the equivalent of 100.0.Jerboa
@trentcl It is completely legit theoretically. But I wonder how the compiler knows what it should be optimized to not knowing what pow is doing in practice... Unless the pow function is hardcoded in the compiler itself.Bioscope
@trentcl Aha! Now you have nailed it. I would give you +5 if I could :)Bioscope
Wait a second - does this mean that when you call pow(10, 2) the compiler understands that it can be optimized away and replaced with the result (so far so good), but then to calculate the result it uses its own implementation of pow, which is built into it, and might be different from the one in the standard library which is called at run time for pow(10, n), as in this case? It sounds crazily dangerous! I hope the fact that the 2 libraries aren't the same in that version of gcc is recognized as a bug!Sapient
@FabioTurati Apparently, the danger of such an optimization is exactly demonstrated by the question :)Bioscope
[...] because squaring 10 should not invoke floating-point rounding errors. Why not? pow(x, y) is typically implemented as exp(log(x) * y) where exp is e^x and log is the natural logarithm.Alurta
+1 @trentcl excellent answer. I've added a comment to the OP regarding a bug fix to glibc, feel free to include in your answer.Mok
@Alurta A fair point. I meant "should" in a more general sense. The C standard technically allows 2.0 + 2.0 to equal 4.000001, but surely it shouldn't! However, even in IEEE 754, it's not required for pow to have correct rounding as for + - * / sqrt, so you are right. (I'm still pretty sure this is a bug though.)Jerboa
@trentcl Can You See The Updated Question what does this mean std::pow<int, int>(int, int)Copaiba
@trentcl: The fact that the exponent 2.0 "is precisely representable as a double" is totally irrelevant. 0.5 is also precisely representable as a double, but I think we would all be surprised if pow(2.0, 0.5) returned the precisely correct value. :) While some standard libraries do check to see if the exponent is a small integer and attempt to use that information to minimize inaccuracy (even making the result precise), pow is under no obligation to do so. A good implementation will probably do better than pow(a, b) => exp(a, log(b)), but that would be legal, too....Antenatal
... You are 100% correct that the compiler is allowed to compute standard library functions at compile time if it can deduce the arguments, and the standard explicitly allows compile-time computations to be different from run-time computations, even for ordinary floating-point arithmetic. The difference should be minor, maybe only an ULP, but it can happen. The real problem with this code is that the floor operation (and implicit double->integer conversion) introduces numerical instability, and should be avoided...Antenatal
... floating-point arithmetic is underspecified in the C standard, but IEEE-749 does a much better job and the requirements it provides are reasonable and sufficient for programmers to analyze code and deduce consequences. If you start with the hypothesis that floating point arithmetic is inherently unpredictable, you will be much less likely to learn how to predict it. :)Antenatal
Floating point numbers are predictable. If they weren't then they would be completely useless.Reformer
@trentcl Does your answer means the same https://mcmap.net/q/162756/-i-got-different-results-with-pow-10-2-and-pow-10-j-j-2 , it is talking about constant folding ?Copaiba
@CuriosGuy: "constant folding" is another term for compile-time computation, so that is the same thing. The standard library is part of the standard, and users are not allowed to define functions with the same name as standard library functions if they have #include'd the standard headers. So the compiler is allowed to compute standard math functions at compile time.Antenatal
@Antenatal So the compiler has another implementation of pow, the answer we are commenting on is talking about built in pow function.Copaiba
What does this line mean " GCC built-in functions are always expanded inline and thus do not have corresponding entry points and their address cannot be obtained. Attempting to use them in an expression other than a function call results in a compile-time error."Copaiba
@Copaiba You should ask that as a separate question!Jerboa
@Antenatal From a C perspective, floating-point math is inherently unpredictable, since it may do different things on different processors (theoretically, even without recompilation). For a given stackup, assuming IEEE 754 and only using + - * / sqrt, sure, that's predictable.Jerboa
@trentcl You have this line in your answer , and i did not understand that "GCC built-in functions are always expanded inline and thus do not have corresponding entry points and their address cannot be obtained. Attempting to use them in an expression other than a function call results in a compile-time error."Copaiba
@curiosguy: i thought I answered that but my comment either vanished or never was registered. Anyway, it means what it says: you cannot take the address of an inline function because it doesn't exist as an entity. This has nothing to do with your question and tentcl's quote of that line is also irrelevant.Antenatal
H
2

Why is the output coming out to be different. ? (in the updated appended code)

We do not know the values are that different.

When comparing the textual out of int/double, be sure to print the double with sufficient precision to see if it is 100.000000 or just near 100.000000 or in hex to remove all doubt.

printf("%d %d\n" , x , y);
// printf("%f %f\n" , k , l);
// Is it the FP number just less than 100?
printf("%.17e %.17e\n" , k , l);  // maybe 9.99999999999999858e+01
printf("%a %a\n" , k , l);        // maybe 0x1.8ffffffffffff0000p+6

Why is the output coming out to be different. ? (in the original code)

C does not specify the accuracy of most <math.h> functions. The following are all compliant results.

// Higher quality functions return 100.0
pow(10,2) --> 100.0   
// Lower quality and/or faster one may return nearby results
pow(10,2) --> 100.0000000000000142...
pow(10,2) --> 99.9999999999999857...

Assigning a floating point (FP) number to an int simple drops the fraction regardless of how close the fraction is to 1.0

When converting FP to an integer, better to control the conversion and round to cope with minor computational differences.

// long int lround(double x);
long i = lround(pow(10.0,2.0));
Horseplay answered 10/2, 2017 at 17:43 Comment(11)
@Copaiba In CODE2 what is the text output with changing printf("%f %f\n" , k , l); to printf("%.17e %.17e\n" , k , l); for you?Horseplay
1.00000000000000000e+002 1.00000000000000000e+002Copaiba
@Copaiba Hmmm. To eliminate that even higher precision might be used "under the hood" for intermediate values, try printf("%d %d\n", l == 100.0, pow(10,n) == 100.0) to see if they is truly 100.0.Horseplay
1.0 is coming iCopaiba
It is coming to be 1.0Copaiba
Sir if it is coming out to be 1 , so what is happening there?Copaiba
@Copaiba Hmmm, result of 1 means a new explanation is needed. Consider y = pow(10,n); double l = pow(10,n); --> changed to l = pow(10,n); y = l; printf("FLT_EVAL_METHOD %d\n", FLT_EVAL_METHOD);. What is your result? C allows intermediate FP to evaluate at higher precision depending on FLT_EVAL_METHOD. If it is 2, than pow(10,n); may be evaluating to long double 99.99999... When converted to int, that is 99 (fraction truncation), when converted to double, that is exactly 100.0 (round to nearest double). In this comment's code, the value is saved as double then to intHorseplay
@Copaiba int y; as in the post.Horseplay
It is saying FLT_EVAL_METHOD not available i used float.hCopaiba
What version of C and what compiler are you using?Horseplay
I am using mingw gcc version 6.2Copaiba
P
0

You're not the first to find this. Here's a discussion form 2013: pow() cast to integer, unexpected result

I'm speculating that the assembly code produced by the tcc guys is causing the second value to be rounded down after calculating a result that is REALLY close to 100. Like mikijov said in that historic post, looks like the bug has been fixed.

Puparium answered 10/2, 2017 at 17:3 Comment(4)
So, it's something with your compiler or your standard library. I could tr to recreate the problem on my university's cse machines... would that help?Puparium
He says it's GCC 4.9, which is close to 3 years old now.Kendre
I could not reproduce the error, I did not modify your code. I am using gcc version 4.9.4 you gotta update. the programmers who make gcc made a bit of a mistake.Puparium
Please see the question again.Copaiba
V
0

As others have mentioned, Code 2 returns 99 due to floating point truncation. The reason why Code 1 returns a different and correct answer is because of a libc optimization.

When the power is a small positive integer, it is more efficient to perform the operation as repeated multiplication. The simpler path removes roundoff. Since this is inlined you don't see function calls being made.

Vicarage answered 11/2, 2017 at 9:37 Comment(0)
M
-1

You've fooled it into thinking that the inputs are real and so it gives an approximate answer, which happens to be slightly under 100, e.g. 99.999999 that is then truncated to 99.

Maguire answered 10/2, 2017 at 17:1 Comment(3)
Please read question again , I know about floating inaccuracy thing it is different question .Copaiba
If you give pow an integer constant, it obviously special-cases it and inlines or somehow calculates with exact arithmetic. If you pass an integer variable, it will be promoted to a double, and must go down the approximation route.Maguire
still your answer is short of many things.. update it to contain all you know about topic including this comment it tells a lot.. and was this question bad... 3 downvotsCopaiba

© 2022 - 2024 — McMap. All rights reserved.