A standard float64 (an IEEE-754 double) carries about 15.95 decimal digits. That is plenty for most work, but ill-conditioned problems — subtracting nearly-equal numbers, accumulating long sums, evaluating a function near a root — can lose all of those digits to rounding. When you need more headroom, Equana offers two extended-precision types:

• float128 — a real number with roughly 31 decimal digits (~106 significant bits).
• complex256 — a complex number whose real and imaginary parts are each a float128.

This tutorial shows how to use them, verifies the extra precision with cells that fail loudly if the accuracy isn't real, and then explains how it works.

Honest marketing. Every claim below is checked by an assert in the same cell. If float128 were secretly computing in float64, these cells would raise an AssertionError instead of running clean — so what you see is what you get.

Creating extended-precision values

Use the float128 and complex256 constructors, exactly like float64 or complex128. The capitalized names Float128 and Complex256 are the type values (handy for eps, typeof, and dispatch).

Machine epsilon: how much precision?

eps(T) returns the machine epsilon of a type — the gap between 1 and the next representable number. finfo(T) prints the full machine profile (à la NumPy). The smaller the epsilon, the more digits you keep.

Note precision = 31 and eps ≈ 4.93e-32 (which is 2^-104) — roughly double the ~53-bit mantissa of a float64.

Verification 1 — retaining a bit that float64 drops

2^-80 is far smaller than the float64 gap at 1.0 (which is eps ≈ 2^-52). So in float64, 1 + 2^-80 rounds straight back to 1, and subtracting 1 gives exactly 0 — the added bit vanished.

In float128 there is room to keep it:

The same computation in float64 collapses to zero:

Verification 2 — catastrophic cancellation

1e16 needs 54 bits to represent 1e16 + 1, one more than a float64 mantissa holds. So 1e16 + 1 rounds to 1e16 and the 1 is lost. Adding a large number and then subtracting it back is a classic way to destroy small quantities:

float128 recovers the exact answer; float64 loses it entirely. This is exactly the failure mode that wrecks ill-conditioned sums, variance formulas, and residuals near a root.

Verification 3 — the digits are really there

Divide 1 by 3 in each precision and look at how many correct threes come back:

The float64 result stops after ~16 digits; the float128 result carries ~31. The same is true for irrational constants — float128(pi) uses a ~32-digit expansion of π rather than round-tripping through a 64-bit pi:

complex256 — quad precision in the complex plane

complex256 applies the same treatment to both components. Here the real part is annihilated by cancellation in complex128 but preserved in complex256:

Multiplication, division, conj, abs, and the component accessors real / imag all run in double-double:

How it works

Why a JS number can't hold it

Equana runs on JavaScript, whose only number type is the IEEE-754 double — a single 64-bit float with a 53-bit significand. There is no 128-bit primitive to lean on. A naive float128 that stored its value in one JS number would silently round every literal and every operation back to 53 bits: quad precision in name only.

Double-double representation

Instead, a float128 value is stored as a double-double: an unevaluated sum of two doubles,

hi is the usual 53-bit approximation; lo captures the part that fell off the end. Together they give ~106 significant bits (~31 digits). A complex256 is simply two of these — a double-double real part and a double-double imaginary part. This is the same scheme used by the classic QD library and by NumPy's longdouble on platforms without hardware binary128.

The storage layout already matches: a float128 element is 16 bytes (two float64s) and a complex256 element is 32 bytes (four float64s), so hi and lo sit side by side.

Avoiding the float64 round-trip

The precision lives in the arithmetic, not the literal. Each operation is built from error-free transformations — small sequences of + and * that compute both the rounded result and its rounding error exactly:

• twoSum(a, b) returns (s, e) where s = fl(a+b) and a + b = s + e exactly.
• twoProd(a, b) does the same for multiplication (via Dekker's splitting).

float128 addition, subtraction, multiplication, division, sqrt, integer powers, and comparisons are assembled from these primitives so no intermediate is ever flattened to a single double. That is why float128(1) + 2^-80 keeps the extra bit: twoSum routes it into lo instead of discarding it.

Decimal literals are still parsed as float64 first, so float128(0.1) carries float64's rounding of 0.1. To obtain a value to full quad precision, compute it in quad — e.g. float128(1) / float128(10) — or use the built-in constants (float128(pi), float128(e)), which are seeded from ~32-digit decimal strings.

Cost versus float64

Double-double is software precision. A float128 add is roughly a dozen float64 flops and a multiply a couple dozen, so expect it to be ~10–30× slower than native float64, with twice the memory per value. Reach for it where accuracy dominates — conditioning-sensitive residuals, reference results, summation of long series — and keep hot inner loops in float64 where 16 digits already suffice.

Scope. Extended precision is implemented for scalar float128 / complex256 arithmetic (the operations above). Element-wise float128 NDArray kernels currently evaluate in float64; build quad-precision reductions from scalar operations when you need them.

Recap

• float128 (~31 digits) and complex256 are backed by genuine double-double arithmetic, not a relabeled float64.
• eps(Float128) ≈ 4.93e-32, ~16 orders of magnitude below float64.
• Extended precision rescues computations that float64 destroys through cancellation — every claim here is enforced by an in-cell assert.
• It costs ~10–30× the time and 2× the memory of float64, so use it deliberately.