IEEE 754, Errors, and Algorithms

1. Motivation

Understand floating-point representation and rounding error, and compare the accuracy of summation algorithms.

1.1 IEEE 754 representation

IEEE 754 covers number formats, sign/exponent/fraction fields, rounding, overflow/underflow, NaN, and infinity handling.

\[x = (-1)^s \times (1.f) \times 2^{e - \text{bias}}\]

Component	Role
Sign bit	0 for positive values, 1 for negative values
Exponent	biased exponent
Mantissa	Fraction bits after the binary point

1.2 Float vs Double

Type	Total bits	Sign	Exponent	Mantissa	Bias
float	32	1	8	23	127
double	64	1	11	52	1023

2. Rounding Errors

2.1 Machine epsilon, relative error, absolute error

Let $fl(x)$ denote the value representing $x$ in the computer system.

\[\epsilon_{relative} = \frac{|x - fl(x)|}{|x|}, \qquad \epsilon_{absolute} = |x - fl(x)|.\]

2.2 Error types and exceptions

• truncation error: error caused by cutting off an infinite decimal expansion.
• rounding-off error: TODO
• overflow / underflow: TODO
• error magnification: TODO
• computation error: error arising in an algorithm or iteration.

Absorption

Original claim: if $|x-y|>2^{52}$ and $y$ is smaller, then $fl(x+y)=fl(x)$.

Loss of significance

Subtracting nearby numbers can cause loss of significance. The original two expressions are

\[f_1(x)=\sqrt{x}(\sqrt{x+1}-\sqrt{x}),\qquad f_2(x)=\frac{\sqrt{x}}{\sqrt{x+1}-\sqrt{x}}.\]

The original text claims that these expressions are mathematically equal.

import decimal
from decimal import Decimal
import math 

decimal.getcontext().prec = 50


f1 = lambda x : Decimal(x).sqrt()*((Decimal(x+1)).sqrt()-(Decimal(x)).sqrt())
f2 = lambda x : Decimal(x).sqrt()/((Decimal(x+1)).sqrt()+(Decimal(x)).sqrt())


for i in range(0,54):
    print(f"at {pow(10,i):d}, \t f1(x) = {f1(pow(10,i))} \t f2(x) = {f2(pow(10,i))}")

# output at 
at 100000000000000000000000000000000000000000000000000000,      f1(x) = 0E-46   f2(x) = 0.5

2.3 Non-associativity

Rounding makes floating-point addition nonassociative.

# Example in Python
a = 1e20
b = -1e20
c = 3.14

print((a + b) + c)  # prints 3.14
print(a + (b + c))  # prints 0.0

3. Precision Algorithms

Kahan Two Sum

Original error claim: summing $n$ terms guarantees $O(n\epsilon)$.

def KahanTwoSum(x:[list:float]):
  c = 0.0 # current error 
  t = 0.0 # global sum 
  for x_ele in x:
    y = x_ele - c
    total = t + y # an error might exists + e1 
    c = (total - t) - y
    t = total 

Neumaier (1974)

This draft aims to improve on Kahan summation. The original error claim is $O(n\epsilon)$.

def Neumaier(x:[list:float]):
  c = 0.0
  t = 0.0
  for x_ele in x:
    total = t + x_ele # f(t + x_ele) = t + x_ele + $\epsilon$
    if abs(total) >= abs(x_ele):
      c += (t - total) - x_ele  
    else:
      c += (x_ele - total) + t
    t = total
  return t + c

Superaccumulator

import struct

def float_to_bits(x:float):
    
    bits = struct.unpack('>Q', struct.pack('>d',x))[0]
    signed_bit = ((bits >> 63) & 1)
    exponent = ((bits >>52) & (0x7FF) )
    mantissa = (bits & ( (1<<52) -1))
    return signed_bit, exponent, mantissa

def superacc(x_list:list[float], bins = 2048):

    acc_neg = [0.0] * bins
    acc_pos = [0.0] * bins
    offset = bins//2
    bias = 1023

    for x in x_list:
        sign_bit, expo, mantissa = float_to_bits(x)

        if expo  == 2047:
            if mantissa == 0:
                return float('inf')
            else:
                return float('nan')


        if expo == 0: #subnormal
            precision = 52
            man = mantissa
            e = 1 - bias
        # this place for NaN and Infty
        else: # normal
            precision = 53
            e = expo - bias
            man = (1<< 52) | mantissa

        for i in range(precision):
            idx = e + offset - i
            if (( man >> (precision-i-1)) & 1) == 1:
                if 0 <= idx < bins:
                    if sign_bit == 0: # positive
                        acc_pos[idx] +=1
                    else:
                        acc_neg[idx] +=1

    # in this case, the exponents are already aligned,
    # so normalization is not required.

    result = 0.0
    for i in range(bins):
        num = acc_pos[i] - acc_neg[i]
        e = i - offset
        result += (num)*2**(e)
        
    print(result)

    

Review references

• Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic: review of representation, rounding, cancellation, and compensated summation.
• Ogita–Rump–Oishi, Accurate Sum and Dot Product: review of the distinction between TwoSum, error-free transformations, and summation error.