Skip to content

Evaluate the incomplete gamma function's series faster #2079

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Oct 22, 2026
+68 −20
Merged

Evaluate the incomplete gamma function's series faster #2079

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
88 changes: 68 additions & 20 deletions Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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));
}

/// <summary>
/// <c>Gamma(a)</c> at <paramref name="work"/>'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 <c>a</c> for every
/// <c>z</c>.
/// </summary>
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();

/// <summary>The power of ten of a nonzero finite decimal's leading digit: 2 for 345.6, -3 for 0.00123.</summary>
private static int DecimalOrder(EDecimal value)
=> value.Exponent.Add(value.Precision()).Subtract(EInteger.One).ToInt32Checked();
Expand All @@ -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));
}

/// <summary><c>z^a = e^(a ln z)</c>, principal.</summary>
Expand Down
Loading