{ } JustinGeoSci

Expansion arithmetic

Shewchuk’s error-free primitives, the building blocks of the adaptive stages, in the public arithmetica::expansion module.

Exactness is conditional here. These are Shewchuk’s error-free primitives, and they are exact only when no operation overflows and, for products, when each product is no smaller than about 2−970 (roughly 10−292). Below that, a result can be finite and silently inexact. The predicates establish these conditions themselves through the exactness window; a caller using the primitives directly takes on that duty. If you cannot establish it for your inputs, use the predicates instead.

expansion — Error-free transformations

Operations that return a rounded result together with its exact rounding error, as a two-term expansion (x, y) with x the larger.

two_sum

pub fn two_sum(a: f64, b: f64) -> (f64, f64)

The rounded sum and its exact rounding error, for any two values. Exact whenever a + b does not overflow; underflow is harmless, because a subnormal sum is always exact. Shewchuk’s Algorithm 2, six floating-point operations.

TwoSumx = \mathrm{fl}(a+b), \qquad x + y = a + b \text{ exactly}
use arithmetica::expansion::two_sum;

// 1e-16 is less than half an ulp of 1.0, so the rounded sum loses it.
let (x, y) = two_sum(1.0, 1e-16);
assert_eq!(x, 1.0);   // the rounded sum
assert_eq!(y, 1e-16); // exactly what rounding discarded

fast_two_sum

pub fn fast_two_sum(a: f64, b: f64) -> (f64, f64)

The same result as two_sum in four operations, under the precondition |a| ≥ |b|. The precondition is checked with debug_assert! in debug builds; in release builds a violation produces a wrong error term without any signal. Shewchuk’s Algorithm 1.

use arithmetica::expansion::fast_two_sum;

// Requires |a| >= |b|.
let (x, y) = fast_two_sum(1.0, 1e-16);
assert_eq!((x, y), (1.0, 1e-16));

two_product

pub fn two_product(a: f64, b: f64) -> (f64, f64)

The rounded product and its exact rounding error, by Veltkamp–Dekker splitting, so it does not depend on a hardware fused multiply-add. Exact when |a| and |b| are below about 2995, the product is finite, and the product is at least about 2−970. Shewchuk’s Algorithm 4.

TwoProductx = \mathrm{fl}(a\,b), \qquad x + y = a\,b \text{ exactly}
use arithmetica::expansion::two_product;

// (2^27 + 1)^2 = 2^54 + 2^28 + 1 needs 55 bits; the final + 1 does not fit.
let a = 134_217_729.0; // 2^27 + 1
let (x, y) = two_product(a, a);
assert_eq!(x, 18_014_398_777_917_440.0); // 2^54 + 2^28, the rounded product
assert_eq!(y, 1.0);                      // the bit rounding dropped

expansion — Expansion operations

An expansion is a slice of f64 components, ordered from smallest to largest magnitude, whose significands do not overlap; its value is the exact sum of the components. The caller provides the output buffer, so nothing is allocated.

grow_expansion

pub fn grow_expansion(e: &[f64], b: f64, h: &mut [f64]) -> usize

Adds a single value to an expansion. Writes the result into h, which must hold at least e.len() + 1 values, and returns how many components it wrote. Zero components are dropped.

Exact sum\sum_i h_i = \sum_i e_i + b
use arithmetica::expansion::grow_expansion;

let tiny = 2f64.powi(-60);
let e = [1.0];              // a one-term expansion
let mut h = [0.0; 2];       // room for e.len() + 1 components
let n = grow_expansion(&e, tiny, &mut h);

// 1 + 2^-60 needs more than 53 bits, so it is kept as two components,
// smallest first, that sum to it exactly.
assert_eq!(&h[..n], &[tiny, 1.0]);

scale_expansion

pub fn scale_expansion(e: &[f64], b: f64, h: &mut [f64]) -> usize

Multiplies an expansion by a single value. h must hold at least 2 * e.len() values. This is the operation most likely to lose exactness silently: the smallest component of an expansion can be many orders of magnitude below its largest, so the trailing products can fall below 2−970 while the leading one looks comfortably scaled.

Exact product\sum_i h_i = b\sum_i e_i
use arithmetica::expansion::scale_expansion;

let tiny = 2f64.powi(-60);
let e = [tiny, 1.0];        // 1 + 2^-60
let mut h = [0.0; 4];       // room for 2 * e.len() components
let n = scale_expansion(&e, 3.0, &mut h);

// 3 + 3 * 2^-60, exactly.
assert_eq!(&h[..n], &[3.0 * tiny, 3.0]);

expansion_sum

pub fn expansion_sum(e: &[f64], f: &[f64], h: &mut [f64]) -> usize

Adds two expansions. h must hold at least e.len() + f.len() values. Both inputs must already be valid expansions; the result is one too.

Exact sum\sum_i h_i = \sum_i e_i + \sum_j f_j
use arithmetica::expansion::expansion_sum;

let tiny = 2f64.powi(-60);
let e = [tiny, 1.0];        // 1 + 2^-60
let f = [-1.0];             // -1
let mut h = [0.0; 3];       // room for e.len() + f.len() components
let n = expansion_sum(&e, &f, &mut h);

// The 1 and -1 cancel exactly; only 2^-60 remains.
assert_eq!(&h[..n], &[tiny]);