2

I am wondering what might explain a floating-point difference that I observed today. Is it a bug, or tripping some undefined behavior? Below is the code. My hope is to understand the behavior so that I might make it consistent across compilers, processors and platforms.

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

#define BUG_ON_GCC
#ifdef BUG_ON_GCC
__attribute__((noinline))
#else
inline __attribute__((always_inline))
#endif
double safe_pow(
  double x,
  double y)
{
  printf("x = %llx\n", *((unsigned long long *)&x));
  printf("y = %llx\n", *((unsigned long long *)&y));
  double result = pow(x, y);
  printf("r = %llx\n", *((unsigned long long *)&result));
  return result;
}

int main() {
  printf("%ld\n", sizeof(pow(15.034465284692086, 3.466120406090667)));
  printf("%.24f\n", pow(15.034465284692086, 3.466120406090667));
  printf("%ld\n", sizeof(safe_pow(15.034465284692086, 3.466120406090667)));
  printf("%.24f\n", safe_pow(15.034465284692086, 3.466120406090667));
  return 0;
}

And the results that I see using Clang in a Linux VM:

8
12020.670425990641888347454369
8
x = 402e11a56f0d331e
y = 400bba9d55e142e0
r = 40c77a55d084d419
12020.670425990641888347454369

Using GCC in the same VM:

8
12020.670425990643707336857915
8
x = 402e11a56f0d331e
y = 400bba9d55e142e0
r = 40c77a55d084d419
12020.670425990641888347454369

Using Clang on macOS on the same machine:

8
12020.670425990643707336857915
8
x = 402e11a56f0d331e
y = 400bba9d55e142e0
r = 40c77a55d084d41a
12020.670425990643707336857915
John Wiegley
  • 6,972
  • 1
  • 16
  • 20
  • 3
    Not all implementations of `pow()` are the same. Please see [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) and [Why Are Floating Point Numbers Inaccurate?](https://stackoverflow.com/questions/21895756/why-are-floating-point-numbers-inaccurate) Work from the assumption that floating point *isn't* accurate. Also, you are asking for 29 significant digits, but `double` can only hold 16-17. – Weather Vane Aug 26 '22 at 06:24
  • 1
    OT: You can use the `%a` format specifier to inspect the exact representation of a floating point value. Also note that 40c77a55d084d41a = 40c77a55d084d419 + 1. – Bob__ Aug 26 '22 at 06:33
  • Perhaps the floating point rounding mode is differently set. The results differ by `1` in the l.s. bit. This is hardly cause for concern, though, seeing as a value can be 1 bit in error anyway. – Weather Vane Aug 26 '22 at 06:47
  • I whole-heartedly agree that floating-point is the wrong thing to use, and my next fix will be to remove it entirely. My question is rather why exactly the answers disagree. – John Wiegley Aug 26 '22 at 07:19

2 Answers2

3

Exhibit A.

As you can see, there are no calls to pow, inline or otherwise. gcc and clang both compute it at compile time, and they compute it differently. The difference is in the least significant bit.

Neither of your operands is exactly representable as an IEEE double. An infinite-precision calculator suggests that GCC uses the exact operands and calculates the answer with greater precision than double offers, and then converts the result to double, while Clang converts the operands to double losing some precision, and calculates with that. There is no right or wrong behaviour, both are right according to the standard.

n. m. could be an AI
  • 112,515
  • 14
  • 128
  • 243
  • Thank you, this sums up the problem nicely. – John Wiegley Aug 26 '22 at 07:50
  • 1
    Calculating x^y with Maple either with the decimal numerals shown for x and y or with the results of converting those numerals to `double` (IEEE-754 binary64) with round-to-nearest-ties-to-even, shows the Clang result, not the GCC result, is closer to the real-number-arithmetic result of exponentiation. – Eric Postpischil Aug 26 '22 at 11:55
  • 1
    Further, I doubt that GCC is using “the exact operands” or that it is using the literals in the source text with more precision than `double`. While the compiler may use extra precision in some circumstances, such as using Intel’s 80-bit floating-point on older processors, I do not believe that extends to translating literals; they are produced in their nominal type… – Eric Postpischil Aug 26 '22 at 11:59
  • 1
    … The difference may be explained by Clang using system’s for both the compile-time exponentiation it computes with constants and the run-time exponentiation it uses for `safe_pow`. On macOS, this is a faithfully rounded `pow`, and we see good results in both the compile-time result and the run-time result. On Linux, the system has a different `pow`, and we see imperfect results in both the compile-time result and the run-time result. Then GCC uses a different `pow` (perhaps a `long double` version) for compile-time computation than for the run-time computation, resulting in different results. – Eric Postpischil Aug 26 '22 at 12:02
1

My hope is to understand the behavior so that I might make it consistent across compilers, processors and platforms.

Simply speaking, consistency for xy obliges a tolerance more than ±0.5000... ULP.

What are the pow(x,y) operand values?

Neither 15.034465284692086 nor 3.466120406090667 are exactly representable as double, instead nearby values are use for x, y

//                          15.0344652846920 86
x = 0x1.e11a56f0d331ep+3    15.0344652846920 85873530231765471398830413818359375
y = 0x1.bba9d55e1420ep+1    3.46612040609066 69610631070099771022796630859375
//                          3.46612040609066 7

What is pow(x,y)?

Mathematically, 15.0344652846920 858735...3.46612040609066 696106... is:

0x1.77a55d084d419 802...p+13  12020.67042599064 2798519464085... extended math

and pow(x,y) is:

0x1.77a55d084d419       p+13  12020.67042599064 1888347454369... Clang Linux
0x1.77a55d084d419       p+13  12020.67042599064 1888347454369... GCC Linux
0x1.77a55d084d41a       p+13  12020.67042599064 3707336857915... Clang MAC

0x1.77a55d084d419 80c   p+13  12020.67042599064 2803171226660... powl(x, y) double args
0x1.77a55d084d41a 640...p+13  12020.67042599064 4417576265463... extended math with code's constants

The two different pow() answers are 1.0 ULP apart.
The math answer is about at the 50.05% point between the 2 possible double results.
With wider powl(), the result if still just over 50% ULP if the double result.
The larger of the two double, (Clang MAC), is the better answer in this case.
Even using the original code constants, Clang MAC is better.

Why does calling pow() in different contexts sometime yield different results?

Review Table-maker's dilemma

Transcendental calculations like pow() (which seeks to approximate math xy) run into a difficult issue: how much precision must internal calculations be to result in the best answer? For xy, the usual approach is to perform the calculation with some extra, but not too many, as that slows the performance yet only marginally increases occurrence rate of the best answer with each extra digit.

It is simply that some pow() calculations perform better in the nearly-half-way-cases than the other due to using more precision, better algorithms, or select vales (the other approach may more often provide a better answer - this is just 1 case.)

Consistency

What should be expected with good math.h library are results that are within 1.0 ULP of the best answer. Here, both answers were within 0.5005 ULP of the best answer. To always expect an answer within 0.5000... ULP is an ever increasing and challenging problem for xy and should not yet be expected.

chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256