Ereal Maximum Limb Count: Mathematical Derivation

Document: Technical analysis of the 19-limb constraint for ereal multi-component arithmetic Date: 2025-01-04 Author: Analysis based on Shewchuk’s expansion arithmetic theory

---

Executive Summary

ereal<maxlimbs, FpType> is limited in how many limbs it may have, by a subtle but critical constraint in Shewchuk’s expansion arithmetic: error terms from two-sum operations on the smallest limb must remain representable as normal values of the limb type.

Key Finding: The constraint applies to all possible starting magnitudes (including values near 1.0), not just the largest representable one. For double limbs that reduces the theoretical limit from ~38 limbs to 19.

The general bound

The derivation below is written for double, because double is the default limb and the concrete numbers make it readable. The result is not specific to double. Each limb is digits bits below the one above it, so starting from a magnitude near 1.0 the n-th limb sits at 2^(-digits * n), and it must stay at or above the smallest normal 2^(min_exponent - 1):

static constexpr unsigned max_safe_limbs =
    static_cast<unsigned>(-(std::numeric_limits<FpType>::min_exponent - 1)
                          / std::numeric_limits<FpType>::digits);

limb type	digits	min_exponent	max_safe_limbs	~decimal digits	where
float	24	-125	5	~36	everywhere
double	53	-1021	19	~303	everywhere
long double = x87 extended	64	-16381	255	~4913	x86-64, MinGW
long double = binary128	113	-16381	144	~4898	aarch64, riscv64, ppc64le with -mabi=ieeelongdouble
long double = IBM double-double	106	-968	–	–	rejected, see below
long double = double	53	-1021	19	~303	MSVC, Apple ARM64

(Measured on the real toolchains under QEMU and wine, not inferred from the standard.)

Which types may be limbs at all

is_expansion_limb_v<FpType> requires a binary, IEC-559 floating-point type whose exponent range is symmetric in the IEEE sense, min_exponent == 3 - max_exponent. That last clause is what does the work. libstdc++ reports is_iec559 == true for IBM extended double-double – the default long double on ppc64le – so an is_iec559 check alone does not reject it. Its exponent range gives it away: min_exponent is -968 against a max_exponent of 1024, where IEEE-754 requires min_exponent == 3 - max_exponent, i.e. -1021.

IBM double-double must be rejected, and not on a technicality: it is itself a two-component unevaluated sum, so using it as a limb would nest an expansion inside an expansion. Its components are not at a fixed relative scale, so “the limbs are non-overlapping and descending” – the invariant every algorithm here depends on – is not something the outer expansion can establish.

On ppc64le this means ereal<n, long double> does not compile, with a static_assert that says why. That is deliberate. Building the whole toolchain with -mabi=ieeelongdouble makes long double binary128 there, and it then works and is error-free (measured); but that is a whole-program ABI choice affecting every library on the system, not something a header can opt into, so ereal does not try to detect or accommodate it. Use ereal<n, double> on ppc64le, or build the toolchain for binary128.

A wider limb is for reach, not speed

binary128 on aarch64 and riscv64 has no hardware: every operation is a libgcc soft-float call, and it is slower than double by a large factor – typically 10x or more – for about twice the significand. x87 is hardware, but its 80-bit loads and stores are not free either. The wide limb is what takes ereal past double’s ~303-digit ceiling; it is not a faster way to get 100 digits. benchmark/accuracy/adaptive/characterize.cpp sweeps accuracy against time across limb types so the trade stays visible.

---

Table of Contents

1. Background: Shewchuk’s Expansion Arithmetic
2. Maximum number of limbs
3. Why 19 libms Is Actually Correct
4. Mathematical Derivation
5. Practical Examples
6. Implementation Constraints
7. References

---

Background: Shewchuk’s Expansion Arithmetic

What is an Expansion?

An expansion is a sequence of non-overlapping floating-point numbers that together represent a high-precision value:

value = limb[0] + limb[1] + limb[2] + ... + limb[n]

Key properties:

1. Non-overlapping: Adjacent limbs differ by at least 53 bits (one mantissa width)
2. Ordered: |limb[0]| ≥ |limb[1]| ≥ |limb[2]| ≥ … ≥ |limb[n]|
3. Error-free: Each limb captures rounding errors from the previous limb

Two-Sum Algorithm

The foundation of expansion arithmetic is the two-sum algorithm, which computes the exact sum of two doubles:

// Fast-Two-Sum (Dekker, 1971):
s = a + b;           // Rounded sum (what hardware gives you)
e = (a - s) + b;     // Exact error (what was lost to rounding)

Critical requirement: Both s and e must be representable as normal IEEE-754 doubles.

If e underflows below DBL_MIN (≈ 2^-1022), the error is lost and the expansion arithmetic breaks down.

---

Maximum number of limbs

Range-based Calculation

IEEE-754 double-precision has:

• Largest normal: DBL_MAX ≈ 2^1023 ≈ 10^308
• Smallest normal: DBL_MIN ≈ 2^-1022 ≈ 10^-308
• Total range: 2^(1023-(-1022)) = 2^2045 ≈ 10^616

Each limb adds approximately 53 bits of precision (the mantissa width), so:

Total range in bits: 2045 bits
Bits per limb: 53 bits
Maximum limbs: 2045 / 53 ≈ 38.6 limbs

This suggests ~38 limbs should work!

The Flawed Assumption

This calculation assumes you can directly use the full range of doubles, spacing limbs from 2^1023 down to 2^-1022:

limb[0]:  2^1023      (top of range)
limb[1]:  2^(1023-54) (54 bits smaller, arbitrary spacing)
limb[2]:  2^(1023-108)
...
limb[37]: 2^-1022     (bottom of range)

This is not how Shewchuk’s expansion arithmetic works!

Why 19 libms Is Actually Correct

How Expansions Are Actually Built

Expansions are constructed through repeated error extraction, not by directly choosing exponents:

Process:

1. Start with a value x at some magnitude (could be near 1.0, 10^8, or any value)
2. Round it to create limb[0]
3. Compute the rounding error using two-sum → this becomes limb[1]
4. Repeat: the error from limb[1] becomes limb[2], and so on

Each successive limb is approximately 53 bits smaller because:

• Rounding error is at most 1 ULP (unit in last place)
• For a value at exponent E, 1 ULP = 2^(E-52)
• Therefore, error magnitude ≈ 2^(E-53)

The Critical Insight

The expansion must work for ANY starting magnitude, not just DBL_MAX!

Consider a value near 1.0:

Starting value: 1.0 = 2^0

limb[0]:  2^0       (magnitude ~1)
limb[1]:  2^-53     (error from limb[0])
limb[2]:  2^-106    (error from limb[1])
limb[3]:  2^-159
...
limb[18]: 2^-954
limb[19]: 2^-1007   ✓ Still above DBL_MIN (2^-1022)
limb[20]: 2^-1060   ✗ BELOW DBL_MIN! UNDERFLOW!

At limb[20], the value underflows below DBL_MIN and:

• Cannot be represented as a normalized double
• May become denormalized (loses precision)
• Two-sum operations produce incorrect error terms
• The expansion has reached its limit

Mathematical Derivation {#mathematical-derivation}

General Limb Magnitude Formula

For an expansion starting at magnitude M (where M = 2^E for some exponent E):

limb[0]: M × 2^0        = M
limb[1]: M × 2^-53      (error from limb[0])
limb[2]: M × 2^-106     (error from limb[1])
...
limb[n]: M × 2^(-53n)   (error from limb[n-1])

Constraint for Normal Representation

For limb[n] to be representable as a normal double:

M × 2^(-53n) >= DBL_MIN
M × 2^(-53n) >= 2^-1022

Case 1: Starting from DBL_MAX (M = 2^1023)

2^1023 × 2^(-53n) >= 2^-1022
2^(1023-53n) >= 2^-1022
1023 - 53n >= -1022
53n <= 2045
n <= 38.6

Result: Up to 38 limbs would work if expansions always started from DBL_MAX.

Case 2: Starting from 1.0 (M = 2^0)

2^0 × 2^(-53n) >= 2^-1022
2^(-53n) >= 2^-1022
-53n >= -1022
53n <= 1022
n <= 19.28

Result: Only 19 limbs work for values near 1.0.

Case 3: General Case (M = 2^E for arbitrary E)

For the expansion to work for any starting magnitude:

2^E × 2^(-53n) >= 2^-1022

For all possible E (where -1022 ≤ E ≤ 1023):
Worst case is E = 0 (magnitude ~1.0)

2^(-53n) >= 2^-1022
n <= 19.28

Conclusion: To guarantee correctness for all possible starting values, we must limit to n ≤ 19 limbs.

Why Not Adapt Based on Magnitude?

One might ask: “Why not allow 38 limbs for large values and 19 for small values?”

Answer: Type safety and implementation complexity.

• ereal<maxlimbs> is a compile-time fixed type
• Operations between ereal<19> and ereal<38> would be undefined
• Arithmetic must preserve the expansion property for all intermediate results
• A multiplication of two large values could produce a small value, requiring the smaller limit
• The type must work correctly for its entire value range

---

Practical Examples {#practical-examples}

Example 1: Value Near 1.0 with 19 Limbs

ereal<19> x(1.0);  // Start with 1.0

// Limb magnitudes after successive operations:
limb[0]:  1.0 × 2^0    = 2^0     = 1.0
limb[1]:  1.0 × 2^-53  = 2^-53   ≈ 1.11e-16
limb[2]:  1.0 × 2^-106 = 2^-106  ≈ 1.23e-32
...
limb[18]: 1.0 × 2^-954 = 2^-954  ≈ 1.89e-287
limb[19]: 1.0 × 2^-1007 = 2^-1007 ≈ 4.75e-304

DBL_MIN = 2^-1022 ≈ 2.23e-308

limb[19] = 2^-1007 > 2^-1022 ✓ SAFE

All limbs remain representable as normal doubles. ✓

Example 2: Value Near 1.0 with 20 Limbs

ereal<20> x(1.0);  // Hypothetical - would fail at compile time

limb[19]: 1.0 × 2^-1007 = 2^-1007 ≈ 4.75e-304  ✓ OK
limb[20]: 1.0 × 2^-1060 = 2^-1060 ≈ 9.21e-320  ✗ BELOW DBL_MIN!

DBL_MIN = 2^-1022 ≈ 2.23e-308

limb[20] = 2^-1060 < 2^-1022 ✗ UNDERFLOW!

limb[20] underflows to zero or denormal, breaking two-sum. ✗

Example 3: Value Near DBL_MAX with 19 Limbs

ereal<19> x(1e308);  // Near DBL_MAX ≈ 2^1023

limb[0]:  2^1023
limb[1]:  2^970
limb[2]:  2^917
...
limb[18]: 2^69    ≈ 5.9e20
limb[19]: 2^16    ≈ 65536

DBL_MIN = 2^-1022 ≈ 2.23e-308

limb[19] = 2^16 >> 2^-1022 ✓ SAFE (huge margin!)

Even for maximum values, 19 limbs works perfectly. ✓

Example 4: Why Range-Based Thinking Fails

Incorrect reasoning:

“If I start with 2^1023, I can go down to 2^-1022, that’s 38 limbs worth of range!”

Why it fails:

ereal<38> x(1e308);  // Starts near 2^1023
ereal<38> y(1.0);    // Starts near 2^0

// These are the SAME TYPE!
// But 'y' would have limbs below DBL_MIN
// Type must work for BOTH values

ereal<38> z = x * 0.0001;  // Now z is small
// z's small limbs would underflow

The type must handle its entire value range, not just large values.

---

The Limb Budget: How Many Limbs a Result Gets

The derivation above bounds how many limbs a value may store. It says nothing about how many limbs an arithmetic result comes back with, and until #1572 nothing did.

The expansion algorithms return whatever limb count the operands produce. The only thing that ever pruned a result was underflow: a component below the limb type’s smallest normal becomes zero, and renormalization drops it. With double limbs that accident looks like a cap – ereal<8>’s 1/3 settles at 16 limbs and stays there. A wide-exponent limb has no such floor. On x87, whose exponent reaches 2^-16382, every quotient grew until it spanned the whole range:

type	1/3	limbs	Newton sqrt(2), 8 steps
ereal<8>	260 digits	16 (maxlimbs is 8)	0.03 s
ereal<8, long double> (x87)		250	41 s
ereal<24, long double> (x87)		250	71 s

250 is not a coincidence: 250 * 64 bits is x87’s exponent range. Every long double quotient reached it whatever maxlimbs said, so a single division cost seconds and the wider limb was unusable in practice – the opposite of what parameterizing the limb was for.

The rule

static constexpr unsigned limb_budget =
    (2 * maxlimbs < max_safe_limbs) ? 2 * maxlimbs : max_safe_limbs;

Truncating is always safe: an expansion in Priest normal form has descending, non-overlapping components, so a prefix is the leading-order value and what is dropped lies below the precision retained. parse() has truncated to maxlimbs on exactly this reasoning since #1006.

Why 2 * maxlimbs and not maxlimbs. The tail is not waste, it is the guard digits an iterative algorithm needs, and the existing suites were tuned against results that carry them – ereal<8>’s 1/3 at 16 limbs is twice the 8 that were asked for. Capping at maxlimbs would have cut that from 260 digits to ~130 and moved accuracy numbers throughout the library. At 2 * maxlimbs no measured double digit count changes at all.

Why the clamp. A limb past max_safe_limbs is subnormal, which is not a valid expansion component – that is what the rest of this document derives.

Where it is applied

Capping only the result is not enough: the Newton reciprocal squares its iterate’s limb count every step, so a division still built the 250-limb intermediate and still cost seconds before the result was trimmed. The budget is passed into expansion_reciprocal, which truncates the iterate each step. Newton is self-correcting – each step recomputes the residual 2 - e*r from scratch rather than accumulating it – so an iterate truncated to b limbs still converges to b limbs of accuracy.

The intermediates e*r_n and 2 - e*r_n are deliberately not truncated. Their low components cancel exactly against 2, which is what lets the residual collapse to a couple of limbs; truncating them leaves a spurious residual at the truncation level, and the next step then multiplies two full-length expansions instead of a long one by a short one. That made double division about 2.7x slower.

The bound also tightens as the iteration proceeds. Newton doubles its correct digits each step, so from a one-limb seed the iterate after step i is accurate to 2^(i+1) limbs; everything past that is noise, and it is that noise which multiplies e on the next step. Without this, a full-width divisor made every step pay a full budget-by-budget product: one ereal<24, long double> division cost 0.46 s, against 18 ms for a full-width multiply.

Truncating each step to the precision it has reached makes the widths grow geometrically and saturate at the budget. step_budget reaches budget once 2^(i+1) + 2 >= budget, so the last step or two run at full width – how many depends on how far 2^iterations overshoots the budget, which is why it is one step for ereal<24> (widths 4, 6, 10, 18, 34, 48 over six iterations) and can be two when maxlimbs is not near a power of two. The total is a geometric sum: about twice the final width, rather than iterations * budget. For ereal<24, long double> that is 120 limb-widths of work instead of 288.

The + 2 guard limbs in that bound are not defensive padding. The doubling is the asymptotic rate and rounding within a step eats into it, so trimming to exactly 2^(i+1) cost the last two digits – ereal<8>’s Newton sqrt(2) went from 257 digits to 255.

Newton sqrt(2), 8 steps	uncapped	result capped only	with the per-step bound
ereal<8> (double)	0.0072 s	0.0067 s	0.0036 s
ereal<19> (double)	0.0100 s	0.0101 s	0.0063 s
ereal<8, long double> (x87)	40.9 s	0.065 s	0.012 s
ereal<24, long double> (x87)	70.9 s	3.65 s	0.35 s

double ends up about twice as fast as it was before any of this, not merely unregressed.

What it costs

Results that used to exceed max_safe_limbs lose the digits those extra limbs were carrying. Measured:

	before	after
ereal<24, long double> sqrt(2)	1235 digits (250 limbs)	953 digits (48 limbs)
ereal<19> sqrt(2)	323 digits (20 limbs)	314 digits (19 limbs)
ereal<5, float> sqrt(2)	45 digits (6 limbs)	38 digits (5 limbs)
ereal<8> 1/3	260 digits (16 limbs)	260 digits (16 limbs)

The 20th double limb and the 6th float limb are past max_safe_limbs: those digits were never ones the type was entitled to.

A consequence worth stating plainly: multiplication is error-free only while its result fits the budget. An exact product of two n-limb expansions can need more components than max_safe_limbs allows in the first place, so exactness is a property of results that fit, not of the operation. Past the budget the guarantee is the truncation contract – what was dropped lies below the last limb that was kept – which elastic/ereal/arithmetic/exact_value_oracle.cpp checks against an exact dyadic oracle.

---

Implementation Constraints {#implementation-constraints}

Static Assertion in Code

From ereal_impl.hpp. There are two, and the first is the one people hit:

static_assert(is_expansion_limb_v<FpType>,
    "ereal<maxlimbs, FpType>: FpType must be a p-bit IEEE-754 binary type (float, double, "
    "x87 extended or binary128 long double). IBM extended double-double -- the default "
    "long double on ppc64le -- is itself a two-component expansion and cannot serve as a "
    "limb; build with -mabi=ieeelongdouble for a binary128 long double there.");

static_assert(maxlimbs <= max_safe_limbs,
    "ereal<maxlimbs, FpType>: maxlimbs must be <= max_safe_limbs = -(min_exponent - 1) / digits "
    "of the limb type (5 for float, 19 for double, 255 for x87, 144 for binary128). More limbs "
    "push the last one below the smallest normal, violating the non-overlapping property "
    "Shewchuk's expansion arithmetic requires, and two_sum/two_product silently lose bits.");

Note that the bound is max_safe_limbs, derived from the limb type, rather than the literal 19 it was when double was the only limb.

Precision vs Limb Count

Limbs	Mantissa Bits	Decimal Digits	Status
4	212	~64	✓ Safe
8	424	~127	✓ Safe
12	636	~191	✓ Safe
16	848	~255	✓ Safe
19	1007	~303	✓ Safe (maximum)
20	1060	~319	✗ FAILS for values near 1.0
38	2014	~606	✗ FAILS for values < 2^969

Error Manifestations with maxlimbs > 19

If the limit were violated (without the static assertion):

1. Silent Arithmetic Errors:

   • Two-sum returns incorrect error terms (lost to underflow)
   • Non-overlapping property violated
   • Results appear valid but are mathematically wrong

2. Precision Loss:

   • Small limbs become denormalized (54-bit precision → fewer bits)
   • Gradual degradation rather than catastrophic failure
   • Hard to detect without rigorous testing

3. Inconsistent Behavior:

   • Works correctly for large values (near DBL_MAX)
   • Fails silently for small values (near 1.0)
   • Breaks the principle of uniform type behavior

---

Why This Matters {#why-this-matters}

Type Safety

// This should work uniformly for all values:
ereal<n> x = ...;  // any value
ereal<n> y = ...;  // any value
ereal<n> z = x + y;  // must be correct regardless of magnitudes

The type cannot have “safe regions” and “unsafe regions” - it must work correctly for its entire value range.

Mathematical Correctness

Shewchuk’s algorithms provide exact arithmetic (within the limits of floating-point representation). This property must be preserved:

// These are GUARANTEED to be exact (error-free):
a + b = (sum, error)  // two_sum
a * b = (prod, error) // two_product

If limbs underflow, these guarantees are violated, leading to:

• Incorrect geometric predicates (orient2d, orient3d)
• Wrong results in high-precision calculations
• Unreliable numerical algorithms

Performance Considerations

While 19 limbs might seem limiting, it provides:

• ~303 decimal digits of precision (more than quadruple-double’s ~64 digits)
• Algorithmic correctness for all operations
• Predictable behavior across the entire value range

For applications needing more precision, consider:

• Arbitrary-precision libraries (MPFR, GMP)
• Symbolic computation (exact rational arithmetic)
• Extended exponent range formats (quadruple precision, custom formats)

---

References {#references}

Primary Sources

1. Shewchuk, J. R. (1997). Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates. Discrete & Computational Geometry, 18(3), 305-363.

   • Defines expansion arithmetic and two-sum/two-product algorithms
   • Establishes correctness requirements for multi-component arithmetic
   • Available: https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf

2. Dekker, T. J. (1971). A Floating-Point Technique for Extending the Available Precision. Numerische Mathematik, 18(3), 224-242.

   • Original two-sum algorithm (Fast-Two-Sum)
   • Foundation for expansion arithmetic

3. Priest, D. M. (1991). Algorithms for Arbitrary Precision Floating Point Arithmetic. Proceedings of the 10th Symposium on Computer Arithmetic.

   • Early work on multi-component arithmetic
   • Discusses precision limits and error propagation

IEEE-754 Standard

4. IEEE Standard 754-2008 for Floating-Point Arithmetic
   • Defines normal vs subnormal numbers
   • DBL_MIN = 2^-1022 (smallest normal double)
   • Rounding modes and error behavior

Related Documentation

5. Universal Numbers Library - ereal implementation
   • include/sw/universal/number/ereal/ereal_impl.hpp (lines 34-61)
   • Static assertion enforcing maxlimbs ≤ 19
   • Implementation of Shewchuk’s algorithms

---

Appendix: Quick Reference {#quick-reference}

Key Formulas

Maximum limbs for starting magnitude M = 2^E:

M × 2^(-53n) >= 2^-1022
n <= (1022 + E) / 53

For E = 0 (magnitude ~1.0):

n <= 1022 / 53 ≈ 19.28  →  n_max = 19

For E = 1023 (magnitude ~DBL_MAX):

n <= 2045 / 53 ≈ 38.6  →  n_max = 38 (but type must work for E=0!)

Decision Tree

Q: How many limbs can I use?
├─ Will values ever be < 2^969?
│  ├─ Yes → Use maxlimbs ≤ 19
│  └─ No  → Could use up to ~38 (but type must be uniform!)
└─ Need uniform type behavior?
   └─ Yes → Use maxlimbs ≤ 19 (always)

Common Misconceptions

❌ Wrong: “The range is 2^2045, so I can use 38 limbs” ✓ Right: “Each limb is 53 bits smaller, and the smallest must be ≥ DBL_MIN”

❌ Wrong: “I only work with large values, so I can use more limbs” ✓ Right: “The type must work for all values, including results of operations”

❌ Wrong: “I can check the magnitude and use variable limb counts” ✓ Right: “maxlimbs is a compile-time constant; the type is fixed”

---

Conclusion

The 19-limb limit for ereal arises from a fundamental requirement of Shewchuk’s expansion arithmetic: all components and error terms must remain representable as normal IEEE-754 doubles.

While naive analysis of the double-precision range suggests ~38 limbs might work, the actual constraint is more subtle: the expansion must work correctly for any starting magnitude, not just the maximum representable value.

For values near 1.0 (which are common in many applications), the 20th limb would have magnitude ~2^-1060, which is below DBL_MIN and therefore cannot be represented as a normal double. This violates the fundamental assumptions of the two-sum algorithm, leading to incorrect arithmetic.

The limit of 19 limbs provides:

• ✓ ~303 decimal digits of precision (mantissa basis)
• ✓ Algorithmic correctness for all value magnitudes
• ✓ Predictable behavior across the type’s entire range
• ✓ Error-free arithmetic properties of expansion arithmetic

This is not a limitation of the implementation, but rather a fundamental mathematical constraint of multi-component floating-point arithmetic using IEEE-754 doubles.

---

Document Status: Ready for review Last Updated: 2025-01-04 Location: docs/ereal_limb_limit_derivation.md