https://gcc.gnu.org/bugzilla/show_bug.cgi?id=116319
Bug ID: 116319
Summary: std::fma should compute as if to infinite precision
and rounded only once to fit the result type.
Product: gcc
Version: 14.1.0
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: c++
Assignee: unassigned at gcc dot gnu.org
Reporter: shihyente at hotmail dot com
Target Milestone: ---
Hi,
Let's use std::bfloat16_t in <stdfloat> to create an example.
Consider 2.375^2 + 2^-133, given there is a small subnormal offset (2^-133), if
we make the output round to the nearest, then it should be 5.65625 in
std::bfloat16_t, but the result is 5.625.
We know 2.375 * 2.375 = 5.640625, which is the middle point of 5.625 and
5.65625 as (5.625 + 5.65625) / 2 = 5.640625.
A small positive offset should make it round to 5.65625 even if we make halfway
cases to even because it is not the true halfway due to infinite precision.
The command and code to reproduce it
g++ -std=c++23 -Wall -Wextra -Werror -O0 test.cpp
#include <cfenv>
#include <cstdint>
#include <cmath>
#include <bit>
#include <stdfloat>
#include <iostream>
#include <iomanip>
using namespace std;
int main() {
fesetround(FE_TONEAREST); // rounding to nearest (halfway cases to even)
const auto b = bit_cast<bfloat16_t>((uint16_t)0x4018); // 2.375
const auto c = bit_cast<bfloat16_t>((uint16_t)0x1); // 2^-133
const auto d = fma(b, b, c);
cout << static_cast<double>(d) << "\n";
cout << "0x" << hex << setfill('0') << bit_cast<uint16_t>(d) << "\n";
}
Best,
John
