Skip to content

Sum the incomplete gamma function's lower series as DLMF 8.7.3 in the left half-plane #2088

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 23, 2026
Merged

Sum the incomplete gamma function's lower series as DLMF 8.7.3 in the left half-plane #2088

Show file tree
Hide file tree
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
93 changes: 79 additions & 14 deletions Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,11 @@ partial record Number
/// algorithm. Near the negative real axis at a moderate <c>|z|</c> it stops short of the value.
/// </para>
/// <para>
/// Short of that it is <c>Gamma(a) - gamma(a, z)</c>, the lower one summed as
/// <c>z^a e^(-z) sum_k z^k / (a (a + 1) ... (a + k))</c> (DLMF 8.7.1). The terms are about <c>e^|z|</c> and the
/// value can be as small as <c>e^(-|z|)</c>, so twice <c>|z|/ln 10</c> digits are carried as guard digits,
/// and the sum is taken again with more where the difference shows it cancelled more.
/// Short of that it is <c>Gamma(a) - gamma(a, z)</c>, the lower one summed in the right half-plane as
/// <c>z^a e^(-z) sum_k z^k / (a (a + 1) ... (a + k))</c> (DLMF 8.7.1), whose terms there are about <c>e^|z|</c>
/// and of one sign on the positive axis, with twice <c>|z|/ln 10</c> guard digits; and in the left half-plane as
/// <c>z^a sum_k (-z)^k / (k! (a + k))</c> (DLMF 8.7.3), whose terms are of one sign on the negative axis, with
/// <c>(|z| + Re z)/ln 10</c>. Each is taken again with more where the difference shows it cancelled more.
/// </para>
/// <para>
/// At a whole <c>a = -n &lt;= 0</c>, where <c>Gamma(a)</c> has a pole, it is
Expand Down Expand Up @@ -76,6 +77,23 @@ public static Complex IncompleteGamma(Complex a, Complex z)
var (re, im) = useTheFraction
? ByContinuedFraction((ar, ai), (zr, zi), Working(context, 10))
: BySeries((ar, ai), (zr, zi), modulus, context, turn);
// On the cut, for a real a, the fraction's coefficients are real and it converges to e^(i pi a) times a real
// function, the mean of the value from above and its conjugate from below, short of the value from above by
// -i e^(i pi a) sin(pi a) Gamma(a). That is below the precision of a value e^|z| in size, but not of a part it
// is the whole of: Gamma(a) at a half-whole a, the real part of Gamma(1/2, -300), which is sqrt(pi), and
// -i (-1)^n pi/n! at a = -n, the imaginary part of Gamma(-3, -300), pi/6. Without it, both are 0.
if (useTheFraction && turn is null && ai.IsZero && zi.IsZero && zr.Sign < 0)
{
var work = Working(context, 10);
if (!ar.IsInteger() && ar.Multiply(EDecimal.FromInt32(2)).IsInteger())
re = re.Add(GammaOf((ar, ai), work).Re, work);
else if (ar.IsInteger() && ar.Sign <= 0)
{
var n = ar.ToEInteger().Negate();
var piOverFactorial = InternalAMExtensions.ConstantCache.Lookup(work).Pi.Divide(EDecimal.FromEInteger(Factorial(n.ToInt32Checked())), work);
im = n.IsEven ? im.Subtract(piOverFactorial, work) : im.Add(piOverFactorial, work);
}
}
// The continued fraction is used on the negative axis only past (P + 5) ln 10, where the value is e^|z| in
// size and (1 - rho) Gamma(a) + rho times the principal value cancels nothing.
if (turn is { } fractionsTurn && useTheFraction)
Expand Down Expand Up @@ -149,13 +167,21 @@ private static (EDecimal Re, EDecimal Im) Turned((EInteger Turns, EInteger Q) tu
private static (EDecimal Re, EDecimal Im) BySeries((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, double modulus, EContext context,
(EInteger Turns, EInteger Q)? turn)
{
var guard = (int)Math.Ceiling(2 * modulus / Math.Log(10)) + 10;
// In the left half-plane the lower function is summed as z^a sum_k (-z)^k/(k! (a + k)), as E_1 is, whose terms
// of size e^|z| sum to e^(-Re z): it loses (|z| + Re z)/ln 10 digits, none on the negative axis, where z^a
// e^(-z) sum_k z^k/(a)_(k + 1) has terms of alternating sign that sum to e^(-|z|) and lost twice |z|/ln 10.
// Carrying those digits asked Gamma(a) at a precision of its own, where it is not remembered: Gamma(1/3, -143)
// took 9.6 seconds, and Gamma(-5/2, -60) 2.4.
var leftHalf = z.Re.Sign < 0;
var guard = leftHalf
? (int)Math.Ceiling(Math.Max(0, modulus + z.Re.ToDouble()) / Math.Log(10)) + 12
: (int)Math.Ceiling(2 * modulus / Math.Log(10)) + 10;
for (var attempt = 0; ; attempt++)
{
var work = Working(context, guard);
var (value, scale) = a.Im.IsZero && a.Re.IsInteger() && a.Re.Sign <= 0
? AtANonPositiveWholeOrder(a.Re.ToEInteger().Negate().ToInt32Checked(), z, work)
: AsGammaLessTheLower(a, z, work, turn);
: AsGammaLessTheLower(a, z, leftHalf, work, turn);
// How many digits the value lost to the cancellation, against the largest thing that cancelled.
var size = Magnitude(value);
var lost = size.IsZero ? guard : DecimalOrder(scale) - DecimalOrder(size) + 1;
Expand All @@ -166,13 +192,16 @@ private static (EDecimal Re, EDecimal Im) BySeries((EDecimal Re, EDecimal Im) a,
}

/// <summary>
/// <c>Gamma(a) - z^a e^(-z) sum_k z^k / (a)_(k + 1)</c>, and the size of the larger of the two terms.
/// <c>Gamma(a) - z^a e^(-z) sum_k z^k / (a)_(k + 1)</c>, or with <paramref name="leftHalf"/>
/// <c>Gamma(a) - z^a sum_k (-z)^k / (k! (a + k))</c>, and the size of the larger of the two terms.
/// </summary>
private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessTheLower((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, EContext work,
(EInteger Turns, EInteger Q)? turn)
private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessTheLower((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, bool leftHalf,
EContext work, (EInteger Turns, EInteger Q)? turn)
{
var epsilon = EDecimal.Create(EInteger.One, EInteger.FromInt32(-work.Precision.ToInt32Checked()));
// term_k = z^k / (a (a + 1) ... (a + k)), starting at 1/a
// term_k = z^k / (a (a + 1) ... (a + k)), or (-z)^k / (k! (a + k)), each starting at 1/a, and each term the
// last times z/(a + k), or times -z (a + k - 1)/(k (a + k))
var minusZ = (Re: z.Re.Negate(), Im: z.Im.Negate());
(EDecimal Re, EDecimal Im) sum;
EDecimal peak;
if (a.Im.IsZero && z.Im.IsZero)
Expand All @@ -185,8 +214,11 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessThe
var aPlusK = a.Re;
for (var k = 1; ; k++)
{
var previous = aPlusK;
aPlusK = aPlusK.Add(EDecimal.One, work);
term = term.Multiply(z.Re, work).Divide(aPlusK, work);
term = leftHalf
? term.Multiply(minusZ.Re, work).Multiply(previous, work).Divide(aPlusK.Multiply(k, work), work)
: term.Multiply(z.Re, work).Divide(aPlusK, work);
real = real.Add(term, work);
var size = term.Abs();
if (size.CompareTo(peak) > 0)
Expand All @@ -204,8 +236,11 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessThe
var aPlusK = a;
for (var k = 1; ; k++)
{
var previous = aPlusK;
aPlusK = (aPlusK.Re.Add(EDecimal.One, work), aPlusK.Im);
term = Multiply(Multiply(term, z, work), Reciprocal(aPlusK, work), work);
term = leftHalf
? Multiply(Multiply(term, minusZ, work), Multiply(previous, Reciprocal((aPlusK.Re.Multiply(k, work), aPlusK.Im.Multiply(k, work)), work), work), work)
: 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)
Expand All @@ -217,7 +252,7 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessThe
var power = ComplexPower(z, a, work);
if (turn is { } realRoot)
power = Multiply(Turned(realRoot, work), power, work);
var prefactor = Multiply(power, ComplexExp((z.Re.Negate(), z.Im.Negate()), work), work);
var prefactor = leftHalf ? power : Multiply(power, ComplexExp(minusZ, work), work);
var lower = Multiply(prefactor, sum, work);
var gamma = GammaOf(a, work);
var scale = EDecimal.Max(Magnitude(gamma), Magnitude(prefactor).Multiply(peak, work));
Expand Down Expand Up @@ -328,7 +363,37 @@ private static (EDecimal Re, EDecimal Im) ComplexLn((EDecimal Re, EDecimal Im) z

/// <summary><c>z^a = e^(a ln z)</c>, principal.</summary>
private static (EDecimal Re, EDecimal Im) ComplexPower((EDecimal Re, EDecimal Im) z, (EDecimal Re, EDecimal Im) a, EContext work)
=> ComplexExp(Multiply(a, ComplexLn(z, work), work), work);
{
// On the cut, for a real a, |z|^a (cos pi a + i sin pi a) with each part accurate to its own size. From pi/2
// rounded to W digits, cos(pi/2) is 10^(-W) where it is 0, and times the lower function's e^143 that was the
// real part of Gamma(1/2, -143), which is sqrt(pi).
if (z.Im.IsZero && z.Re.Sign < 0 && a.Im.IsZero)
{
var modulus = z.Re.Negate().Log(work).Multiply(a.Re, work).Exponential(work);
var (cos, sin) = CosAndSinOfPiTimes(a.Re, work);
return (modulus.Multiply(cos, work), modulus.Multiply(sin, work));
}
return ComplexExp(Multiply(a, ComplexLn(z, work), work), work);
}

/// <summary>
/// <c>(cos pi t, sin pi t)</c>, each accurate to its own size: with <c>n</c> the whole number nearest <c>t</c> and
/// <c>s = 1/2 - |t - n|</c>, both exact, they are <c>(-1)^n sin pi s</c> and <c>(-1)^n sgn(t - n) cos pi s</c>, and
/// <c>cos pi t</c> is exactly 0 at a half-whole <c>t</c>.
/// </summary>
private static (EDecimal Cos, EDecimal Sin) CosAndSinOfPiTimes(EDecimal t, EContext work)
{
var n = t.RoundToExponent(EInteger.Zero, ERounding.HalfEven);
var offset = t.Subtract(n);
var s = EDecimal.Create(5, -1).Subtract(offset.Abs());
var sign = n.ToEInteger().IsEven ? EDecimal.One : EDecimal.One.Negate();
if (s.IsZero)
return (EDecimal.Zero, offset.Sign > 0 ? sign : sign.Negate());
var angle = InternalAMExtensions.ConstantCache.Lookup(work).Pi.Multiply(s, work);
var cos = angle.Sin(work).Multiply(sign, work);
var sin = offset.IsZero ? EDecimal.Zero : angle.Cos(work).Multiply(offset.Sign > 0 ? sign : sign.Negate(), work);
return (cos, sin);
}
}
}
}
13 changes: 12 additions & 1 deletion Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.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 @@ -45,6 +45,16 @@ private static void AssertClose(string reference, EDecimal actual)
// On the cut, from above.
[InlineData("gamma(1/2, -3)", "1.77245385090551602729816748334114518279754945612238712821381", "-14.6261713840190916850829290647365438829669286131782239586917")]
[InlineData("gamma(-1, -20)", "1.35739289356707492202577927702104389515074891049797988869358E+6", "3.14159265358979323846264338327950288419716939937510582097494")]
// Far into the left half-plane, by the lower function's series of one sign there: each part to its own size,
// the real part of the first being Gamma(1/2).
[InlineData("gamma(1/2, -143)", "1.772453850905516027298167483341145182797549456122387128213808", "-1.066534514082282609180752885388079901236847457667946132456796E+61")]
[InlineData("gamma(-5/2, -60)", "-0.9453087204829418812256893244486107641586930432652731350473642", "7.256677574613207107369164321252421746598199910922310065217093E+19")]
[InlineData("gamma(1/2, -100 + 5i)", "2.607638471329857822055492550981639715386578693808752929552738E+42", "-7.00326380520744724063748887540346529570189148404061756494028E+41")]
// Past the continued fraction's threshold on the cut, with the part that Gamma(a) is the whole of.
[InlineData("gamma(1/2, -300)", "1.772453850905516027298167483341145182797549456122387128213808", "-1.123338927148883734198084176745306173925451665800872666587958E+129")]
[InlineData("gamma(3/2, -250)", "0.8862269254527580136490837416705725913987747280611935641069039", "5.911793698136623273673640330606776777031254613823953153197099E+109")]
[InlineData("gamma(-3, -300)", "2.430575190291251015731676197614951916724433231125758635038235E+120", "0.5235987755982988730771072305465838140328615665625176368291574")]
[InlineData("gamma(0, -300)", "-6.496482508088665789025691894934247269238089700566375975103098E+127", "-3.141592653589793238462643383279502884197169399375105820974945")]
[InlineData("gamma(100, 10)", "9.3326215443944152681699238856266700490715968264381621468593E+155", "0")]
[InlineData("gamma(1/2 + i, 2 + 3i)", "-0.0106186868505695828645674618919283599789635032416382778377351", "-0.026574511726313175577690285230774989836262216835255934709977")]
public void IsItsValue(string expression, string real, string imaginary)
Expand All @@ -63,6 +73,7 @@ public void IsItsValue(string expression, string real, string imaginary)
[InlineData("gamma(1/3, -2)", "10.212552064243320672871765046103178972863292927441879189946367")]
[InlineData("gamma(2/3, -4)", "-36.897696051078160810409733567406007863692044903007209488356707")]
[InlineData("gamma(2/5, -17/15)", "-1.6376459092577785167485577680520445343763746015778210809826328")]
[InlineData("gamma(1/3, -143)", "4.669432864644246857974117704263949442598436311750677119332775E+60")]
[InlineData("gamma(-1/3, -300)", "9.71532351293767903599760645811774844613945799010165863950178E+126")]
public void OnTheCutItTakesThePowersRealRoot(string expression, string real)
{
Expand Down
Loading