"""Proof: the rounding error is catchable, exactly. 3 facts, each
checked against exact rational arithmetic (fractions.Fraction):
1. Fast2Sum: for |a| >= |b|, 3 operations recover the rounding
   error of a + b perfectly: s = a+b; z = s-a; e = b-z. Then
   s + e equals the true sum of the stored a and b, exactly.
   Run on the page's opening example, 0.1 + 0.2.
2. Sterbenz: when x and y are within a factor of 2 of each
   other, x - y is computed exactly: cancellation commits no
   error of its own; it only exposes error the inputs carried.
3. 2MultFMA: 1 fused multiply-add recovers a product's whole
   rounding error: e = fma(a, b, -a*b).
"""
from fractions import Fraction as F
import ctypes, ctypes.util

# ---- 1. Fast2Sum on 0.1 + 0.2
a, b = 0.2, 0.1
s = a + b
z = s - a
e = b - z
print(f"s = a + b        = {s!r}")
print(f"e = (b - (s - a)) = {e!r}")
assert F(a) + F(b) == F(s) + F(e)
print("checked: stored 0.2 + stored 0.1 == s + e, EXACTLY.")
print("the 3 float operations caught the rounding error whole.\n")

# ---- 2. Sterbenz: close subtraction is exact
x, y = 1.0000001, 1.0
d = x - y
assert F(x) - F(y) == F(d)
print(f"x - y = {d!r}")
print("checked: equals the exact difference of the stored inputs.")
print("the 19% error of the cancellation section was already in")
print("the inputs; the subtraction added nothing.\n")

# ---- 3. 2MultFMA: a product's error in one instruction
libm = ctypes.CDLL(ctypes.util.find_library("m"))
libm.fma.restype = ctypes.c_double
libm.fma.argtypes = [ctypes.c_double] * 3
u = v = 1.0 + 2.0 ** -27
p = u * v
e2 = libm.fma(u, v, -p)
assert F(u) * F(v) == F(p) + F(e2)
print(f"p = u * v       = {p!r}")
print(f"e = fma(u,v,-p) = {e2!r}")
print("checked: stored u times stored v == p + e, exactly.")
print("\nKahan's compensated sum is fact 1 run in a loop; the")
print("double-word arithmetic that stretches precision in software")
print("is facts 1 and 3 kept instead of thrown away.")