From 2c25febac61859edf5d2dfbd9664ee71dcb5da5e Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 15 Sep 2026 00:35:42 +0000 Subject: [PATCH] The sine, the cosine and the arctangent are series in the same fixed point as the logarithm and the exponential #1342 made the trigonometric functions argument reduction and short series, in EDecimal; #1346 found that a series in EDecimal pays two to three microseconds a term for the alignment in each addition and moved the logarithm and the exponential to fixed point. The sine, the cosine and the arctangent are there too now: the reduction to [-pi, pi] is an integer division by 2 pi in fixed point, the folds and halvings are subtractions and shifts, the square root in the arctangent's halving x/(1 + sqrt(1 + x^2)) is the integer square root of the value shifted up, exactly floored, and the halvings come back as a shift. Pi, 2 pi and pi/2 join the constant cache in fixed point. The arcsine and arccosine are the arctangent's, and the tangent the quotient, as before. The hundred-digit values are the ones pinned before, to the last digit. EvalTrig 552 us / 1,192,089 B -> 153 us / 260,728 B EvalTrigPrecise 13.8 ms / 7,935,226 B -> 2.16 ms / 2,525,707 B EvalTranscendentalFresh 245 us / 423,032 B -> 83 us / 128,664 B SolveEasy 2.88 ms / 4,425,194 B -> 1.19 ms / 1,675,487 B every other entry within 1.5% Alone, at a hundred digits: sin(x) is 45 us where it was 176, arctan(x) 46 where it was 202, arcsin(x) 59 where it was 238; at five hundred digits arctan(x) is 1.2 ms where it was 9.3. Against #1341, before any of this, EvalTrig is 709 -> 153 us and EvalTrigPrecise 22.5 -> 2.2 ms. The gate's baseline is from this run; the performance log has the entry. Four suites green (11,252 / 134 / 18 / 41). Part of #1338. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../WhatsNew/version_performance_control.md | 26 ++++ .../Functions/InternalAMExtensions.cs | 129 +++++++++--------- .../DotnetBenchmark/performance-baseline.json | 64 ++++----- 3 files changed, 122 insertions(+), 97 deletions(-) diff --git a/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md b/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md index 834cbfa55..e1b417ce4 100644 --- a/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md +++ b/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md @@ -256,6 +256,32 @@ cent — its measured run-to-run spread on mean time is up to 51.8%. A move of 8 outside that band by a wide margin and agrees in sign and rough size with the allocation column beside it. The small rows from the same run are still not worth reading, and are not quoted. +## The 2055th: the sine, cosine and arctangent in the same fixed point + +The 2050th made the trigonometric functions argument reduction and short series, in +`EDecimal`; the 2054th found that a series in `EDecimal` pays two to three microseconds a term +for the alignment in each addition and moved the logarithm and the exponential to fixed point. +The sine, the cosine and the arctangent are there too now: the reduction to `[-pi, pi]` is an +integer division by `2 pi` in fixed point, the folds and halvings are subtractions and shifts, +the square root in the arctangent's halving `x/(1 + sqrt(1 + x^2))` is the integer square root +of the value shifted up, exactly floored, and the halvings come back as a shift. The arcsine +and arccosine are the arctangent's, and the tangent the quotient. The hundred-digit values +are the ones pinned before, to the last digit. + +| benchmark | 2054th | 2055th | allocation | time | +|---|--:|--:|--:|--:| +| `EvalTrig` | 1,192,089 | **260,728** | **−78.1%** | 552 → **153 µs** | +| `EvalTrigPrecise` | 7,935,226 | **2,525,707** | **−68.2%** | 13.8 → **2.16 ms** | +| `EvalTranscendentalFresh` | 423,032 | **128,664** | **−69.6%** | 245 → **83 µs** | +| `SolveEasy` | 4,425,194 | **1,675,487** | **−62.1%** | 2.88 → **1.19 ms** | +| every other entry | | | within 1.5% | | + +Bytes allocated per call, same machine, both columns by the gate in one session; the gate's +baseline was taken from this run. Alone, at a hundred digits: `sin(x)` is 45 µs where it was +176, `arctan(x)` 46 where it was 202, `arcsin(x)` 59 where it was 238; at five hundred digits +`arctan(x)` is 1.2 ms where it was 9.3. Against the 2049th, before any of this, `EvalTrig` is +709 → 153 µs and `EvalTrigPrecise` 22.5 → 2.2 ms. + ## The 2054th: the logarithm and the exponential as series in fixed point `EvalTranscendentalFresh` -- `ln(1 + x^2) arctan(x)` at a fresh point -- was 1.3 ms after the diff --git a/Sources/AngouriMath/Functions/InternalAMExtensions.cs b/Sources/AngouriMath/Functions/InternalAMExtensions.cs index b294118fe..a2716adbd 100644 --- a/Sources/AngouriMath/Functions/InternalAMExtensions.cs +++ b/Sources/AngouriMath/Functions/InternalAMExtensions.cs @@ -205,6 +205,9 @@ public static ConstantCache Lookup(EContext context) FixedOne = EInteger.One.ShiftLeft(FixedBits); FiveToFixedBits = EInteger.FromInt32(5).Pow(FixedBits); FixedSqrt2 = ToFixed(EDecimal.FromString("1.41421356237309504880168872420969807856967187537694"), FixedBits); + FixedPi = ToFixed(Pi, FixedBits); + FixedTwoPi = FixedPi.ShiftLeft(1); + FixedHalfPi = FixedPi.ShiftRight(1); // ln 2 = 2 artanh(1/3), and ln 10 = ln 8 + ln 5/4 = 3 ln 2 + 2 artanh(1/9). FixedLn2 = TwiceArtanh(FixedOne.Divide(3), FixedBits); FixedLn10 = FixedLn2.Multiply(3).Add(TwiceArtanh(FixedOne.Divide(9), FixedBits)); @@ -219,6 +222,12 @@ public static ConstantCache Lookup(EContext context) public EInteger FiveToFixedBits { get; } /// The square root of 2, in fixed point, to fifty digits public EInteger FixedSqrt2 { get; } + /// Pi in fixed point, to the context's digits + public EInteger FixedPi { get; } + /// 2 pi in fixed point + public EInteger FixedTwoPi { get; } + /// Pi / 2 in fixed point + public EInteger FixedHalfPi { get; } /// The natural logarithm of 2, in fixed point public EInteger FixedLn2 { get; } /// The natural logarithm of 10, in fixed point @@ -406,69 +415,68 @@ internal static (EDecimal Sin, EDecimal Cos) SinAndCos(EDecimal x, EContext cont { var working = WithGuardDigits(context, 8); var consts = ConstantCache.Lookup(working); + var bits = consts.FixedBits; + var one = consts.FixedOne; - //truncating to [-2*PI;2*PI] - TruncateToPeriodicInterval(ref x, consts, working); - // now x in (-2pi, 2pi); into [-pi, pi], with the signs the shift costs + // In fixed point from here: the reduction is integer arithmetic on pi's digits, + // and a halving is a shift. + var fixedX = ToFixed(x, bits); + // Into [-pi, pi]: the nearest multiple of 2 pi off. + var shifted = fixedX.Add(consts.FixedPi); + var turns = shifted.Divide(consts.FixedTwoPi); + if (shifted.Sign < 0 && !shifted.Remainder(consts.FixedTwoPi).IsZero) + turns = turns.Subtract(1); // the division truncates, and this floors + if (!turns.IsZero) + fixedX = fixedX.Subtract(consts.FixedTwoPi.Multiply(turns)); + // Into [-pi/2, pi/2], with the sign the fold costs the cosine. var negateCos = false; - if (x.GreaterThan(consts.Pi)) - { - x = x.Subtract(consts.TwoPi, working); - } - else if (x.LessThan(-consts.Pi)) - { - x = x.Add(consts.TwoPi, working); - } - if (x.GreaterThan(consts.HalfPi)) + if (fixedX.CompareTo(consts.FixedHalfPi) > 0) { // sin(x) = sin(pi - x), cos(x) = -cos(pi - x) - x = consts.Pi.Subtract(x, working); + fixedX = consts.FixedPi.Subtract(fixedX); negateCos = true; } - else if (x.LessThan(-consts.HalfPi)) + else if (fixedX.CompareTo(consts.FixedHalfPi.Negate()) < 0) { // sin(x) = sin(-pi - x), cos(x) = -cos(-pi - x) - x = consts.Pi.Negate().Subtract(x, working); + fixedX = consts.FixedPi.Negate().Subtract(fixedX); negateCos = true; } // Halve until |x| < 1/20. var halvings = 0; - var twentieth = EDecimal.Create(5, -2); - while (x.Abs().GreaterThan(twentieth)) + var twentieth = one.Divide(20); + while (fixedX.Abs().CompareTo(twentieth) > 0) { - x = x.Multiply(consts.Half, working); + fixedX = fixedX.Sign < 0 ? fixedX.Negate().ShiftRight(1).Negate() : fixedX.ShiftRight(1); halvings++; } - var xx = x.Multiply(x, working); + var square = MultiplyFixed(fixedX, fixedX, bits); // sin: x - x^3/3! + ...; cos: 1 - x^2/2! + ... - var sinTerm = x; - var sin = x; - var cosTerm = EDecimal.One; - var cos = EDecimal.One; - for (var i = 1; i < 200; i++) + var sinTerm = fixedX; + var sin = fixedX; + var cosTerm = one; + var cos = one; + for (var i = 1; i < 100000; i++) { - cosTerm = -cosTerm.Multiply(xx, working).Divide((2 * i - 1) * (2 * i), working); - sinTerm = -sinTerm.Multiply(xx, working).Divide((2 * i) * (2 * i + 1), working); - var nextCos = cos.Add(cosTerm, working); - var nextSin = sin.Add(sinTerm, working); - var settled = nextCos.Equals(cos) && nextSin.Equals(sin); - cos = nextCos; - sin = nextSin; - if (settled) + cosTerm = MultiplyFixed(cosTerm, square, bits).Negate().Divide((2 * i - 1) * (2 * i)); + sinTerm = MultiplyFixed(sinTerm, square, bits).Negate().Divide((2 * i) * (2 * i + 1)); + if (cosTerm.IsZero && sinTerm.IsZero) break; + cos = cos.Add(cosTerm); + sin = sin.Add(sinTerm); } for (var i = 0; i < halvings; i++) { // sin(2y) = 2 sin(y) cos(y), cos(2y) = 2 cos(y)^2 - 1 - var doubledSin = sin.Multiply(cos, working).Multiply(2, working); - cos = cos.Multiply(cos, working).Multiply(2, working).Decrement(); + var doubledSin = MultiplyFixed(sin, cos, bits).ShiftLeft(1); + cos = MultiplyFixed(cos, cos, bits).ShiftLeft(1).Subtract(one); sin = doubledSin; } if (negateCos) - cos = -cos; - return (sin.RoundToPrecision(context), cos.RoundToPrecision(context)); + cos = cos.Negate(); + return (FromFixed(sin, consts, working).RoundToPrecision(context), FromFixed(cos, consts, working).RoundToPrecision(context)); } /// Analogy of @@ -533,21 +541,6 @@ public static EDecimal Arccosecant(this EDecimal x, EContext context) } /// Truncates to [-2*, 2*] - private static void TruncateToPeriodicInterval(ref EDecimal x, ConstantCache consts, EContext context) - { - while (x.GreaterThanOrEquals(consts.TwoPi)) - { - var divide = x.Divide(consts.TwoPi, context).Floor().Abs(); - x = divide.MultiplyAndAdd(-consts.TwoPi, x, context); - } - - while (x.LessThanOrEquals(-consts.TwoPi)) - { - var divide = x.Divide(consts.TwoPi, context).Floor().Abs(); - x = divide.MultiplyAndAdd(consts.TwoPi, x, context); - } - } - /// Analogy of /// /// arcsin(x) = arctan(x/sqrt(1 - x^2)), with 1 - x^2 as (1 - x)(1 + x) @@ -598,27 +591,33 @@ public static EDecimal Arctan(this EDecimal x, EContext context) if (x.GreaterThan(EDecimal.One)) return ConstantCache.Lookup(working).HalfPi.Subtract(Arctan(EDecimal.One.Divide(x, working), working), working).RoundToPrecision(context); + // In fixed point from here: the square root of a fixed-point number is the + // integer square root of it shifted up, exactly floored, and the halvings come + // back as a shift. + var workingConsts = ConstantCache.Lookup(working); + var bits = workingConsts.FixedBits; + var one = workingConsts.FixedOne; + var fixedX = ToFixed(x, bits); var halvings = 0; - var twentieth = EDecimal.Create(5, -2); - while (x.GreaterThan(twentieth)) + var twentieth = one.Divide(20); + while (fixedX.CompareTo(twentieth) > 0) { - x = x.Divide(EDecimal.One.Add(x.MultiplyAndAdd(x, EDecimal.One, working).Sqrt(working), working), working); + var root = one.Add(MultiplyFixed(fixedX, fixedX, bits)).ShiftLeft(bits).Sqrt(); + fixedX = fixedX.ShiftLeft(bits).Divide(one.Add(root)); halvings++; } - var xx = x.Multiply(x, working); - var power = x; - var sum = x; - for (var i = 1; i < 400; i++) + var square = MultiplyFixed(fixedX, fixedX, bits); + var power = fixedX; + var sum = fixedX; + for (var i = 1; i < 100000; i++) { - power = -power.Multiply(xx, working); - var next = sum.Add(power.Divide(2 * i + 1, working), working); - if (next.Equals(sum)) + power = MultiplyFixed(power, square, bits).Negate(); + var term = power.Divide(2 * i + 1); + if (term.IsZero) break; - sum = next; + sum = sum.Add(term); } - if (halvings > 0) - sum = sum.Multiply(EDecimal.FromInt32(1 << halvings), working); - return sum.RoundToPrecision(context); + return FromFixed(sum.ShiftLeft(halvings), workingConsts, working).RoundToPrecision(context); } /// Analogy of public static EDecimal Acos(this EDecimal x, EContext context) diff --git a/Sources/Tests/DotnetBenchmark/performance-baseline.json b/Sources/Tests/DotnetBenchmark/performance-baseline.json index f363dda35..978b86796 100644 --- a/Sources/Tests/DotnetBenchmark/performance-baseline.json +++ b/Sources/Tests/DotnetBenchmark/performance-baseline.json @@ -1,94 +1,94 @@ { "comment": "Allocated bytes and mean nanoseconds per operation for the popular use cases of https://github.com/asc-community/AngouriMath/issues/746. Checked by PerformanceGate.cs; see Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md for when updating it is legitimate.", "benchmark": "CommonFunctionsInterVersion", - "commit": "bfb60e9cdcdc855c13b5ac58001db0ec98e3f3b2", + "commit": "95082e2f814874b99f152fbb1072fc977735337f", "measuredOn": "2026-09-15", "runtime": ".NET 10.0.10", "machine": "Ubuntu 26.04 LTS, X64, 8 logical cores", "cases": { "CompileEasy": { - "allocatedBytes": 11023, - "meanNanoseconds": 184521.07338867188 + "allocatedBytes": 11031, + "meanNanoseconds": 184131.37203543526 }, "CompileHard": { - "allocatedBytes": 20737, - "meanNanoseconds": 309526.83076171874 + "allocatedBytes": 20434, + "meanNanoseconds": 310381.04733072914 }, "Derivate": { - "allocatedBytes": 53099, - "meanNanoseconds": 14289.474793570382 + "allocatedBytes": 53111, + "meanNanoseconds": 14270.049594116212 }, "EvalEasy": { "allocatedBytes": 0, - "meanNanoseconds": 1.8900332285889558 + "meanNanoseconds": 1.883551713079214 }, "EvalPolynomialFresh": { "allocatedBytes": 12328, - "meanNanoseconds": 4151.610092798869 + "meanNanoseconds": 4164.772733835073 }, "EvalPolynomialFresh15Digits": { "allocatedBytes": 11656, - "meanNanoseconds": 4141.208742959158 + "meanNanoseconds": 4157.980944497244 }, "EvalTranscendentalFresh": { - "allocatedBytes": 423032, - "meanNanoseconds": 245377.50509207588 + "allocatedBytes": 128664, + "meanNanoseconds": 83463.06102643695 }, "EvalTrig": { - "allocatedBytes": 1192089, - "meanNanoseconds": 551794.5797293527 + "allocatedBytes": 260728, + "meanNanoseconds": 153426.1548339844 }, "EvalTrigPrecise": { - "allocatedBytes": 7935226, - "meanNanoseconds": 13781485.40625 + "allocatedBytes": 2525707, + "meanNanoseconds": 2164306.404296875 }, "ParseEasy": { - "allocatedBytes": 18125, - "meanNanoseconds": 5981.867161090558 + "allocatedBytes": 18104, + "meanNanoseconds": 5932.770784231333 }, "ParseHard": { "allocatedBytes": 3555057, - "meanNanoseconds": 1367037.3726283482 + "meanNanoseconds": 1392708.7341145833 }, "RunEasy": { "allocatedBytes": 0, - "meanNanoseconds": 20.056596242464504 + "meanNanoseconds": 20.095689982175827 }, "RunHard": { "allocatedBytes": 0, - "meanNanoseconds": 295.51566811970304 + "meanNanoseconds": 295.40826177597046 }, "RunMedium": { "allocatedBytes": 0, - "meanNanoseconds": 163.17946696281433 + "meanNanoseconds": 162.8344373519604 }, "SimplifyEasy": { "allocatedBytes": 78578, - "meanNanoseconds": 45092.335445149736 + "meanNanoseconds": 44253.979728190105 }, "SimplifyHard": { - "allocatedBytes": 189726120, - "meanNanoseconds": 109370614.94736843 + "allocatedBytes": 189604392, + "meanNanoseconds": 111491250.2413793 }, "SolveEasy": { - "allocatedBytes": 4425194, - "meanNanoseconds": 2884867.8782552085 + "allocatedBytes": 1675487, + "meanNanoseconds": 1185599.8670479911 }, "SolveEasyMedium": { "allocatedBytes": 65184, - "meanNanoseconds": 18597.25531768799 + "meanNanoseconds": 18840.902629307337 }, "SolveHard": { "allocatedBytes": 11706656, - "meanNanoseconds": 64677752.72 + "meanNanoseconds": 63240162.666666664 }, "SolveMedium": { - "allocatedBytes": 446445, - "meanNanoseconds": 266309.50540597097 + "allocatedBytes": 446414, + "meanNanoseconds": 259625.15890066963 }, "SolveMediumHard": { "allocatedBytes": 1367798, - "meanNanoseconds": 8561109.407291668 + "meanNanoseconds": 8776522.877083333 } }, "ungated": {