From 9f03081a9bf9fa9767022984f105872a34d17409 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 09:52:47 +0000 Subject: [PATCH] Sum the incomplete gamma function's lower series as DLMF 8.7.3 in the left half-plane Gamma(a) - gamma(a, z) summed gamma(a, z) as z^a e^(-z) sum_k z^k/(a)_(k + 1) (DLMF 8.7.1) everywhere. On the negative axis its terms alternate, are e^|z| in size and sum to e^(-|z|), so twice |z|/ln 10 guard digits were carried, and at that precision Gamma(a) was computed afresh: Gamma(1/3, -143) took 9.6 s, Gamma(1/2, -100 + 5i) 5.6 s, Gamma(-5/2, -60) 2.4 s. In the left half-plane it is now z^a sum_k (-z)^k/(k! (a + k)) (DLMF 8.7.3), as E_1 already was, whose terms are of one sign on that axis and which loses (|z| + Re z)/ln 10 digits: 6 ms, 9 ms and 6 ms. The large guard had been hiding two defects on the cut, where each part of the value should be accurate to its own size: - z^a was e^(a ln|z| + i a pi) with a pi rounded, so cos(pi/2) was 10^(-W) where it is 0. Times e^143, that was the real part of Gamma(1/2, -143), which is sqrt(pi). The angle is reduced exactly before pi multiplies it, and is exactly 0 at a half-whole a. - Past (P + 5) ln 10 the continued fraction, whose coefficients are real there, gives e^(i pi a) times a real function, the mean of the two sides. It falls short of the value from above by -i e^(i pi a) sin(pi a) Gamma(a). That is below the precision of the modulus, but it is the whole of the real part at a half-whole a and of the imaginary part at a = -n. On master Gamma(1/2, -300)'s real part was 3e13 times sqrt(pi), Gamma(11/2, -300)'s 1.6e25 times Gamma(11/2), and Gamma(-3, -300)'s imaginary part 9e4 times pi/6. Checked part by part against mpmath at 300-400 digits on 320 points: 16 orders, half-whole, whole, complex and odd-denominator ones among them, at 20 points of the left half-plane including both sides of the threshold. Every part agrees to 1e-97 of its own size. On master, 18 of the 80 cases past the threshold had a part wrong, by up to 1e113 times its size. The grid took 48.4 s on master and 3.4 s here. Part of #2075. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Continuous/Number/IncompleteGamma.cs | 93 ++++++++++++++++--- .../Convenience/IncompleteGammaTest.cs | 13 ++- 2 files changed, 91 insertions(+), 15 deletions(-) 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) {