"""Proof: the fused multiply-add really is one rounding instead of
2, and the difference is visible. Choose a, b and c so that the
exact product a*b carries a tail the separate path must round away:
  a = b = 1 + 2^-27, so a*b = 1 + 2^-26 + 2^-54 exactly
  c = -(1 + 2^-26)
The separate path rounds a*b to 1 + 2^-26 (the 2^-54 tail is under
half an ulp), then adds c and returns exactly 0.0: it reports that
a*b + c is zero when it is not. The fused path computes the exact
product, adds c, rounds once, and returns the true answer, 2^-54.
Both answers are legal for the source line "a*b + c", which is why
a compiler's fusion choice can change a program's bits.
"""
import math
import torch

a = b = 1.0 + 2.0 ** -27
c = -(1.0 + 2.0 ** -26)

separate = a * b + c
if hasattr(math, "fma"):
    fused = math.fma(a, b, c)          # Python 3.13+
    how = "math.fma"
else:
    import ctypes, ctypes.util
    libm = ctypes.CDLL(ctypes.util.find_library("m"))
    libm.fma.restype = ctypes.c_double
    libm.fma.argtypes = [ctypes.c_double] * 3
    fused = libm.fma(a, b, c)          # the C library's fma
    how = "libm fma"

print(f"a = b = 1 + 2^-27, c = -(1 + 2^-26)")
print(f"separate (a*b) + c : {separate}")
print(f"fused fma(a, b, c) : {fused}   (via {how})")
print(f"true answer        : {2.0 ** -54}")
assert separate == 0.0
assert fused == 2.0 ** -54

t = (torch.tensor(a, dtype=torch.float64) *
     torch.tensor(b, dtype=torch.float64) +
     torch.tensor(c, dtype=torch.float64)).item()
print(f"torch eager a*b+c  : {t}   (separate ops: also 0.0)")
print("\nthe separate path says the answer is exactly 0 when it is")
print("not; the fused path keeps it. both are legal readings of the")
print("same source line, which is why fusion decisions (a compiler")
print("flag, an inlined function, torch.compile) can change a")
print("program's last bits without anyone touching the math.")