Skip to content

Decimo v0.14.0

Latest

Choose a tag to compare

@forfudan forfudan released this 01 Sep 23:19
ffc697b

20260901 (v0.14.0)

Decimo v0.14.0 is a Python and optimization release, on top of a base change in both integer types. Seventeen bugs go with it, two of them in versions already released.

The Python package becomes a near drop-in for the standard library's decimal, and a superset of it. Every method decimal.Decimal has, decimo.Decimal now has: a full Context, every rounding mode but ROUND_05UP, and the whole specification surface. On top of that are pi(), e(), a rounding= argument on sqrt, exp, ln and log10 that decimal has no equivalent of, and Decimal128, which brings trigonometry, cbrt, root and the IEEE 754 interchange bytes. What it refuses rather than answers differently: NaN and infinity, ROUND_05UP, signals and traps, and one context per process rather than per thread. Wheels for macOS and Linux.

BigUInt moves from base 10^9 to base 10^18 and BigInt from base 2^32 to base 2^64, both gain a number-theoretic transform, and small values now live inside the struct. Against v0.13.0, on the same machine:

speedup
addition and subtraction 2-3x
multiplication 2-3x
division 1.2-2x
reading and writing text 2-7x
small BigDecimal operations about 2x
BigDecimal.pi() four orders of magnitude

Ratios rather than times, since the times belong to whatever machine took them. The absolute numbers, on one machine and one commit, are in docs/benchmarks.md. BigInt is ahead of CPython's int on every large operation, and is measured against GMP now, which is the comparison it should be held to.

Every transcendental also decides its rounding rather than assuming it: it takes the interval its own error bound allows and checks that the whole of it rounds to one answer. The trigonometric functions size their argument reduction to the argument instead of to a flat ninety-nine digits, which they were silently wrong past.

⭐️ New in v0.14.0

The Python package:

  1. Wheels for macOS and Linux, CPython 3.13 and 3.14 — macOS arm64 (11 and later) and Linux x86_64 and arm64 (glibc 2.35 and later). The Mojo runtime libraries travel inside the wheel, so nothing else is needed. release_python.yaml builds the set and uploads through PyPI trusted publishing. The library is compiled and tested on Linux now too.

  2. Every rounding mode. getcontext().rounding accepts all seven but ROUND_05UP, and arithmetic, quantize, round(x, n) and to_integral_value follow it, checked against decimal digit for digit. Context is a value until installed, localcontext() takes keyword overrides, and the Mojo arithmetic methods take a rounding_mode.

  3. The rest of decimal's surface. Keyword arguments where decimal takes them; Decimal((sign, digits, exponent)), closing the as_tuple() round trip; pow(x, y, modulus); the sixty Context methods that compute without disturbing the current context; and the specification methods from remainder_near through the four logical_ ones to number_class. The arithmetic lives in Mojo, in the new decimo.bigdecimal.spec.

  4. decimo.pi() and decimo.e(), to the context precision or to the digits asked for. decimal has neither.

  5. sqrt, exp, ln and log10 are always half to even, as in decimal; ** follows the context mode, also as in decimal. All four additionally take decimo's own rounding=, and are correctly rounded under whichever mode applies.

  6. Every operation applies the context, as in decimal. abs(), max, min, normalize, scaleb, fma and the remainder from % and divmod were returning an exact value where decimal returns a rounded one. fma is exact now where it went through * and +, which round at each step.

  7. Decimal128 in Python, as decimo.Decimal128 (Dec128 for short): 96 bits of coefficient and a scale from 0 to 28, in sixteen bytes that own nothing. It does arithmetic, compares, hashes, rounds, copies, pickles and formats like Decimal, takes an int, a float or a str on either side of an operator, and carries the methods and the mathematics, the IEEE 754 bytes included. Its hash agrees with int, float, decimal.Decimal and Decimal, so the four are interchangeable as dictionary keys, and a mixed expression settles in the wider type. Its results never allocate, which makes it quicker than Decimal on every operation and quicker than decimal on all of them but str.

  8. Four conversions no longer go through a string. hash() is a tp_hash slot reducing the coefficient over its own words, Decimal(x) copies the struct when x is already one, int() uses PyLong_FromSsize_t for anything that fits a machine word, and float() no longer imports builtins on every call. Between 2x and 15x each.

Decimal128:

  1. The IEEE 754 decimal128 interchange format. decimo.ieee754 encodes and decodes the sixteen bytes in the binary integer decimal layout — what MongoDB's BSON decimal128 and Intel's library store — with Decimal128.to_ieee754(), from_ieee754() and the BigDecimal equivalents on top. It is a codec and brings no IEEE arithmetic with it: trailing zeros are part of the encoding and are kept, an infinity or a NaN is refused by name, and densely packed decimal is not read here.

  2. Trigonometry. sin, cos, tan, cot, sec and csc, correctly rounded across the whole range the type holds. The work is the argument reduction: Decimal128 reaches 7.9E+28, so subtracting k * (pi/2) with the 28-digit quarter turn the type itself holds answers a different question than the caller asked. The quarter turn is kept here in four exact pieces of 38 digits. 881 checks against CPython's decimal: none wrong.

Rounding that is decided rather than assumed:

  1. exp, ln, log10, **, sin, cos and tan computed a fixed number of digits beyond what was asked and rounded once, which is right whenever the discarded tail misses a boundary and silently wrong when it does not — for half to even as much as for a directional mode, since a tie is just another boundary. They now take the interval their own error bound allows and check that the whole of it rounds to one answer, widening when a boundary falls inside. sqrt needs no loop, being algebraic. A second attempt is needed on about eight calls in a thousand. Checked over seventeen thousand cases against decimal.

  2. arctan, cot, csc and sec decide theirs too. They were the last ones adding nine guard digits and rounding once. An argument built to catch that — a value whose 29th digit is a 5, put through the inverse function — had arctan come back one unit low and log10 one unit high. Both are correct now, and the arguments are pinned as tests.

  3. The trigonometric reduction budget is measured, not guessed. An argument close to a multiple of pi/2 loses digits to the subtraction there, by as much as it is close, and no constant can cover that: sin of pi taken to 250 digits came back with no digit of the answer in it. budget_for() now measures the distance to the nearest multiple, at a narrow width first and widening when the measurement comes back at its own noise floor. The measurement is the reduction the function was going to do anyway, so an argument nowhere near a multiple pays nothing.

Integers:

  1. BigInt keeps small values inside the struct. WordList, written for BigUInt, gains an inline-capacity parameter and moves to decimo.wordlist; BigInt uses it under the name Magnitude, at seven words after the move to base 2^64 — the sum of two hundred-digit values. Addition at a hundred digits is 3.6x, division 1.4x. Above a thousand digits nothing moves.

  2. BigInt is measured against GMP, timed in C, which is the comparison a big-integer library should be held to; CPython's int can only be reached through the interpreter and loses on call overhead before the arithmetic starts. GMP wins most rows. The exception is small values, where an mpz_t still goes to the heap and decimo no longer does.

  3. BigInt.sqrt() gains Zimmermann's recursion. Above 64 words it uses the Karatsuba square root of INRIA RR-3805 instead of CPython's precision-doubling: the division at the last step is half the width, and the remainder falls out of the recursion. Roughly 1.8x from ten thousand digits up. Below the crossover the older path still wins and stays, as the recursion's base case.

  4. BigUInt multiplication gains a number-theoretic transform. Toom-3 was the largest algorithm available to it, so BigDecimal multiplication was stuck at O(n^1.465) while libmpdec switches to a transform. decimo.biguint.ntt supplies that tier, reusing the field arithmetic in decimo.bigint.ntt — only the packing differs, since a decimal magnitude can only be cut at a power of ten.

Conversions (PR #269):

  1. Rational conversions. Constructors from an integral scalar, a String and a BigDecimal; the factories from_string(), from_integral_scalar(), from_float_scalar() and from_bigdecimal(); and the outbound __int__(), __float__(), to_int(), to_integer(), to_float() and to_bigdecimal().

  2. from_integral_scalar() and from_float_scalar() on BigInt, BigDecimal, Decimal128 and Rational — one pair of names across the library, constrained with where clauses rather than comptime assert, so the wrong scalar kind is an overload mismatch and not an assertion. Plus BigInt.from_biguint() / to_biguint() and the BigInt10 pair.

🦋 Changed in v0.14.0

Allocation, and the small operations it dominates:

  1. Heap blocks are reused instead of freed. alloc and dealloc cost about the same whatever the size, so every value too long to sit inside the struct paid the same toll — a six-word addition spent more than half its time there. A released block now goes on a small stack sorted by size, and the next request of that size takes it back. The pool is process-wide, holds at most eight blocks per size, and is shared under an atomic flag.

  2. Small operations are about twice as fast. At these sizes an operation's speed is very nearly its allocation count, and several were allocating for nothing: add() and subtract() scaled both coefficients when only one ever needs it, a buffer was sized exactly and then grown, and three debug_assert calls built their message with + String(n), which allocates even with assertions compiled out (modular/modular#6439).

  3. BigDecimal construction no longer copies its coefficient. The component constructor took it borrowed and then copied, so all 27 call sites already written as coefficient=coef^ paid a heap allocation for a move they had explicitly asked for.

Multiplication, division and roots:

  1. Multiplication is 2-3x faster. Toom-3 comes to BigInt, which had stopped at Karatsuba. The schoolbook base case is rewritten to product scanning, so each column stays in registers, and then packed into base-2^64 limbs, which quarters the word pairs. Both crossovers were re-swept afterwards, worth another 20% on their own.

  2. BigInt multiplication gains a number-theoretic transform. Above Toom-3 the product is a cyclic convolution modulo the Goldilocks prime 2^64 - 2^32 + 1, bringing the exponent from n^1.465 to n log n. One prime, so no CRT step. The chunk width is left free rather than fixed at 16 bits, since the transform length must be a power of two and the rounding up is waste, and the dispatch compares fitted cost models rather than a word count, because the transform's cost steps at powers of two while Toom-3's climbs smoothly.

  3. Addition and subtraction are 2-3x faster, and the recursive algorithms inherit most of it. BigInt's little-endian words are already base-2^64 limbs in pairs, so one 64-bit add does the work of two. A base-10^9 BigUInt word carries by comparison against BASE, which put the carry on the loop-carried chain; both answers are computed off that chain now and the incoming carry only selects. subtract_simd() and add_slices_simd() become subtract_carry_select() and add_slices_carry_select(), since neither is vectorized any more, and normalize_borrows() is gone.

  4. Division is 1.2x to 2x faster. Knuth D was spending six of its own multiplications on work a multiplication does in one. Its multiply-subtract runs two words at a time now, and it no longer allocates a fifth word list for its remainder: the working dividend is the remainder by then.

  5. BigUInt schoolbook division is 1.3-2x faster. Each quotient word built q * y as a fresh BigUInt, shifted it, compared it against the whole remainder and subtracted it — four passes plus an allocation, per word. It is one fused multiply-subtract over the n + 1 word window now.

  6. Burnikel-Ziegler starts where it begins to pay, and pads to j * 2^k words. The cutoff was a 24-word divisor where the recursion does not earn its keep until about 48. And the recursion falls back to schoolbook on the first odd block, so the padding has to keep the block size even the whole way down; BigInt rounded up only to even, so a large divisor lost the algorithm's asymptotics at the first step.

  7. sqrt_via_reciprocal_iteration() is 1.6x faster at high precision. Written around the residual — r + r * (1 - x * r^2) / 2 rather than r * (3 - x * r^2) / 2, algebraically the same — the correction is tiny, which a BigDecimal keeps in the scale rather than the coefficient, so the multiply is half-width by half-width instead of half by full.

pi(), and the constants:

  1. BigDecimal.pi() is four orders of magnitude faster. Five changes that compound: the Chudnovsky binary splitting uses the P/Q/T recurrence, so each leaf is O(1); the leaves are built in machine arithmetic; the term count is sized to the precision rather than to a flat margin; the square root of 10005 is no longer taken exactly; and the pipeline stays binary to the end, so the irrational factor enters as a reciprocal square root, which Newton reaches without a division. A hundred thousand digits was out of practical reach before. Digits are unchanged, exact against MPFR from 1 to 100 000, and everything that range-reduces against pi inherits the gain.

  2. ln picks its series by how small the argument is, not how long it is. The choice between the Taylor series and the atanh identity read the digit count of z = x - 1, where what decides the term count is z's magnitude. ln(2.3456789) has eight digits but leaves z = 0.34, so it took the Taylor path and paid twice over. About 2x, which turns the one benchmark row where the logarithm lost against libmpdec into a win.

Text:

  1. Reading a BigUInt from text is up to seven times faster. The parser normalized every string first, writing one digit per byte and reading it back into words. A string that is nothing but digits goes straight into the words now, which is nearly all of them.

  2. A plain decimal string is parsed straight into the words. The same for BigDecimal; anything with an exponent or a separator still takes the general parser. from_string and the Parsable trait take a StringSlice now, so a substring or a foreign buffer costs no allocation — which is what lets the Python binding read CPython's own string rather than copy it.

  3. Writing a BigDecimal out is two to four times cheaper. The text was three allocations and two copies for what is one row of digits. The digits go straight from the coefficient's words into one buffer now, on the stack when the value is short enough, and BigUInt emits two digits per division against a table of pairs. tp_str is a real slot for both Python types rather than a __str__ CPython has to dispatch to.

  4. String(BigInt) is 1.3x to 1.9x faster above about 600 digits. Its two divide-and-conquer thresholds were derived rather than measured, and the derivation was wrong: what D&C buys is the balanced split, not large enough internal divisions, and that pays long before any division inside it reaches Burnikel-Ziegler.

  5. BigInt.to_biguint() no longer detours through a decimal string. The conversion splits on powers of 10^18 instead of powers of 10, so each half lands on a word boundary and its words go straight into the result.

Decimal128 internals:

  1. Division is one wide division rather than a walk. The quotient was built a digit at a time and could run out of digits before reaching the position the rounding needed — 1 in 300 random pairs was wrong in the last place. The numerator is raised until the integer quotient is about thirty digits, divided once, and rounded from the remainder. That needed a 256-by-128 divider, udiv_u256_by_u128, Knuth D over 64-bit limbs. Division is 3.2x and correct on all 300 pairs.

  2. Wide is written once and used at two widths. WideValue[DIGITS] carries the mantissa the series run on; Wide is 38 digits of it and Extended 75. The constants exist at both widths, and a test narrows each wide one to check it gives the narrow one exactly. The reciprocal divider reaches 10^48 now, so the second width no longer falls back to software division — which made every second attempt in exp, ln and power about ten times cheaper.

  3. An exact integer power is known to be exact, and trailing zeros come off in one step. A base of d digits raised to the n-th has at most d * n digits, so below the working width nothing is rounded and the answer needs no second look: 1.05^12 is one pass again rather than three. The zeros at the end were stripped one at a time; halving finds them in five steps.

Other:

  1. Rational.__add__ and __sub__ cancel before they multiply, by Algorithm A of Knuth 4.5.1, so the result is in lowest terms without a second gcd over full-width operands. Summing 1/k^2 to 1 200 terms is 40x. This pays because gcd() balances its operands before entering Stein's binary loop, which makes about one bit of progress per full-width subtraction and so is quadratic when one operand dwarfs the other.

  2. BigInt10 is no longer used by any other module. Bridging goes through BigUInt; BigInt10 keeps its own conversions for code that wants them (PR #269).

  3. List._data is no longer used anywhere. All 63 sites moved to unsafe_ptr(), which returns the same address with an origin attached, so the library no longer depends on a private field of the standard library. Restoring origins exposed twelve deliberate in-place aliases, which say so through alias_as_immutable_source() now.

  4. tests/test.sh runs in parallel by default. DECIMO_TEST_JOBS defaults to the machine's logical CPU count. CI stays at 1, since each suite is already its own job there.

  5. The out-of-range message for power_of_10_unsafe[uint256] says 0..77. It still said 0..58 after the table was extended.

🩹 Fixed in v0.14.0

  1. sin, cos and tan were wrong for a large argument, silently. The reduction x mod 2*pi cancels everything above the remainder, so 10^k spends k digits of pi before the remainder starts. The budget was a flat ninety-nine digits: right up to 10^99, half wrong by 10^105, and at 10^150 nothing in the answer was correct, with nothing to say so. reduction_digits() sizes the budget to the argument now, and cot, csc and sec inherit it. Pinned against an independent computation to 10^300.

  2. BigInt.sqrt() never returned for values at the top of a word. The one- and two-word paths refined an estimate with while (guess + 1) * (guess + 1) <= value. Near the top of the range that square wraps, so the test reads true forever. sqrt(2^32 - 1) hung, as did 131 071 other single-word values and about 2^33 two-word ones. Both paths go through isqrt_uint64() now, which clamps the estimate first. Present since v0.13.0, so released versions are affected.

  3. BigUInt division crashed, or silently lost a factor of 10^9, for three-word dividends. floor_divide_by_uint128() consumes the dividend four words at a time, and when the count is not a multiple of four the leading group is short — its quotient was discarded. A seven-word dividend over a three-word divisor came back 10^9 too small, and a three-word one produced a BigUInt with no words at all, which faults the next operation that reads it. 80 wrong results in an 1 824-case sweep. Present since PR #111 (2025-07-23), so released versions are affected.

  4. BigUInt in-place subtraction returned garbage for equal operands, and with it long division for a whole class of dividends. subtract_inplace() handled x == y by shortening x to one zero word and then fell through into the general path, reading and writing past the end. The reach went past -=: Burnikel-Ziegler's base case computes its remainder as a_slice -= q * b_slice, and a block that divides exactly makes those equal, so // and % were wrong whenever the recursion met one — 676 wrong results in a 5 148-case sweep. The fix is a return.

  5. Toom-3 wrote one word past the end of its result buffer for lopsided operands. A three-way split sizes its limbs by the longer operand, so a short second operand can leave 4 * m, where the w4 coefficient is recomposed, past the end of the result. w4 is zero in exactly those cases, so it wrote a zero out of bounds and no value comparison could see it. _add_at_offset_inplace() states the precondition and asserts it now.

  6. sqrt_via_reciprocal_iteration() returned fewer correct digits than asked for, from two causes: the iteration schedule credited the seed with more digits than it carries, and the seed itself was x ** -0.5, good for some inputs to only about ten digits. A 1500-digit request came back correct to 1248, with the full count returned and nothing raised. isqrt_via_reciprocal_seed() shared both and hid them behind full-size corrective divisions, so sqrt_exact() is about 30% faster now.

  7. Decimal128's exp, ln, log10 and sqrt were wrong in the last digits. All four summed their series in Decimal128 arithmetic, which rounds to 28 digits after every term, so the answer inherited every one of those roundings — ln was out on 108 of 200 random arguments. The series run in a fixed-width accumulator ten digits wider now, rounded once at the end: none wrong in the same 720 checks. sqrt no longer refines a floating estimate at all. Decimal128 still imports nothing from BigDecimal.

  8. Decimal128's powers and roots were wrong in the last digits. x^y went through exp(y * ln(x)) with the logarithm and the product each rounded to 28 digits on the way, and an absolute error in y * ln(x) is a relative error in the answer: 51 of 60 random arguments were wrong. The integer path was wrong on 22 of 40 and root on a quarter of those tried. All three compute at 38 digits and round once now, running again at 75 when the digits below the answer do not settle it. 540 checks: none wrong.

  9. Decimal128's logarithms and exponentials decide their rounding. When the value sits on a boundary within the computation's own error, the whole thing runs again at 75 digits, which has forty-six digits below the answer instead of nine. Two such arguments were found by searching three million. It used to abort rather than answer, because a 75-digit mantissa asks for powers of ten the reciprocal divider's table did not reach.

  10. ln of a value close to one lost most of its digits. The reduction wrote x as m * 2^p * 10^q and added p * ln(2) + q * ln(10) back at the end. For x just under one those terms are each about two while their sum is tiny, so most of the digits carried went into cancelling them out. Arguments already in [0.5, 2) go straight to the series now.

  11. A computed value whose dropped digits were all zeros claimed to be exact. to_decimal_decided refused to round when the digits below the answer sat near a boundary, but treated a remainder of exactly zero as settled however much room the computation had asked for. Zero is the boundary: with room to be wrong the true value may sit either side of the multiple, and the claim decides whether the trailing zeros are dropped.

  12. A value too large for Decimal128 came back with a scale of four billion. When more digits had to be dropped than there were places after the point, the scale went negative and wrapped. Nothing reached it before: exp refuses its argument above 66.54 and a logarithm is small. tan of an angle a hair past a pole reaches it, and raises OverflowError now.

  13. A value below the smallest scale returned zero instead of rounding. ln(1.0000000000000000000000000001) has every digit below the 1E-28 that Decimal128 stops at. Rounding them says 1E-28; the conversion returned zero whenever the digits being dropped were all of them.

  14. The digit count stopped at 58 and returned 59 for anything larger. number_of_digits covered the 58-digit product of two Decimal128 coefficients and answered wrongly, rather than refusing, above that. It covers both types to the top now, and is about sixty times faster: the old binary search compared against 10 ** k, which is not folded for 128- and 256-bit scalars and so was built at run time by repeated multiplication.

  15. Burnikel-Ziegler's "add one more block" guard tested the wrong length. It compared len(a.words) against t * n while t counts blocks of the normalized dividend, which the normalization has usually lengthened, so the guard fired more or less at random. It tests the normalized length now, as BigInt's copy of the algorithm always has.

  16. The in-place single-word divisions left a BigUInt with no words at all when the quotient was zero, and floor_divide_by_word_inplace() read its loop bound from the already-shortened list, skipping a word whenever the leading word was smaller than the divisor. Neither is called anywhere in the library today; the out-of-place versions that // uses were correct.

  17. BigUInt.is_two() could not return True for any value. It asked for a two-word value and then for the second word to be zero, which the no-leading-zero invariant forbids; words[0] was never compared against 2 at all. It has no callers, which is why nothing caught it (issue #312).

🗑️ Deprecated in v0.14.0

  1. BigDecimal.from_float() and Decimal128.from_float() are deprecated in favour of from_float_scalar(), the name that lines up with from_integral_scalar(). They forward unchanged (PR #269).

💥 Breaking in v0.14.0

  1. BigUInt's words are UInt64, and its base is 10^18. It held nine decimal digits in a UInt32 and now holds eighteen in a UInt64 — the same digits per byte, half the words. BigUInt.words, Coefficient and BigUInt(raw_words=...) all follow, so raw_words= takes a List[UInt64] and a list literal becomes [UInt64(1)]. Write BigUInt.Word for the type of a coefficient word and BigUInt.DIGITS_PER_WORD for how many digits it holds, rather than a literal UInt32 or 9. Values, strings and every arithmetic result are unchanged; only the representation is.

  2. BigInt's words are UInt64, and its base is 2^64. BigInt.words, Magnitude and BInt(raw_words=..., sign=...) all follow. Code that reads or builds the magnitude directly has to change: a list literal becomes [UInt64(1)], shifting by 32 becomes shifting by 64, masking with 0xFFFF_FFFF becomes masking with 0xFFFF_FFFF_FFFF_FFFF or dropping the mask, and a word count derived as (bits + 31) // 32 becomes (bits + 63) // 64. Values, strings and every arithmetic result are unchanged; only the representation is. It held its words in base 2^32 while doing all its arithmetic in 64-bit registers, so schoolbook multiplication and Knuth D both made twice the passes they needed to.

  3. BInt(raw_words=..., sign=...) takes a Magnitude, not a List[UInt32]. That is the inline word storage BigInt moved to, and the constructor moves into it rather than copying. A list literal still works unchanged; an existing List goes in as BInt(raw_words=Magnitude(words^), sign=False). Magnitude is exported from decimo.

  4. The Integer alias for BigInt is removed. BInt remains, and matches BDec and Dec128 in shape. Integer named a general concept rather than one concrete type, and collided with the ordinary English word used throughout the documentation. Replace Integer with BInt or BigInt.

  5. gcd() and BigInt.gcd() are now raises. They take a remainder on unbalanced operands, and BigInt division raises. Callers already inside a raises function need no change.

  6. product_range() caps the number of factors, not the size of the bounds. Its old bound, high <= 2^32 - 1, was there because each factor is cast to a word; every non-negative Int fits a word now. The cap is FACTORIAL_MAX_INPUT, the same one factorial() and permutation() already answer to.

  7. BigInt.from_bigint10() and BigInt.to_bigint10() are removed. Use BigInt10.from_bigint() and BigInt10.to_bigint() (PR #269).