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()