Every Digit It Is Given

Arbitrary-precision arithmetic isn't one technique. It's a stack of independent decisions, and most confusion about it comes from conflating the layers — arguing about algorithms when the disagreement is really about what kind of number is being represented.

1 · What kind of number are you representing?

Five families, and they answer different questions.

Big integers — GMP's mpz, Java's BigInteger, Python's int — are the foundation everything else is built on. Floats (MPFR) are what scientific work wants. Scaled decimal (Java BigDecimal, Python decimal, IBM's decNumber) is what you want when the decimal representation is the ground truth: money, invoices, anything a person will check by hand.

Rationals are seductive and usually a trap outside computer-algebra work — a few dozen operations and the denominator is forty thousand digits. Exact reals (iRRAM, Boehm's calculator in Android) genuinely give you π to whatever precision you ask for later, at the price that x == y can loop forever.

A sixth model sits orthogonal to all of them. Ball arithmetic — Arb is the state of the art — carries a midpoint plus a radius, so every result reports how many of its digits are actually true. It is the only model that answers "is this answer good?" rather than assuming it.

2 · What radix?

Binary limbs (2³², 2⁶⁴) win on speed: the hardware adder, the hardware multiplier and shifts-for-scaling are all free. Decimal limbs — BCD, or 10⁹ packed per uint32 — win on fidelity: 0.1 is exact, decimal rounding rules apply directly, and there is no conversion error at the I/O boundary.

The cost of decimal is real but bounded. Carry propagation needs a compare-and-subtract instead of falling out of the word width, and roughly 7% of each word goes unused. The cost of binary in a decimal-facing application is that base conversion moves onto the hot path — and naive base conversion is O(n²), which makes it the single most under-measured cost in bignum libraries.

IEEE 754-2008 standardised decimal floating point for exactly this reason, with two encodings: BID, used by Intel's library, and DPD, used by IBM hardware on z and POWER.

3 · Who owns the precision?

Three answers, and the choice shapes the whole API.

• The value. MPFR's model: each variable carries its own precision, fixed at construction, and operations round to the destination's. Explicit, slightly tedious, impossible to get wrong by accident.
• The operation. A working precision passed into every call.
• An ambient context. Python's decimal, Java's MathContext, the IEEE decimal spec — a thread-local carrying digits, rounding mode and trap flags. Convenient, and the source of every "why did it round there" bug ever filed.

A related question with no default answer: for operations whose exact result is unbounded — division, and every transcendental — something has to say when to stop. Multiplication is the easy case: an n-digit by an m-digit product is exactly n+m digits and needs no policy at all.

4 · Rounding, and the dilemma underneath it

Correctly rounded means the result is the representable number nearest the true value — which requires knowing the true value well enough to break the tie. That is the Table Maker's Dilemma: for some inputs the exact result sits arbitrarily close to a halfway point, so no fixed working precision suffices to decide the rounding.

The practical answer is Ziv's strategy: compute at modest excess precision, test whether the rounding is decidable, and if it isn't, widen and retry. MPFR guarantees correct rounding this way for every function it exports. Most other libraries offer faithful rounding — within one ulp — and don't say so loudly.

Guard digits are the cheap version of the same idea, and the thing to know about them is that they do not compose. Two operations each carrying five guard digits do not hand you a five-guard-digit result; error growth through a series is closer to multiplicative, which is why long summations want a running error bound rather than a fixed pad.

5 · The algorithm ladder

Addition and subtraction are linear and boring. Everything else is a ladder, and each rung is only worth climbing past a machine-tuned crossover — GMP ships a tune program to find them empirically per CPU.

Multiplication goes schoolbook O(n²) → Karatsuba O(n^1.585), crossing over around 20–40 limbs → Toom-Cook 3/4 → FFT-based: Schönhage–Strassen, a floating-point FFT, or a number-theoretic transform modulo one or more primes with CRT recombination. Harvey–van der Hoeven's O(n log n) is a theoretical result, not something you would implement.

Division starts at Knuth's Algorithm D, with normalisation and a two-limb quotient estimate; then Burnikel–Ziegler recursive division; then Newton–Raphson on the reciprocal, which reduces division to a small constant times the cost of a multiplication. There is also exact division — Jebelean, or Hensel lifting — for the common case where you already know the remainder is zero; it runs low-to-high and is roughly twice as fast.

Square root is Newton iteration, or Zimmermann's Karatsuba square root. Base conversion is the one people forget: naive is O(n²) in both directions, and the fix is divide-and-conquer against a precomputed table of powers of the target base.

Transcendentals reduce the argument into a small range and then run a series. Binary splitting is the key trick — evaluating a rational hypergeometric series as a balanced product tree instead of a linear accumulation turns O(n²) into something fast multiplication can attack. At very high precision the AGM takes over for log and π, with quadratic convergence.

6 · How much of this do you actually need?

Most arbitrary-precision work lives between twenty and a thousand digits, and at that size the asymptotically fast algorithms are irrelevant. Karatsuba over schoolbook, a decent division, tight allocation and a base conversion that isn't quadratic will beat a library with an FFT and a sloppy memory model. Save the ladder for when a measurement asks for it.

C++ has the tools to make a bignum feel like a built-in type, and enough rope to make it four times slower than it needs to be. The specifics divide into a structural decision, a set of ergonomic ones, and a short list of traps that are genuinely language-specific.

The landscape, if you aren't building one

Boost's cpp_int with a fixed-size backend — number<cpp_int_backend<256,256,unsigned_magnitude,unchecked,void>> — is stack-allocated and often the right answer when you know the bound in advance.

The structural decision: split the algorithms from the values

This is GMP's mpn/mpz separation, and it is the one design choice worth copying wholesale. The algorithm layer works on borrowed limb ranges and never allocates:

// Returns the carry out. dst may alias a or b.
uint64_t add( std::span<uint64_t> dst,
              std::span<const uint64_t> a,
              std::span<const uint64_t> b ) noexcept;

// dst.size() >= a.size() + b.size(); scratch from mul_scratch_size().
void mul( std::span<uint64_t> dst,
          std::span<const uint64_t> a,
          std::span<const uint64_t> b,
          std::span<uint64_t> scratch ) noexcept;

The value layer owns memory, handles signs, normalisation and growth, and calls into that. The payoff is that the hard, performance-critical code has no allocation, no exceptions, no ownership and no class invariants in it — you can fuzz it directly, and the recursive algorithms carve their temporaries out of one caller-supplied scratch buffer instead of allocating per level.

std::span is what makes this pleasant. Before C++20 it was pointer-and-length pairs, which is why GMP's signatures look the way they do.

Value semantics, and the move payoff

The classic complaint about operator overloading for bignums is temporaries. C++11 largely fixed it, if the operators are written the right way:

// Take by value, mutate, move out — the temporary's buffer is reused.
inline BigInt operator+( BigInt lhs, const BigInt &rhs )
{
    lhs += rhs;
    return lhs;              // NRVO or move, never a copy
}

In a = b + c + d the first + allocates once and the second reuses that same buffer. That recovers most of what expression templates were invented for, at no cost in complexity.

Making the move constructor noexcept is not cosmetic — without it, std::vector<BigInt> copies on reallocation instead of moving:

BigInt( BigInt &&other ) noexcept;
BigInt &operator=( BigInt &&other ) noexcept;
std::strong_ordering operator<=>( const BigInt & ) const noexcept;

Expression templates still earn their place in one case: fusing a = b * c + d into a single pass with no temporary at all. Boost does this, and it genuinely helps for small fixed-size backends. For heap-backed arbitrary sizes the gain over move semantics is modest, and the cost is that auto x = a + b; captures an expression node holding references to dead temporaries. The lower-effort alternative is to expose GMP-style three-operand functions and let hot code call those directly — mul( dst, a, b ) creates no temporary to begin with.

Small-buffer optimisation

At the sizes people actually use, allocation dominates. A union between an inline array and a heap pointer removes it entirely for small values:

class BigInt {
    static constexpr std::size_t kInlineLimbs = 3;   // 192 bits
    union {
        uint64_t  m_inline[ kInlineLimbs ];
        uint64_t *m_heap;
    };
    uint32_t m_size;        // limbs in use
    uint32_t m_capacity;    // <= kInlineLimbs means inline
    bool     m_negative;
};

Note what it rules out: no stable pointers into the limbs across a move, and every accessor takes the inline/heap branch. Worth doing, but decide early — retrofitting means auditing every data() call.

If you would rather not, std::vector<uint64_t> is fine with one caveat: resize() value-initialises, so growing a working buffer zeroes memory you are about to overwrite. Either reserve and construct manually, or use a no-init allocator. And std::pmr::monotonic_buffer_resource is a natural fit for Karatsuba scratch — one stack buffer at the top of the recursion, allocations that are pointer bumps, nothing to free.

The limb primitives

Carry chains, portably. Compilers recognise this shape reasonably well:

inline uint64_t addc( uint64_t a, uint64_t b, uint64_t &carry ) noexcept
{
    uint64_t s  = a + b;
    uint64_t c1 = ( s < a );
    uint64_t t  = s + carry;
    carry = c1 | ( t < s );
    return t;
}

Where they don't, _addcarry_u64 (MSVC and Intel), __builtin_addcll (Clang) or the type-generic __builtin_add_overflow (GCC and Clang) get you the real ADC instruction. Wide multiply needs a platform split:

inline void mul_wide( uint64_t a, uint64_t b,
                      uint64_t &hi, uint64_t &lo ) noexcept
{
#if defined( __SIZEOF_INT128__ )
    unsigned __int128 p = (unsigned __int128)a * b;
    lo = (uint64_t)p;  hi = (uint64_t)( p >> 64 );
#elif defined( _MSC_VER ) && defined( _M_X64 )
    lo = _umul128( a, b, &hi );
#else
    /* 32-bit split fallback */
#endif
}

For a decimal engine, base 10⁹ in uint32_t with uint64_t intermediates needs none of this and is genuinely portable. Base 10¹⁸ in uint64_t needs __int128 plus a division-by-constant trick for the carry — which the compiler will do for you if the base is a literal constant, an argument for making it one rather than a runtime parameter.

Error handling in an embeddable engine

Three options, and mixing them is the failure mode. Exceptions are natural for std::overflow_error on divide-by-zero and std::bad_alloc on growth, but rule you out of -fno-exceptions hosts. Status returns from every operation make operators unusable. Sticky flags in a context object — the IEEE model, used by decNumber and Python's decimal — keep operators clean and let the caller check at a boundary of its own choosing, at the cost of a context parameter or a thread-local.

For a decimal engine specifically the flag model composes best, because the spec already defines the set: inexact, rounded, division-by-zero, overflow, underflow, invalid-operation. Callers arriving from other decimal libraries will expect exactly those.

constexpr

C++20 allows heap allocation during constant evaluation as long as it doesn't escape into the runtime object, so a constexpr bignum can compute a compile-time answer — but it cannot be a compile-time constant if it holds a pointer. For persistent compile-time values, back them with std::array at a fixed size. Marking the limb layer constexpr costs nothing and lets you write static_assert tests that run at build time.

The traps that are specifically C++

• Shift by the full width is undefined. x >> 64 on a uint64_t is UB — and it is exactly what a limb-shift loop computes when the shift amount is zero. Special-case it. This bug is close to universal in first drafts and shows on only one compiler.
• Signed overflow is UB, so limbs must be unsigned. Related: uint16_t * uint16_t integer-promotes to int and can overflow signed — a real hazard in BCD or 16-bit-limb designs.
• Aliasing. a = a * a and mul( x, x, x ) must work. Either document that the destination may not alias and assert it, or always compute into scratch. Silent wrong answers otherwise.
• Self-assignment in operator=, when you free the old buffer first.
• Comparisons that allocate. a < b + 1 should not touch the heap. Provide small-integer overloads, or you pay for one in every loop condition.
• Header-only inlining that vanishes in debug builds, quietly moving your profile somewhere it will never be in release.

Testing

The language-specific leverage here is that a differential oracle is cheap: link GMP into the test binary only, and property-test every operation against it over randomly sized random inputs, plus a boundary corpus — zero, one limb, limb boundaries ±1, powers of the base, values that carry all the way up. For a decimal engine, Python's decimal makes a good oracle over a pipe, since it implements the same spec.

Then run the same suite under -fsanitize=address,undefined. The shift-by-width and signed-overflow traps above are precisely what UBSan exists to catch, and they will not surface any other way.

A survey is easier to read against something concrete. Addaptic is an arbitrary-precision decimal arithmetic engine in C++, and it makes a set of choices that land it somewhere most libraries are not — which makes it a useful worked example of how the axes interact.

It is a BigDecimal, not a bignum

The scale is held as a count of fractional digits, taken from the right. That is exactly the unscaled-value-plus-scale representation of Java's BigDecimal, Python's decimal and IBM's decNumber — a decimal floating-point number whose exponent happens to be stored as a negated digit count. Not a big integer, not a rational, not an exact real, not a ball.

So the third and fifth families from Part I are out of scope by design, and correctly so. You do not want lazy reals or exploding denominators in a calculator engine.

Two choices that are genuinely unusual

A decimal engine with a complete transcendental library. This combination barely exists. decNumber and BigDecimal give you arithmetic and stop. Python's decimal adds exp, ln, log10, sqrt and power — five functions, and the documentation is explicit that it stops there. Everything carrying the full trigonometric, inverse, hyperbolic and gamma set — MPFR, mpmath, Arb — is binary. Addaptic has that set in base ten, every function of it correctly rounded, which puts it in a niche I cannot name another occupant of.

A runtime-variable base. Practically every engine fixes its radix at compile time. Addaptic carries the base per number — any base from two to sixteen — with live conversion between them, derives its rounding thresholds from the base so that rounding is base-generic, and converts a mixed-base pair into the left operand's base before operating on it, rather than silently reading one value's digits as though they were the other's. Base twelve as a first-class option is something I have not seen elsewhere at all, and the odd bases are harder still: an odd base has no half digit, so a half never terminates there — 0.333… in base seven — and every rounding decision has to be made without one.

Both choices cost the packed-limb representation as a stored form. A number addressable digit by digit, in a base that can change under it, cannot also be limbs of 10⁹: base conversion, rounding at a digit position and every digit-level accessor would have to unpack first. What they do not cost is packing inside an operation, where the base has only to be respected at the boundaries — this engine's long division, its multiplication — schoolbook and Karatsuba — and its transform tier all pack several base-b digits into each element, in every base, base twelve and the odd bases included, and unpack on the way out. Base twelve does give up the shift-and-mask convenience of a power-of-two limb, but so does base ten; dividing by the limb radix costs a multiply. That is a coherent trade — but it is worth naming what it buys and what it spends.

What one digit per word actually costs

Eight bytes holding a value from zero to nine. Against a packed 10⁹-per-uint32 scheme that is roughly ten times the memory and nine times the iterations in every linear operation — which is the figure below.

The consequence that is easy to miss is what it does to the thresholds. A Karatsuba crossover around thirty elements is completely normal; GMP's is in the same range. But thirty elements here is thirty decimal digits, where thirty limbs in GMP is about 580 decimal digits. So a digit-array engine reaches for Karatsuba at values a packed engine is still handling in two limbs of schoolbook. The ladder is correctly built — it simply starts much lower down the value range than the same numbers would suggest elsewhere.

This engine did exactly that until it began packing inside its multiplication and division. Its packed schoolbook then outran its digit-level Karatsuba at every size up to where its transform takes over, and the Karatsuba was removed: a split whose additions still go a digit at a time costs more than the products it saves. It came back splitting the packed limbs instead, where a carry is a comparison rather than a division, and now takes over at under a hundred limbs — some 850 digits in base ten, close to where a packed binary engine climbs the same rung. It is quick enough there that the transform above it now waits, for operands of similar length, until past a hundred thousand digits in every base — about four times further out, in base ten, than it began. Finding where showed the crossing to be a staircase rather than a point: a transform's cost is flat between power-of-two lengths and doubles at each, so just past every step Karatsuba wins again.

The rest of the ladder

Part I's other rungs are climbed too, each from a crossover measured for every base:

• Division takes Newton's reciprocal and one product when the divisor is long, chosen by the lengths of both the quotient and the divisor.
• The exponential, sine and cosine go by binary splitting from a few hundred digits, the tangent and cotangent from one run of it.
• The logarithm goes by Newton on the exponential, and by the AGM — in most bases from under a hundred digits — once its constants are on hand.
• The inverse circular functions go by Newton on the sine and cosine.
• The square root uses the AGM's own core — Newton on the reciprocal root, multiplications only — and the cube root goes by its reciprocal, both climbing through widths that double.
• Base conversion goes by products past a crossover, writing the value in the new base by halves joined with one product — the divide-and-conquer Part I calls the cost people forget, and hundreds of times faster than the quadratic loops it replaced, at ten thousand places.

How it rounds

Every answer is correctly rounded at the working precision, in its own base and rounding mode — Part I's strongest guarantee, reached by nearly the same route as MPFR. Each function works past its answer with guard digits sized for it. When the final rounding lands within a hair of a half, the function is worked again wider, the extra width doubling each round until the answer is clear of it: Ziv's strategy. A value exactly on a half never clears, however wide the retry, so wherever one can arise — roots, powers to a fraction, logarithms that are fractions — it is decided exactly, with whole numbers.

The one difference is a cap. The retry stops at one more answer's width of digits, and a true value nearer a half than that is rounded as the widest attempt stands — rare in ordinary use, though such values can be built on purpose. An audit against an independent decimal implementation, across eighty orders of magnitude and run under sanitizers on every change, checks the result rather than asserting it.

What the placement is worth

Reading an engine against the axes is more useful than benchmarking it, because it tells you what the engine is for. This one is not competing with GMP and would lose if it tried; the packed representation and the assembly cores are not recoverable from a base-generic decimal design. What it has instead is the thing GMP does not: exact decimal behaviour, a base the caller chooses at runtime, and a function library that goes all the way up to gamma, correctly rounded, without ever leaving base ten.

That is the general lesson, and it survives the specific example. The axes are not a scorecard. Nothing sits well on all of them at once, and an engine that looks compromised on one axis is usually paying for a position on another that its callers actually need.