diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 2cad5a998..3064be998 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -134,11 +134,12 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa **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. +- **Branch.** It takes the branch of `z^a` the power takes: the principal one, with the cut on the negative real axis, + and on that axis the real root for an exact `a` of odd denominator below 10, as `(-8)^(1/3)` is `-2` + ([#204](https://github.com/asc-community/AngouriMath/issues/204)). So `gamma(1/3, -2)` is real, where Mathematica's + `Gamma[1/3, -2]` is not, and the derivative below is the derivative of the value there too. - **Evaluation.** It is evaluated to the working precision for any complex `a` and `z`. -- **Derivative.** In `z` it is `-z^(a - 1) e^(-z)`. At a negative `z` that power takes the real root for an exact - `a` of small odd denominator, so there it is not the derivative of the value - ([#204](https://github.com/asc-community/AngouriMath/issues/204)). +- **Derivative.** In `z` it is `-z^(a - 1) e^(-z)`. - **Exact forms.** A whole `a >= 1` and `a = 1/2` are elementary. - **Integrals.** Of a linear `u`, beside a polynomial, or beside any power of `x` but `1/x` where `u` is a multiple of `x`, by parts, in gamma functions of `u`. @@ -149,6 +150,7 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `"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)` | +| `"gamma(1/3, -2)".EvalNumerical()` | a `FunctionArgumentCountException` | `10.212552064243320672…`, real | | `"x^2*gamma(0, a*x)".ToEntity().Integrate("x")` | a `FunctionArgumentCountException` | `x ^ 3 / 3 * gamma(0, a * x) - 1 / (3 * a ^ 3) * 2 * e ^ (-a * x) * (1 + a * x + (a * x) ^ 2 / 2) + C` | | `"x^m*gamma(n, b*x)".ToEntity().Integrate("x")` | a `FunctionArgumentCountException` | `(x ^ (m + 1) * gamma(n, b * x) - x ^ m * (b * x) ^ (-m) * gamma(n + m + 1, b * x) / b) / (m + 1) + C` | diff --git a/Sources/AngouriMath/Convenience/MathS.cs b/Sources/AngouriMath/Convenience/MathS.cs index c588e75ea..732f7adf6 100644 --- a/Sources/AngouriMath/Convenience/MathS.cs +++ b/Sources/AngouriMath/Convenience/MathS.cs @@ -1479,11 +1479,10 @@ public static Entity VectorLength(Matrix v) /// 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). At a negative - /// z the power there takes the real root for an exact a of small odd denominator, not the principal - /// branch, so for such an a it is not the derivative of the value - /// (https://github.com/asc-community/AngouriMath/issues/204). It is written + /// z^a is the power's: its principal branch, so the cut is the negative real axis, which keeps the value from + /// above, and on that axis the real root where the power takes it, for an exact a of odd denominator below 10, + /// where the value is real (https://github.com/asc-community/AngouriMath/issues/204). + /// 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 /// diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs index 2be14ef1e..a4f4511c4 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/IncompleteGamma.cs @@ -16,8 +16,9 @@ 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. + /// precision, for complex a and z, with the branch of z^a the power takes: the principal + /// one, whose cut is the negative real axis and which keeps the value from above, and on that axis the real root + /// where the power takes it, for an exact a of odd denominator below 10. /// /// /// @@ -65,12 +66,30 @@ public static Complex IncompleteGamma(Complex a, Complex z) // 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)); + // Gamma(a, z) is Gamma(a) - z^a g(a, z) with g entire, and its z^a is the power's: on the negative axis, for + // an exact a of odd denominator q below 10, the power takes the real root, rho = (-1)^p e^(-i pi p/q) times + // the principal one. With the principal branch there, -z^(a - 1) e^(-z), whose power takes the same root, was + // not the derivative of the value, and an antiderivative written in the function had the values of one + // turned by rho: the integral of e^(x^3) from 1 to 2 came out -137.8 + 238.6 i, where it is 275.5. + // SimplificationContract.md section 6; https://github.com/asc-community/AngouriMath/issues/204 + var turn = TheRealRootsTurn(a, zr, zi); 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)) + : BySeries((ar, ai), (zr, zi), modulus, context, turn); + // 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) + { + var work = Working(context, 10); + var rho = Turned(fractionsTurn, work); + var gamma = GammaOf((ar, ai), work); + var turned = Multiply(rho, (re, im), work); + (re, im) = (gamma.Re.Subtract(gamma.Re.Multiply(rho.Re, work), work).Add(turned.Re, work), turned.Im.Subtract(gamma.Re.Multiply(rho.Im, work), work)); + } + // Real where it is real, exactly: for a real a and a positive z, for a positive whole a and any real z, + // where Gamma(n, z) is (n - 1)! e^(-z) times a polynomial and z^a leaves a residue of e^(i pi n), and where + // z^a is the real root, each factor of Gamma(a) - z^a g(a, z) being real. + if (ai.IsZero && zi.IsZero && (zr.Sign > 0 || ar.IsInteger() && ar.Sign > 0 || turn is not null)) 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. @@ -104,7 +123,31 @@ private static (EDecimal Re, EDecimal Im) ByContinuedFraction((EDecimal Re, EDec 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) + /// + /// Where the power takes the real root of z^a, the turn from the principal value to it, e^(i pi k/q) + /// as (k, q): an exact of odd denominator q below 10 and a negative real + /// z, the case searches for a real root in. elsewhere. + /// + private static (EInteger Turns, EInteger Q)? TheRealRootsTurn(Complex a, EDecimal zr, EDecimal zi) + { + if (a is not Rational { ERational: var exponent } || a is Integer || !zi.IsZero || zr.Sign >= 0) + return null; + var q = exponent.Denominator.Abs(); + if (q.IsEven || q.CompareTo(EInteger.FromInt32(10)) >= 0) + return null; + // (-1)^p e^(-i pi p/q) = e^(i pi p (q - 1)/q), the angle taken modulo 2 pi. + var p = exponent.Denominator.Sign < 0 ? exponent.Numerator.Negate() : exponent.Numerator; + return (p.Multiply(q.Subtract(EInteger.One)).Mod(q.Multiply(EInteger.FromInt32(2))), q); + } + + private static (EDecimal Re, EDecimal Im) Turned((EInteger Turns, EInteger Q) turn, EContext work) + { + var angle = InternalAMExtensions.ConstantCache.Lookup(work).Pi.Multiply(EDecimal.FromEInteger(turn.Turns), work).Divide(EDecimal.FromEInteger(turn.Q), work); + return (angle.Cos(work), angle.Sin(work)); + } + + 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; for (var attempt = 0; ; attempt++) @@ -112,7 +155,7 @@ private static (EDecimal Re, EDecimal Im) BySeries((EDecimal Re, EDecimal Im) a, 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); + : AsGammaLessTheLower(a, z, 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; @@ -125,7 +168,8 @@ 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. /// - private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessTheLower((EDecimal Re, EDecimal Im) a, (EDecimal Re, EDecimal Im) z, EContext work) + 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) { var epsilon = EDecimal.Create(EInteger.One, EInteger.FromInt32(-work.Precision.ToInt32Checked())); // term_k = z^k / (a (a + 1) ... (a + k)), starting at 1/a @@ -170,7 +214,10 @@ private static ((EDecimal Re, EDecimal Im) Value, EDecimal Scale) AsGammaLessThe break; } } - var prefactor = Multiply(ComplexPower(z, a, work), ComplexExp((z.Re.Negate(), z.Im.Negate()), work), work); + 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 lower = Multiply(prefactor, sum, work); var gamma = GammaOf(a, work); var scale = EDecimal.Max(Magnitude(gamma), Magnitude(prefactor).Multiply(peak, work)); diff --git a/Sources/AngouriMath/Docs/Contributing/SimplificationContract.md b/Sources/AngouriMath/Docs/Contributing/SimplificationContract.md index ce3a5d071..8da3532ba 100644 --- a/Sources/AngouriMath/Docs/Contributing/SimplificationContract.md +++ b/Sources/AngouriMath/Docs/Contributing/SimplificationContract.md @@ -1,4 +1,4 @@ -# What a simplification rule may assume +# What a simplification rule may assume A rewrite rule in this library is a claim: *this expression and that one are two ways of writing the same thing.* This file says what has to be true for that claim to hold, so that a new rule can be @@ -256,6 +256,8 @@ release touches a branch. A rewrite that changes any of them changes answers. | `(-8)^(1/5)` | `-1.5157...` | the real root | | `(-8)^(1/4)` | `1.1892 + 1.1892i` | principal — there is no real fourth root | | `(-8)^0.3333333333` | `1.0 + 1.7320i` | **principal**, not the real root | +| `(-8)^(1/11)` | `1.1592 + 0.3404i` | principal: the real root is searched for below a denominator of 10 | +| `gamma(1/3, -2)` | `10.2125...` | the power's real root, which `z^a` in `Gamma(a) - z^a g(a, z)` is | Read the last three rows together, because they are a trap. For a negative base and an exact `Rational` exponent with an **odd denominator**, this library takes the real root. Write the same diff --git a/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs index 319d40338..e99269e85 100644 --- a/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs +++ b/Sources/Tests/UnitTests/Convenience/IncompleteGammaTest.cs @@ -54,6 +54,40 @@ public void IsItsValue(string expression, string real, string imaginary) AssertClose(imaginary, value.ImaginaryPart.EDecimal); } + /// + /// On the cut, for an exact order of odd denominator below 10, with the real root the power takes there: + /// Gamma(a) - z^a g(a, z) with g(a, z) = sum_k (-z)^k/(k! (a + k)) (DLMF 8.7.3) and z^a real. The + /// references are that sum at 90 digits; the last is past the continued fraction's threshold. + /// + [Theory] + [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, -300)", "9.71532351293767903599760645811774844613945799010165863950178E+126")] + public void OnTheCutItTakesThePowersRealRoot(string expression, string real) + { + var value = Assert.IsAssignableFrom(expression.EvalNumerical()); + AssertClose(real, value.RealPart.EDecimal); + Assert.True(value.ImaginaryPart.EDecimal.IsZero, $"{value}"); + } + + /// + /// And so its derivative, whose power takes the same root, is the derivative of its value there: with the principal + /// branch it was the derivative turned by e^(-2 i pi/3). + /// + [Fact] + public void OnTheCutTheDerivativeIsTheDerivativeOfTheValue() + { + var at = Rational.Create(-2, 1); + var h = EDecimal.FromString("1E-30"); + var above = "gamma(1/3, z)".ToEntity().Substitute("z", Real.Create(EDecimal.FromInt32(-2).Add(h))).EvalNumerical(); + var below = "gamma(1/3, z)".ToEntity().Substitute("z", Real.Create(EDecimal.FromInt32(-2).Subtract(h))).EvalNumerical(); + var difference = Assert.IsType(above - below).EDecimal.Divide(h.Multiply(2), EContext.ForPrecision(100)); + var derivative = Assert.IsType("gamma(1/3, z)".ToEntity().Differentiate("z").Substitute("z", at).EvalNumerical()).EDecimal; + var error = difference.Subtract(derivative, EContext.ForPrecision(100)).Abs().Divide(derivative.Abs(), EContext.ForPrecision(10)); + Assert.True(error.LessThan(EDecimal.FromString("1E-25")), $"the value changes at {difference} where the derivative is {derivative}"); + } + [Theory] [InlineData("gamma(1, x)", "e ^ (-x)")] [InlineData("gamma(3, x)", "2 * e ^ (-x) * (1 + x + x ^ 2 / 2)")]