From ff5cb0617824c868324fb42be7741675c1e8f586 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 22:11:05 +0000 Subject: [PATCH] Evaluate the incomplete gamma function's series faster Gamma(a, z) by its series is Gamma(a) - gamma(a, z), and Gamma(a) was computed by Spouge's approximation on every evaluation: six of the nine milliseconds at a hundred digits. It is remembered now, per order and precision, since a function the integrator samples at points is asked it at one a for every z. A real a and z sum the series in real arithmetic, one product and one quotient a term where the complex arithmetic took nineteen operations, and the logarithm on the cut takes pi from the cache rather than computing it afresh. gamma(2.91, -3.3) 9.52 -> 0.92 ms, gamma(1/3, 2) 23.0 -> 2.2 ms, gamma(1/2 + i, 2 + 3i) 11.8 -> 2.7 ms; the continued fraction and the exponential integral's series are unchanged. Found by #2075's next integration rule, whose answers one declined integral sampled three thousand times. Part of #2075. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Continuous/Number/IncompleteGamma.cs | 88 ++++++++++++++----- 1 file changed, 68 insertions(+), 20 deletions(-) diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs index b870d7dc4..2be14ef1e 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs @@ -1,4 +1,4 @@ -// +// // Copyright (c) 2019-2026 Angouri. // AngouriMath is licensed under MIT. // Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. @@ -129,28 +129,50 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessThe { var epsilon = EDecimal.Create(EInteger.One, EInteger.FromInt32(-work.Precision.ToInt32Checked())); // term_k = z^k / (a (a + 1) ... (a + k)), starting at 1/a - var term = Reciprocal(a, work); - var sum = term; - var peak = Magnitude(term); - var aPlusK = a; - for (var k = 1; ; k++) + (EDecimal Re, EDecimal Im) sum; + EDecimal peak; + if (a.Im.IsZero && z.Im.IsZero) { - aPlusK = (aPlusK.Re.Add(EDecimal.One, work), aPlusK.Im); - term = Multiply(Multiply(term, z, work), Reciprocal(aPlusK, work), work); - sum = (sum.Re.Add(term.Re, work), sum.Im.Add(term.Im, work)); - var size = Magnitude(term); - if (size.CompareTo(peak) > 0) - peak = size; - if (k > 4 && size.CompareTo(Magnitude(sum).Multiply(epsilon, work)) < 0 && size.CompareTo(peak) < 0) - break; + // Real, as a sampled function's arguments are: one product and one quotient a term, where the complex + // arithmetic took nineteen operations, and an evaluation four milliseconds. + var term = EDecimal.One.Divide(a.Re, work); + var real = term; + peak = term.Abs(); + var aPlusK = a.Re; + for (var k = 1; ; k++) + { + aPlusK = aPlusK.Add(EDecimal.One, work); + term = term.Multiply(z.Re, work).Divide(aPlusK, work); + real = real.Add(term, work); + var size = term.Abs(); + if (size.CompareTo(peak) > 0) + peak = size; + if (k > 4 && size.CompareTo(real.Abs().Multiply(epsilon, work)) < 0 && size.CompareTo(peak) < 0) + break; + } + sum = (real, EDecimal.Zero); + } + else + { + var term = Reciprocal(a, work); + sum = term; + peak = Magnitude(term); + var aPlusK = a; + for (var k = 1; ; k++) + { + aPlusK = (aPlusK.Re.Add(EDecimal.One, work), aPlusK.Im); + term = Multiply(Multiply(term, z, work), Reciprocal(aPlusK, work), work); + sum = (sum.Re.Add(term.Re, work), sum.Im.Add(term.Im, work)); + var size = Magnitude(term); + if (size.CompareTo(peak) > 0) + peak = size; + if (k > 4 && size.CompareTo(Magnitude(sum).Multiply(epsilon, work)) < 0 && size.CompareTo(peak) < 0) + break; + } } var prefactor = Multiply(ComplexPower(z, a, work), ComplexExp((z.Re.Negate(), z.Im.Negate()), work), work); var lower = Multiply(prefactor, sum, work); - Complex gammaOfA; - using (MathS.Settings.DecimalPrecisionContext.Set(work)) - using (MathS.Settings.DowncastingEnabled.Set(false)) - gammaOfA = Factorial(Complex.Create(a.Re.Subtract(EDecimal.One, work), a.Im) as Complex); - var gamma = (Re: gammaOfA.RealPart.EDecimal, Im: gammaOfA.ImaginaryPart.EDecimal); + var gamma = GammaOf(a, work); var scale = EDecimal.Max(Magnitude(gamma), Magnitude(prefactor).Multiply(peak, work)); return ((gamma.Re.Subtract(lower.Re, work), gamma.Im.Subtract(lower.Im, work)), scale); } @@ -210,6 +232,29 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AtANonPositive return ((difference.Re.Multiply(coefficient, work), difference.Im.Multiply(coefficient, work)), scale.Multiply(coefficient.Abs(), work)); } + /// + /// Gamma(a) at 's digits, remembered: it is most of the cost of an evaluation, six + /// milliseconds at a hundred digits, and a function sampled at points is asked it at one a for every + /// z. + /// + private static (EDecimal Re, EDecimal Im) GammaOf((EDecimal Re, EDecimal Im) a, EContext work) + { + var key = (a.Re, a.Im, work.Precision.ToInt32Checked()); + if (gammaOfTheOrder.TryGetValue(key, out var known)) + return known; + Complex value; + using (MathS.Settings.DecimalPrecisionContext.Set(work)) + using (MathS.Settings.DowncastingEnabled.Set(false)) + value = Factorial(Complex.Create(a.Re.Subtract(EDecimal.One, work), a.Im) as Complex); + var gamma = (value.RealPart.EDecimal, value.ImaginaryPart.EDecimal); + // Emptied rather than grown without end, as the integrator's answers are. + if (gammaOfTheOrder.Count >= 256) + gammaOfTheOrder.Clear(); + gammaOfTheOrder[key] = gamma; + return gamma; + } + [ConcurrentField] private static readonly System.Collections.Concurrent.ConcurrentDictionary<(EDecimal, EDecimal, int), (EDecimal, EDecimal)> gammaOfTheOrder = new(); + /// The power of ten of a nonzero finite decimal's leading digit: 2 for 345.6, -3 for 0.00123. private static int DecimalOrder(EDecimal value) => value.Exponent.Add(value.Precision()).Subtract(EInteger.One).ToInt32Checked(); @@ -227,8 +272,11 @@ private static (EDecimal Re, EDecimal Im) ComplexLn((EDecimal Re, EDecimal Im) z { if (z.Im.IsZero && z.Re.Sign > 0) return (z.Re.Log(work), EDecimal.Zero); + // On the cut, ln |z| + i pi, with pi from the cache: computed afresh it was a third of an evaluation. + if (z.Im.IsZero) + return (z.Re.Negate().Log(work), InternalAMExtensions.ConstantCache.Lookup(work).Pi); var modulus = z.Re.Multiply(z.Re, work).Add(z.Im.Multiply(z.Im, work), work).Sqrt(work); - return (modulus.Log(work), z.Im.IsZero ? EDecimal.PI(work) : z.Im.Arctan2(z.Re, work)); + return (modulus.Log(work), z.Im.Arctan2(z.Re, work)); } /// z^a = e^(a ln z), principal.