diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs
index a4f4511c4..2f7889784 100644
--- a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs
+++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs
@@ -28,10 +28,11 @@ partial record Number
/// algorithm. Near the negative real axis at a moderate |z| it stops short of the value.
///
///
- /// Short of that it is Gamma(a) - gamma(a, z), the lower one summed as
- /// z^a e^(-z) sum_k z^k / (a (a + 1) ... (a + k)) (DLMF 8.7.1). The terms are about e^|z| and the
- /// value can be as small as e^(-|z|), so twice |z|/ln 10 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 Gamma(a) - gamma(a, z), the lower one summed in the right half-plane as
+ /// z^a e^(-z) sum_k z^k / (a (a + 1) ... (a + k)) (DLMF 8.7.1), whose terms there are about e^|z|
+ /// and of one sign on the positive axis, with twice |z|/ln 10 guard digits; and in the left half-plane as
+ /// z^a sum_k (-z)^k / (k! (a + k)) (DLMF 8.7.3), whose terms are of one sign on the negative axis, with
+ /// (|z| + Re z)/ln 10. Each is taken again with more where the difference shows it cancelled more.
///
///
/// At a whole a = -n <= 0, where Gamma(a) has a pole, it is
@@ -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)
@@ -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;
@@ -166,13 +192,16 @@ private static (EDecimal Re, EDecimal Im) BySeries((EDecimal Re, EDecimal Im) a,
}
///
- /// Gamma(a) - z^a e^(-z) sum_k z^k / (a)_(k + 1), and the size of the larger of the two terms.
+ /// Gamma(a) - z^a e^(-z) sum_k z^k / (a)_(k + 1), or with
+ /// Gamma(a) - z^a sum_k (-z)^k / (k! (a + k)), and the size of the larger of the two terms.
///
- 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)
@@ -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)
@@ -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)
@@ -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));
@@ -328,7 +363,37 @@ private static (EDecimal Re, EDecimal Im) ComplexLn((EDecimal Re, EDecimal Im) z
/// z^a = e^(a ln z), principal.
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);
+ }
+
+ ///
+ /// (cos pi t, sin pi t), each accurate to its own size: with n the whole number nearest t and
+ /// s = 1/2 - |t - n|, both exact, they are (-1)^n sin pi s and (-1)^n sgn(t - n) cos pi s, and
+ /// cos pi t is exactly 0 at a half-whole t.
+ ///
+ 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);
+ }
}
}
}
diff --git a/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs
index e99269e85..3d74fac35 100644
--- a/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs
+++ b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.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.
@@ -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)
@@ -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)
{