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

CodePudding user response：

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.

CodePudding user response：

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

Simply speaking, consistency for x^y 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 x^y) run into a difficult issue: how much precision must internal calculations be to result in the best answer? For x^y, 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 x^y and should not yet be expected.