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.

Caveat: the Veltkamp split overflows for |x| > ~1e300; mixle’s log-densities are far from that.

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.

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).

Parameters:

  • hi (Any)

  • lo (Any)

hi
lo
classmethod from_float(x)[source]
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.

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 bit-for-bit identical); otherwise each product is split error-free by two_prod() and the products + errors are summed by dd_sum(). Defeats the cancellation that wrecks a naive float64 dot.

Parameters:
Return type:

DoubleDouble