From 80aad7f7e3f5868e9732a60cbfbab8828e1f5854 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 02:03:53 +0000 Subject: [PATCH] Integrate a power beside an exponential of a power as an incomplete gamma function int x^m F^(a + b x^n) dx is -F^a x^(m + 1) z^(-s) gamma(s, z)/n with z = -b ln(F) x^n and s = (m + 1)/n, and int (c + d x)^m F^(a + b x) dx the same with n = 1 under t = c + d x: x^(m + 1) z^(-s) is constant between the cuts, and gamma(s, z) differentiates to -z^(s - 1) e^(-z) z'. Only where s is neither whole nor half a whole number, whose integrals are elementary, exponential integrals or error functions and keep those answers, or where it is past 30, where parts would take s steps: x^100 e^(-a x) is -a^(-101) gamma(101, a x). Where the constant factor is one number on both half-lines it is written so, and the answer has a value at 0: e^(x^3) is gamma(1/3, -x^3)/3. Two costs this exposed, each fixed here. By parts integrated (a + b ln(c x^n))^p into a gamma function of a logarithm and handed its remainder on, which nothing integrates and the searches sampled thousands of times: three by-parts paths now decline a gamma function of an argument that is not linear. And a definite integral of an answer undefined at a bound, x^2 e^(x^4) from 0, took a limit through l'Hopital's rule for over a minute; it is left as written there. gamma's own limit is taken by continuity where its value is finite. The order is read through Simplify where its symbols cancel: (-1 + n + 1)/n is 1, and (p x)^(n - 1) e^(x^n) keeps its elementary answer. Rubi's suite on the harness that reads GAMMA, against master e62a56d6: 48,502 -> 48,634 of 49,052, 0 wrong, 0 rows lost, timeouts 5 -> 2; the gains are 68 in 2.3, 24 in 3.1.2, 12 in 8.6 and the logarithm files. Two rows' allocation grew in the run and is identical alone on both builds. Four tests that pinned a decline answer now, each checked by value. Part of #2075; answers #180's item 22. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 8 + .../Integration/IndefiniteIntegralSolver.cs | 31 ++- .../Integration/IntegralPatterns.cs | 206 ++++++++++++++++++ .../Integration/Integration.Definition.cs | 27 ++- .../Limits/Solvers/Limit.Classes.cs | 19 ++ ...BesideAnExponentialOfAPowerIntegralTest.cs | 108 +++++++++ .../APowerTimesAFunctionOfItsPowerTest.cs | 13 +- .../Calculus/BreakpointIntegrationTest.cs | 16 +- .../Calculus/ExponentialAnsatzTest.cs | 5 +- .../Calculus/ProductOfPowersDerivativeTest.cs | 21 +- 10 files changed, 439 insertions(+), 15 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 777cb0d88..3d83a5225 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -145,6 +145,10 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa of `x`, by parts, in gamma functions of `u`. And of a whole `a <= 0` over a power of another linear, by parts, where what is left is an exponential over linears. +- **Integrals in it.** `x^m F^(a + b x^n)` and `(c + d x)^m F^(a + b x)` are `-F^a x^(m + 1) z^(-s) gamma(s, z)/n`, with + `z = -b ln(F) x^n` and `s = (m + 1)/n`, where `s` is neither whole nor half a whole number, or is past 30, as + `x^31 e^(-x)` is, whose sum of 32 terms is `-gamma(32, x)` now. A definite integral of such an answer that is + undefined at a bound is left as written. ([#2075](https://github.com/asc-community/AngouriMath/issues/2075)) | Input | Was (2.6.0) | Now | @@ -156,6 +160,10 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `"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` | | `"gamma(0, a + b*x)/(c + d*x)^2".ToEntity().Integrate("x")` | a `FunctionArgumentCountException` | `-gamma(0, a + b x)/(d (c + d x))` and two exponential integrals, of `a + b x` and of `b (c + d x)/d` | +| `"x^(1/3)*e^(-x)".ToEntity().Integrate("x")` | `integral(x ^ (1/3) * e ^ (-x), x)` | `-gamma(4/3, x) + C` | +| `"(c + d*x)^m*e^(a + b*x)".ToEntity().Integrate("x")` | `integral((c + d * x) ^ m * e ^ (a + b * x), x)` | `e ^ (a - b * c / d) * (c + d * x) ^ m * (-b * (c + d * x) / d) ^ (-m) * gamma(m + 1, -b * (c + d * x) / d) / b + C` | +| `"x^31*e^(-x)".ToEntity().Integrate("x")` | a sum of 32 terms, 6,172 characters | `-gamma(32, x) + C` | +| `"integral(e^(x^3), x, 0, 1)".ToEntity().Simplify()` | `integral(e ^ x ^ 3, x, 0, 1)` | `(gamma(1/3, -1) - (-2/3)!) / 3`, which is `1.3419044179774…` | ### Numerical equality is to half the evaluation's digits diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 6a8256eee..4b4ca4295 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -260,8 +260,11 @@ private static bool OnlyInsideTranscendentalFunctions(Entity expr, Entity.Variab // And one with a polynomial in x beside something in its antiderivative: what parts // leaves then has the power it took away, and asking it does not descend -- // `x erfc(a + b x)^2`, whose cofactor's antiderivative has `x erfc(a + b x)^2` in it. + // Nor one whose antiderivative is an incomplete gamma function of an argument that is not linear, which + // the remainder would hold: see TryIntegrateByPartsOnce. `(g + h x)^2 (a + b ln(c (d (e + f x)^p)^q))^n` + // went from a decline in three seconds to one in two minutes. if (Integration.ComputeAsAQuestionOfItsOwn(rest, x, integrateByParts: true) is not { } antiderivative - || HasAPolynomialFactorIn(antiderivative, x)) + || HasAPolynomialFactorIn(antiderivative, x) || IntroducesAnIncompleteGammaFunctionOfANonlinear(antiderivative, rest, x)) { remainder.Add(term); continue; @@ -5166,6 +5169,9 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO ?? Integration.ComputeIndefiniteIntegral(toIntegrate, x, false); if (integral is null) return null; var differential = polynomialToDifferentiate.Differentiate(x); + // And the next round integrates what this one returned, which nothing does where that is an incomplete + // gamma function of an argument that is not linear: see TryIntegrateByPartsOnce. + if (!TreeAnalyzer.IsZero(differential) && IntroducesAnIncompleteGammaFunctionOfANonlinear(integral, toIntegrate, x)) return null; var result = IntegrateByPartsPolynomial(differential, integral, x, currentRecursion + 1); return (result is null) ? null : polynomialToDifferentiate * integral - result; } @@ -5194,6 +5200,12 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO byThePolynomialsParts = integralOfU is not null; } if (integralOfU is null) return null; + // Not where the integral of u is an incomplete gamma function of an argument that is not linear, which u + // did not hold: the remainder holds it then, and nothing here integrates it. (a + b ln(c x^n))^p beside + // d + e ln(f x^r) integrates the power to gamma(p + 1, -(a + b ln(c x^n))/(b n)), and the search on its + // remainder sampled it three thousand times before declining, thirty-five seconds where it had been one. + // https://github.com/asc-community/AngouriMath/issues/2075 + if (IntroducesAnIncompleteGammaFunctionOfANonlinear(integralOfU, u, x)) return null; // Differentiate v var derivativeOfV = v.Differentiate(x).InnerSimplified; @@ -50047,6 +50059,23 @@ Entity Term(Entity? sum, Entity coefficient, int power) /// [System.ThreadStatic] private static bool askingTheTermsOfARemainder; + /// + /// Whether holds an incomplete gamma function of an argument not linear in + /// where holds none. + /// + private static bool IntroducesAnIncompleteGammaFunctionOfANonlinear(Entity integral, Entity integrand, Variable x) + { + foreach (var node in integral.Nodes) + if (node is IncompleteGammaf { LowerLimit: var argument } && !TreeAnalyzer.TryGetPolyLinear(argument, x, out _, out _)) + { + foreach (var written in integrand.Nodes) + if (written is IncompleteGammaf) + return false; + return true; + } + return false; + } + private static bool IsASpecialFunction(Entity node) => node is Entity.Erff or Entity.Erfcf or Entity.Erfif or Entity.Eif or Entity.Lif or Entity.Sif or Entity.Cif or Entity.Shif or Entity.Chif; diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index bf99e98a6..25dfb935c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -347,6 +347,205 @@ private static bool TakeMomentFactor(Entity factor, bool inverted, Entity.Variab return coefficient * (integrated * gamma - quotient * shifted / slope) / mPlusOne; } + /// + /// int x^m F^(a + b x^n) dx, and int (c + d x)^m F^(a + b x) dx, where the integral is a gamma + /// function of s = (m + 1)/n, for an s neither whole nor half a whole number -- where it is + /// elementary, an exponential integral or an error function, and those are answered as such -- or not a + /// number at all. With z = -b ln(F) x^n, x^(m + 1) z^(-s) is constant between the cuts and + /// Gamma(s, z) differentiates to -z^(s - 1) e^(-z) z', so the first is + /// -F^a x^(m + 1) z^(-s) Gamma(s, z)/n; the second is the first with n = 1 under + /// t = c + d x. Rubi's 2.3. https://github.com/asc-community/AngouriMath/issues/2075 + /// + private static Entity? APowerBesideAnExponentialOfAPower(Entity expr, Entity.Variable x) + { + // Asked of every integrand that reaches the table, so a type test comes first, and + // nothing here captures a local. + if (expr is not (Entity.Mulf or Entity.Divf or Entity.Powf)) + return null; + var (numerator, denominator) = expr is Entity.Divf(var n, var d) ? (n, d) : (expr, (Entity)Entity.Number.Integer.One); + Entity.Powf? exponential = null; + var exponentialBelow = false; + foreach (var factor in Entity.Mulf.LinearChildren(numerator)) + if (factor is Entity.Powf(var exponentialBase, var power) candidate && !exponentialBase.ContainsNode(x) && power.ContainsNode(x)) + { + if (exponential is not null) + return null; + exponential = candidate; + } + foreach (var factor in Entity.Mulf.LinearChildren(denominator)) + if (factor is Entity.Powf(var exponentialBase, var power) candidate && !exponentialBase.ContainsNode(x) && power.ContainsNode(x)) + { + if (exponential is not null) + return null; + exponential = candidate; + exponentialBelow = true; + } + if (exponential is null) + return null; + + // What is beside it: constants, and powers of multiples of x, (d x)^m, whose product P has + // P'/P = M/x with M the sum of their exponents -- or one power of a linear with an offset. + Entity coefficient = Entity.Number.Integer.One; + Entity? written = null; + Entity exponents = Entity.Number.Integer.Zero; + var factorsInX = 0; + var allOfX = true; + var allMultiples = true; + Entity? single = null; + var singleBelow = false; + Entity singleBase = x; + Entity singleScale = Entity.Number.Integer.One; + for (var side = 0; side < 2; side++) + foreach (var factor in Entity.Mulf.LinearChildren(side == 0 ? numerator : denominator)) + { + if (ReferenceEquals(factor, exponential)) + continue; + if (!factor.ContainsNode(x)) + { + coefficient = side == 0 ? coefficient * factor : coefficient / factor; + continue; + } + factorsInX++; + single = factor; + singleBelow = side == 1; + if (TryReadAPowerOfAMultiple(factor, x, out var @base, out var scale, out var m)) + { + written = written is null ? (side == 0 ? factor : 1 / factor) : side == 0 ? written * factor : written / factor; + exponents = side == 0 ? exponents + m : exponents - m; + allOfX &= @base == x; + singleBase = @base; + singleScale = scale; + } + else + allMultiples = false; + } + var exponent = exponentialBelow ? -exponential.Exponent : exponential.Exponent; + var logarithm = exponential.Base == MathS.e ? Entity.Number.Integer.One : MathS.Ln(exponential.Base); + + // The exponent as a + b x^n, the term in x as written. + Entity constant = Entity.Number.Integer.Zero; + Entity? term = null; + Entity? degree = null; + foreach (var summand in Entity.Sumf.LinearChildren(exponent)) + if (!summand.ContainsNode(x)) + constant += summand; + else if (term is null && ThePowerOfX(summand, x) is { } read) + { + term = summand; + degree = read; + } + else + return null; + if (term is null || degree is null || TreeAnalyzer.IsZero(degree)) + return null; + + if (allMultiples) + { + var mPlusOne = (exponents + 1).InnerSimplified; + var s = TheOrder(mPlusOne / degree); + if (WholeOrHalfWhole(s)) + return null; + // z = k x for a linear exponent, with the rate read as the number it may be. + var linearRate = degree.Equals(Entity.Number.Integer.One) && TreeAnalyzer.TryGetPolyLinear(term, x, out var rateOfTerm, out _) + ? (-logarithm * rateOfTerm).InnerSimplified + : null; + var z = linearRate is null ? (-logarithm * term).InnerSimplified : (linearRate * x).InnerSimplified; + var leading = (-coefficient / degree).InnerSimplified * Exponential(exponential.Base, constant); + // For a positive k, x^(m + 1) z^(-s) is k^(-s) on every branch, and for a whole s with any k. + if (linearRate is not null && allOfX && (linearRate is Entity.Number.Real { IsPositive: true } || s is Entity.Number.Integer)) + return leading * MathS.Pow(linearRate, (-s).InnerSimplified).InnerSimplified * MathS.IncompleteGamma(s, z); + // P x: x^(M + 1) for powers of x, and (d x)^(m + 1)/d for one power of a multiple. + var integrated = allOfX ? Power(x, mPlusOne) + : factorsInX == 1 ? MathS.Pow(singleBase, mPlusOne) / singleScale + : written! * x; + var constantFactor = integrated * MathS.Pow(z, (-s).InnerSimplified); + // x^(m + 1) z^(-s), constant on each half-line, is K^(-s) for x > 0 with z = K x^n, and its value at -1 for + // x < 0. Where the two are one number it is that number on the whole line, and the answer has a value at 0 + // rather than 0 times an infinity there: e^(x^3) is gamma(1/3, -x^3)/3, and its integral from 0 to 1 was NaN. + if (allOfX && mPlusOne is Entity.Number && degree is Entity.Number + && (-logarithm * term.Substitute(x, Entity.Number.Integer.One)).InnerSimplified is Entity.Number theRate + && MathS.Pow(theRate, (-s).InnerSimplified).InnerSimplified is Entity.Number onTheRight + && constantFactor.Substitute(x, Entity.Number.Integer.MinusOne).InnerSimplified is Entity.Number onTheLeft + && onTheRight == onTheLeft) + return leading * onTheRight * MathS.IncompleteGamma(s, z); + return leading * constantFactor * MathS.IncompleteGamma(s, z); + } + + // A power of a linear with an offset, beside an exponential of a linear. + if (factorsInX != 1 || !degree.Equals(Entity.Number.Integer.One) + || single is not Entity.Powf(var linear, var exponentOfLinear) || exponentOfLinear.ContainsNode(x) + || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var offset) + || !TreeAnalyzer.TryGetPolyLinear(term, x, out var rate, out _)) + return null; + var exponentOfPower = singleBelow ? (-exponentOfLinear).InnerSimplified : exponentOfLinear; + if (WholeOrHalfWhole(TheOrder(exponentOfPower + 1))) + return null; + var under = (-logarithm * rate * linear / slope).InnerSimplified; + return coefficient * Exponential(exponential.Base, (constant - rate * offset / slope).InnerSimplified) * (singleBelow ? 1 / single : single) + * MathS.Pow(under, -exponentOfPower) * MathS.IncompleteGamma((exponentOfPower + 1).InnerSimplified, under) / (logarithm * rate); + } + + /// + /// The gamma function's order , as the number it is where its symbols cancel: (-1 + n + 1)/n + /// is 1, which inner simplification leaves as written, and (p x)^(n - 1) e^(x^n), whose integral is + /// e^(x^n)/n, was answered as gamma((-1 + n + 1)/n, -x^n). + /// + private static Entity TheOrder(Entity s) + { + var inner = s.InnerSimplified; + if (inner is Entity.Number || !inner.Vars.Any()) + return inner; + // 1 provided not n = 0, as the integrator answers the generic case. + return inner.Simplify() switch + { + Entity.Number number => number, + Entity.Providedf(Entity.Number number, _) => number, + _ => inner + }; + } + + /// + /// Whether is a number that is whole or half of a whole number, and no more than 30 in size: + /// where the integral is elementary, an exponential integral or an error function, and is answered as such. Past + /// that it is the gamma function, which stays one, where parts or the Gaussian's moments would take as many + /// steps as s and write as many terms: x^100 e^(-a x) is -Gamma(101, a x)/a^101. + /// + private static bool WholeOrHalfWhole(Entity s) + => s is Entity.Number.Rational { ERational: var exact } + ? exact.Denominator.CompareTo(EInteger.FromInt32(2)) <= 0 && exact.Abs().CompareTo(ERational.FromInt32(30)) <= 0 + : s is Entity.Number.Real { EDecimal: var approximate } && approximate.IsFinite && approximate.Multiply(2).IsInteger() + && approximate.Abs().CompareTo(EDecimal.FromInt32(30)) <= 0; + + /// to , without a zeroth power. + private static Entity Exponential(Entity @base, Entity power) + => TreeAnalyzer.IsZero(power) ? Entity.Number.Integer.One : MathS.Pow(@base, power); + + /// + /// The n of as b x^n, with b and n free of + /// , or : each factor in x is x, a power of it, or a whole + /// power of one, which is how a quotient's divisor is read. + /// + private static Entity? ThePowerOfX(Entity term, Entity.Variable x) + { + Entity? degree = null; + foreach (var factor in Entity.Mulf.LinearChildren(term)) + { + if (!factor.ContainsNode(x)) + continue; + Entity? read = factor == x ? Entity.Number.Integer.One + : factor is Entity.Powf(var @base, var power) && !power.ContainsNode(x) + ? (@base == x ? power + : @base is Entity.Powf(var inner, var innerPower) && inner == x && !innerPower.ContainsNode(x) && power is Entity.Number.Integer + ? innerPower * power + : null) + : null; + if (read is null) + return null; + degree = degree is null ? read : degree + read; + } + return degree?.InnerSimplified; + } + /// /// 's coefficients by power of , with /// read as 1, or where it is not a polynomial @@ -429,6 +628,9 @@ private static Entity BesideAPolynomialWithoutAnOffset(Dictionary MathS.Pow(@base, power) }; + private static Entity Power(Entity @base, Entity power) + => power is Entity.Number.Integer { EInteger: var whole } && whole.CanFitInInt32() ? Power(@base, whole.ToInt32Checked()) : MathS.Pow(@base, power); + /// /// read as (d x)^m, with d and m free of /// : x itself, a power of it, or a power of a multiple of it. @@ -584,6 +786,10 @@ _ when GaussianMoment(expr, x) is { } moment => moment, // polynomial. https://github.com/asc-community/AngouriMath/issues/2075 _ when IncompleteGammaBesideAPower(expr, x) is { } gamma => gamma, + // And a power beside an exponential of a power, where that is a gamma function and + // neither elementary nor an error function. https://github.com/asc-community/AngouriMath/issues/2075 + _ when APowerBesideAnExponentialOfAPower(expr, x) is { } gammaForm => gammaForm, + Entity.Absf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => // ∫ |ax + b| dx = sgn(ax + b) * (ax + b)^2 / (2a) MathS.Signum(arg) * MathS.Pow(arg, 2) / (2 * a), diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 62ca25bb0..5330f7a75 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -80,10 +80,15 @@ public Entity Integrate(Variable x, Entity from, Entity to) => || (lower < upper ? AntiderivativeBreaks.Inside(antiderivative, x, lower, upper) : AntiderivativeBreaks.Inside(antiderivative, x, upper, lower)) is not { Count: > 0 } breaks) - return AtTheBound(antiderivative, x, to, from) - AtTheBound(antiderivative, x, from, to); + return AtTheBound(antiderivative, x, to, from) is { } atTheTop && AtTheBound(antiderivative, x, from, to) is { } atTheBottom + ? atTheTop - atTheBottom + : null; if (lower > upper) breaks.Reverse(); - var total = -AtTheBound(antiderivative, x, from, to); + // Not across a break of an antiderivative in the incomplete gamma function: see AtTheBound. + if (HoldsAnIncompleteGammaFunction(antiderivative) || AtTheBound(antiderivative, x, from, to) is not { } atTheStart) + return null; + var total = -atTheStart; foreach (var (at, _) in breaks) { // Approached from the side the range comes from, then left on the other. @@ -94,7 +99,15 @@ public Entity Integrate(Variable x, Entity from, Entity to) => return null; total += before - after; } - return total + AtTheBound(antiderivative, x, to, from); + return AtTheBound(antiderivative, x, to, from) is { } atTheEnd ? total + atTheEnd : null; + } + + private static bool HoldsAnIncompleteGammaFunction(Entity expr) + { + foreach (var node in expr.Nodes) + if (node is IncompleteGammaf) + return true; + return false; } /// @@ -112,7 +125,7 @@ public Entity Integrate(Variable x, Entity from, Entity to) => /// decided, the value stays undefined, as it is for sin(x) from 0 to +oo, /// whose antiderivative has no limit there. /// - private static Entity AtTheBound(Entity antiderivative, Variable x, Entity bound, Entity otherBound) + private static Entity? AtTheBound(Entity antiderivative, Variable x, Entity bound, Entity otherBound) { var value = antiderivative.Substitute(x, bound); if (!value.InnerSimplified.IsNaN) @@ -128,6 +141,12 @@ private static Entity AtTheBound(Entity antiderivative, Variable x, Entity bound }; if (side is not { } inside) return value; + // An antiderivative in the incomplete gamma function undefined at a bound is left as written rather than + // taken to its limit there: x^3 (-x^4)^(-3/4) gamma(3/4, -x^4), the integral of x^2 e^(x^4), is 0 times an + // infinity at 0, and its limit went through l'Hopital's rule on the whole product, differentiating the function + // and evaluating Gamma(3/4) at every step, for over a minute. https://github.com/asc-community/AngouriMath/issues/2075 + if (HoldsAnIncompleteGammaFunction(antiderivative)) + return null; var limit = antiderivative.Limit(x, bound, inside); return limit.IsNaN || limit.Nodes.Any(node => node is Limitf) ? value : limit; } diff --git a/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Classes.cs b/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Classes.cs index c5a0ac6f5..d4dba0e85 100644 --- a/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Classes.cs +++ b/Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Limit.Classes.cs @@ -347,6 +347,25 @@ partial record Factorialf x, dist, side); } + // The upper incomplete gamma function, continuous in z away from 0, and at 0 for Re a > 0, where it is Gamma(a): the + // limit is the function at its argument's limit wherever that value is a finite number. An infinite one is not taken, + // since it can differ by side: gamma(-1, z) is about 1/z. An answer in it from 0 is 0 times an infinity there, which a + // definite integral takes the limit of. https://github.com/asc-community/AngouriMath/issues/2075 + partial record IncompleteGammaf + { + internal override Entity? ComputeLimitDivideEtImpera(Variable x, Entity dist, ApproachFrom side) + { + if (Parameter.ContainsNode(x)) + return null; + // By continuity alone: the value at a point where it is infinite is not the limit there -- gamma(-1, 0) is +oo, + // the integral diverging, where gamma(-1, z) is about 1/z and tends to -oo from the left. + return LowerLimit.ComputeLimitDivideEtImpera(x, dist, side) is { IsFinite: true } limit && !limit.ContainsNode(x) + && New(Parameter, limit) is var value && value.Evaled is Number.Complex { IsFinite: true, IsNaN: false } + ? value + : null; + } + } + partial record Derivativef { internal override Entity? ComputeLimitDivideEtImpera(Variable x, Entity dist, ApproachFrom side) => diff --git a/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs new file mode 100644 index 000000000..4a6c6f022 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs @@ -0,0 +1,108 @@ +// +// 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 System.Linq; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// x^m F^(a + b x^n) and (c + d x)^m F^(a + b x), whose integral is a gamma function of + /// s = (m + 1)/n where s is neither whole nor half a whole number. Rubi's 2.3. + /// #2075 + /// + [Trait("Area", "Calculus")] + public sealed class APowerBesideAnExponentialOfAPowerIntegralTest + { + private static readonly (string Name, string Value)[] Pins = + { ("a", "1/3"), ("b", "-2/3"), ("c", "5/4"), ("d", "1/2"), ("f", "3"), ("m", "3/7"), ("n", "5/2") }; + + /// Where every integrand here is real. + private static readonly double[] Points = { 0.35, 0.9, 1.45, 2.3 }; + + private static Entity DifferentiatesBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => Pins.Aggregate(e, (current, pin) => current.Substitute(pin.Name, pin.Value.ToEntity())); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in Points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart)); + Assert.True(difference / scale < 1e-9, $"d/dx of {integral} is {got} at x = {at}, where {integrand} is {want}"); + } + return integral; + } + + [Theory] + [InlineData("x^m*e^(b*x)")] + [InlineData("(d*x)^m*e^(b*x)")] + [InlineData("x^(1/3)*e^(-x)")] + [InlineData("f^(a + b*x^3)")] + [InlineData("f^(a + b*x^3)*x")] + [InlineData("e^(b*x^3)/x^2")] + [InlineData("x^m*f^(a + b*x^n)")] + [InlineData("f^(a + b/x)*x^m")] + [InlineData("f^(a + b/x^2)*x^m")] + [InlineData("(d*x)^m*x^2*e^(b*x)")] + [InlineData("x^(1/3)*sqrt(x)*e^(-x)")] + // A whole order is elementary before it is integrated, which leaves exactly these: Rubi's 8.6. + [InlineData("(d*x)^m*gamma(3, b*x)")] + public void APowerOfXBesideAnExponentialOfAPowerOfX(string integrand) => DifferentiatesBack(integrand); + + /// + /// Past 30, a whole s as the gamma function, where parts would write s terms. + /// + [Fact] + public void AHighWholePowerIsOneGammaFunction() + => Assert.Contains("gamma(101, a * x)", DifferentiatesBack("x^100*e^(-a*x)").Stringize()); + + [Theory] + [InlineData("(c + d*x)^m*e^(a + b*x)")] + [InlineData("(c + d*x)^(1/3)*f^(a + b*x)")] + [InlineData("e^(a + b*x)/(c + d*x)^(4/3)")] + public void APowerOfALinearBesideAnExponentialOfAnother(string integrand) => DifferentiatesBack(integrand); + + /// + /// Where x^(m + 1) z^(-s) is one number on both half-lines, it is written as that number, and the answer has a + /// value at 0: from 0 to 1, by quadrature, 1.34190441797742 and 0.807511182139671. + /// + [Theory] + [InlineData("integral(e^(x^3), x, 0, 1)", 1.34190441797742)] + [InlineData("integral(e^(-x^3), x, 0, 1)", 0.807511182139671)] + public void ADefiniteIntegralFromZero(string integral, double value) + { + var result = integral.ToEntity().Simplify().EvalNumerical(); + Assert.True(Math.Abs((double)result.ImaginaryPart) < 1e-20 && Math.Abs((double)result.RealPart - value) < 1e-12, $"{result}"); + } + + /// + /// With a positive rate, x^(m + 1) (k x)^(-s) is k^(-s), and is written so. + /// + [Fact] + public void APositiveRateIsTakenOutOfThePower() + => Assert.Equal("-gamma(4/3, x)", DifferentiatesBack("x^(1/3)*e^(-x)").Substitute("C", 0).InnerSimplified.Stringize()); + + /// + /// Not where the integral is elementary, an exponential integral or an error function: those keep the form + /// they had. + /// + [Theory] + [InlineData("x^2*e^(b*x^3)")] + [InlineData("sqrt(x)*e^(-x)")] + [InlineData("e^(b*x^2)")] + [InlineData("e^(b*x)/x")] + public void NotWhereTheIntegralIsElementaryOrAnErrorFunction(string integrand) + => Assert.DoesNotContain("gamma(", DifferentiatesBack(integrand).Stringize()); + } +} diff --git a/Sources/Tests/UnitTests/Calculus/APowerTimesAFunctionOfItsPowerTest.cs b/Sources/Tests/UnitTests/Calculus/APowerTimesAFunctionOfItsPowerTest.cs index 961430160..1ce6d02db 100644 --- a/Sources/Tests/UnitTests/Calculus/APowerTimesAFunctionOfItsPowerTest.cs +++ b/Sources/Tests/UnitTests/Calculus/APowerTimesAFunctionOfItsPowerTest.cs @@ -78,9 +78,16 @@ public void TheScaleIsSaidToBePositive() Assert.Equal("e^(x^n)/n".ToEntity().Simplify(), "x^(n-1)*e^(x^n)".ToEntity().Integrate("x").Substitute("C", 0).Simplify()); } - /// A power in front that is not one less than a multiple of the power inside is left alone. + /// + /// A power in front that is not one less than a multiple of the power inside is not substituted for: its integral is + /// a gamma function of -x^n, of order (n + 1)/n. + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// [Fact] - public void APowerThatIsNotTheDerivativeIsLeftAsWritten() - => Assert.Contains("integral(", "x^n*e^(x^n)".ToEntity().Integrate("x").Stringize()); + public void APowerThatIsNotTheDerivativeIsAGammaFunction() + { + DifferentiatesBack("x^n*e^(x^n)"); + Assert.Contains("gamma(", "x^n*e^(x^n)".ToEntity().Integrate("x").Stringize()); + } } } diff --git a/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs b/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs index 3f2d0bcb2..85e8febd5 100644 --- a/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs +++ b/Sources/Tests/UnitTests/Calculus/BreakpointIntegrationTest.cs @@ -5,6 +5,7 @@ // Website: https://am.angouri.org. // +using System; using AngouriMath; using AngouriMath.Extensions; using Xunit; @@ -119,10 +120,23 @@ public void TheSumTheSplitLeavesIsWhatAnswersTheIntegral() // is not a break this reads, and an integrand the unit interval cannot integrate stays. [Theory] [InlineData("integral(floor(x^2), x, 0, 2)")] - [InlineData("integral(floor(x) * e^(x^3) , x, 0, 2)")] + [InlineData("integral(floor(x) * e^(x^3 + x) , x, 0, 2)")] public void WhatCannotBeSplitIsLeftAsWritten(string integral) => Assert.IsType(integral.ToEntity().Simplify()); + /// + /// And split, where the unit interval integrates it: e^(x^3) from 1 to 2 as a gamma function of -x^3, + /// whose values on the cut take the power's real root. Quadrature gives 275.510983763311601498257886996. + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// + [Fact] + public void AGammaFunctionOnItsCutHasTheValuesOfAnAntiderivative() + { + var value = "integral(floor(x) * e^(x^3), x, 0, 2)".ToEntity().Simplify().EvalNumerical(); + Assert.True(Math.Abs((double)value.ImaginaryPart) < 1e-20, $"{value}"); + Assert.True(Math.Abs((double)value.RealPart - 275.510983763311601498257886996) < 1e-12, $"{value}"); + } + // ---- what does not break is unaffected -------------------------------------------------- [Theory] diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs index 585daf23d..85d43eb9f 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs @@ -180,10 +180,11 @@ public void AnExponentialTimesARootOfAQuadraticIsDeclinedWhereThereIsNone() /// /// What has no antiderivative in the library's functions is declined, not answered: the /// ansatz finds no N, which by Liouville is the proof that none is elementary. - /// e^x/x and e^x/(1 + x) are not elementary either, but they are Ei. + /// e^x/x and e^x/(1 + x) are not elementary either, but they are Ei, and + /// e^(x^3) is a gamma function of -x^3; e^(x^3 + x) is neither. /// [Theory] - [InlineData("e^(x^3)")] + [InlineData("e^(x^3 + x)")] public void ANonElementaryOneIsDeclined(string integrand) => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); diff --git a/Sources/Tests/UnitTests/Calculus/ProductOfPowersDerivativeTest.cs b/Sources/Tests/UnitTests/Calculus/ProductOfPowersDerivativeTest.cs index 425e3cb77..ce7f4ab31 100644 --- a/Sources/Tests/UnitTests/Calculus/ProductOfPowersDerivativeTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ProductOfPowersDerivativeTest.cs @@ -98,12 +98,25 @@ public void TheConstantIsRead() "p*x^m*ln(x)^n*(1 + n + (1 + m)*ln(x))".ToEntity().Integrate("x").Stringize()); /// - /// A sum that is not the derivative's is declined, not answered: x^m ln(x)^n (1 + ln x) - /// with symbols for m and n is not elementary. + /// A sum that is not the derivative's is not answered as one: x^m ln(x)^n (1 + ln x) with symbols for + /// m and n is not elementary, and under t = ln x it is two gamma functions of + /// -(m + 1) ln x. https://github.com/asc-community/AngouriMath/issues/2075 /// [Fact] - public void WhatIsNotADerivativeIsDeclined() - => Assert.Contains("integral(", "x^m*ln(x)^n*(1 + ln(x))".ToEntity().Integrate("x").Stringize()); + public void WhatIsNotADerivativeIsNotAnsweredAsOne() + { + var integrand = "x^m*ln(x)^n*(1 + ln(x))".ToEntity(); + var integral = integrand.Integrate("x"); + Assert.Contains("gamma(", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("m", "1/2".ToEntity()).Substitute("n", "3/4".ToEntity()).Substitute("C", 0); + var derivative = Pinned(integral).Differentiate("x"); + foreach (var at in new[] { "7/5", "23/10" }) + { + var got = derivative.Substitute("x", at.ToEntity()).EvalNumerical(); + var want = Pinned(integrand).Substitute("x", at.ToEntity()).EvalNumerical(); + Assert.True(((double)(got - want).Abs()) < 1e-9 * Math.Max(1, (double)want.Abs()), $"d/dx {integral} is {got} at {at}, where {want}"); + } + } /// /// Rubi's 1.3.2, two symbolic powers of a quadratic and a cubic beside a power of x: