mixle.engines.extended module

Fast extended precision via error-free transformations – the “big boy math” end of mixle’s numeric spectrum, with no mpmath/gmpy2 in the hot path.

mpmath and gmpy2 are per-object and non-vectorized (~1000x slower than float64 at ~100-bit precision); they belong in tests as a correctness oracle, never in compute. The fast way to exceed float64 is error-free transformations (Dekker/Knuth TwoSum / TwoProd): represent a number as an unevaluated sum of float64 components (hi + lo) and carry the rounding error explicitly. Every operation is a handful of float64 ops that vectorize over numpy arrays, so a double-double (~106-bit mantissa, “fp128”) costs ~5-25x a float64 op – versus ~1000x for mpmath at the same precision.

This module provides the double-double primitives and the two reductions that matter for mixle’s EM hot paths – an accurate sum and dot product – where catastrophic cancellation (E[x^2]-E[x]^2, log-sum-exp of near-equal terms) otherwise eats precision. Beyond double-double, quad-double / multi-limb take over (a Cython/C job); this file is the part that needs only numpy.

The Veltkamp split is power-of-two scaled for large operands, so every finite representable product either produces a finite error-free pair or fails closed when the exact result cannot be represented.

two_sum(a, b)[source]

Error-free transformation of a sum: returns (s, e) with a + b == s + e exactly.

Knuth’s TwoSum – no assumption on the relative magnitudes of a and b. Vectorized.

Parameters:
Return type:

tuple[Any, Any]

quick_two_sum(a, b)[source]

Error-free sum assuming |a| >= |b| (Dekker). One fewer op than two_sum().

Parameters:
Return type:

tuple[Any, Any]

two_prod(a, b)[source]

Error-free transformation of a product: (p, e) with a * b == p + e exactly (Dekker).

Uses the Veltkamp split because numpy exposes no fused-multiply-add. Vectorized.

A product of two non-zero operands that underflows all the way to zero fails closed: this transform’s contract is that p + e reproduces a * b exactly, and (0, 0) does not. A reduction over many such products is a different question – see dd_dot().

Parameters:
Return type:

tuple[Any, Any]

class DoubleDouble(hi, lo=0.0)[source]

Bases: object

An (almost) ~106-bit float as two non-overlapping float64 arrays hi + lo.

Scalars or numpy arrays; operations broadcast. The invariant is |lo| <= 0.5 * ulp(hi). A component may be +inf or -inf – exactly representable in float64, carried with lo == 0 – but never NaN, and never an overflow of finite operands.

Parameters:

  • hi (Any)

  • lo (Any)

classmethod from_float(x)[source]

Create a double-double value from float64 data.

Parameters:

x (Any)

Return type:

DoubleDouble

to_float()[source]

Collapse back to the nearest float64 (hi rounded with lo).

Return type:

ndarray

dd_sum(x)[source]

Accurate sum of a float64 array in double-double precision – vectorized, no Python loop over elements.

Pairwise tree reduction with an error-free two_sum() combine at every node: O(n) work in O(log n) vectorized passes, accumulating the rounding error into the lo component. The result is correct to ~106 bits even for catastrophically cancelling inputs that float64 sums get wrong.

An infinite term short-circuits to the IEEE answer (-inf for a zero-probability log-density) instead of raising; NaN, and +inf mixed with -inf, still fail closed – see _infinite_sum().

Parameters:

x (Any)

Return type:

DoubleDouble

dd_dot(a, b)[source]

Accurate dot product sum(a_i * b_i) in double-double precision.

Uses the compiled hardware-FMA kernel when available (one fma per element, ~3x faster than the pure-numpy Veltkamp-split path, and exactly as accurate – each product’s error is exact either way, though the two paths accumulate in different orders, sequential versus pairwise tree, so their results agree to double-double precision rather than bit for bit). Otherwise each product is split error-free by _product_and_residual() and the products + errors are summed by dd_sum(). Defeats the cancellation that wrecks a naive float64 dot.

Every check that decides whether an input is acceptable is made once, before dispatch, so the two implementations cannot drift apart on what they accept: equal length (see _require_equal_length(), ValueError), finite operands (ValueError), and a result that is representable at all (OverflowError). A single product that underflows to zero is not one of them: it contributes less than one smallest-subnormal to the sum, so it is summed as the zero it rounds to rather than failing the whole reduction – dd_dot([1e-200, 1.0], [1e-200, 1.0]) is 1.0, which is the right answer, not an error. (two_prod() on its own still fails closed there, because a lone product has nothing to be negligible against.)

Parameters:
Return type:

DoubleDouble