Floating Point Arithmetic

Floating Point Arithmetic

Definitions

Formulas

IEEE 754 Deep Dive

Double Precision (float64) — 64 bits total

Exponent special values:

• e = 0, m = 0: ±0
• e = 0, m ≠ 0: Denormalized numbers
• e = 2047, m = 0: ±Infinity
• e = 2047, m ≠ 0: NaN (Not a Number)

Range of float64:

• Minimum positive normal: 2^(−1022) ≈ 2.2 × 10^(−308)
• Maximum finite: (2 − 2^(−52)) × 2^1023 ≈ 1.8 × 10^308

Half Precision (float16) — Used in ML Training

• Range: ±65504
• Smallest positive normal: 6.1 × 10^(−5)
• Machine epsilon: 9.77 × 10^(−4)
• Used extensively in deep learning for memory savings and speed

Single Precision (float32)

• Machine epsilon: ≈ 1.19 × 10^(−7)
• Range: ±3.4 × 10^38
• Common default for ML inference

Examples

Theorems

Python Implementation

import numpy as np
import struct
import sys


def inspect_float64(value):
    """Break down a float64 into its IEEE 754 components."""
    # Pack as binary (8 bytes) and unpack as unsigned int
    bits = struct.unpack('Q', struct.pack('d', value))[0]

    sign = (bits >> 63) & 1
    exponent = ((bits >> 52) & 0x7FF) - 1023  # Remove bias
    mantissa = bits & 0xFFFFFFFFFFFFF          # Lower 52 bits

    print(f"Value: {value}")
    print(f"  Binary: {bits:064b}")
    print(f"  Sign:     {sign} ({'negative' if sign else 'positive'})")
    print(f"  Exponent: {exponent + 1023} (unbiased: {exponent})")
    print(f"  Mantissa: {mantissa:052b}")
    print(f"  Implicit leading 1: {1 + mantissa / (1 << 52)}")
    print()


def demonstrate_epsilon():
    """Show machine epsilon and its effects."""
    eps = np.finfo(float).eps
    print(f"Machine epsilon (float64): {eps}")
    print(f"  = 2^(-52) = {2**-52}")

    # Find the next representable float after 1.0
    one_plus = np.nextafter(1.0, 2.0)
    print(f"Next float after 1.0: {one_plus}")
    print(f"  = 1 + {one_plus - 1.0}")
    print(f"  = 1 + 2^(-52) = 1 + {2**-52}")
    print()


def demonstrate_rounding():
    """Show how different rounding modes affect results."""
    # 0.1 in binary is a repeating fraction
    print("Rounding of 0.1 in float64:")
    print(f"  fl(0.1) = {np.float64(0.1):.20f}")
    print(f"  Error:   {np.float64(0.1) - 0.1:.2e}")
    print()

    # Show precision loss in subtraction
    a = 1.0
    b = 1.0 + 1e-15
    print(f"Precision loss demonstration:")
    print(f"  1.0 + 1e-15 = {b:.20f}")
    print(f"  (1.0 + 1e-15) - 1.0 = {b - a:.20f}")
    print(f"  Expected: 1e-15 = {1e-15:.20f}")
    print()


def demonstrate_accumulation_error():
    """Show error accumulation in repeated operations."""
    # Add 0.1 ten times
    result = 0.0
    for i in range(10):
        result += 0.1

    print(f"Summing 0.1 ten times: {result}")
    print(f"  Expected: 1.0")
    print(f"  Error: {result - 1.0:.2e}")
    print(f"  Direct: 10 * 0.1 = {10 * 0.1}")
    print()


def demonstrate_denormals():
    """Show denormalized number behavior."""
    tiny = np.finfo(np.float64).tiny
    print(f"Smallest normal: {tiny}")
    print(f"Smallest subnormal: {np.finfo(np.float64).smallest_subnormal}")
    print()

    # Gradual underflow
    x = tiny
    step = 0
    while x > 0 and step < 60:
        x /= 2
        step += 1
        if step <= 5 or x == 0:
            print(f"  Divide by 2 (step {step}): {x}")
    print()


def demonstrate_special_values():
    """Show IEEE 754 special values."""
    print("Special values:")
    print(f"  +inf: {np.inf}")
    print(f"  -inf: {-np.inf}")
    print(f"  NaN: {np.nan}")
    print(f"  inf + inf = {np.inf + np.inf}")
    print(f"  inf - inf = {np.inf - np.inf}")  # NaN
    print(f"  0/0 = {0.0 / 0.0}")              # NaN
    print(f"  1/0 = {1.0 / 0.0}")              # inf
    print(f"  NaN == NaN: {np.nan == np.nan}")  # False!
    print(f"  isnan(NaN): {np.isnan(np.nan)}")  # True
    print()


def kahan_summation(values):
    """Kahan compensated summation algorithm."""
    total = 0.0
    compensation = 0.0

    for value in values:
        y = value - compensation
        temp = total + y
        compensation = (temp - total) - y
        total = temp

    return total


if __name__ == "__main__":
    print("=" * 60)
    print("IEEE 754 FLOATING POINT ANALYSIS")
    print("=" * 60)

    # Inspect some values
    print("\n--- Value Inspection ---")
    inspect_float64(0.1)
    inspect_float64(0.2)
    inspect_float64(0.1 + 0.2)

    # Epsilon
    print("\n--- Machine Epsilon ---")
    demonstrate_epsilon()

    # Rounding
    print("\n--- Rounding Effects ---")
    demonstrate_rounding()

    # Accumulation
    print("\n--- Accumulation Error ---")
    demonstrate_accumulation_error()

    # Denormals
    print("\n--- Denormalized Numbers ---")
    demonstrate_denormals()

    # Special values
    print("\n--- Special Values ---")
    demonstrate_special_values()

    # Kahan summation
    print("\n--- Kahan Summation ---")
    np.random.seed(42)
    values = np.random.uniform(1e10, 1e10 + 1, 1_000_000)
    naive = sum(values)
    kahan = kahan_summation(values)
    print(f"Naive:  {naive:.6f}")
    print(f"Kahan:  {kahan:.6f}")
    print(f"Difference: {abs(naive - kahan):.6f}")

Applications in AI/ML

Training Considerations

Loss Scaling

When using float16, gradients often fall below the smallest representable normal (6.1 × 10^(-5)). Loss scaling multiplies the loss by a large factor before backpropagation, shifting gradients into the representable range, then divides after the optimizer step.

Gradient Clipping

To prevent float16 overflow during gradient accumulation, gradients are clipped to a maximum norm before the optimizer step. Without clipping, exploding gradients in float16 cause NaN propagation that corrupts the entire model.

Quantization

Model quantization converts float32 weights to int8 or float4 representations. Understanding the rounding behavior of float32 is essential for designing quantization schemes that minimize accuracy loss.

Common Pitfalls in ML

• Log-sum-exp overflow: Computing log(Σexp(x_i)) where x_i are large causes exp to overflow. Use the identity log(Σexp(x_i)) = max(x) + log(Σexp(x_i − max(x))).
• Softmax instability: If inputs differ by more than ~30 in float32, the softmax denominator overflows. Always subtract the max value before exponentiating.
• Batch norm with float16: Running statistics computed in float16 lose precision. Always compute batch norm in float32.

Common Mistakes

Interview Questions

Q1: Why is 0.1 + 0.2 ≠ 0.3 in floating point? A: 0.1 and 0.2 are not exactly representable in binary. Their float64 representations are slightly off from the true decimal values. When added, the result is close to 0.3 but not equal to the float64 representation of 0.3. The difference is approximately 5.55 × 10^(−17), which is within one ULP (unit in the last place) of the true sum.

Q2: What is machine epsilon and why does it matter? A: Machine epsilon is the smallest ε such that fl(1 + ε) ≠ 1. For float64 it is 2^(−52) ≈ 2.22 × 10^(−16). It bounds the relative rounding error of every floating point operation. Algorithms with more than ~1/ε operations may lose all significant digits.

Q3: How would you implement a numerically stable softmax? A: Subtract the maximum value from all inputs before exponentiation: softmax(x_i) = exp(x_i − max(x)) / Σ exp(x_j − max(x)). This prevents overflow since all exponent arguments are ≤ 0. The result is mathematically identical but numerically stable.

Q4: What is the difference between float16 and bfloat16? A: float16 has 5 exponent bits and 10 mantissa bits (range ±65504, ~3 digits precision). bfloat16 has 8 exponent bits and 7 mantissa bits (range ±3.4 × 10^38, ~2 digits precision). bfloat16 is preferred for ML because its float32-like dynamic range prevents overflow, even though it has slightly less precision.

Q5: Explain Kahan summation and why it works. A: Kahan summation tracks a compensation variable that captures the low-order bits lost when adding a value to the running sum. Each iteration: (1) subtract compensation from the input, (2) add to sum, (3) recover the lost part using algebraic identity, (4) store it as new compensation. This achieves O(1) error regardless of n, versus O(n) for naive summation.

Q6: When would you use float32 vs float64 in ML? A: float64 is rarely needed in ML — it doubles memory and compute with no accuracy benefit for most models. float32 is the standard for training. float16/bfloat16 are used for inference and mixed-precision training to save memory and leverage tensor cores. Use float64 only for numerical algorithms requiring high precision (e.g., solving linear systems, computing condition numbers).

Practice Problems

Quick Reference

Cross-References

• 087-Root-Finding — Newton's method depends on floating point arithmetic; understanding rounding is essential for choosing convergence tolerances
• 088-Integration — Quadrature error analysis requires understanding truncation vs roundoff error tradeoffs
• 089-Interpolation — Polynomial interpolation accuracy is limited by floating point precision; Runge's phenomenon interacts with roundoff
• 090-Error Analysis — Condition numbers and stability analysis build directly on IEEE 754 properties
• Linear Algebra — Matrix operations (LU decomposition, eigenvalue computation) are highly sensitive to floating point rounding; pivoting strategies address this
• Optimization — Gradient descent convergence depends on learning rates being representable; underflow in gradient updates causes silent failures