diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a7c295ac9..b4dbdeba9 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -679,6 +679,26 @@ back after. Each answer is in three cases by the sign of `4 a/c - (b/c)^2` | `"x^5/(a + b*x + c*x^2)^5".ToEntity().Integrate("x")` | `integral(...)` | 4,066 characters | | `"x^7/(a + b*x + c*x^2)^5".ToEntity().Integrate("x")` | `integral(...)` | 5,602 characters | +### A power of x below -1 beside a power of a binomial is raised by the reduction that takes out its pole + +**Answers where there were none.** `1/(x^3 (1 - x^3)^(2/3) (1 + x^3))`, Rubi's 1.1.3.4, was declined: split over +`1 + x^3` its terms are not elementary apart. So were `1/(x^6 (a + b x^2)^10)` and `1/(x^7 (a + b x^2)^10)`, Rubi's +1.1.2.2, which since 2.5.0 have run out of time in the partial fractions instead. The derivative of +`x^(m + 1) (a + b x^n)^(p + 1)` is `x^m (a + b x^n)^p ((m + 1) a + b (m + 1 + n (p + 1)) x^n)`, so a multiple of it +takes the pole at zero out of `x^m (a + b x^n)^p R(x^n)` and raises the power of x by n; what is left at `x^(-1)` or +above is asked on its own. Beside a root it is asked last, at the question asked only; over a whole power, with a +symbol in the binomial and past the twenty-fourth degree below the bar only, where it answers +`1/(x^5 (a + b x^2)^10)` in 1,860 characters and the partial fractions have since given 67,808 +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/(x^3*(1-x^3)^(2/3)*(1+x^3))".ToEntity().Integrate("x")` | `integral(...)` | 715 characters | +| `"1/(x^3*(1+x^3)^(2/3)*(1-x^3))".ToEntity().Integrate("x")` | `integral(...)` | 577 characters | +| `"1/(x^6*(a+b*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 2,286 characters | +| `"1/(x^7*(a+b*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 1,892 characters | +| `"1/(x^5*(a+b*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 1,860 characters | + ### A power of x below the bar beside a root of a binomial is rationalised by dividing by the root **Answers where there were none.** `(a + b x^3)^(1/3)/(x^2 (c + d x^3))` was declined, and so were sixteen more of Rubi's @@ -2924,6 +2944,21 @@ Rubi's `x^m (a + b x^n)^p` and `(a + b x^n)^p (c + d x^n)^q` | `"1/(2+3/x^2)^3".Integrate("x")`, `1/(a + b/x^3)` | left unevaluated | an antiderivative | | `"(a+b*x^n)*(c+d*x^n)^3".Integrate("x")` | left unevaluated | written out, eight powers of `x` integrated: `a c^3 x + ... + b d^3 x^(4 n + 1)/(4 n + 1)` | +### A product of powers with symbolic exponents is integrated where it is a derivative whose constant the brackets hide + +**Answers where there were none.** `(a - b x^(n/2))^p (a + b x^(n/2))^p (c + d x^n)^q` with +`c = -a^2 d n (p + 1)/(b^2 (1 + n + n p))` and `q = -(1 + 2 n + n p)/n`, Rubi's 1.1.3.3, was declined. It is +`-b^2 (1 + n + n p)/(a^4 d n (p + 1))` times the derivative of `x (a - b x^(n/2))^(p + 1) (a + b x^(n/2))^(p + 1) (c + d x^n)^(q + 1)`, +and the rule that reads such derivatives found the two brackets proportional at sampled points but could not show their +quotient free of x, which takes `x x^(n - 1) = x^n` and the symbolic coefficients of `x^n` and `x^(2 n)` cancelling. +Asked once more last in the chain, at the question asked only, it reads the constant at `x = 1` +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a - b*x^(n/2))^p*(a + b*x^(n/2))^p*(a^2*d*(1 + p)/(b^2*(1 + (-1 - 2*n - n*p)/n)) + d*x^n)^((-1 - 2*n - n*p)/n)".ToEntity().Integrate("x")` | `integral(...)` | 247 characters | +| `"(a^2 - b^2*x^n)^p*(-a^2*d*n*(1 + p)/(b^2*(1 + n + n*p)) + d*x^n)^((-1 - 2*n - n*p)/n)".ToEntity().Integrate("x")` | `integral(...)` | 200 characters | + ### The derivative of a power of `x` over a binomial is integrated as an arctangent over one plus a multiple of its square **Answers where there were none.** `(k - 4 f x^3)/(k^2 + 4 d f x^2 + 4 k f x^3 + 4 f^2 x^6)` is @@ -3447,6 +3482,23 @@ improper fraction is declined before the first division rather than after the la | `"(1 - b*x^2)^3/(c*(1 - b*x^2) + a*d*x^2)^3".ToEntity().Integrate("x")` | no answer within a minute | the antiderivative | | `"(a + b*x)^(5/2)/(c + d*x)^4".ToEntity().Integrate("x")` | `integral(...)` | the antiderivative | +### A root over the root of a quadratic, or over two roots of linears with one slope, is integrated in the variable that rationalises them + +**Answers where there were none.** `sqrt(1 - x^2 + x sqrt(x^2 - 1))` and `sqrt(-x + sqrt(x) sqrt(1 + x))/sqrt(1 + x)`, +Rubi's 1.3.2, were declined, and Rubi gives no answer for either. Euler's substitution `t = x + sqrt(x^2 - 1)` makes +the root inside rational and the radicand outside the rational function `(t^2 - 1)/(2 t^2)`; for two roots of linears +with one slope `v = sqrt(x) + sqrt(1 + x)` does the same, and the integrand in the new variable is answered. The +answer is checked where the integrand is real: for the first that is `x >= 1`, and where `x <= -1`, where the +integrand is imaginary, its derivative is the integrand's negative +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(1 - x^2 + x*sqrt(x^2 - 1))".ToEntity().Integrate("x")` | `integral(...)` | 199 characters | +| `"sqrt(-x + sqrt(x)*sqrt(1 + x))/sqrt(1 + x)".ToEntity().Integrate("x")` | `integral(...)` | 197 characters | +| `"sqrt(1 + x^2 + x*sqrt(1 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | 230 characters | +| `"sqrt(x^2 + x*sqrt(x^2 - 1))".ToEntity().Integrate("x")` | `integral(...)` | 228 characters | + ### Four nested roots with a closed form are integrated **Answers where there were none.** Rubi's 1.3.3 gives four nested roots a closed form, and they were diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a1d92a989..065a34212 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -8897,6 +8897,227 @@ static Entity ToThe(Entity @base, int power) => power == 0 ? Number.Integer.One : MathS.Pow(@base, Number.Integer.Create(power)); } + /// + /// A whole power of x below -1 beside a power of a binomial B = a + b x^n that is + /// not whole and a rational function R of x^n with a value at zero, its pole + /// at zero taken out by the reduction that raises the power of x by n: + /// + /// int x^m B^p R(x^n) = C x^(m + 1) B^(p + 1) + int x^(m + n) B^p S(x^n), + /// C = R(0)/(a (m + 1)), S(t) = (R(t) - C ((m + 1) a + b (m + 1 + n (p + 1)) t))/t, + /// + /// repeated until the power of x is -1 or more, and the rest asked on its own. + /// + /// + /// + /// The derivative of x^(m + 1) B^(p + 1) is x^m B^p ((m + 1) a + b (m + 1 + n (p + 1)) x^n) + /// on every branch, so the step is exact, and S is a rational function since what + /// it divides by t vanishes at zero. Rubi's 1.1.3.4 1/(x^3 (1 - x^3)^(2/3) (1 + x^3)) + /// was declined: split over 1 + x^3 its terms are not elementary apart, and the + /// cube-root ansatz stops short of the ninth degree its denominator has. One step leaves + /// -(1 - x^3)^(1/3)/(2 (1 + x^3)), which is answered. + /// + /// + /// Last in the chain, and at the question asked only: asked before the scaling of the + /// variable and by parts, or below them, it answered six solved rows of Rubi's 1.1.3.4, + /// 1/(x^6 (a + b x^3)^(1/3) (c + d x^3)) among them, in other spellings, three far + /// shorter and one twice as long. Not for the binomial alone, which Chebyshev's + /// theorem decides. A factor B left in the rest is taken into the power. The answer + /// is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveANegativePowerOfXBesideARootOfABinomialByRaisingThePowerOfX(Entity expr, Entity.Variable x, bool integrateByParts) + => SolveANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfX(expr, x, integrateByParts, wholePower: false); + + /// + /// The reduction of + /// for x^m B^p alone with a whole p <= -2 and a symbol in B: Rubi's + /// 1.1.2.2 1/(x^6 (a + b x^2)^10) is three steps from 1/(a + b x^2)^10. + /// + /// + /// The partial fractions take the binomial apart over its symbolic roots, and past the + /// twenty-fourth degree below the bar they ran out the budget: 1/(x^4 (a + b x^2)^10) + /// took a second, 1/(x^6 (a + b x^2)^10) and 1/(x^7 (a + b x^2)^10) more than + /// twenty, and 1/(x^5 (a + b x^2)^10) about two, for 67,808 characters where this gives + /// 1,860. Asked before the reduction over a power of a symbolic quadratic, which reads x as + /// a linear beside it, and past the twenty-fourth degree only, so that what the partial + /// fractions answered in time is answered as before. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveANegativePowerOfXOverAWholePowerOfABinomialByRaisingThePowerOfX(Entity expr, Entity.Variable x, bool integrateByParts) + => SolveANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfX(expr, x, integrateByParts, wholePower: true); + + private static Entity? SolveANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfX(Entity expr, Entity.Variable x, bool integrateByParts, bool wholePower) + { + if (!wholePower && !Integration.AnsweringTheQuestionAsked) + return null; + Entity constant = Number.Integer.One; + var m = 0; + Entity? binomial = null; + var p = ERational.Zero; + var above = new List(); + var below = new List(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, exponent) = factor is Powf(var raised, Number.Rational power) ? (raised, power.ERational) : (factor, ERational.One); + if (underneath) + exponent = exponent.Negate(); + if (@base == x) + { + if (!exponent.Denominator.Equals(EInteger.One) || !exponent.Numerator.CanFitInInt32()) + return null; + m += exponent.Numerator.ToInt32Checked(); + continue; + } + if (!exponent.Denominator.Equals(EInteger.One) || wholePower && @base is Sumf or Minusf) + { + if (binomial is not null && binomial != @base) + return null; + binomial = @base; + p = p.Add(exponent); + continue; + } + if (wholePower || !exponent.Numerator.CanFitInInt32() || exponent.Numerator.Abs().CompareTo(EInteger.FromInt32(MostWrittenFactorsRaised)) > 0) + return null; + var times = exponent.Numerator.ToInt32Checked(); + for (var i = 0; i < System.Math.Abs(times); i++) + (times > 0 ? above : below).Add(@base); + } + if (binomial is null || m > -2 || p.Denominator.Equals(EInteger.One) != wholePower + || !wholePower && above.Count + below.Count == 0 + || wholePower && (p.CompareTo(ERational.FromInt32(-2)) > 0 || !binomial.Vars.Any(symbol => symbol != x)) + || !TryReadAsABinomialIn(binomial, x, out var a, out var b, out var n) + || wholePower && -m + n * p.Numerator.Negate().ToInt32Checked() <= MostDegreeThePartialFractionsTake + || VanishesIdentically(a) || VanishesIdentically(b) + || (-1 - m + n - 1) / n > MostPowersOfXRaised) + return null; + // Beside a root, a power of x one less than a multiple of n is `x^(n - 1)` times a function + // of `t = x^n`, and the root a root of a linear in t, which the substitutions read whole: + // raised, Rubi's `1/(x^5 (1 - x^2)^(1/3) (3 + x^2)^2)` asked what was left as a question of + // its own, declined it in twelve seconds and twenty-two gigabytes, and was answered as before. + if (!wholePower && (m + 1) % n == 0) + return null; + + // The rest in t = x^n: each written factor a polynomial in it. + Dictionary? InT(Entity polynomial) + { + if (!TreeAnalyzer.TryGetPolynomial(polynomial, x, out var monomials)) + return null; + var read = new Dictionary(); + foreach (var pair in monomials) + { + if (pair.Value.ContainsNode(x) || pair.Key.Sign < 0 || !pair.Key.CanFitInInt32() || pair.Key.ToInt32Checked() % n != 0) + return null; + read[pair.Key.ToInt32Checked() / n] = pair.Value; + } + return read; + } + Dictionary? Product(List factors) + { + var product = new Dictionary { [0] = Number.Integer.One }; + foreach (var factor in factors) + { + if (InT(factor) is not { } read) + return null; + product = Multiplied(product, read); + if (product.Keys.Max() > MostDegreeInThePower) + return null; + } + return product; + } + if (Product(above) is not { } numerator || Product(below) is not { } denominator + || VanishesIdentically(At(numerator, 0)) || VanishesIdentically(At(denominator, 0))) + return null; + + Entity raisedPower = Exponent(p.Add(ERational.One)); + Entity sum = Number.Integer.Zero; + while (m <= -2) + { + // C = R(0)/(a (m + 1)), and the numerator of R - C L over the denominator of R + // divided by t, for L = (m + 1) a + b (m + 1 + n (p + 1)) t. + var coefficient = LowestOverTheSymbols(At(numerator, 0) / (At(denominator, 0) * Number.Integer.Create(m + 1) * a)); + var linear = new Dictionary + { + [0] = Number.Integer.Create(m + 1) * a, + [1] = b * Number.Rational.Create(ERational.FromInt32(m + 1).Add(ERational.FromInt32(n).Multiply(p.Add(ERational.One)))), + }; + var taken = Multiplied(denominator, linear); + var next = new Dictionary(); + foreach (var degree in numerator.Keys.Union(taken.Keys)) + next[degree] = LowestOverTheSymbols(At(numerator, degree) - coefficient * At(taken, degree)); + if (!VanishesIdentically(At(next, 0))) + return null; + numerator = next.Where(pair => pair.Key > 0).ToDictionary(pair => pair.Key - 1, pair => pair.Value); + if (!VanishesIdentically(coefficient)) + sum += coefficient * (m + 1 == -1 ? Number.Integer.One / x : MathS.Pow(x, Number.Integer.Create(m + 1))) * MathS.Pow(binomial, raisedPower); + m += n; + } + + // A factor B of what is left above is taken into its power. + var restPower = p; + for (var divided = 0; divided < MostWrittenFactorsRaised && numerator.Count > 0 && numerator.Keys.Max() >= 1; divided++) + { + var top = numerator.Keys.Max(); + var quotient = new Dictionary(); + for (var k = top; k >= 1; k--) + quotient[k - 1] = LowestOverTheSymbols((At(numerator, k) - a * At(quotient, k)) / b); + if (!VanishesIdentically(At(numerator, 0) - a * At(quotient, 0))) + break; + numerator = quotient; + restPower = restPower.Add(ERational.One); + } + Entity above0 = Number.Integer.Zero; + foreach (var pair in numerator.OrderBy(pair => pair.Key)) + if (!VanishesIdentically(pair.Value)) + above0 += pair.Key == 0 ? pair.Value : pair.Value * MathS.Pow(x, Number.Integer.Create(pair.Key * n)); + // Nothing left where the integrand is the derivative of the terms. + Entity integral = Number.Integer.Zero; + if (above0 != Number.Integer.Zero) + { + Entity rest = constant * above0 * MathS.Pow(binomial, Exponent(restPower)); + if (m != 0) + rest *= MathS.Pow(x, Number.Integer.Create(m)); + foreach (var factor in below) + rest /= factor; + if (Integration.ComputeAsAQuestionOfItsOwn(Functions.PartialFractions.Bare(rest), x, integrateByParts) is not { } asked) + return null; + integral = asked; + } + var answer = (constant * sum + integral).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + + static Entity At(Dictionary polynomial, int degree) + => polynomial.TryGetValue(degree, out var coefficient) ? coefficient : Number.Integer.Zero; + static Entity Exponent(ERational power) + => power.Denominator.Equals(EInteger.One) ? Number.Integer.Create(power.Numerator) : Number.Rational.Create(power); + static Dictionary Multiplied(Dictionary left, Dictionary right) + { + var product = new Dictionary(); + foreach (var l in left) + foreach (var r in right) + product[l.Key + r.Key] = At(product, l.Key + r.Key) + l.Value * r.Value; + return product.ToDictionary(pair => pair.Key, pair => LowestOverTheSymbols(pair.Value)); + } + } + + /// The most steps takes. + private const int MostPowersOfXRaised = 8; + + /// The most times one written factor stands in what reads, and the most factors B it takes into the power. + private const int MostWrittenFactorsRaised = 6; + + /// The largest degree in x^n of the rational function reads. + private const int MostDegreeInThePower = 12; + + /// The degree below the bar up to which leaves x^m B^p to the partial fractions. + private const int MostDegreeThePartialFractionsTake = 24; + /// /// Reads as a multiple of x^m (a + b x^n)^(p/q), with /// m and n rational, q above one, n neither zero nor one, and @@ -18599,6 +18820,30 @@ static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.R /// https://github.com/asc-community/AngouriMath/issues/718 /// internal static Entity? SolveAsTheDerivativeOfAProductOfPowers(Entity expr, Entity.Variable x, bool integrateByParts) + => AsTheDerivativeOfAProductOfPowers(expr, x, integrateByParts, readTheRatioAtAPoint: false); + + /// + /// with the + /// constant read at a point where the brackets are proportional at the sampled points and + /// their quotient is not read free of x: Rubi's 1.1.3.3 + /// (a - b x^(n/2))^p (a + b x^(n/2))^p (c + d x^n)^q at c = -a^2 d n (p + 1)/(b^2 (1 + n + n p)) + /// and q = -(1 + 2 n + n p)/n is the derivative of + /// x (a - b x^(n/2))^(p + 1) (a + b x^(n/2))^(p + 1) (c + d x^n)^(q + 1) times + /// -b^2 (1 + n + n p)/(a^4 d n (p + 1)), and was declined. + /// + /// + /// Its bracket is a constant only once x x^(n - 1) is x^n and the symbolic + /// coefficients of x^n and x^(2 n) cancel, which neither the reading off one + /// monomial nor the simplification of the quotient sees; at x = 1 the quotient is the + /// constant. Last in the chain and at the question asked only, so that what another rule + /// answers is answered as before; the answer is differentiated back at sampled points, as + /// the first's is. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveAsTheDerivativeOfAProductOfPowersReadAtAPoint(Entity expr, Entity.Variable x, bool integrateByParts) + => Integration.AnsweringTheQuestionAsked ? AsTheDerivativeOfAProductOfPowers(expr, x, integrateByParts, readTheRatioAtAPoint: true) : null; + + private static Entity? AsTheDerivativeOfAProductOfPowers(Entity expr, Entity.Variable x, bool integrateByParts, bool readTheRatioAtAPoint) { if (!Integration.AnsweringTheQuestionAskedOrOneBelow) return null; @@ -18609,13 +18854,13 @@ static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.R // elementary. A base with whole exponents may be missing from some terms: Rubi's // `(1 + e^x) x/sqrt(e^x + x) + 2 sqrt(e^x + x)` is `(2 x sqrt(e^x + x))'`. if (expr is Sumf or Minusf) - return WithTheCommonPowersOut(expr, x) is { } factored ? SolveAsTheDerivativeOfAProductOfPowers(factored, x, integrateByParts) : null; + return WithTheCommonPowersOut(expr, x) is { } factored ? AsTheDerivativeOfAProductOfPowers(factored, x, integrateByParts, readTheRatioAtAPoint) : null; // A bracket that is a polynomial with a power of x in every term has that power taken out, to // join the powers it multiplies: `x^(2n) (c + d x)^n (2 c x + 3 d x^2)` is // `x^(2n + 1) (c + d x)^n (2 c + 3 d x)`, the derivative of `x^(2n + 2) (c + d x)^(n + 1)/(n + 1)`, // and with the x in the bracket nothing raised x by enough. if (WithAPowerOfXOutOfTheBracket(expr, x) is { } outOfTheBracket) - return SolveAsTheDerivativeOfAProductOfPowers(outOfTheBracket, x, integrateByParts); + return AsTheDerivativeOfAProductOfPowers(outOfTheBracket, x, integrateByParts, readTheRatioAtAPoint); Entity constant = Number.Integer.One; Entity? bracket = null; var raisable = new List<(Entity Base, Entity Exponent)>(); @@ -18788,6 +19033,10 @@ static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.R if (k.ContainsNode(x) && ASmallRationalRatioAtAPoint(bracket, bracketOfTheDerivative, x) is { } ratio) k = ratio; k = Functions.PartialFractions.Bare((constant * k).InnerSimplified); + // Proportional at the sampled points, and the quotient not read free of x: the + // constant is then the quotient at a point, asked for last in the chain. + if (readTheRatioAtAPoint && k.ContainsNode(x) && TheRatioOfProportionalBracketsAtAPoint(bracket, bracketOfTheDerivative, x) is { } atAPoint) + k = Functions.PartialFractions.Bare((constant * atAPoint).InnerSimplified); if (k.ContainsNode(x) || k.Nodes.Any(node => node == MathS.NaN)) continue; var answer = (k * candidate).InnerSimplified; @@ -19536,6 +19785,27 @@ static EInteger Binomial(int n, int r) /// private const int MaximumLinearPowerSum = 12; + /// + /// The constant over is, for two brackets + /// already shown proportional at sampled points, read at x = 1, or at x = 2 + /// where the right one vanishes at the first; null where it has no value at either. + /// + private static Entity? TheRatioOfProportionalBracketsAtAPoint(Entity left, Entity right, Entity.Variable x) + { + foreach (var at in new[] { 1, 2 }) + { + var below = right.Substitute(x, at).InnerSimplified; + if (VanishesIdentically(below)) + continue; + var ratio = Functions.PartialFractions.Bare((left.Substitute(x, at) / below).InnerSimplified); + if (ratio.ContainsNode(x) || ratio.Nodes.Any(node => node == MathS.NaN)) + continue; + var simplified = Functions.PartialFractions.Bare(Functions.SingleQuotient.Combine(ratio).Simplify()); + return simplified.Complexity < ratio.Complexity && !simplified.Nodes.Any(node => node == MathS.NaN) ? simplified : ratio; + } + return null; + } + /// /// The constant k with equal to k times /// , read off one monomial with both written as polynomials @@ -28386,6 +28656,226 @@ static Entity Lowest(Entity constant) return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; } + /// + /// A root whose radicand holds roots one substitution rationalises -- the root of a + /// quadratic, by Euler's, or the roots of two linears with one slope, by their sum -- + /// integrated in the substitution's variable, where the radicand is a rational function: + /// sqrt(1 - x^2 + x sqrt(x^2 - 1)) is sqrt((t^2 - 1)/(2 t^2)) (t^2 - 1)/(2 t^2) at + /// t = x + sqrt(x^2 - 1), and sqrt(-x + sqrt(x) sqrt(1 + x))/sqrt(1 + x) is + /// sqrt((v^2 - 1)/(2 v^2)) (v^2 - 1)/v^2 at v = sqrt(x) + sqrt(1 + x). + /// + /// + /// + /// For Q = a x^2 + b x + c, Euler's first substitution t = sqrt(Q) + sqrt(a) x has + /// x = (t^2 - c)/(2 sqrt(a) t + b) and sqrt(Q) = t - sqrt(a) x, and his second, + /// t = (sqrt(Q) - sqrt(c))/x, has x = (2 sqrt(c) t - b)/(a - t^2) and + /// sqrt(Q) = x t + sqrt(c). For L1 - L2 = k with one slope m, + /// v = sqrt(L1) + sqrt(L2) has sqrt(L1) = (v + k/v)/2, sqrt(L2) = (v - k/v)/2 + /// and dx = (v^4 - k^2)/(2 m v^3) dv. Each is an identity in the values of the roots + /// wherever they are defined, whatever their branches, so the radicand written in the new + /// variable is the radicand, and the integrand there a rational function beside roots of + /// rational functions, which the chain is asked. + /// + /// + /// Rubi's 1.3.2 has both above without an answer; the Euler substitution and the sum of two + /// roots read only a rational function of the roots, and declined. Numbers only, roots + /// nested two deep, at the question asked or one below it, and late; the answer is + /// differentiated back at the sampled points where the integrand is real. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveANestedRootByRationalisingTheRootsInIt(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAskedOrOneBelow || expr.Vars.Any(symbol => symbol != x)) + return null; + static bool IsARootOf(Entity node, Entity.Variable x) + => node is Powf(var @base, Number.Rational power) && power is not Number.Integer && @base.ContainsNode(x); + // The roots inside, of x alone, and at least one root over them. + var radicands = new List(); + var nested = false; + foreach (var node in expr.Nodes) + { + if (!node.ContainsNode(x) || node is Entity.Variable or Sumf or Minusf or Mulf or Divf or Powf(_, Number.Integer)) + continue; + if (!IsARootOf(node, x)) + return null; + var radical = (Powf)node; + var @base = radical.Base; + var power = ((Number.Rational)radical.Exponent).ERational; + if (@base.Nodes.Any(inner => IsARootOf(inner, x))) + { + if (@base.Nodes.Any(inner => IsARootOf(inner, x) && ((Powf)inner).Base.Nodes.Any(deeper => IsARootOf(deeper, x)))) + return null; + // A radicand that is a polynomial in x and the roots, as Rubi's are: the tangent's + // substitution writes `(1 + coth(x))^(7/2)` as a power of `(u + 1)/u + sqrt(1 - u^2)/u`, + // which this rationalised in u at ten times the allocation and declined, or answered + // a level down in 1,661 characters where the rules after it give 287. + if (@base.Nodes.Any(inner => inner.ContainsNode(x) && (inner is Divf || inner is Powf(_, Number.Rational { IsNegative: true })))) + return null; + nested = true; + continue; + } + if (!power.Denominator.Equals(EInteger.FromInt32(2)) || !power.Numerator.CanFitInInt32()) + return null; + if (!radicands.Contains(@base)) + radicands.Add(@base); + } + if (!nested) + return null; + // Roots and nothing else of x, as Rubi's are: beside a rational factor, the tangent's + // substitution of `1/sqrt(1 - coth(x)^2)` was rationalised in u to be declined, at three + // times the allocation of the rules after it. + if (FactorsOfTheIntegrand(expr).Any(pair => pair.Factor.ContainsNode(x) && !IsARootOf(pair.Factor, x))) + return null; + + var t = Variable.CreateUnique(expr, "t_nested"); + Entity xInT, dxdt, back; + Entity[] rootsInT; + if (radicands.Count == 1 && TryReadAQuadraticOverAConstant(radicands[0], x, out var a, out var b, out var c)) + { + if (a.Evaled is not Number.Real { IsZero: false } aValue || b.Evaled is not Number.Real || c.Evaled is not Number.Real cValue + || IsAPerfectSquareDiscriminant(a, b, c)) + return null; + var root = MathS.Sqrt(radicands[0]); + if (aValue.IsPositive) + { + var sqrtA = MathS.Sqrt(a).InnerSimplified; + xInT = (MathS.Sqr(t) - c) / (2 * sqrtA * t + b); + rootsInT = new[] { t - sqrtA * xInT }; + back = root + sqrtA * x; + } + else if (cValue.IsPositive) + { + var sqrtC = MathS.Sqrt(c).InnerSimplified; + xInT = (2 * sqrtC * t - b) / (a - MathS.Sqr(t)); + rootsInT = new[] { xInT * t + sqrtC }; + back = (root - sqrtC) / x; + } + else + return null; + dxdt = xInT.Differentiate(t); + } + else if (radicands.Count == 2 + && TreeAnalyzer.TryGetPolyLinear(radicands[0], x, out var m1, out var n1) + && TreeAnalyzer.TryGetPolyLinear(radicands[1], x, out var m2, out var n2) + && !VanishesIdentically(m1) && VanishesIdentically(m1 - m2) && !VanishesIdentically(n1 - n2)) + { + var k = (n1 - n2).InnerSimplified; + rootsInT = new[] { (t + k / t) / 2, (t - k / t) / 2 }; + xInT = (MathS.Sqr(rootsInT[1]) - n2) / m1; + dxdt = (MathS.Pow(t, 4) - MathS.Sqr(k)) / (2 * m1 * MathS.Pow(t, 3)); + back = MathS.Sqrt(radicands[0]) + MathS.Sqrt(radicands[1]); + } + else + return null; + + // The roots inside in t, x in t, and every radicand left a rational function in + // lowest terms, as is the rest. + var inT = expr.Replace(node => IsARootOf(node, x) && radicands.IndexOf(((Powf)node).Base) is var index and >= 0 + ? MathS.Pow(rootsInT[index], Number.Integer.Create(((Number.Rational)((Powf)node).Exponent).ERational.Numerator)) + : node) + .Substitute(x, xInT) * dxdt; + Entity roots = Number.Integer.One; + Entity rational = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(inT)) + { + var piece = factor; + if (IsARootOf(factor, t)) + { + var power = (Number.Rational)((Powf)factor).Exponent; + if (InLowestTerms(((Powf)factor).Base, t) is not { } radicand) + return null; + piece = MathS.Pow(radicand, power); + roots = underneath ? roots / piece : roots * piece; + } + else + rational = underneath ? rational / piece : rational * piece; + } + if (InLowestTerms(rational, t) is not { } rest) + return null; + inT = rest * roots; + if (inT.ContainsNode(x) || inT.Nodes.Any(node => node == MathS.NaN) + || inT.Nodes.Any(node => IsARootOf(node, t) && ((Powf)node).Base.Nodes.Any(inner => IsARootOf(inner, t)))) + return null; + if (Integration.ComputeAsAQuestionOfItsOwn(inT, t, integrateByParts) is not { } inTermsOfT) + return null; + + // Checked where the integrand is real, at two points at least, with the sign of t + // held constant in the derivative, as it is between the zeros of t: the derivative + // of `sgn(sqrt(x) + sqrt(1 + x))` is not answered, since its argument is not shown real. + // First with that sign one and `|t|` as t, as they are where the roots inside are real + // and t positive, which is the integrand's real domain for both above. + using var _ = MathS.Settings.DowncastingEnabled.Set(false); + var real = NestedRootSamplePoints.Where(point => + expr.Substitute(x, Number.Real.Create(EDecimal.FromString(point))).Evaled is Number.Complex value && !value.IsNaN + && value.ImaginaryPart.EDecimal.Abs().CompareTo(EDecimal.FromString("1e-12")) <= 0).ToArray(); + if (real.Length < 2) + return null; + var sign = Variable.CreateUnique(expr * inTermsOfT, "s_nested"); + var signed = inTermsOfT.Replace(node => node switch + { + Signumf(var argument) when argument == t => sign, + Absf(var argument) when argument == t => sign * t, + _ => node, + }); + foreach (var value in new Entity[] { Number.Integer.One, MathS.Signum(back) }) + { + var answer = signed.Substitute(sign, value).Substitute(t, back).InnerSimplified; + if (answer.Nodes.Any(node => node == MathS.NaN)) + return null; + var derivative = Functions.PartialFractions.Bare(signed.Substitute(t, back)).Differentiate(x).Substitute(sign, value); + if (Functions.PartialFractions.HoldsAtSampledPoints(derivative, expr, x, real)) + return answer; + } + return null; + + // A rational function of t as one quotient of polynomials, cancelled, the content of + // numbers taken out of the two sides. + static Entity? InLowestTerms(Entity e, Entity.Variable t) + { + var (above, below) = Functions.SingleQuotient.Of(e.InnerSimplified); + if (!TreeAnalyzer.TryGetPolynomial(above, t, out var top) || !TreeAnalyzer.TryGetPolynomial(below, t, out var bottom) + || top.Keys.Concat(bottom.Keys).Any(degree => degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(MostNestedRootDegree)) > 0)) + return null; + var coefficients = top.Values.Concat(bottom.Values).Select(value => value.InnerSimplified.Evaled).ToList(); + var content = ERational.One; + if (coefficients.All(value => value is Number.Rational)) + { + var rationals = coefficients.Select(value => ((Number.Rational)value).ERational).Where(value => !value.IsZero).ToList(); + var numerators = rationals.Aggregate(EInteger.Zero, (gcd, value) => gcd.Gcd(value.Numerator.Abs())); + var denominators = rationals.Aggregate(EInteger.One, (lcm, value) => lcm * value.Denominator / lcm.Gcd(value.Denominator)); + if (!numerators.IsZero) + content = ERational.Create(numerators, denominators); + } + Entity Rebuilt(Dictionary read) + { + Entity built = Number.Integer.Zero; + foreach (var pair in read.OrderBy(pair => pair.Key)) + { + var coefficient = (pair.Value / Number.Rational.Create(content)).InnerSimplified; + if (coefficient.Evaled is Number.Complex { IsZero: true }) + continue; + built += pair.Key.IsZero ? coefficient : coefficient * (pair.Key.Equals(EInteger.One) ? t : MathS.Pow(t, Number.Integer.Create(pair.Key))); + } + return built.InnerSimplified; + } + var numerator = Rebuilt(top); + var denominator = Rebuilt(bottom); + if (denominator == Number.Integer.Zero) + return null; + return Functions.PolynomialGcd.TryCancel(numerator, denominator, out var cancelled) + ? Functions.PartialFractions.Bare(cancelled) + : (numerator / denominator).InnerSimplified; + } + } + + /// The largest degree a radicand or the rest of has in t. + private const int MostNestedRootDegree = 16; + + /// Where looks for the integrand to be real: both sides of zero and of one. + [ConstantField] + private static readonly string[] NestedRootSamplePoints = { "-3.17", "-1.43", "-0.61", "0.29", "0.83", "1.43", "2.3", "3.17" }; + /// /// A square root of a perfect square in is the modulus: /// sqrt(x^2) is |x|, which for a real x is sgn(x) x, and diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index a35a41750..42b3790b4 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -773,6 +773,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of a linear that is not whole, times a linear, over a whole power of a quadratic, // the power of the quadratic lowered one at a time. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearOverAPowerOfAQuadratic(expr, x, integrateByParts)) is { }) return answer; + // A power of x below -1 over a high whole power of a symbolic binomial, the power of x + // raised to -1 or more: before the next, which reads x as a linear beside the quadratic, + // and the partial fractions, which ran out the budget on `1/(x^6 (a + b x^2)^10)`. + if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfXOverAWholePowerOfABinomialByRaisingThePowerOfX(expr, x, integrateByParts)) is { }) return answer; // A rational function over a whole power of a symbolic quadratic, the power lowered one // at a time, before the partial fractions take the quadratic apart over its roots. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverAPowerOfASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer; @@ -1461,6 +1465,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // apart where its residue, the radicand at its root, is zero. if ((answer = IndefiniteIntegralSolver.SolveByWritingALinearBelowTheBarOverTheRadicand(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByEulerSubstitution(expr, x)) is { }) return answer; + // And a root over the root of a quadratic, or over the roots of two linears with one + // slope, in the variable that rationalises the roots inside: `sqrt(1 - x^2 + x sqrt(x^2 - 1))`. + if ((answer = IndefiniteIntegralSolver.SolveANestedRootByRationalisingTheRootsInIt(expr, x, integrateByParts)) is { }) return answer; // And a root of a palindromic quartic, which is a root of a quadratic in x -+ 1/x: // Charlwood's `(1 + x^2)/((1 - x^2) sqrt(1 + x^4))` is `-du/(u sqrt(u^2 + 2))`. if ((answer = IndefiniteIntegralSolver.SolveByReciprocalSubstitution(expr, x)) is { }) return answer; @@ -1606,6 +1613,14 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A rational function of x and ln(x) over distinct linears in the logarithm, by partial // fractions in it, where nothing has answered it whole: `(2 L^2 - L - x^2)/(L^3 - x^2 L)`. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOfTheLogarithmByPartialFractionsInIt(expr, x, integrateByParts)) is { }) return answer; + // A power of x below -1 beside a root of a binomial and a rational function of its + // power, the pole at zero taken out by the reduction raising the power of x. Last: + // the rules above answer the shapes they share with it in other spellings, and the + // cube-root ansatz stops at the ninth degree of `1/(x^3 (1 - x^3)^(2/3) (1 + x^3))`. + if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfXBesideARootOfABinomialByRaisingThePowerOfX(expr, x, integrateByParts)) is { }) return answer; + // And a product of powers that is the derivative of one with symbolic exponents, its + // constant read at a point where the quotient of the brackets does not simplify free of x. + if ((answer = IndefiniteIntegralSolver.SolveAsTheDerivativeOfAProductOfPowersReadAtAPoint(expr, x, integrateByParts)) is { }) return answer; // An arcoth of a linear as its artanh plus a constant, last of all, so that every arcoth // answered as written keeps that answer: `acoth(x)^2/(1 - x^2)^2` is // `(atanh(x) + K)^2/(1 - x^2)^2` with `K = acoth(x) - atanh(x)`, constant where it is diff --git a/Sources/Tests/UnitTests/Calculus/ANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfXIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfXIntegralTest.cs new file mode 100644 index 000000000..e0eea1b01 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfXIntegralTest.cs @@ -0,0 +1,62 @@ +// +// 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 AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A power of x below -1 beside a power of a binomial a + b x^n, its pole at zero taken + /// out by the reduction that raises the power of x by n. Beside a root of the binomial and a + /// rational function of x^n, Rubi's 1.1.3.4 1/(x^3 (1 - x^3)^(2/3) (1 + x^3)), whose + /// terms apart are not elementary; over a high whole power of a symbolic one, 1.1.2.2 + /// 1/(x^6 (a + b x^2)^10), which ran past the budget in the partial fractions. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ANegativePowerOfXBesideAPowerOfABinomialByRaisingThePowerOfXIntegralTest + { + [Theory] + [InlineData("1/(x^3*(1-x^3)^(2/3)*(1+x^3))")] + [InlineData("1/(x^3*(1+x^3)^(2/3)*(1-x^3))")] + public void BesideARootOfTheBinomial(string integrand) + => DifferentiatesBack(integrand, new[] { -1.7, -0.6, 0.3, 0.7, 1.4, 2.2 }); + + [Theory] + [InlineData("1/(x^6*(a + b*x^2)^10)")] + [InlineData("1/(x^7*(a + b*x^2)^10)")] + [InlineData("1/(x^5*(a + b*x^2)^10)")] + public void OverAHighWholePowerOfTheBinomial(string integrand) + => DifferentiatesBack(integrand, new[] { -1.9, -0.8, 0.3, 0.7, 1.4, 2.2 }); + + private static void DifferentiatesBack(string integrand, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in points) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/ANestedRootByRationalisingTheRootsInItIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ANestedRootByRationalisingTheRootsInItIntegralTest.cs new file mode 100644 index 000000000..1f4cd42e5 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ANestedRootByRationalisingTheRootsInItIntegralTest.cs @@ -0,0 +1,55 @@ +// +// 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 AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A root over the root of a quadratic, or over the roots of two linears with one slope, + /// integrated in the variable that rationalises the roots inside it: Rubi's 1.3.2 + /// sqrt(1 - x^2 + x sqrt(x^2 - 1)) under t = x + sqrt(x^2 - 1), and + /// sqrt(-x + sqrt(x) sqrt(1 + x))/sqrt(1 + x), which Rubi does not answer. Compared where + /// the integrand is real, which is where the answer is claimed. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ANestedRootByRationalisingTheRootsInItIntegralTest + { + [Theory] + [InlineData("sqrt(1 - x^2 + x*sqrt(x^2 - 1))", new[] { 1.1, 1.43, 1.9, 2.3, 3.17, 4.6 })] + [InlineData("sqrt(-x + sqrt(x)*sqrt(1 + x))/sqrt(1 + x)", new[] { 0.1, 0.29, 0.83, 1.43, 3.17, 4.6 })] + [InlineData("sqrt(1 + x^2 + x*sqrt(1 + x^2))", new[] { -3.17, -1.43, -0.61, 0.29, 1.43, 3.17 })] + [InlineData("1/sqrt(1 + x^2 + x*sqrt(1 + x^2))", new[] { -3.17, -1.43, -0.61, 0.29, 1.43, 3.17 })] + [InlineData("sqrt(x^2 + x*sqrt(x^2 - 1))", new[] { -3.17, -2.2, -1.43, 1.43, 2.2, 3.17 })] + [InlineData("sqrt(x + sqrt(x)*sqrt(1 + x))", new[] { 0.1, 0.29, 0.83, 1.43, 3.17, 4.6 })] + public void InTheVariableThatRationalisesTheRootsInside(string integrand, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + var derivative = integral.Substitute("C", 0).Differentiate("x"); + var original = integrand.ToEntity(); + var compared = 0; + foreach (var at in points) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN || Math.Abs((double)want.ImaginaryPart) > 1e-12) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/AsTheDerivativeOfAProductOfPowersReadAtAPointIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AsTheDerivativeOfAProductOfPowersReadAtAPointIntegralTest.cs new file mode 100644 index 000000000..3b0cbfb84 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AsTheDerivativeOfAProductOfPowersReadAtAPointIntegralTest.cs @@ -0,0 +1,54 @@ +// +// 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 AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A product of powers with symbolic exponents that is the derivative of another, where the + /// quotient of the two brackets is a constant only once the coefficients of x^n and + /// x^(2 n) cancel: Rubi's 1.1.3.3 (a - b x^(n/2))^p (a + b x^(n/2))^p (c + d x^n)^q + /// at c = -a^2 d n (p + 1)/(b^2 (1 + n + n p)) and q = -(1 + 2 n + n p)/n, which is + /// (x (a - b x^(n/2))^(p + 1) (a + b x^(n/2))^(p + 1) (c + d x^n)^(q + 1))' up to a constant. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AsTheDerivativeOfAProductOfPowersReadAtAPointIntegralTest + { + [Theory] + [InlineData("(a - b*x^(n/2))^p*(a + b*x^(n/2))^p*(a^2*d*(1 + p)/(b^2*(1 + (-1 - 2*n - n*p)/n)) + d*x^n)^((-1 - 2*n - n*p)/n)")] + [InlineData("(a - b*x^(n/2))^p*(a + b*x^(n/2))^p*(-a^2*d*n*(1 + p)/(b^2*(1 + n + n*p)) + d*x^n)^((-1 - 2*n - n*p)/n)")] + [InlineData("(a^2 - b^2*x^n)^p*(-a^2*d*n*(1 + p)/(b^2*(1 + n + n*p)) + d*x^n)^((-1 - 2*n - n*p)/n)")] + public void TheConstantReadAtAPoint(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 1000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.1).Substitute("b", 0.6).Substitute("d", 0.9) + .Substitute("n", 2.3).Substitute("p", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.7, 1.1, 1.6, 2.3, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}