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: