Repository navigation
Optimized numerical evaluation #1338
Description
Activity
Some measurements to anchor the first two questions, taken on master today (Release, net10.0, one core, 2,000 evaluations each;
expr.Substitute(x, v).EvalNumerical()with a fresh value each time, andexpr.Compile<double, double>(x)for the double column):shape 100 digits (default) 15 digits double, compiled 100/15 100/double x^3 + 2x^2 - x + 1500 µs 117 µs 0.073 µs 4.3× ~7,000× sin(x)^2 + cos(x)^2752 µs 114 µs 0.032 µs 6.6× ~24,000× e^(-x^2/2)/sqrt(2 pi)1,062 µs 428 µs 0.017 µs 2.5× ~63,000× ln(1 + x^2) arctan(x)3,547 µs 782 µs 0.019 µs 4.5× ~190,000× (x + 1)/(x^2 + 1) + sqrt(x^2 + 3)624 µs 164 µs 0.020 µs 3.8× ~32,000× Two things the table says:
- Precision is not where most of the time is. Going from 100 digits to 15 buys 2.5–6.6×, not the ~40× the digit count would suggest, so the bulk of an
EvalNumericalis per-node overhead — theEContextarithmetic, the series for the transcendental functions, and allocation — rather than the width of the numbers.Substitutealone is 28 µs and a secondEvalNumericalof the same entity is 1.6 µs (it is memoised on the node), so the cost is the evaluation itself. - A double path is four to five orders of magnitude faster, and the compiled delegate is already that path. So for the "would double suffice at 10 digits" question the answer is numerically yes for the arithmetic — a double carries 15–16 significant digits — but not unconditionally:
x - ywithxandyagreeing to 12 digits leaves 3, and the tree does not know it lost them. That is the case the current precision guards against, and it is the one an auto-tuned default would have to detect (an interval or double-double evaluation alongside the double one gives a bound cheaply; a plain double does not).
On the PeterO question: I do not know of a managed arbitrary-precision decimal library that is both maintained and faster; .NET itself has
BigIntegerand nothing above it, and the fast libraries (MPFR, Arb) are native. Since the measured cost is overhead rather than digits, the cheaper win is probably inside the evaluator — the per-node overhead and the series behindsin/cosinInternalAMExtensions, which take the context but were not written for speed — before the number type is changed at all. And for a caller who asks for 15 digits or fewer, evaluating the tree indoublethe way the compiler already does, and only falling back toEDecimalwhere a cancellation is detected, would give the table's last column toEvalNumericalfor most inputs.Happy to turn any of these rows into a tracked benchmark if that is the shape you want the frontier measured in.
- Precision is not where most of the time is. Going from 100 digits to 15 buys 2.5–6.6×, not the ~40× the digit count would suggest, so the bulk of an
Sure! It looks like you have identified specific opportunities for optimization. Proceed with that; see how fast arbitrary precision decimal evaluation can be against native doubles.
Will do. I'll start with a profile of where an
EvalNumericalspends its microseconds at 100 digits and at 15, then the evaluator's own overhead, then a double fast path with a cancellation check, each as its own PR against a tracked benchmark with these five rows. (If you want this tracked as an accepted goal, anAcceptedorAgentic goallabel on the issue is the signal I read; I will proceed either way.)- added 4 commits that reference this issue
on Sep 14, 2026 #1342 is merged. What it does, measured against the three rows of #1341 and against the primitives alone (bytes per call, both columns by the gate in one session on one machine):
benchmark before after allocation time EvalPolynomialFresh(100 digits)650,888 14,000 −97.8% 349 → 9.4 µs EvalPolynomialFresh15Digits265,816 12,280 −95.4% 83.7 → 9.0 µs EvalTranscendentalFresh(ln(1+x²)·arctan(x))5,668,395 2,247,769 −60.3% 3.25 → 1.32 ms EvalTrig1,341,457 1,190,217 −11.3% 709 → 555 µs EvalTrigPrecise12,742,285 7,933,690 −37.7% 22.5 → 13.8 ms SolveEasy8,865,872 6,280,444 −29.2% 4.83 → 3.54 ms SimplifyHard331,683,360 314,496,408 −5.2% 190 → 178 ms Where the time went, and what changed:
- The rational search (
Real.Create→Rational.FindRational, a 15-level continued fraction in 100-digit arithmetic on every arithmetic result, ~20 µs) was most of a node. It is now decided first in a double read off the leading mantissa bits (~0.3 µs), then inSystem.Decimalwhere the double's error has grown past deciding, and only a value neither rules out reaches the exact search. Nothing that was a rational stops being one; the tests pin both directions and the reading was checked against the exact search on 12,000 rationals and 12,000 near-rationals. - Powers: the binary power computed its half twice; a whole power of an inexact decimal is one
EDecimal.Pow, a half power isSqrt(Newton, ~70× faster than exp(log/2)). - Sin/cos: argument reduction to [−π/2, π/2], halving below 1/20, two short series, doubling back — instead of sixty Taylor terms and
sin = sqrt(1 − cos²), which lost half the digits near zero. - Arctan/arcsin: the arctangent had gone through the arcsine series (2 ms at x = 1/3); it is its own series after reciprocal and halving, and arcsin is arctan of x/√((1−x)(1+x)).
Against native doubles, for the question you asked: the 15-digit polynomial evaluation is 9 µs per call for a cubic, against roughly 10 ns for the same in doubles — three orders, of which the numbers themselves are now a small part. What remains is the evaluator's own overhead:
EvalNumericalisInnerSimplifyWithCheckper node, with its domain check and the lazy caches, ~2–3 µs a node; and PeterO'sLog/Expat ~300 µs a call, which is why the transcendental row is still 1.3 ms. Those two are the next steps: a dedicated numeric walker, and a double fast path for contexts of ≤ 15 digits with a cancellation check that falls back to the decimal when the double loses digits. I'll report each with the same rows.The last digits of sin, cos, tan, arcsin and arctan changed in the last 1–3 of 100 places (twenty pinned values in
NumericDigits); checked against mpmath at 130 digits, sin/cos/tan(1) are exact now where they were 1–3 ulps off, the complex trig functions went from ~20 ulps to ≤ 2. Recorded in BREAKING-CHANGES.md and in the performance log.- The rational search (
12 remaining items
- added a commit that references this issue
on Sep 16, 2026 Step 2 is in (#1366), and the pre-check fix before it (#1365). The gate rows, same machine, one session each, pinned digits unmoved throughout:
gate row #1364 #1365 #1366 against #1341 EvalTrig83 µs 37 32 µs 709 → 32 EvalTrigPrecise(500 digits)1.05 ms 0.62 0.41 ms 22.5 → 0.41 EvalTranscendentalFresh33 µs 33 16 µs 245 → 16 SolveEasy0.49 ms 0.49 0.39 ms 2.88 → 0.39 The functions alone at 100 digits, against mpmath on the same machine:
#1347 now mpmath ln52.4 µs 8.4 4.5 arctan48.8 7.7 3.9 exp32.3 10.0 4.3 sin46.9 13.1 5.6 At 500 digits
ln3,921 → 241 µs,sin1,273 → 141, andsqrt41 → 22 against PeterO's (441 → 174 at 2,000). What did it: the logarithm's mantissa divided by the nearest1 + j/64with those logarithms cached per context (80 → 25 terms), the arctangent by the nearestj/64likewise (47 → 30 terms, no square roots), the cosine assqrt(1 − sin²)at the reduced argument, and a precision-doubling integer square root with a sticky digit so the decimal root is still correctly rounded in every mode. The remaining factor of two is constant overhead around the series — theEDecimalat the boundary,Real.Create, the entity — which is the tower of step 2 in the plan above rather than more arithmetic.One more thing found on the way, filed as #1367 and fixed in #1368:
pi,eand any held entity's cachedEvaledkept the digits of the first precision context they were evaluated under, sopiat 300 digits after a 100-digit evaluation was 100 digits ande^xat 300 was correct to 100. The constant is now looked up at evaluation, and the precision setting advances an epoch that every entity checks against its cache — oneintper entity, +0.1–1.0% on the gate's allocation rows, no baseline change.Yes, queue this for v3 (do you have a queue for it yet?)
Queued. The docket is #1019 (Pending breakages, the one you asked for on #1009), and this is its item 8 now: the 48 published members that take or return PeterO's types, listed, with the three steps — the series on
BigInteger(done), the internal tower behind the current properties (next, measurable on the whole gate with no signature change), and the type change itself at v3 with the migration inBREAKING-CHANGES.md. Nothing before v3 touches a signature.Reacted by Hadrian Tang- added a commit that references this issue
on Sep 16, 2026 - added this to the milestone
on Sep 18, 2026 - removed this from the milestone
on Sep 22, 2026 Per #1486 (comment) some of the new architecture may want to be piloted before v3 to support optimizing internal evaluations like for the integrator.
Evaledwould then be a deprecated property to be removed on v3.Agreed. I've replied on #1486 with a scope: an internal interval evaluator that takes its accuracy as an argument and reads no global settings, with the integrator's sampled checks and screens as its first clients. This issue can hold the evaluator itself as it grows.
A first measurement of the pilot. The double-interval tier is in
AngouriMath.Numericsand reads no settings. It covers arithmetic, powers, exp and ln, the trigonometric functions and their inverses,abs,sgn,providedandpiecewise, with conditions decided three ways.I audited it against the old evaluation inside
HoldsAtSampledPoints, the integrator's differentiate-back check. At every call it computed both verdicts and returned the old one, so no answer moved. That covered the 1774-problem suite, Rubi 4.7.2 + 4.3.10 + 4.4.10, and the family 1 and 4 samples:checks same verdict intervals alone short disagree suite 659 659 0 0 trig files 264 249 15 (too wide) 0 family 1 651 628 20 too wide, 3 undefined 0 family 4 1205 1204 1 undefined 0 all 2779 2740 (98.6%) 39 0 Doubles never disagree with a hundred digits, and in 98.6% of checks they reach the same verdict. The 39 checks they can't settle are cancellation, plus four points where a double lands on a cut. Those are what the precision escalation is for, so that's next: the same evaluation at more digits, with the precision passed per call. Once that settles the 39,
HoldsAtSampledPointsstops callingEvalNumericalaltogether.The escalation is in: the same evaluation at 40 and then 80 digits, with the precision passed as an argument. On the same four sets:
checks same verdict intervals short intervals hold, old evaluation fails suite 659 659 0 0 trig files 266 263 0 3 family 1 651 649 0 2 family 4 1203 1202 0 1 Every point is now settled: 100% on the suite, the trig files and family 1. There's no check where the old evaluation holds and intervals alone wouldn't.
I logged the six disagreements, and all six are the old evaluation failing. In three of them it returns
NaNfor the Hermite ansatz'srationalPart' + logPart, a large rational expression, at every point. In the other two it returns a wrong value: 689.94 where the value is 690.168…, and exactly zero where it is −0.02575…. The intervals had both sides equal, and exact rational arithmetic in Python confirms it at every point, with the difference exactly 0. So each of the six was a correct antiderivative rejected by the old check. I'll make a small reproducer for that evaluation defect and file it.Next:
HoldsAtSampledPointsmoves onto intervals completely, with noEvalNumericaland noDowncastingEnabled. It will be measured over the families, the suite and the gate, since accepting those six answers may change a few rows.Two corrections to my last comment, then where this stands.
- I wrote that I had logged all six disagreements. I had logged five. The sixth, in family 4, came from a run made before the logging went in, and I described it without having looked at it. I've now rerun that family with logging on and checked it. It's the old evaluation failing too: it evaluated the answer's derivative to 10.55 at the first point, where 80-digit mpmath gives 3.0216 for both sides, which is what the intervals said.
- The old evaluation's failures aren't an arithmetic bug. They come from the limits of its decimals.
DecimalPrecisionContextbounds exponents to [-100, 1000], so with downcasting off a decimal below about 1e-199 becomes zero and one above 1e1000 becomes infinite, and a divisor below 1e-16 counts as zero. In the five I bisected, a factor near 1e-316 times one near 1e+316 came out 0, and a factor flushed to zero times one saturated to infinity came outNaN. Filed as With downcasting off, EvalNumerical loses small divisors and both ends of the exponent range #1498, with a one-line reproducer for each limit.
#1497 moves
HoldsAtSampledPointsandIsZeroAtPinnedSymbolsonto the intervals. The comparison makes noEvalNumericalcall and reads noDowncastingEnabled. Measuring every family found two things the four-set audit hadn't:- Conditions
… in RRinside the pinned answer's derivative. The evaluator didn't decide them, and Rubi 7.1.4/7.1.5 lost five problems until it did. - The decimal tier's cost. On
1/(a + b csch(c + d x)^2)^3its checks took 3.9 s where the old evaluation took 0.33 s. That lasted until the contexts were cached per precision and the work was shared across a check's points.
With both fixed, against master on the same machine, interleaved:
problems solved verdicts changed answers changed time allocation Rubi 4.7.2, 4.3.10, 4.4.10 312 269 → 269 0 0 −18.2% −32.2% the 1774-problem suite 1774 1719 → 1719 0 0 −0.2% −1.3% Rubi 7.1.4 – 7.2.5 376 345 → 345 0 0 −2.9% −1.8% families 1, 4, 5, 6, 7 (samples) 1594 1356 → 1355 2 1 −2.3% −0.3% The changed answer and both verdicts are explained in the PR: a correct answer that master's check turned away, and two timeouts at the 5 s budget. The suite passes, and so does the allocation gate on all 19 gated benchmarks.
- added a commit that references this issue
on Sep 29, 2026
The default precision today enables precise computation without necessarily merging terms just because they differ in floating point computation.
However, for some applications we may opt for less precision, because if we can use hardware floating point types for example, then the speedups will be great.
The same applies for some workloads on some computers that may benefit from specialized instructions or GPU delegation.
This is separate from function compilation because we are optimizing on node manipulation here.
So here are the following items to be determined:
Use this issue to track improvements; only close when all visible improvements are done.