Speculative Arbitrary Precision Speed

Arbitrary-precision speed

Speculative. A proposal, not adopted.

The nightly bench puts our arbitrary-precision special functions at 2.5 to 14 times mpmath, and the greatest common divisor of two 10000-bit integers at 12 to 77 times the other systems. This page says where that time goes, whose code it is in, what would cut it, and what to do first.

"Measured" here means on one laptop, one process, warmed, with compute-engine 0.187.0: good for ratios and for counting operations, not for absolute times. The bench's ratios are the nightly's. "Prototype" means a short standalone script written to size an option, not a patch.

1. Whose code it is

Every special-function case on the bench runs in compute-engine's own kernels (src/compute-engine/numerics/special-functions.ts). Our packages wrap those heads only for other inputs: a complex operand, a closed form, an interval. For a real rational argument under N() they fall through to the native code.

Head	Native kernel	Algorithm
Zeta(s)	zetaCore	Cohen, Villegas and Zagier's acceleration of the eta series: about 1.32 terms per working digit, one (k+1)^(-s) per term
HurwitzZeta(s, a)	bigHurwitzZetaEM	Euler-Maclaurin, with a planner that picks the split and the working digits
Gamma(x)	gammaCore	Shift the argument up by a running product, then Stirling's series, then exp
PolyGamma(1, x)	trigammaCore	Shift by adding 1/w² per step, then the asymptotic series
PolyLog(s, z)	bigPolyLog over lerchPhiSeries	The direct series, one power per term
LerchPhi	bigLerchPhi, and our patch where the native code declines	The direct series
GCD	Euclid with % on bigint	Ours replaces it for integers: declareFastGcd in number-theory, over lehmerGcd in residues

So the special-function options below are changes to compute-engine. They go upstream; a patch in packages/ce-patches is how one is carried until it lands. Only the GCD is ours to fix directly.

2. Where the time goes

The primitives are not it

compute-engine's BigDecimal is a bigint significand with a decimal exponent. Its primitives, beside mpmath's (pure Python backend) on the same operands, in microseconds:

Digits		add	multiply	divide	sqrt	ln	exp	37^(3/7)
30	BigDecimal	0.10	0.08	0.35	2.1	3.3	2.7	3.2
30	mpmath	0.97	0.85	1.16	2.3	2.9	3.5	5.6
100	BigDecimal	0.35	0.20	0.63	3.8	6.4	5.4	6.7
100	mpmath	1.16	1.21	1.96	4.3	6.9	10.3	12.9

Arithmetic is three to ten times faster than mpmath's and the transcendentals are level or better. A faster BigDecimal would not close the gap. The gap is in how many operations a kernel does, and at what precision.

The working precision is the asked precision plus 34 digits

Two guards stack. An exact operand such as Rational(3, 7) is read at ten digits more than asked (EXACT_OPERAND_GUARD_DIGITS), and the kernel then adds 24 to what it was handed (SPECIAL_FN_GUARD). Observed inside the kernels:

Asked	Working digits
Zeta at 30	64, and 67 inside each power
HurwitzZeta at 30	77 to 81
Gamma at 50	84

mpmath computes at the precision asked plus a few guard bits of its own. More working digits cost twice: each operation is dearer, and the series that scale with precision run longer. Zeta at 30 digits sums 88 terms; at 40 working digits it would sum 56.

Per kernel

Times are for one function value. Counts are BigDecimal calls per value.

Case	Ours	mpmath	What the time is
Zeta(3/7) at 30	1.18 ms	0.18 ms	88 non-integer powers, each a full ln and exp at 67 digits and 11 to 14 µs: about two thirds. 180 divisions: about a fifth.
HurwitzZeta at 30	0.77 ms	0.47 ms	25 non-integer powers at 81 digits, about half; a rounding after every operation.
PolyLog at 30	0.80 ms	0.31 ms	About 270 terms of the direct series, each an integer power, an add and a rounding.
Gamma at 50	0.19 to 0.25 ms	0.03 ms	About 115 rounded multiplications in the argument shift at 84 digits, then 36 divisions in Stirling's series, three ln and an exp.
PolyGamma(1, x) at 30	0.37 ms	0.21 ms	About 70 divisions: one per shift step and one per series term, where each Bernoulli number is rebuilt from a string.

Three patterns account for most of it:

• A power per term. A non-integer pow is exp(s · ln k). The logarithm of each small integer is recomputed on every call, nothing is shared between k and its multiples, and the product handed to exp is exact, 107 digits long at a working precision of 67.
• The guard. Above.
• Per-term setup. A Bernoulli number is turned into a string and back into a BigDecimal, then divided, for every term of every call, in Gamma and PolyGamma. The Hurwitz kernel caches the same coefficients once per precision.

Boxing and conversion are not it. A prototype of the Zeta series written straight on BigDecimal, with no engine around it, takes 1.10 ms against the engine's 1.18. In a CPU profile our own packages are a few percent in total. The installed compute-engine is minified, so its share is not split further here.

3. Options

	Option	Expected gain	Whose code	Evidence
A	Powers at primes. In Zeta's series, compute p^(-s) only for primes and build each composite's by one rounded multiplication	2.1× on Zeta at non-integer s	compute-engine	Prototype: 1.10 ms to 0.51 ms, the same 30 digits
B	One guard, sized by the kernel. The kernel's guard covers the operand's, so 30 digits asked is about 40 to 44 working digits, not 64	1.6× to 2× on Zeta; something like it on every kernel	compute-engine	Prototype: 0.51 ms to 0.33 ms at 44 digits, 0.28 at 40, still 30 correct digits on that input. The accuracy grid is not run.
C	Cache ln k for small integers per working precision, as mpmath does	Up to 2× on what powers remain after A; small in total	compute-engine	Counted, not prototyped
D	Round the exponent before exp inside non-integer pow	10 to 20% of a power	compute-engine	Read from the operand sizes
E	Bernoulli coefficients as BigDecimal, once per precision, in Gamma and PolyGamma, as the Hurwitz kernel has them	1.2× to 1.5× on those two	compute-engine	Counted, not prototyped
F	Binary splitting for the argument shift in Gamma at a rational: the product x(x+1)…(x+m−1) as one exact integer product and one division	1.5× to 2× on the shift, which is the larger half of Gamma	compute-engine	Not prototyped
G	The sieve for PolyLog and LerchPhi at integer s, as A	At most 1.5×: an integer power is already cheap	compute-engine	Counted
H	Faster BigDecimal primitives	None at these sizes	compute-engine	The table in §2
I	Less boxing and conversion	Under 7%	compute-engine and ours	The standalone prototype
J	Kernels of our own in place of the native ones	Whatever A to F give, sooner	ours	Against the direction: mathematics goes upstream, and a second kernel is a second thing to keep correct

A and B together take Zeta(3/7) at 30 digits from 1.10 ms to about 0.3 ms in the prototype, which is under twice mpmath on the same machine where it is now more than six times.

4. Lean

1. A first. It is measured, it is small, and it has no accuracy question: the same series, the same terms and the same working precision, with each (k+1)^(-s) reached by a handful of rounded multiplications where the guard has 34 digits to spare.
2. B second, as a measured study and not a constant. The guards were set against hard cases, and the rule is to decline before answering wrong. So the work is a grid against mpmath per kernel, and the outcome is a guard each kernel sizes from its own cancellation, as the Hurwitz planner already does. It needs a build box.
3. E and F for Gamma, which has the largest ratio and the least of its time in powers.
4. Not H, I or J.

Each goes to compute-engine as an issue with its pull request, through the session that handles upstream posts. A is identical in result to what ships, so it needs no local patch while it waits.

5. Plan for the first step

For a builder working in a compute-engine checkout.

1. src/compute-engine/numerics/special-functions.ts, zetaCore: when s is not an integer, fill a table of (k+1)^(-s) for the n terms before the loop. A smallest-prime-factor sieve over 1..n; a prime's entry is pow, a composite's is the rounded product of two earlier entries. The loop reads the table. Integer s keeps the ladder it has.
2. Tests, as that repository writes them: the existing Zeta values unchanged; a grid of non-integer s (rationals on both sides of 1, large s, negative s through the reflection) at 30, 50 and 100 digits giving the digits the current kernel gives.
3. The pull request states the before and after for Zeta(3/7) at 30 and 100 digits, and for one large and one negative s.
4. Nothing changes here until the release that carries it. The nightly bench's Zeta cases are the confirmation; they are N() cases, so they measure.

6. The greatest common divisor

This one is ours, and a profile changes what the fix is.

What runs. compute-engine's GCD is Euclid with %. declareFastGcd (packages/symbols/arithmetic/number-theory/src/declare-fast-gcd.ts) replaces it for integers with hybridGcd, which calls lehmerGcd (packages/symbols/arithmetic/residues/src/arith.ts) at 512 bits and above. So the kernel is already Lehmer's algorithm.

Measured, eight pairs of 10000-bit integers, each evaluation on a freshly boxed expression:

	ms
Plain bigint Euclid with %	61
Bare compute-engine	64
Our engine (lehmerGcd)	34
Plain bigint binary GCD	17 to 21
Prototype Lehmer: a 60-bit window, bit length tracked in 64-bit steps	4.1

Our Lehmer is slower than the plain binary GCD written for this comparison. (The comment on lehmerGcd records it as faster than the repository's own binary gcd; that one was not timed here.) Two things in it, read from the code:

• A 26-bit window. Each outer step removes about 24 bits, so a 10000-bit pair takes over 400 steps where a 60-bit window takes about 170.
• The bit length by string. Every outer step calls x.toString(2).length on the full number, which is linear in its size and builds a 10000-character string.

The fix. In lehmerGcd: a window of 60 bits with bigint cofactors (or 53 with doubles), and the bit length carried from step to step and walked down in 64-bit chunks. The threshold is measured again: the prototype is slower than Euclid at 1000 bits and eighteen times faster at 10000. Tests hold it to Euclid on random pairs around the threshold, with and without a large common factor.

Expected. About eight times: 34 ms to 4 or 5 on the laptop, so roughly 63 ms to 8 on the bench's runner. That is about 3.5 times SymPy, 1.5 times Rust's num, and 10 times Julia. Julia's is GMP's subquadratic algorithm in assembly; a half-GCD in JavaScript pays only well past 10^5 bits and is not proposed.

Then upstream. Euclid on bigint is what every compute-engine user gets. Once ours is fixed and measured it is a candidate for compute-engine's own gcd, and the override here is retired when it lands.

7. To measure on a build box

• The accuracy grid for B: each kernel against mpmath at 30, 50 and 100 digits, over the bench's generators and the hard cases (near a pole, large arguments, the reflection), at each candidate guard.
• A CPU profile on an unminified compute-engine, per kernel, to split what this page could only count.
• Gamma: the shift against the series, to size E and F.
• The GCD threshold on the bench's runner.
• The bench's own inputs: some random rationals land on exact paths (an integer or half-integer argument, an integer a), so a few samples time something other than the kernel. How many, per case.