From 8ccb118da3a0a21865e9bbe0471f88c6f59d37ed Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 22 Oct 2026 20:08:46 +0000 Subject: [PATCH] The upper incomplete gamma function gamma(a, z) with a second argument is the upper incomplete gamma function Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt, as Mathematica's Gamma[a, z] and mpmath's gammainc(a, z) are; with one, gamma is still Gamma(x), and a second argument was a parse error. z^a takes its principal branch, so the cut is the negative real axis, which keeps the value from above. It is printed so, in LaTeX as \Gamma(a, z) and in SymPy as uppergamma, and its derivative in z is -z^(a - 1) e^(-z); a parameter that varies declines. A whole a >= 1 is (n - 1)! e^(-z) times a polynomial (DLMF 8.4.8), a half is sqrt(pi) erfc(sqrt(z)) (8.4.6), and Gamma(a, 0) is Gamma(a). It is evaluated to the working precision for any complex a and z. Where |z| > 2|a| + 10, in the right half-plane or past (P + 5) ln 10, it is Legendre's continued fraction (8.9.2) by the modified Lentz algorithm; elsewhere it is Gamma(a) less the lower one's series (8.7.1), with twice |z|/ln 10 guard digits and more where the difference shows it cancelled more; and at a whole a <= 0, where Gamma(a) has a pole, E_1(z) by its series and the recurrence (8.4.15). It is real exactly where it is real: for a positive z, and for a whole a >= 1 at any real z. It agrees with mpmath at 168 points to 1e-90 at a hundred digits and at 54 to 1e-290 at three hundred, near the poles and on the cut. Part of #2075. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 17 ++ Sources/.editorconfig | 6 + Sources/AngouriMath/Convenience/MathS.cs | 13 + Sources/AngouriMath/Core/Antlr/AngouriMath.g | 4 +- .../Core/Antlr/AngouriMathParser.cs | 2 +- Sources/AngouriMath/Core/Domains.Classes.cs | 8 + ...tity.Continuous.IncompleteGamma.Classes.cs | 36 +++ .../Continuous/Number/IncompleteGamma.cs | 239 ++++++++++++++++++ .../Serialization/Entity.Serialization.cs | 1 + .../Transformations/Matching/MatchPattern.cs | 1 + .../Docs/Contributing/WritingARule.md | 2 +- Sources/AngouriMath/Docs/Usage/Syntax.md | 8 + .../Functions/Continuous/Differentiation.cs | 11 + .../EquationSolver/InvertNode.Classes.cs | 7 + ...aluation.Continuous.Arithmetics.Classes.cs | 38 +++ .../Output/Latex/Latex.Arithmetics.Classes.cs | 7 + .../ToString/ToString.Arithmetics.Classes.cs | 8 + .../ToSympy/ToSympy.Arithmetics.Classes.cs | 6 + Sources/AngouriMath/Functions/Substitute.cs | 7 + .../Functions/TreeAnalyzer/Sort.Classes.cs | 6 + Sources/MCP/AngouriMath.Mcp/Resources.cs | 3 +- Sources/Tests/UnitTests/Common/PublicApi.txt | 25 ++ .../Convenience/IncompleteGammaTest.cs | 90 +++++++ .../Transformations/RuleAuthoringGuideTest.cs | 2 +- 24 files changed, 542 insertions(+), 5 deletions(-) create mode 100644 Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.IncompleteGamma.Classes.cs create mode 100644 Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs create mode 100644 Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 726403fd2..a6acbfd03 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -129,6 +129,23 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `"polygamma(1, 3)".ToEntity().InnerSimplified` | an `UnhandledParseException` | `1/6 * pi ^ 2 - 5/4` | | `"digamma(x)".ToEntity().Integrate("x")` | `digamma * x ^ 2 / 2 + C`, of the product | `ln((x - 1)!) + C` | +### `gamma(a, z)` is the upper incomplete gamma function + +**Not silent: what was a parse error is a function.** `gamma` took one argument, and a second was a +`FunctionArgumentCountException`. With two it is now the upper incomplete gamma function, +`Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt`, as Mathematica's `Gamma[a, z]` is. +- **Branch.** It takes the principal branch of `z^a`, with the cut on the negative real axis. +- **Evaluation.** It is evaluated to the working precision for any complex `a` and `z`. +- **Derivative.** In `z` it is `-z^(a - 1) e^(-z)`. +- **Exact forms.** A whole `a >= 1` and `a = 1/2` are elementary. +([#2075](https://github.com/asc-community/AngouriMath/issues/2075)) + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"gamma(1/3, x)".ToEntity()` | a `FunctionArgumentCountException` | `gamma(1/3, x)` | +| `"gamma(1/3, 2)".EvalNumerical()` | a `FunctionArgumentCountException` | `0.068136414441459136…` | +| `"gamma(3, x)".ToEntity().InnerSimplified` | a `FunctionArgumentCountException` | `2 * e ^ (-x) * (1 + x + x ^ 2 / 2)` | + ### Numerical equality is to half the evaluation's digits **Silent: wrong answers.** `MathS.UnsafeAndInternal.AreEqualNumerically` evaluated both sides at its check points and diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 46679c092..43723eecc 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -699,6 +699,12 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [AngouriMath/Core/Entity/Continuous/Number/Polygamma.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[AngouriMath/Core/Entity/Continuous/Entity.Continuous.IncompleteGamma.Classes.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [AngouriMath/Core/Entity/Continuous/Entity.Continuous.Polygamma.Classes.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index 3b01c00dc..5e937153c 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -1474,6 +1474,19 @@ public static Entity VectorLength(Matrix v) [MethodImpl(MethodImplOptions.AggressiveInlining)] public static Entity Digamma(Entity argument) => new Polygammaf(Integer.Zero, argument); + /// The upper incomplete gamma function, Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt + /// The parameter a, any complex number + /// The lower limit z, any complex number; 0 only where Re a > 0 + /// The node + /// + /// z^a takes its principal branch, so the cut is the negative real axis, which keeps the value from above. + /// Gamma(a, 0) is Gamma(a), and the derivative in z is -z^(a - 1) e^(-z). It is written + /// gamma(a, z), as Mathematica's Gamma[a, z] and mpmath's gammainc(a, z) are. + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static Entity IncompleteGamma(Entity parameter, Entity lowerLimit) => new IncompleteGammaf(parameter, lowerLimit); + /// The sine integral, Si(z) = int_0^z sin(t)/t dt /// The argument, any complex number /// The node diff --git a/Sources/AngouriMath/Core/Antlr/AngouriMath.g b/Sources/AngouriMath/Core/Antlr/AngouriMath.g index efcf3b848..c3831dd1a 100644 --- a/Sources/AngouriMath/Core/Antlr/AngouriMath.g +++ b/Sources/AngouriMath/Core/Antlr/AngouriMath.g @@ -543,7 +543,9 @@ atom returns[Entity value] | 'factorial(' args = function_arguments ')' { Assert("factorial", 1, $args.list.Count); $value = MathS.Factorial($args.list[0]); } - | 'gamma(' args = function_arguments ')' { Assert("gamma", 1, $args.list.Count); $value = MathS.Gamma($args.list[0]); } + /* With a second argument, the upper incomplete gamma function Gamma(a, z), as Mathematica's Gamma[a, z]. + https://github.com/asc-community/AngouriMath/issues/2075 */ + | 'gamma(' args = function_arguments ')' { if ($args.list.Count == 2) $value = MathS.IncompleteGamma($args.list[0], $args.list[1]); else { Assert("gamma", 1, $args.list.Count); $value = MathS.Gamma($args.list[0]); } } | 'derivative(' args = function_arguments ')' { if (Assert("derivative", (3, 2), $args.list.Count)) diff --git a/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs b/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs index 0afd596d3..95dfd7743 100644 --- a/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs +++ b/Sources/AngouriMath/Core/Antlr/AngouriMathParser.cs @@ -4369,7 +4369,7 @@ public AtomContext atom() { _localctx.args = function_arguments(); State = 886; Match(T__30); - Assert("gamma", 1, _localctx.args.list.Count); _localctx.value = MathS.Gamma(_localctx.args.list[0]); + if (_localctx.args.list.Count == 2) _localctx.value = MathS.IncompleteGamma(_localctx.args.list[0], _localctx.args.list[1]); else { Assert("gamma", 1, _localctx.args.list.Count); _localctx.value = MathS.Gamma(_localctx.args.list[0]); } } break; case 94: diff --git a/Sources/AngouriMath/Core/Domains.Classes.cs b/Sources/AngouriMath/Core/Domains.Classes.cs index 08a75c9c7..fb56d1790 100644 --- a/Sources/AngouriMath/Core/Domains.Classes.cs +++ b/Sources/AngouriMath/Core/Domains.Classes.cs @@ -384,6 +384,14 @@ partial record Polygammaf internal override Domain DefaultCodomain => Domain.Complex; } + partial record IncompleteGammaf + { + /// + public override Domain Codomain { get; protected init; } = Domain.Complex; + /// + internal override Domain DefaultCodomain => Domain.Complex; + } + partial record Sif { /// diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.IncompleteGamma.Classes.cs b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.IncompleteGamma.Classes.cs new file mode 100644 index 000000000..ad1489784 --- /dev/null +++ b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.IncompleteGamma.Classes.cs @@ -0,0 +1,36 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; + +namespace AngouriMath +{ + partial record Entity + { + /// + /// The upper incomplete gamma function, Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt, with the principal branch + /// of z^a: the cut is the negative real axis. + /// + /// https://github.com/asc-community/AngouriMath/issues/2075 + public sealed partial record IncompleteGammaf(Entity Parameter, Entity LowerLimit) : Function, IBinaryNode + { + private IncompleteGammaf New(Entity parameter, Entity lowerLimit) => + ReferenceEquals(Parameter, parameter) && ReferenceEquals(LowerLimit, lowerLimit) ? this : new(parameter, lowerLimit) { Codomain = Codomain }; + + /// + public Entity NodeFirstChild => Parameter; + + /// + public Entity NodeSecondChild => LowerLimit; + + /// + public override Entity Replace(Func func) => func(New(Parameter.Replace(func), LowerLimit.Replace(func))); + /// + protected override Entity[] InitDirectChildren() => new[] { Parameter, LowerLimit }; + } + } +} diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs new file mode 100644 index 000000000..b870d7dc4 --- /dev/null +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs @@ -0,0 +1,239 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using PeterO.Numbers; + +namespace AngouriMath +{ + partial record Entity + { + partial record Number + { + /// + /// The upper incomplete gamma function, Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt, at the working + /// precision, for complex a and z, with the principal branch of z^a: the cut is the + /// negative real axis, which keeps the value from above. + /// + /// + /// + /// Where |z| > 2 |a| + 10, in the right half-plane or anywhere past |z| = (P + 5) ln 10 with P + /// digits of working precision, it is Legendre's continued fraction (DLMF 8.9.2), + /// e^(-z) z^a / (z + 1 - a - 1 (1 - a) / (z + 3 - a - 2 (2 - a) / ...)), by the modified Lentz + /// 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. + /// + /// + /// At a whole a = -n <= 0, where Gamma(a) has a pole, it is + /// Gamma(-n, z) = (-1)^n/n! (E_1(z) - e^(-z) sum_(k < n) (-1)^k k!/z^(k + 1)) (DLMF 8.4.15), with + /// E_1(z) = Gamma(0, z) = -gamma - ln z - sum_(k >= 1) (-z)^k/(k k!) (DLMF 6.6.2). + /// + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// + public static Complex IncompleteGamma(Complex a, Complex z) + { + var context = MathS.Settings.DecimalPrecisionContext.Value; + var (ar, ai) = (a.RealPart.EDecimal, a.ImaginaryPart.EDecimal); + var (zr, zi) = (z.RealPart.EDecimal, z.ImaginaryPart.EDecimal); + if (ar.IsNaN() || ai.IsNaN() || zr.IsNaN() || zi.IsNaN() || !ar.IsFinite || !ai.IsFinite) + return Real.NaN; + if (zr.IsInfinity() || zi.IsInfinity()) + return zi.IsZero && zr.IsPositiveInfinity() ? Integer.Zero : Real.NaN; + if (zr.IsZero && zi.IsZero) + { + // Gamma(a, 0) is Gamma(a) where the integral converges at 0, and grows without bound where it + // does not, which for a real a is +oo. + if (ar.Sign > 0) + return Factorial(a - Integer.One); + return ai.IsZero ? Real.PositiveInfinity : Real.NaN; + } + var digits = context.Precision.ToInt32Checked(); + var modulus = Math.Sqrt(Math.Pow(zr.ToDouble(), 2) + Math.Pow(zi.ToDouble(), 2)); + // The continued fraction where it converges quickly and to the principal branch: where |z| is large against + // |a|, in the right half-plane, or anywhere past (P + 5) ln 10. Near the negative real axis at a moderate + // |z| it stops short -- Gamma(1/2, -10) came out wrong in its fourth digit -- and with |a| large against |z| + // it converges slowly -- Gamma(100, 1/10 - 10 i) had the wrong imaginary part -- and there the series is + // summed instead. + var order = Math.Sqrt(Math.Pow(ar.ToDouble(), 2) + Math.Pow(ai.ToDouble(), 2)); + var useTheFraction = modulus > 2 * order + 10 && (zr.Sign >= 0 || modulus > (digits + 5) * Math.Log(10)); + var (re, im) = useTheFraction + ? ByContinuedFraction((ar, ai), (zr, zi), Working(context, 10)) + : BySeries((ar, ai), (zr, zi), modulus, context); + // Real where it is real, exactly: for a real a and a positive z, and for a positive whole a and any real + // z, where Gamma(n, z) is (n - 1)! e^(-z) times a polynomial; z^a leaves a residue of e^(i pi n) there. + if (ai.IsZero && zi.IsZero && (zr.Sign > 0 || ar.IsInteger() && ar.Sign > 0)) + im = EDecimal.Zero; + // Otherwise each part is accurate to its own size, whatever that is, and is never rounded onto zero by + // size: Gamma(1/2, 300) is 3e-132, and Gamma(100, -5 + i/10) has an imaginary part 1e-86 of its real one. + return Complex.CreateAtScales(re.RoundToPrecision(context), im.RoundToPrecision(context), EDecimal.Zero, EDecimal.Zero); + } + + private static (EDecimal Re, EDecimal Im) ByContinuedFraction((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, EContext work) + { + var epsilon = EDecimal.Create(EInteger.One, EInteger.FromInt32(-work.Precision.ToInt32Checked())); + var tiny = EDecimal.Create(EInteger.One, EInteger.FromInt32(-10 * work.Precision.ToInt32Checked())); + (EDecimal Re, EDecimal Im) Guard((EDecimal Re, EDecimal Im) v) => Magnitude(v).CompareTo(tiny) < 0 ? (tiny, EDecimal.Zero) : v; + // b = z + 1 - a + var b = (Re: z.Re.Add(EDecimal.One, work).Subtract(a.Re, work), Im: z.Im.Subtract(a.Im, work)); + var c = (Re: EDecimal.One.Divide(tiny, work), Im: EDecimal.Zero); + var d = Reciprocal(Guard(b), work); + var h = d; + for (var i = 1; i <= 50 * work.Precision.ToInt32Checked() + 1000; i++) + { + // a_i = -i (i - a) + var an = (Re: EDecimal.FromInt32(-i).Multiply(EDecimal.FromInt32(i).Subtract(a.Re, work), work), Im: EDecimal.FromInt32(i).Multiply(a.Im, work)); + b = (b.Re.Add(EDecimal.FromInt32(2), work), b.Im); + var ad = Multiply(an, d, work); + d = Reciprocal(Guard((ad.Re.Add(b.Re, work), ad.Im.Add(b.Im, work))), work); + var anOverC = Multiply(an, Reciprocal(c, work), work); + c = Guard((b.Re.Add(anOverC.Re, work), b.Im.Add(anOverC.Im, work))); + var delta = Multiply(d, c, work); + h = Multiply(h, delta, work); + if (Magnitude((delta.Re.Subtract(EDecimal.One, work), delta.Im)).CompareTo(epsilon) < 0) + break; + } + return Multiply(Multiply(ComplexExp((z.Re.Negate(), z.Im.Negate()), work), ComplexPower(z, a, work), work), h, work); + } + + private static (EDecimal Re, EDecimal Im) BySeries((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, double modulus, EContext context) + { + var guard = (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); + // 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; + if (lost <= guard - 10 || attempt >= 3) + return value; + guard = lost + 20; + } + } + + /// + /// Gamma(a) - z^a e^(-z) sum_k z^k / (a)_(k + 1), 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) + { + 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++) + { + 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 scale = EDecimal.Max(Magnitude(gamma), Magnitude(prefactor).Multiply(peak, work)); + return ((gamma.Re.Subtract(lower.Re, work), gamma.Im.Subtract(lower.Im, work)), scale); + } + + /// + /// Gamma(-n, z) = (-1)^n/n! (E_1(z) - e^(-z) sum_(k < n) (-1)^k k!/z^(k + 1)), and the size of what + /// cancelled. + /// + private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AtANonPositiveWholeOrder(int n, (EDecimal Re, EDecimal Im) z, EContext work) + { + var epsilon = EDecimal.Create(EInteger.One, EInteger.FromInt32(-work.Precision.ToInt32Checked())); + // E_1(z) = -gamma - ln z - sum_(k >= 1) (-z)^k/(k k!) + var minusZ = (Re: z.Re.Negate(), Im: z.Im.Negate()); + var power = (Re: EDecimal.One, Im: EDecimal.Zero); + var sum = (Re: EDecimal.Zero, Im: EDecimal.Zero); + var peak = EDecimal.Zero; + for (var k = 1; ; k++) + { + // (-z)^k / k! + power = Multiply(power, minusZ, work); + power = (power.Re.Divide(k, work), power.Im.Divide(k, work)); + var term = (Re: power.Re.Divide(k, work), Im: power.Im.Divide(k, 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 ln = ComplexLn(z, work); + var e1 = (Re: EulerGamma(work).Add(ln.Re, work).Add(sum.Re, work).Negate(), Im: ln.Im.Add(sum.Im, work).Negate()); + var scale = EDecimal.Max(peak, Magnitude(ln)); + if (n == 0) + return (e1, scale); + // e^(-z) sum_(k < n) (-1)^k k!/z^(k + 1) + var inverse = Reciprocal(z, work); + var inversePower = inverse; + var tail = (Re: EDecimal.Zero, Im: EDecimal.Zero); + var factorial = EDecimal.One; + for (var k = 0; k < n; k++) + { + if (k > 0) + { + factorial = factorial.Multiply(k, work); + inversePower = Multiply(inversePower, inverse, work); + } + var term = (Re: inversePower.Re.Multiply(factorial, work), Im: inversePower.Im.Multiply(factorial, work)); + tail = k % 2 == 0 ? (tail.Re.Add(term.Re, work), tail.Im.Add(term.Im, work)) : (tail.Re.Subtract(term.Re, work), tail.Im.Subtract(term.Im, work)); + } + var exponentialTail = Multiply(ComplexExp(minusZ, work), tail, work); + scale = EDecimal.Max(scale, Magnitude(exponentialTail)); + var difference = (Re: e1.Re.Subtract(exponentialTail.Re, work), Im: e1.Im.Subtract(exponentialTail.Im, work)); + // (-1)^n / n! + var nFactorial = EDecimal.FromEInteger(Factorial(n)); + var sign = n % 2 == 0 ? EDecimal.One : EDecimal.One.Negate(); + var coefficient = sign.Divide(nFactorial, work); + return ((difference.Re.Multiply(coefficient, work), difference.Im.Multiply(coefficient, work)), scale.Multiply(coefficient.Abs(), work)); + } + + /// 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(); + + private static (EDecimal Re, EDecimal Im) ComplexExp((EDecimal Re, EDecimal Im) z, EContext work) + { + var modulus = z.Re.Exponential(work); + if (z.Im.IsZero) + return (modulus, EDecimal.Zero); + return (modulus.Multiply(z.Im.Cos(work), work), modulus.Multiply(z.Im.Sin(work), work)); + } + + /// The principal logarithm: ln |z| + i arg z, with arg in (-pi, pi]. + private static (EDecimal Re, EDecimal Im) ComplexLn((EDecimal Re, EDecimal Im) z, EContext work) + { + if (z.Im.IsZero && z.Re.Sign > 0) + return (z.Re.Log(work), EDecimal.Zero); + 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)); + } + + /// 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); + } + } +} diff --git a/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs b/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs index e9ff31a3a..9c9e3e340 100644 --- a/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs +++ b/Sources/AngouriMath/Core/Serialization/Entity.Serialization.cs @@ -127,6 +127,7 @@ [EntityJsonConverter] partial record SpecialSet [EntityJsonConverter] partial record Eif; [EntityJsonConverter] partial record Lif; [EntityJsonConverter] partial record Polygammaf; + [EntityJsonConverter] partial record IncompleteGammaf; [EntityJsonConverter] partial record Sif; [EntityJsonConverter] partial record Cif; [EntityJsonConverter] partial record Shif; diff --git a/Sources/AngouriMath/Core/Transformations/Matching/MatchPattern.cs b/Sources/AngouriMath/Core/Transformations/Matching/MatchPattern.cs index d0fc54f57..a82524217 100644 --- a/Sources/AngouriMath/Core/Transformations/Matching/MatchPattern.cs +++ b/Sources/AngouriMath/Core/Transformations/Matching/MatchPattern.cs @@ -479,6 +479,7 @@ internal abstract bool ETryBuild( [typeof(Entity.Eif)] = (1, static c => new Entity.Eif(c[0])), [typeof(Entity.Lif)] = (1, static c => new Entity.Lif(c[0])), [typeof(Entity.Polygammaf)] = (2, static c => new Entity.Polygammaf(c[0], c[1])), + [typeof(Entity.IncompleteGammaf)] = (2, static c => new Entity.IncompleteGammaf(c[0], c[1])), [typeof(Entity.Sif)] = (1, static c => new Entity.Sif(c[0])), [typeof(Entity.Cif)] = (1, static c => new Entity.Cif(c[0])), [typeof(Entity.Shif)] = (1, static c => new Entity.Shif(c[0])), diff --git a/Sources/AngouriMath/Docs/Contributing/WritingARule.md b/Sources/AngouriMath/Docs/Contributing/WritingARule.md index 57f1ae9e6..1af3fdf87 100644 --- a/Sources/AngouriMath/Docs/Contributing/WritingARule.md +++ b/Sources/AngouriMath/Docs/Contributing/WritingARule.md @@ -110,7 +110,7 @@ itself, and needs no `when:` to say so. It is also strictly more specific than ` the ordering knows about — see below. **Matchable is not the same as buildable.** A pattern can be *matched* against any node type, but a -pattern used as a **replacement** must be one the library can construct: there are **62** such types. +pattern used as a **replacement** must be one the library can construct: there are **63** such types. Two families are deliberately absent, and neither is an oversight: - **binders** — `Lambda`, `Set.ConditionalSet`. They bind a variable, and `DirectChildren` renames it diff --git a/Sources/AngouriMath/Docs/Usage/Syntax.md b/Sources/AngouriMath/Docs/Usage/Syntax.md index c5f869b74..dfd3f14a5 100644 --- a/Sources/AngouriMath/Docs/Usage/Syntax.md +++ b/Sources/AngouriMath/Docs/Usage/Syntax.md @@ -340,6 +340,14 @@ argument is exact, a rational times a power of `pi` plus a rational: `polygamma( `1/6 * pi ^ 2` and `polygamma(1, 3)` is `1/6 * pi ^ 2 - 5/4`. Other special values, such as `digamma(1) = -gamma`, stay as written, since the library has no symbol for Euler's constant. +**Incomplete gamma function** — `gamma(a, z)`. With a second argument, `gamma` is the upper incomplete +gamma function `Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt`, as Mathematica's `Gamma[a, z]` and +mpmath's `gammainc(a, z)` are. `z^a` takes its principal branch, so the cut is the negative real axis, +which keeps the value from above. It is evaluated to the working precision for any complex `a` and `z`, +and `gamma(a, 0)` is `gamma(a)` where `Re a > 0`. Its derivative in `z` is `-z ^ (a - 1) * e ^ (-z)`. +A whole `a >= 1` is elementary, `gamma(3, x)` being `2 * e ^ (-x) * (1 + x + x ^ 2 / 2)`, and +`gamma(1/2, x)` is `sqrt(pi) * erfc(sqrt(x))`. + **Calculus** — `derivative(expr, var, order)`, `integral(expr, var)`, `integral(expr, var, from, to)`, `limit(expr, var, dest)`, `limitleft(...)`, `limitright(...)`; `sum(expr, var in set)` for a sum over a set; `max(expr, var in set)` and `min` for the extremum of diff --git a/Sources/AngouriMath/Functions/Continuous/Differentiation.cs b/Sources/AngouriMath/Functions/Continuous/Differentiation.cs index 8ca962ba9..f392358ce 100644 --- a/Sources/AngouriMath/Functions/Continuous/Differentiation.cs +++ b/Sources/AngouriMath/Functions/Continuous/Differentiation.cs @@ -338,6 +338,17 @@ protected override Entity InnerDifferentiate(Variable variable) => : MathS.Polygamma(Order + 1, Argument) * Argument.DerivativeOver(variable); } + partial record IncompleteGammaf + { + // Gamma(a, u)' = -u^(a - 1) e^(-u) u', the integrand at its lower limit. A parameter that varies wants + // a Meijer G function, and this declines. https://github.com/asc-community/AngouriMath/issues/2075 + /// + protected override Entity InnerDifferentiate(Variable variable) => + Parameter.ContainsNode(variable) + ? base.InnerDifferentiate(variable) + : -MathS.Pow(LowerLimit, Parameter - 1) * MathS.Pow(MathS.e, -LowerLimit) * LowerLimit.DerivativeOver(variable); + } + partial record Sif { // Si(a)' = sin(a)/a a', the integrand of its definition, and so for the other three. diff --git a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs index a43bbdbda..864df3dd2 100644 --- a/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs +++ b/Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/InvertNode.Classes.cs @@ -285,6 +285,13 @@ partial record Polygammaf null; } + // Neither has an inverse among the library's nodes. https://github.com/asc-community/AngouriMath/issues/2075 + partial record IncompleteGammaf + { + private protected override IEnumerable? InvertNode(Entity value, Entity x) => + null; + } + // None of the four has an inverse among the library's nodes, and none of their preimages // is empty: Si(x) = 1 near 1.1654, Ci(x) = 0 near 0.6165. // https://github.com/asc-community/AngouriMath/issues/1501 diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs index e0245a716..72a0609ce 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Continuous/Evaluation.Continuous.Arithmetics.Classes.cs @@ -610,6 +610,44 @@ protected override Entity InnerSimplify(bool isExact) (@this, order, argument) => ((Polygammaf)@this).New(order, argument), isExact); } + public partial record IncompleteGammaf + { + // Defined but at z = 0 where Re a <= 0, where the integral diverges at its lower limit. + // https://github.com/asc-community/AngouriMath/issues/2075 + private protected override Entity IntrinsicCondition => + !LowerLimit.EqualTo(0) | Parameter > 0; + + /// + protected override Entity InnerSimplify(bool isExact) + => ExpandOnTwoArguments(Parameter, LowerLimit, + (a, z) => (a, z) switch + { + (Complex n, Complex x) when !isExact => Number.IncompleteGamma(n, x), + (_, Real { EDecimal: var x }) when x.IsPositiveInfinity() => Integer.Zero, + // Gamma(a, 0) is Gamma(a) where the integral converges at 0. + (Real { IsPositive: true } positive, Integer { IsZero: true }) => MathS.Gamma(positive).InnerSimplified, + // A whole n >= 1 is elementary, (n - 1)! e^(-z) sum_(k < n) z^k/k! (DLMF 8.4.8), and a half is + // sqrt(pi) erfc(sqrt(z)) (DLMF 8.4.6). + (Integer { EInteger: var n }, _) when n.Sign > 0 && n.CompareTo(30) <= 0 => WholeOrder(n.ToInt32Checked(), z), + (Rational { ERational: var h }, _) when h.Numerator.Equals(PeterO.Numbers.EInteger.One) && h.Denominator.Equals(PeterO.Numbers.EInteger.FromInt32(2)) + => MathS.Sqrt(MathS.pi) * MathS.Erfc(MathS.Sqrt(z)), + _ => null + }, + (@this, a, z) => ((IncompleteGammaf)@this).New(a, z), isExact); + + private static Entity WholeOrder(int n, Entity z) + { + Entity sum = Integer.One; + Entity term = Integer.One; + for (var k = 1; k < n; k++) + { + term = term * z / k; + sum += term; + } + return (Number.Factorial(Integer.Create(n - 1)) * MathS.Pow(MathS.e, -z) * sum).InnerSimplified; + } + } + public partial record Sif { // Entire, so defined for every argument. https://github.com/asc-community/AngouriMath/issues/1501 diff --git a/Sources/AngouriMath/Functions/Output/Latex/Latex.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Output/Latex/Latex.Arithmetics.Classes.cs index 4167661f1..d083cdd65 100644 --- a/Sources/AngouriMath/Functions/Output/Latex/Latex.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/Latex/Latex.Arithmetics.Classes.cs @@ -220,6 +220,13 @@ private protected override string LatexizeNode() : $@"\psi^{{({Order.Latexize()})}}\left({Argument.Latexize()}\right)"; } + partial record IncompleteGammaf + { + /// + private protected override string LatexizeNode() + => $@"\Gamma\left({Parameter.Latexize()}, {LowerLimit.Latexize()}\right)"; + } + partial record Sif { /// diff --git a/Sources/AngouriMath/Functions/Output/ToString/ToString.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Output/ToString/ToString.Arithmetics.Classes.cs index 52aaadefa..16433e4bb 100644 --- a/Sources/AngouriMath/Functions/Output/ToString/ToString.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/ToString/ToString.Arithmetics.Classes.cs @@ -199,6 +199,14 @@ private protected override string StringizeNode() public override string ToString() => Stringize(); } + public partial record IncompleteGammaf + { + /// + private protected override string StringizeNode() => $"gamma({Parameter.Stringize()}, {LowerLimit.Stringize()})"; + /// + public override string ToString() => Stringize(); + } + public partial record Sif { /// diff --git a/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Arithmetics.Classes.cs b/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Arithmetics.Classes.cs index 6db3325af..c38bd6299 100644 --- a/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Arithmetics.Classes.cs +++ b/Sources/AngouriMath/Functions/Output/ToSympy/ToSympy.Arithmetics.Classes.cs @@ -133,6 +133,12 @@ internal override string ToSymPy() : $@"sympy.polygamma({Order.ToSymPy()}, {Argument.ToSymPy()})"; } + public partial record IncompleteGammaf + { + internal override string ToSymPy() + => $@"sympy.uppergamma({Parameter.ToSymPy()}, {LowerLimit.ToSymPy()})"; + } + public partial record Sif { internal override string ToSymPy() diff --git a/Sources/AngouriMath/Functions/Substitute.cs b/Sources/AngouriMath/Functions/Substitute.cs index 1c4bd25e7..18ffd769f 100644 --- a/Sources/AngouriMath/Functions/Substitute.cs +++ b/Sources/AngouriMath/Functions/Substitute.cs @@ -237,6 +237,13 @@ public override Entity Substitute(Entity x, Entity value) => this == x ? value : New(Order.Substitute(x, value), Argument.Substitute(x, value)); } + partial record IncompleteGammaf + { + /// + public override Entity Substitute(Entity x, Entity value) + => this == x ? value : New(Parameter.Substitute(x, value), LowerLimit.Substitute(x, value)); + } + partial record Sif { /// diff --git a/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs b/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs index b72c951df..5e64abc51 100644 --- a/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs +++ b/Sources/AngouriMath/Functions/TreeAnalyzer/Sort.Classes.cs @@ -285,6 +285,12 @@ private protected override string SortHashName(SortLevel level) => Choice(level, "polygamma_", "polygamma_", "polygamma_"); } + public partial record IncompleteGammaf + { + private protected override string SortHashName(SortLevel level) + => Choice(level, "incompletegamma_", "incompletegamma_", "incompletegamma_"); + } + public partial record Sif { private protected override string SortHashName(SortLevel level) diff --git a/Sources/MCP/AngouriMath.Mcp/Resources.cs b/Sources/MCP/AngouriMath.Mcp/Resources.cs index 812bbd271..200deeb82 100644 --- a/Sources/MCP/AngouriMath.Mcp/Resources.cs +++ b/Sources/MCP/AngouriMath.Mcp/Resources.cs @@ -95,7 +95,8 @@ that every response returns. - `ln log log2 log10 exp pow sqrt cbrt sqr abs signum sign sgn` - `floor ceil ceiling round min max gcd lcm` - `gamma factorial binomial phi prime valuation` - - special functions: `erf erfc erfi Ei li Si Ci Shi Chi`, and `digamma(z)` and `polygamma(n, z)` + - special functions: `erf erfc erfi Ei li Si Ci Shi Chi`, `digamma(z)` and `polygamma(n, z)`, and `gamma(a, z)`, the + upper incomplete gamma function - sets: `card powerset union intersection complement image preimage domain aleph iverson` - `derivative integral limit limitleft limitright sum product argmax argmin piecewise diff --git a/Sources/Tests/UnitTests/Common/PublicApi.txt b/Sources/Tests/UnitTests/Common/PublicApi.txt index b609d322a..f42b990c5 100644 --- a/Sources/Tests/UnitTests/Common/PublicApi.txt +++ b/Sources/Tests/UnitTests/Common/PublicApi.txt @@ -1212,6 +1212,28 @@ AngouriMath.Entity+Impliesf.ToString() : System.String AngouriMath.Entity+Impliesf.ctor(AngouriMath.Entity, AngouriMath.Entity) AngouriMath.Entity+Impliesf.op_Equality(AngouriMath.Entity+Impliesf, AngouriMath.Entity+Impliesf) : System.Boolean AngouriMath.Entity+Impliesf.op_Inequality(AngouriMath.Entity+Impliesf, AngouriMath.Entity+Impliesf) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.$() : AngouriMath.Entity+IncompleteGammaf +AngouriMath.Entity+IncompleteGammaf.Codomain { } : AngouriMath.Core.Domain +AngouriMath.Entity+IncompleteGammaf.Deconstruct(AngouriMath.Entity&, AngouriMath.Entity&) : System.Void +AngouriMath.Entity+IncompleteGammaf.EqualityContract { } : System.Type +AngouriMath.Entity+IncompleteGammaf.Equals(AngouriMath.Entity+Function) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.Equals(AngouriMath.Entity+IncompleteGammaf) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.Equals(System.Object) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.GetHashCode() : System.Int32 +AngouriMath.Entity+IncompleteGammaf.InitDirectChildren() : AngouriMath.Entity[] +AngouriMath.Entity+IncompleteGammaf.InnerDifferentiate(AngouriMath.Entity+Variable) : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.InnerSimplify(System.Boolean) : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.LowerLimit { } : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.NodeFirstChild { } : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.NodeSecondChild { } : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.Parameter { } : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.PrintMembers(System.Text.StringBuilder) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.Replace(System.Func) : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.Substitute(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity +AngouriMath.Entity+IncompleteGammaf.ToString() : System.String +AngouriMath.Entity+IncompleteGammaf.ctor(AngouriMath.Entity, AngouriMath.Entity) +AngouriMath.Entity+IncompleteGammaf.op_Equality(AngouriMath.Entity+IncompleteGammaf, AngouriMath.Entity+IncompleteGammaf) : System.Boolean +AngouriMath.Entity+IncompleteGammaf.op_Inequality(AngouriMath.Entity+IncompleteGammaf, AngouriMath.Entity+IncompleteGammaf) : System.Boolean AngouriMath.Entity+Integralf.$() : AngouriMath.Entity+Integralf AngouriMath.Entity+Integralf.Codomain { } : AngouriMath.Core.Domain AngouriMath.Entity+Integralf.Deconstruct(AngouriMath.Entity&, AngouriMath.Entity&, System.Nullable>&) : System.Void @@ -1869,6 +1891,7 @@ AngouriMath.Entity+Number.Gamma(AngouriMath.Entity+Number+Complex) : AngouriMath AngouriMath.Entity+Number.GetAllRoots(AngouriMath.Entity+Number+Complex, PeterO.Numbers.EInteger) : System.Collections.Generic.HashSet AngouriMath.Entity+Number.GetAllRootsOf1(PeterO.Numbers.EInteger) : System.Collections.Generic.IEnumerable AngouriMath.Entity+Number.GetHashCode() : System.Int32 +AngouriMath.Entity+Number.IncompleteGamma(AngouriMath.Entity+Number+Complex, AngouriMath.Entity+Number+Complex) : AngouriMath.Entity+Number+Complex AngouriMath.Entity+Number.InitDirectChildren() : AngouriMath.Entity[] AngouriMath.Entity+Number.InnerDifferentiate(AngouriMath.Entity+Variable) : AngouriMath.Entity AngouriMath.Entity+Number.InnerSimplify(System.Boolean) : AngouriMath.Entity @@ -3232,6 +3255,7 @@ AngouriMath.MathS.I_4 : AngouriMath.Entity+Matrix AngouriMath.MathS.IdentityMatrix(System.Int32) : AngouriMath.Entity+Matrix AngouriMath.MathS.IdentityMatrix(System.Int32, System.Int32) : AngouriMath.Entity+Matrix AngouriMath.MathS.Implication(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity +AngouriMath.MathS.IncompleteGamma(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Integral(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Integral(AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity AngouriMath.MathS.Intersection(AngouriMath.Entity, AngouriMath.Entity) : AngouriMath.Entity+Set @@ -3415,6 +3439,7 @@ class AngouriMath.Entity+Gcdf class AngouriMath.Entity+GreaterOrEqualf class AngouriMath.Entity+Greaterf class AngouriMath.Entity+Impliesf +class AngouriMath.Entity+IncompleteGammaf class AngouriMath.Entity+Integralf class AngouriMath.Entity+Iversonf class AngouriMath.Entity+Lambda diff --git a/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs new file mode 100644 index 000000000..319d40338 --- /dev/null +++ b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs @@ -0,0 +1,90 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Tests.Convenience +{ + /// + /// The upper incomplete gamma function Gamma(a, z) = int_z^oo t^(a - 1) e^(-t) dt, written gamma(a, z), + /// with the principal branch of z^a. References from mpmath's gammainc at 130 digits: by the series + /// near the origin, by the continued fraction past it, at the poles of Gamma(a), and on the cut. + /// #2075 + /// + [Trait("Area", "Convenience")] + public sealed class IncompleteGammaTest + { + private static void AssertClose(string reference, EDecimal actual) + { + var expected = EDecimal.FromString(reference); + if (expected.IsZero) + { + Assert.True(actual.IsZero, $"{actual} where it is 0"); + return; + } + var error = actual.Subtract(expected, EContext.ForPrecision(200)).Abs().Divide(expected.Abs(), EContext.ForPrecision(10)); + Assert.True(error.LessThan(EDecimal.FromString("1E-55")), $"{actual} where it is {expected}"); + } + + [Theory] + [InlineData("gamma(1/3, 2)", "0.0681364144414591361659959249099704570952527539499344393170888", "0")] + [InlineData("gamma(-1/2, 3)", "0.00677613600177021229377514699942294840656536635691110578503547", "0")] + // E_1(1), and a pole of Gamma(a) handled by the recurrence. + [InlineData("gamma(0, 1)", "0.219383934395520273677163775460121649031047293406908207577979", "0")] + [InlineData("gamma(-3, 1/2)", "1.32194260686678451978439959176024018549333450108404397738537", "0")] + // By the continued fraction, and far below the floor at zero. + [InlineData("gamma(5/2, 40)", "1.11555920556816836584315244570307218006275734693311260218414E-15", "0")] + [InlineData("gamma(7, 300)", "3.82936658212129710001498193653594875603891571118487941125183E-116", "0")] + // On the cut, from above. + [InlineData("gamma(1/2, -3)", "1.77245385090551602729816748334114518279754945612238712821381", "-14.6261713840190916850829290647365438829669286131782239586917")] + [InlineData("gamma(-1, -20)", "1.35739289356707492202577927702104389515074891049797988869358E+6", "3.14159265358979323846264338327950288419716939937510582097494")] + [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) + { + var value = Assert.IsAssignableFrom(expression.EvalNumerical()); + AssertClose(real, value.RealPart.EDecimal); + AssertClose(imaginary, value.ImaginaryPart.EDecimal); + } + + [Theory] + [InlineData("gamma(1, x)", "e ^ (-x)")] + [InlineData("gamma(3, x)", "2 * e ^ (-x) * (1 + x + x ^ 2 / 2)")] + [InlineData("gamma(1/2, x)", "sqrt(pi) * erfc(sqrt(x))")] + [InlineData("gamma(5/2, 0)", "gamma(5/2)")] + public void IsElementaryWhereItCanBe(string expression, string elementary) + => Assert.Equal(elementary.ToEntity().InnerSimplified, expression.ToEntity().InnerSimplified); + + [Fact] + public void TheDerivativeInTheLowerLimitIsTheIntegrandThere() + { + var derivative = "gamma(1/3, x^2)".ToEntity().Differentiate("x"); + var at = derivative.Substitute("x", Rational.Create(3, 2)).EvalNumerical(); + var expected = "-(x^2)^(-2/3) * e^(-x^2) * 2 * x".ToEntity().Substitute("x", Rational.Create(3, 2)).EvalNumerical(); + AssertClose(expected.Stringize(), Assert.IsType(at).EDecimal); + } + + [Fact] + public void AVaryingParameterHasNoDerivativeHere() + => Assert.IsType("gamma(x, 2)".ToEntity().Differentiate("x").InnerSimplified); + + [Fact] + public void IsParsedAndPrinted() + { + var entity = "gamma(1/3, x)".ToEntity(); + Assert.IsType(entity); + Assert.Equal("gamma(1/3, x)", entity.Stringize()); + Assert.Equal(@"\Gamma\left(\frac{1}{3}, x\right)", entity.Latexize()); + Assert.Contains("sympy.uppergamma(sympy.Rational(1, 3), x)", MathS.ToSympyCode(entity)); + Assert.Equal(entity, entity.Stringize().ToEntity()); + Assert.IsType("gamma(x)".ToEntity()); + } + } +} diff --git a/Sources/Tests/UnitTests/Core/Transformations/RuleAuthoringGuideTest.cs b/Sources/Tests/UnitTests/Core/Transformations/RuleAuthoringGuideTest.cs index 31b3b67e7..ae4db26ef 100644 --- a/Sources/Tests/UnitTests/Core/Transformations/RuleAuthoringGuideTest.cs +++ b/Sources/Tests/UnitTests/Core/Transformations/RuleAuthoringGuideTest.cs @@ -79,7 +79,7 @@ public void TheIdentityIsNotTheName() [Fact] public void ThePatternLanguage() - => Stated(62, MatchPattern.BuildableNodeTypes.Count, "the number of buildable node types"); + => Stated(63, MatchPattern.BuildableNodeTypes.Count, "the number of buildable node types"); [Fact] public void ReplacementAPatternOrCode()