diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 0909ba42e..05fdd6165 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -55,6 +55,21 @@ numbers for every coefficient the solve is now done in exact rationals, to degre | `"1/((3-2*x+x^2)^(11/2)*(1+x+2*x^2)^5)".ToEntity().Integrate("x")` | `integral(...)` | 1,853 characters | | `"1/((3-2*x+x^2)^(21/2)*(1+x+2*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 3,475 characters | +### Two powers of exponentials with one logarithmic derivative are read as multiples of one exponential + +**Answers where there were none.** Rubi's 2.3 `(a + b (F^(e (c + d x)))^n)^p (G^(h (f + g x)))^(d e n ln(F)/(g h ln(G)))` +was declined, and on unreleased master ran past its budget, and so were the power of `G` over a power of the sum and +beside its root. A power of an exponential is `K e^(r x)`, `K` constant wherever it is differentiable and `r` its logarithmic derivative; two whose `r` +is one, decided by value, are written over the one exponential, the question in it is `(a + b K e^(r x))^p e^(r x)`, and +the answer is Rubi's. A power of an exponential was written so only beside a polynomial below the bar before +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a+b*(F^(k*(c+d*x)))^n)^p*(G^(h*(f+g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))".ToEntity().Integrate("x")` | `integral(...)` | `(G^(h (f + g x)))^(...) (a + b (F^(k (c + d x)))^n)^(p + 1)/((p + 1) b (F^(k (c + d x)))^n n d k ln(F))` | +| `"(G^(h*(f+g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))/(a+b*(F^(k*(c+d*x)))^n)^2".ToEntity().Integrate("x")` | `integral(...)` | 466 characters | +| `"(a+b*(F^(k*(c+d*x)))^n)^p*(F^(h*(f+g*x)))^(d*k*n/(g*h))".ToEntity().Integrate("x")` | `integral(...)` | 167 characters | + ### An exponential of something of the reciprocal of a linear is integrated beside another linear **Answers where there were none.** `f^(c/(a + b x)) x^2` was declined, and so were 24 more of Rubi's 2.3: @@ -246,6 +261,24 @@ It is written out again in the precision asked for now | `Rational.Create(128, 12167)`'s `EDecimal`, read at the default precision and then under 300 digits | 100 digits | 300 digits | | `sqrt(46) * 128/12167` evaluated under 300 digits, with that rational | relative error `4.7e-101` | correct to 300 digits | +### A polynomial over a power of a linear beside `a + b ln(c x^n)` is integrated by parts against its antiderivative vanishing at zero + +**Answers where there were none.** Rubi's 3.1.4 `x^5 (a + b ln(c x^n))/(d + e x)^7` was declined, and on unreleased +master ran past forty seconds, and so were `x^6 (a + b ln(c x^n))/(d + e x)^8` and `x^5 (a + b ln(x))/(d + e x)^7`. +By parts against the logarithm, the antiderivative of `P/(d + e x)^q` that vanishes at zero is `x M/(d + e x)^(q - 1)`, its coefficients read off +`N' (d + e x) - (q - 1) e N = P` from `N(0) = 0`, and beside the logarithm's derivative `b n/x` what is left is +`b n M/(d + e x)^(q - 1)`, a polynomial over a power of the linear. The quotient handed on before held every power of the +linear in the antiderivative multiplied together, of degree twenty-two for the first. Where the degree of `P` is at most +`q - 2`, and only where that quotient does not split over its written factors +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^5*(a+b*ln(c*x^n))/(d+k*x)^7".ToEntity().Integrate("x")` | `integral(...)` | 387 characters | +| `"x^6*(a+b*ln(c*x^n))/(d+k*x)^8".ToEntity().Integrate("x")` | `integral(...)` | 438 characters | +| `"x^5*(a+b*ln(x))/(d+k*x)^7".ToEntity().Integrate("x")` | `integral(...)` | 375 characters | +| `"x^4*(a+b*ln(c*x^n))/(d+k*x)^7".ToEntity().Integrate("x")` | `integral(...)` | 426 characters | + ### A power of x times the logarithm of a polynomial plus a root is integrated by parts over the conjugate **Answers where there were none.** `x^3 ln(-1 + 4 x + 4 sqrt((x - 1) x))` was declined, and so were three more of Rubi's 3.5, @@ -2556,6 +2589,22 @@ answered it wrongly ([#1770](https://github.com/asc-community/AngouriMath/issues |---|---|---| | `"1/((x^2 + 3)^2*(x + sqrt(x^2 + 3)))".ToEntity().Integrate("x")` | `integral(...)` | `-2/((sqrt(x^2 + 3) + x)^2 + 3)^2` | +### A root over the root of the quartic inside it, beside a power of a linear, is split over the quartic's factors + +**Answers where there were none.** Rubi's 1.3.2 `sqrt(2 x^2 + sqrt(3 + 4 x^4))/((c + d x) sqrt(3 + 4 x^4))` and the same +over `(c + d x)^2` were declined, and so were Welz's `sqrt(x^2 + sqrt(1 + x^4))/((1 + x) sqrt(1 + x^4))` and its square. +Where `c^2 = b d^2`, `sqrt(c x^2 + d sqrt(a + b x^4))/sqrt(a + b x^4)` is `sqrt(d/(2 sqrt(b)))` times +`1/sqrt(x^2 + i beta) + 1/sqrt(x^2 - i beta)` on the real line, `beta = sqrt(a/b)`, and each over a power of a linear is +an inverse hyperbolic tangent and roots by the recurrence. The sum is real, and continuous on each side of the pole: the +logarithm of a quotient that the table writes for a real root jumps on the real line for an imaginary one. For positive +numbers `a`, `b`, `c` and `d` ([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(2*x^2+sqrt(3+4*x^4))/((c+d*x)*sqrt(3+4*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 608 characters | +| `"sqrt(2*x^2+sqrt(3+4*x^4))/((c+d*x)^2*sqrt(3+4*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 928 characters | +| `"sqrt(x^2+sqrt(1+x^4))/((1+x)*sqrt(1+x^4))".ToEntity().Integrate("x")` | `integral(...)` | 273 characters | + ### A whole power of a sum of two square roots below the bar is rationalised **Answers where there were none.** `x/(sqrt(a + b x) + sqrt(c + b x))^3` was declined, while the @@ -3466,6 +3515,20 @@ arctangents and two inverse hyperbolic tangents for each sign of `c r` | `"x/((b*x^3 + 2*a*(5 - 3*sqrt(3)))*sqrt(a + b*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 1,703 characters, by the sign of `a` | | `"x/((x^3 + 10 + 6*sqrt(3))*sqrt(1 + x^3))".ToEntity().Integrate("x")` | `integral(...)` | 626 characters | +### A linear over the root of a quartic whose symbols are a scale of x is integrated in the scaled variable + +**Answers where there were none.** Rubi's 1.3.2 `x/sqrt(-44375 b^4 + 576000 b^3 c x + 576000 b^2 c^2 x^2 + 5308416 c^4 x^4)` +was declined, and so was Bronstein's `x/sqrt(x^4 + 10 x^2 - 96 x - 71)` with a scale in it, while each was answered with +numbers for the symbols. Under `x = (b/c) t` the quartic is `5308416 b^4` times a monic one in `t` with rational +coefficients, and the logarithmic ansatz answers that: `ln(p + q sqrt(P))/(18432 c^2)` for `p` of degree eight. The +scaling of the variable took one symbol and did not take the square out from under the root +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x/sqrt(-44375*b^4+576000*b^3*c*x+576000*b^2*c^2*x^2+5308416*c^4*x^4)".ToEntity().Integrate("x")` | `integral(...)` | 520 characters | +| `"x/sqrt(x^4+10*a^2*x^2-96*a^3*x-71*a^4)".ToEntity().Integrate("x")` | `integral(...)` | 312 characters | + ### `x^2` over a three-quarter power of a quadratic binomial beside another is integrated where that is elementary **Answers where there were none.** `x^2/((A + B x^2)^(3/4) (C + D x^2))` at `B C - 2 A D = 0`, Rubi's diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 31dc3156e..194821e16 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -7548,6 +7548,127 @@ private static bool TryReadAsABinomialIn(Entity expr, Entity.Variable x, out Ent return !c.ContainsNode(x) && n >= 1; } + /// + /// A linear over the square root of a quartic whose symbols are a scale of the variable, + /// (alpha x + beta)/sqrt(P(x)) with P(s t) = C M(t), M monic with rational + /// coefficients and C a positive square: under x = s t it is a constant times + /// (t + nu)/sqrt(M(t)) for a rational nu, which the logarithmic ansatz answers + /// where it is elementary. Rubi's 1.3.2 + /// x/sqrt(-44375 b^4 + 576000 b^3 c x + 576000 b^2 c^2 x^2 + 5308416 c^4 x^4), which is + /// ln(p + q sqrt(P))/(18432 c^2) for p of degree eight. + /// + /// + /// The ansatz reads rational coefficients, and the scaling of the variable takes one symbol + /// and reads what is left by simplifying, which does not take C out from under the + /// root, so the row was declined with b = c = 1 answered. Here s is the + /// symbols of the ratio of two adjacent coefficients, each coefficient of M is + /// checked to be a number, and sqrt(C) is written as the monomial it is, its powers + /// multiples of four. The question in t has no symbols, so it is not scaled again. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveALinearOverTheRootOfAQuarticWithItsSymbolsScaledOut(Entity expr, Entity.Variable x, bool integrateByParts) + { + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + Entity? radicand = null; + Entity constantBelow = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + constantBelow = constantBelow == Number.Integer.One ? factor : constantBelow * factor; + continue; + } + if (radicand is not null || factor is not Powf(var @base, Number.Rational half) + || !half.ERational.Equals(ERational.Create(1, 2))) + return null; + radicand = @base; + } + if (radicand is null || !radicand.Vars.Any(symbol => symbol != x) + || !TreeAnalyzer.TryGetPolyLinear(numerator, x, out var alpha, out var beta) + || alpha.ContainsNode(x) || beta.ContainsNode(x) || VanishesIdentically(alpha) + || !TreeAnalyzer.TryGetPolynomial(radicand, x, out var read) + || read.Keys.Any(degree => degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(4)) > 0) + || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + var coefficients = Enumerable.Range(0, 5) + .Select(degree => read.TryGetValue(EInteger.FromInt32(degree), out var c) ? c.InnerSimplified : Number.Integer.Zero).ToArray(); + if (VanishesIdentically(coefficients[4])) + return null; + Entity? ratio = null; + for (var j = 0; j < 4 && ratio is null; j++) + if (!VanishesIdentically(coefficients[j]) && !VanishesIdentically(coefficients[j + 1])) + ratio = Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficients[j] / coefficients[j + 1]); + if (ratio is null) + return null; + // Its symbols alone: the number in it would only rescale t. + var (ratioAbove, ratioBelow) = Functions.SingleQuotient.Of(ratio); + Entity Symbolic(Entity side) => Mulf.LinearChildren(side).Where(factor => factor is not Number) + .Aggregate(Number.Integer.One as Entity, (product, factor) => product == Number.Integer.One ? factor : product * factor); + var (symbolsAbove, symbolsBelow) = (Symbolic(ratioAbove), Symbolic(ratioBelow)); + if (symbolsAbove == Number.Integer.One && symbolsBelow == Number.Integer.One) + return null; + var scale = symbolsBelow == Number.Integer.One ? symbolsAbove : symbolsAbove / symbolsBelow; + var leading = Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficients[4] * MathS.Pow(scale, 4)); + if (RootOfAPositiveSquare(leading) is not { } rootOfLeading) + return null; + var t = Variable.CreateUnique(expr, "t_scale"); + Entity monic = MathS.Pow(t, 4); + for (var j = 3; j >= 0; j--) + { + var normalised = Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficients[j] * MathS.Pow(scale, j) / leading); + // The lowest terms leave `b^4/b^4` standing past the gcd's bounds. + if (normalised.Vars.Any()) + normalised = Functions.PartialFractions.Bare(normalised.Simplify()); + if (normalised.Vars.Any() || normalised.Evaled is not Number.Rational rational) + return null; + if (!rational.IsZero) + monic += j == 0 ? rational : rational * (j == 1 ? t : MathS.Pow(t, j)); + } + // (alpha s t + beta) s dt/(K sqrt(C) sqrt(M)) = alpha s^2/(K sqrt(C)) (t + nu) dt/sqrt(M). + var shift = Functions.PartialFractions.InLowestTermsOverTheSymbols(beta / (alpha * scale)); + if (shift.Vars.Any()) + shift = Functions.PartialFractions.Bare(shift.Simplify()); + if (shift.Vars.Any() || shift.Evaled is not Number.Rational) + return null; + var inT = (shift.Evaled is Number.Complex { IsZero: true } ? (Entity)t : t + shift.Evaled) / MathS.Sqrt(monic); + if (Integration.ComputeAsAQuestionOfItsOwn(inT, t, integrateByParts) is not { } integrated + || integrated.Nodes.Any(node => node is Integralf)) + return null; + var multiple = Functions.PartialFractions.InLowestTermsOverTheSymbols(alpha * MathS.Pow(scale, 2) / (constantBelow * rootOfLeading)); + // The root as the integrand's, over the monomial, where the answer kept it as written. + var root = MathS.Sqrt(monic); + var answer = multiple * integrated.Replace(node => node == root ? MathS.Sqrt(radicand) / rootOfLeading : node) + .Substitute(t, symbolsBelow == Number.Integer.One ? x / symbolsAbove : symbolsBelow * x / symbolsAbove); + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + } + + /// + /// The square root of , a positive rational times powers of + /// symbols, each power a multiple of four, as the rational's root times their halves; + /// null for any other value. + /// + private static Entity? RootOfAPositiveSquare(Entity value) + { + var (top, bottom) = Functions.SingleQuotient.Of(value); + Entity root = Number.Integer.One; + var number = ERational.One; + foreach (var (side, above) in new[] { (top, true), (bottom, false) }) + foreach (var factor in Mulf.LinearChildren(side)) + switch (factor) + { + case Number.Rational rational: + number = above ? number.Multiply(rational.ERational) : number.Divide(rational.ERational); + break; + case Powf(Variable symbol, Number.Integer power) when power.EInteger.Remainder(4).IsZero && power.EInteger.Sign > 0: + var half = MathS.Pow(symbol, Number.Integer.Create(power.EInteger.Divide(2))); + root = above ? root * half : root / half; + break; + default: + return null; + } + return number.Sign > 0 ? (MathS.Sqrt(Number.Rational.Create(number)) * root).InnerSimplified : null; + } + /// /// A linear over a linear, or a quadratic over a quadratic, beside the reciprocal of the /// square root of a cubic binomial a + b x^3, at the coefficients where its @@ -26434,11 +26555,22 @@ when IsAnExponential(@base, inner) && !outer.ContainsNode(x) && outer.Evaled is /// https://github.com/asc-community/AngouriMath/issues/718 /// /// + /// /// (F^u)^n is F^(n u) only where F^u is on the principal branch, which /// it is for a positive F and need not be for a symbol; K carries the /// difference, and the answer holds for every F. Only beside a polynomial below the /// bar, where nothing else reads the power: above it, the rules for a polynomial times an /// exponential answer it as written. + /// + /// + /// And beside another power of an exponential with the same rate, decided by value, which + /// nothing reads from the two as written: Rubi's 2.3 + /// (a + b (F^(e (c + d x)))^n)^p (G^(h (f + g x)))^(d e n ln(F)/(g h ln(G))) is + /// K_2 e^(r x) (a + b K_1 e^(r x))^p, which the substitution u = e^(r x) + /// answers, and was a search past its budget. The constants come back beside the + /// exponential they multiplied and are cancelled against it, and the answer is Rubi's. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// /// internal static Entity? SolveByWritingAPowerOfAnExponentialAsAMultipleOfOne(Entity expr, Entity.Variable x, bool integrateByParts) { @@ -26449,8 +26581,32 @@ when IsAnExponential(@base, inner) && !outer.ContainsNode(x) && outer.Evaled is .Distinct().ToList(); if (powers.Count == 0) return null; + var rates = new List(); + foreach (var node in powers) + { + if (node is not Powf(Powf(var @base, var inner), var outer) + || !TreeAnalyzer.TryGetPolyLinear(inner, x, out var slope, out _) || TreeAnalyzer.IsZero(slope)) + return null; + rates.Add(outer.Evaled is Number.Integer ? null + : (outer * slope * (@base == MathS.e ? Number.Integer.One : MathS.Ln(@base))).InnerSimplified); + } + // Two powers written apart whose rates are one, read as multiples of one exponential: + // Rubi's `(a + b (F^(e (c + d x)))^n)^p (G^(h (f + g x)))^(d e n ln(F)/(g h ln(G)))`, + // where the second is a constant multiple of the first on each interval and nothing + // reads that from the two as written. Decided by value, each later rate then written as + // the first that equals it. + var ofOneRate = false; + for (var later = 1; later < rates.Count; later++) + for (var earlier = 0; earlier < later; earlier++) + if (rates[later] is { } laterRate && rates[earlier] is { } earlierRate + && (laterRate == earlierRate || VanishesIdentically(laterRate - earlierRate))) + { + rates[later] = earlierRate; + ofOneRate = true; + break; + } // A polynomial below the bar, or under a negative whole power. - if (!FactorsOfTheIntegrand(expr).Any(pair => + if (!ofOneRate && !FactorsOfTheIntegrand(expr).Any(pair => pair.Factor.ContainsNode(x) && !pair.Factor.Nodes.Any(IsAnExponential) && (pair.Underneath || pair.Factor is Powf(_, Number.Integer { IsNegative: true })) && TreeAnalyzer.TryGetPolynomial(pair.Factor is Powf(var raised, Number.Integer) ? raised : pair.Factor, x, out var read) @@ -26458,18 +26614,17 @@ when IsAnExponential(@base, inner) && !outer.ContainsNode(x) && outer.Evaled is return null; var back = new List<(Entity.Variable Constant, Entity Value)>(); var rewritten = expr; - foreach (var node in powers) + for (var at = 0; at < powers.Count; at++) { - if (node is not Powf(Powf(var @base, var inner), var outer) - || !TreeAnalyzer.TryGetPolyLinear(inner, x, out var slope, out _) || TreeAnalyzer.IsZero(slope)) + var node = powers[at]; + if (node is not Powf(Powf(var @base, var inner), var outer)) return null; Entity replacement; - if (outer.Evaled is Number.Integer) + if (rates[at] is not { } rate) replacement = MathS.Pow(@base, (outer * inner).InnerSimplified); else { var constant = Variable.CreateUnique(rewritten, "k_pow"); - var rate = (outer * slope * (@base == MathS.e ? Number.Integer.One : MathS.Ln(@base))).InnerSimplified; replacement = constant * MathS.Pow(MathS.e, rate * x); back.Add((constant, node * MathS.Pow(MathS.e, -rate * x))); } @@ -26479,9 +26634,46 @@ when IsAnExponential(@base, inner) && !outer.ContainsNode(x) && outer.Evaled is return null; foreach (var (constant, value) in back) answer = answer.Substitute(constant, value); + if (!ofOneRate) + return answer; + // Each constant came back as its power times `e^(-r x)`, beside the `e^(r x)` it + // multiplied: the two are cancelled as written, a product at a time and then across + // the bar, which an exponential, never zero, allows. + foreach (var rate in rates.OfType().Distinct()) + answer = WithOppositeExponentialsCancelled(answer, MathS.Pow(MathS.e, rate * x), MathS.Pow(MathS.e, -rate * x)); return answer; } + /// + /// with and , an + /// exponential and its reciprocal, cancelled where they stand as factors of one product, + /// and then, for a product or quotient, each cancelled where it stands above and below + /// the bar. Exact, since an exponential is never zero. + /// + private static Entity WithOppositeExponentialsCancelled(Entity expr, Entity up, Entity down) + { + var cancelled = expr.Replace(node => + { + if (node is not Mulf) + return node; + var factors = Mulf.LinearChildren(node).ToList(); + var pairs = System.Math.Min(factors.Count(factor => factor == up), factors.Count(factor => factor == down)); + if (pairs == 0) + return node; + for (var removed = 0; removed < pairs; removed++) + { + factors.Remove(up); + factors.Remove(down); + } + return factors.Count == 0 ? Number.Integer.One : factors.Aggregate((product, factor) => product * factor); + }); + if (cancelled is Sumf or Minusf || !cancelled.Nodes.Any(node => node == up || node == down)) + return cancelled; + var (above, below) = Functions.SingleQuotient.Of(cancelled); + var acrossTheBar = CancelCommonFactors(above, below); + return acrossTheBar.Nodes.Any(node => node == up || node == down) ? cancelled : acrossTheBar; + } + /// /// An exponential of a multiple of a logarithm is a power of the argument: /// e^(k ln(q)) is q^k, since e^(k ln q) is the definition of the @@ -26969,7 +27161,7 @@ internal static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) /// with a^2 = b^2 c; sqrt(c x^2 + d sqrt(a + b x^4))/sqrt(a + b x^4) with /// c^2 = b d^2; 1/((a + b x^n) sqrt(c x^2 + d (a + b x^n)^(2/n))); and /// sqrt(a x^2 + b x sqrt(c + d x^2))/(x sqrt(c + d x^2)) with a^2 = b^2 d and - /// b^2 c + a = 0. + /// b^2 c + a = 0. And the second over a power of a linear, for numbers. /// /// /// @@ -26983,6 +27175,19 @@ internal static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) /// Timofeev suites declined all of them. /// https://github.com/asc-community/AngouriMath/issues/718 /// + /// + /// Over a power of a linear L the substitution leaves x as a fourth root in + /// w; there the quartic is split instead. With beta = sqrt(a/b), + /// a + b x^4 = b (x^2 + i beta)(x^2 - i beta), and on the real line the two roots of + /// the factors are conjugates whose sum is sqrt(2 x^2 + 2 sqrt(x^4 + beta^2)), so + /// sqrt(S)/R is sqrt(d/(2 sqrt(b))) times + /// 1/sqrt(x^2 + i beta) + 1/sqrt(x^2 - i beta), each over L^k an inverse + /// hyperbolic tangent and roots, and the two together real. For positive numbers + /// a, b, c, d: Rubi's 1.3.2 + /// sqrt(2 x^2 + sqrt(3 + 4 x^4))/((c + d x) sqrt(3 + 4 x^4)) and its square, and + /// Welz's 47 and 48, declined before. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// /// internal static Entity? SolveANestedRootByItsClosedForm(Entity expr, Entity.Variable x, bool integrateByParts) { @@ -27013,7 +27218,7 @@ internal static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) { 1 => ARootOfAConstantPlusARootOfAQuadratic(factors[0]), 2 => ARootOverTheRootInsideIt(factors) ?? AReciprocalOfTheRootInsideIt(factors), - _ => AProductRootOverTheRootInsideIt(factors), + _ => AProductRootOverTheRootInsideIt(factors) ?? ARootOverTheRootInsideItBesideAPowerOfALinear(factors), }; if (answer is null) return null; @@ -27131,6 +27336,58 @@ internal static bool IsAPerfectSquareDiscriminant(Entity a, Entity b, Entity c) return MathS.Sqrt(2) * b / a * (q is Variable named ? inX.Substitute(named, reciprocal) : inX); } + // sqrt(c x^2 + d sqrt(P))/(sqrt(P) L^k), P = a + b x^4 with numbers, a/b > 0 and + // c = d sqrt(b) > 0: with beta = sqrt(a/b), P = b (x^2 + i beta)(x^2 - i beta), and + // on the real line sqrt(x^2 + i beta) + sqrt(x^2 - i beta) is the root of + // 2 x^2 + 2 sqrt(x^4 + beta^2), so the integrand is + // sqrt(d/(2 sqrt(b))) (1/sqrt(x^2 + i beta) + 1/sqrt(x^2 - i beta))/L^k. + Entity? ARootOverTheRootInsideItBesideAPowerOfALinear(List three) + { + if (three.Find(f => f is Powf(_, var p) && p == half) is not Powf(var sum, _) + || three.Find(f => f is Powf(_, var p) && p == -half) is not Powf(var radicand, _) + || three.Find(f => f is Powf(var l, Number.Integer { EInteger.Sign: < 0 }) && TreeAnalyzer.TryGetPolyLinear(l, x, out _, out _)) + is not Powf(var linear, Number.Integer power) + || power.EInteger.CompareTo(EInteger.FromInt32(-4)) < 0 + || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var offset) + // A power of x alone is the substitution's, which answers it shorter. + || VanishesIdentically(slope) || VanishesIdentically(offset) + || Functions.PartialFractions.HoldsTheImaginaryUnit(slope) || Functions.PartialFractions.HoldsTheImaginaryUnit(offset) + || Coefficients(radicand, 0, 4) is not [var a, var b] + || SquareOfXPlusMultipleOfRoot(sum, radicand, half) is not var (c, d) + || a.Evaled is not Number.Real { IsPositive: true } || b.Evaled is not Number.Real { IsPositive: true } + || c.Evaled is not Number.Real { IsPositive: true } || d.Evaled is not Number.Real { IsPositive: true } + || !Functions.PartialFractions.IsZeroAsAValue(c * c - b * d * d)) + return null; + var beta = MathS.Sqrt(a / b).InnerSimplified; + var multiple = MathS.Sqrt(d / (2 * MathS.Sqrt(b))).InnerSimplified; + var k = (-power.EInteger).ToInt32Checked(); + return multiple * (OverAPowerOfALinearBesideARoot(MathS.i * beta) + OverAPowerOfALinearBesideARoot(-MathS.i * beta)); + + // The integral of 1/(L^k sqrt(x^2 + q)), L = offset + slope x, by + // I_(m + 1) = ((1 - m) I_(m - 1) + (2m - 1) offset I_m - slope sqrt(x^2 + q)/L^m)/(m D), + // D = offset^2 + slope^2 q, from I_0 = ln(x + sqrt(x^2 + q)) and + // I_1 = -artanh((q slope - offset x)/(sqrt(D) sqrt(x^2 + q)))/sqrt(D). For an + // imaginary q these are continuous on each side of the pole: one minus the square + // of the artanh's argument is q L^2/(D (x^2 + q)), which is a negative number nowhere + // on the real line. The table's ln((A + B)/(A - B)) is twice the inverse hyperbolic + // cotangent of the same argument, whose cut is where that is a positive number, and + // for `2 x^2 + sqrt(3 + 4 x^4)` over `0.7 + 1.3 x` it jumped at x = -1.7. + Entity OverAPowerOfALinearBesideARoot(Entity q) + { + var root = MathS.Sqrt(MathS.Pow(x, 2) + q); + var discriminant = (offset * offset + slope * slope * q).InnerSimplified; + var rootOfIt = MathS.Sqrt(discriminant); + Entity before = MathS.Ln(x + root); + Entity current = -MathS.Hyperbolic.Artanh((q * slope - offset * x) / (rootOfIt * root)) / rootOfIt; + for (var m = 1; m < k; m++) + { + var next = ((1 - m) * before + (2 * m - 1) * offset * current - slope * root / MathS.Pow(linear, m)) / (m * discriminant); + (before, current) = (current, next); + } + return current; + } + } + // multiple times the integral of 1/(1 - q w^2) at w = at, with q named while it is asked. Entity? InW(Entity multiple, Entity coefficient, Entity at) { @@ -29084,6 +29341,12 @@ static bool AreEqualAsPolynomials(Entity left, Entity right) /// a constant has been taken out in front; against a constant alone it is by parts /// against one, asked at the top or one below it like the rest of by parts. /// https://github.com/asc-community/AngouriMath/issues/718 + /// Where the antiderivative times the derivative of A + B ln(C x^n) does not split + /// over its written factors and the rational function is a polynomial over a power of one + /// linear, the antiderivative is the one that vanishes at zero, as one quotient: Rubi's 3.1.4 + /// x^5 (a + b ln(c x^n))/(d + e x)^7 left a remainder of degree twenty-two that the + /// substitution search did not finish. + /// https://github.com/asc-community/AngouriMath/issues/233 /// internal static Entity? SolveARationalFunctionTimesAPowerOfALogarithm(Entity expr, Entity.Variable x, bool integrateByParts) { @@ -29208,8 +29471,28 @@ static bool AreEqualAsPolynomials(Entity left, Entity right) // and each on its own beside the power of the logarithm is a dilogarithm. var remainders = new List(); if (power == 1) + { + var splitEach = true; foreach (var term in derivativeTerms) - remainders.Add(AsPartialFractions(term * integralOfTheRational, x)); + { + remainders.Add(AsPartialFractions(term * integralOfTheRational, x, out var split)); + splitEach &= split; + } + // A product that does not split was handed on as one quotient over every + // denominator of the antiderivative's terms multiplied together: against + // `a + b ln(c x^n)`, `x^5/(d + k x)^7` has six, and the remainder was a quotient + // of degree twenty-two with symbols in it, which the substitution search did + // not finish in forty seconds. Where the logarithm's derivative is a constant over + // `x` and the rational function a polynomial over one power of a linear, the + // antiderivative vanishing at zero is `x M/(s + t x)^(q - 1)` with `M` a + // polynomial, and the remainder is that constant times `M/(s + t x)^(q - 1)`. + if (!splitEach && TheConstantOverX(derivativeTerms, x) is { } overX + && TheAntiderivativeVanishingAtZeroOverAPowerOfALinear(rational, x) is var (vanishing, vanishingOverX)) + { + integralOfTheRational = vanishing; + remainders = new List { overX * vanishingOverX }; + } + } else { // Keyed by the part of the denominator holding x; the constants go with the @@ -29219,7 +29502,7 @@ static bool AreEqualAsPolynomials(Entity left, Entity right) // with nothing cancelled. var collected = new List<(Entity Below, Entity Above)>(); foreach (var term in derivativeTerms) - foreach (var fraction in TermsDistributed(AsPartialFractions(term * integralOfTheRational, x))) + foreach (var fraction in TermsDistributed(AsPartialFractions(term * integralOfTheRational, x, out _))) { var (fractionAbove, fractionBelow) = Functions.SingleQuotient.Of(fraction); Entity holdingX = Number.Integer.One, constantBelow = Number.Integer.One; @@ -29289,14 +29572,17 @@ static Entity Distributed(Entity inner, Number.Integer exponent) /// brought over one bar and written as its whole part plus /// its partial fractions over the written linear factors of the denominator, the /// symbolic ones by the derivatives at the roots; over one - /// bar as it came where the denominator is not such a product. + /// bar as it came where the denominator is not such a product, and then + /// is false. /// - private static Entity AsPartialFractions(Entity product, Entity.Variable x) + private static Entity AsPartialFractions(Entity product, Entity.Variable x, out bool split) { + split = true; var combined = Functions.SingleQuotient.Combine(product).InnerSimplified; var (above, below) = Functions.SingleQuotient.Of(combined); if (below == Number.Integer.One || !below.ContainsNode(x)) return combined; + split = false; above = Functions.PartialFractions.Bare(above.Expand().InnerSimplified); if (!TreeAnalyzer.TryGetPolynomial(above, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(below, x, out var belowRead) || aboveRead.Count == 0 || belowRead.Count == 0) @@ -29315,7 +29601,10 @@ private static Entity AsPartialFractions(Entity product, Entity.Variable x) (proper, below) = Functions.SingleQuotient.Of(Functions.PartialFractions.Bare(remainder.InnerSimplified)); } if (proper == Number.Integer.Zero || proper.Evaled is Number.Complex { IsZero: true }) + { + split = true; return whole; + } Entity fractions; if (below == Number.Integer.One) fractions = proper; @@ -29328,7 +29617,8 @@ private static Entity AsPartialFractions(Entity product, Entity.Variable x) : decomposition; } else - fractions = proper / below; + return whole == Number.Integer.Zero ? proper / below : whole + proper / below; + split = true; return whole == Number.Integer.Zero ? fractions : whole + fractions; } @@ -29383,6 +29673,108 @@ private static Entity MatchedToAPoleOf(Entity antiderivative, List deriv return antiderivative; } + /// + /// The constant c where the sum of is + /// c/x, each of them a constant over , as the derivative of + /// A + B ln(C x^n) is; null otherwise. + /// + private static Entity? TheConstantOverX(List derivativeTerms, Entity.Variable x) + { + Entity? sum = null; + foreach (var term in derivativeTerms) + { + var (above, below) = Functions.SingleQuotient.Of(term); + if (above.ContainsNode(x)) + return null; + Entity constant = above; + var seenX = false; + foreach (var factor in Mulf.LinearChildren(below)) + if (factor == x && !seenX) + seenX = true; + else if (factor.ContainsNode(x)) + return null; + else + constant /= factor; + if (!seenX) + return null; + sum = sum is null ? constant : sum + constant; + } + return sum is null || VanishesIdentically(sum) ? null : sum; + } + + /// + /// For P/(K (s + t x)^q) with P a polynomial of + /// degree at most q - 2, K, s and t free of + /// and s not zero: its antiderivative that vanishes at zero, + /// x M/(K (s + t x)^(q - 1)) with M a polynomial, and that over + /// ; null for any other . + /// + /// + /// N = x M is read off N' (s + t x) - (q - 1) t N = P coefficient by + /// coefficient: (j + 1) s n_(j + 1) = p_j + (q - 1 - j) t n_j from n_0 = 0, + /// which is the constant that makes it vanish at zero, and the equation at + /// x^(q - 1) is p_(q - 1) = 0, the degree bound, which says that there is + /// no logarithm in it. Each coefficient is put in lowest terms over the symbols as it is + /// made, so that those of the lower powers that are zero are written as zero: for + /// x^5/(d + k x)^7 only n_6 = 1/(6 d) is not, and the antiderivative is + /// x^6/(6 d (d + k x)^6), where the integral of the rational function, matched to + /// the pole, was six powers of the linear and a constant. + /// + private static (Entity Antiderivative, Entity OverX)? TheAntiderivativeVanishingAtZeroOverAPowerOfALinear(Entity rational, Entity.Variable x) + { + var (above, below) = Functions.SingleQuotient.Of(rational); + Entity? linear = null; + var power = 0; + Entity constantBelow = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(below)) + { + if (!factor.ContainsNode(x)) + { + constantBelow = constantBelow == Number.Integer.One ? factor : constantBelow * factor; + continue; + } + if (linear is not null || factor is not Powf(var @base, Number.Integer { EInteger: var exponent }) + || exponent.CompareTo(EInteger.FromInt32(2)) < 0 || exponent.CompareTo(EInteger.FromInt32(24)) > 0) + return null; + linear = @base; + power = exponent.ToInt32Checked(); + } + if (linear is null || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var offset) + || VanishesIdentically(slope) || VanishesIdentically(offset) + || !TreeAnalyzer.TryGetPolynomial(above, x, out var read) + || read.Keys.Any(degree => degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(power - 2)) > 0) + || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + var coefficients = new Entity[power]; + coefficients[0] = Number.Integer.Zero; + for (var j = 0; j < power - 1; j++) + { + var lower = read.TryGetValue(EInteger.FromInt32(j), out var p) ? p : Number.Integer.Zero; + coefficients[j + 1] = lower == Number.Integer.Zero && coefficients[j] == Number.Integer.Zero + ? Number.Integer.Zero + : Functions.PartialFractions.InLowestTermsOverTheSymbols( + (lower + Number.Integer.Create(power - 1 - j) * slope * coefficients[j]) / (Number.Integer.Create(j + 1) * offset)); + } + Entity polynomial = Number.Integer.Zero; + for (var j = 1; j < power; j++) + { + var coefficient = coefficients[j]; + if (coefficient == Number.Integer.Zero || coefficient.Evaled is Number.Complex { IsZero: true }) + continue; + var term = j == 1 ? coefficient : j == 2 ? coefficient * x : coefficient * MathS.Pow(x, j - 1); + polynomial = polynomial == Number.Integer.Zero ? term : polynomial + term; + } + if (polynomial == Number.Integer.Zero) + return null; + Entity overPower = MathS.Pow(linear, power - 1); + if (constantBelow != Number.Integer.One) + overPower = constantBelow * overPower; + var overX = polynomial / overPower; + var antiderivative = x * polynomial / overPower; + return Functions.PartialFractions.HoldsAtSampledPoints(antiderivative.Differentiate(x), rational, x) + ? (antiderivative, overX) : null; + } + /// /// The terms of with every product distributed over the sums /// among its factors: B (ln(a + b x) - ln(c + d x)) is two terms. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index b771c2129..c4cff24d8 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -836,8 +836,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of an exponential with a positive base is the exponential of the product, // exactly, and only that spelling is one the exponential rules read. if ((answer = IndefiniteIntegralSolver.SolveByFlatteningAPowerOfAnExponential(expr, x, integrateByParts)) is { }) return answer; - // And any other power of an exponential beside a polynomial below the bar, as a constant - // multiple of an exponential wherever it is differentiable. + // And any other power of an exponential beside a polynomial below the bar, or beside + // another of the same rate, as a constant multiple of an exponential wherever it is + // differentiable. if ((answer = IndefiniteIntegralSolver.SolveByWritingAPowerOfAnExponentialAsAMultipleOfOne(expr, x, integrateByParts)) is { }) return answer; // An exponential of a quadratic in 1/L beside a power of L, onto the Gaussian under u = 1/L. if ((answer = IndefiniteIntegralSolver.SolveAGaussianInAReciprocal(expr, x)) is { }) return answer; @@ -1522,6 +1523,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // is wrong with -- it clears a parameter rather than a shape -- so everything that // can answer the problem as written gets to try first. if ((answer = IndefiniteIntegralSolver.SolveByScalingTheVariable(expr, x, integrateByParts)) is { }) return answer; + // And a linear over the root of a quartic whose symbols are a scale of x, where that + // takes them out from under the root, which the simplifier does not: the logarithmic + // ansatz answers the numeric quartic that is left. + if ((answer = IndefiniteIntegralSolver.SolveALinearOverTheRootOfAQuarticWithItsSymbolsScaledOut(expr, x, integrateByParts)) is { }) return answer; // After every rule that can answer a power in its own terms, because expanding one // throws away whatever structure it had: (1 + x^2)^2 is answered as a power and only // wants writing out if that fails. diff --git a/Sources/Tests/UnitTests/Calculus/ALinearOverTheRootOfAQuarticWithItsSymbolsScaledOutIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ALinearOverTheRootOfAQuarticWithItsSymbolsScaledOutIntegralTest.cs new file mode 100644 index 000000000..7d4748426 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ALinearOverTheRootOfAQuarticWithItsSymbolsScaledOutIntegralTest.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 linear over the root of a quartic whose symbols are a scale of the variable, read under + /// x = s t as a numeric quartic that the logarithmic ansatz answers: Rubi's 1.3.2 + /// x/sqrt(-44375 b^4 + 576000 b^3 c x + 576000 b^2 c^2 x^2 + 5308416 c^4 x^4), and + /// Bronstein's x/sqrt(x^4 + 10 x^2 - 96 x - 71) with a scale in it, declined before. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ALinearOverTheRootOfAQuarticWithItsSymbolsScaledOutIntegralTest + { + [Theory] + [InlineData("x/sqrt(-44375*b^4 + 576000*b^3*c*x + 576000*b^2*c^2*x^2 + 5308416*c^4*x^4)")] + [InlineData("x/sqrt(-44375*b^4 + 576000*b^3*x + 576000*b^2*x^2 + 5308416*x^4)")] + [InlineData("x/sqrt(x^4 + 10*b^2*x^2 - 96*b^3*x - 71*b^4)")] + [InlineData("x/sqrt(c^4*x^4 + 10*b^2*c^2*x^2 - 96*b^3*c*x - 71*b^4)")] + public void IsIntegrated(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 1500, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("b", -0.7).Substitute("c", 1.3); + var antiderivative = Pinned(integral.Substitute("C", 0)); + var derivative = antiderivative.Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -4.1, -2.3, -1.6, 1.1, 2.2, 3.7 }) + { + 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) + 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 >= 4, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfALinearTimesALogarithmIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfALinearTimesALogarithmIntegralTest.cs new file mode 100644 index 000000000..50afb1657 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfALinearTimesALogarithmIntegralTest.cs @@ -0,0 +1,56 @@ +// +// 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 polynomial over a power of a linear times A + B ln(C x^n), by parts with the + /// antiderivative of the rational function that vanishes at zero written as + /// x M(x)/(d + k x)^(q - 1): the remainder is then a polynomial over a power of the + /// linear. Rubi's 3.1.4 x^5 (a + b ln(c x^n))/(d + e x)^7 handed the substitution + /// search a quotient of degree twenty-two and ran past forty seconds. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ARationalFunctionOverAPowerOfALinearTimesALogarithmIntegralTest + { + [Theory] + [InlineData("x^5*(a + b*ln(c*x^n))/(d + k*x)^7")] + [InlineData("x^6*(a + b*ln(c*x^n))/(d + k*x)^8")] + [InlineData("x^5*(a + b*ln(x))/(d + k*x)^7")] + [InlineData("x^4*(a + b*ln(c*x^n))/(d + k*x)^7")] + [InlineData("x^5*(a + b*ln(c*x^n))/(d + k*x)^8")] + public void IsIntegrated(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 1500, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.6) + .Substitute("d", 1.1).Substitute("k", -0.9).Substitute("n", 2); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -2.3, -0.7, 0.3, 0.8, 2.9, 4.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}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/ARootOverTheRootInsideItBesideAPowerOfALinearIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARootOverTheRootInsideItBesideAPowerOfALinearIntegralTest.cs new file mode 100644 index 000000000..59c0d92d7 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARootOverTheRootInsideItBesideAPowerOfALinearIntegralTest.cs @@ -0,0 +1,74 @@ +// +// 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 +{ + /// + /// sqrt(c x^2 + d sqrt(a + b x^4))/(sqrt(a + b x^4) (p + q x)^k) with c^2 = b d^2, + /// written as sqrt(d/(2 sqrt(b))) times 1/sqrt(x^2 + i beta) + 1/sqrt(x^2 - i beta) + /// over the power, beta = sqrt(a/b): Rubi's 1.3.2 + /// sqrt(2 x^2 + sqrt(3 + 4 x^4))/((c + d x) sqrt(3 + 4 x^4)) and its square, declined before. + /// #233 + /// + /// + /// The answer is checked by differentiating back and also against Simpson's rule on an + /// interval to the left of the pole: the logarithm of a quotient, which the table writes + /// for a real root, jumps there for an imaginary one, and only its value says so. + /// + [Trait("Area", "Calculus")] + public sealed class ARootOverTheRootInsideItBesideAPowerOfALinearIntegralTest + { + [Theory] + [InlineData("sqrt(2*x^2 + sqrt(3 + 4*x^4))/((c + d*x)*sqrt(3 + 4*x^4))")] + [InlineData("sqrt(2*x^2 + sqrt(3 + 4*x^4))/((c + d*x)^2*sqrt(3 + 4*x^4))")] + [InlineData("sqrt(2*x^2 + sqrt(3 + 4*x^4))/((c + d*x)^3*sqrt(3 + 4*x^4))")] + [InlineData("sqrt(3*x^2 + sqrt(2 + 9*x^4))/((c + d*x)*sqrt(2 + 9*x^4))")] + [InlineData("sqrt(2*x^2 + sqrt(3 + 4*x^4))/((1 + 2*x)*sqrt(3 + 4*x^4))")] + public void IsIntegrated(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 4000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("c", 0.7).Substitute("d", 1.3); + var antiderivative = Pinned(integral.Substitute("C", 0)); + var derivative = antiderivative.Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -3.1, -1.7, -0.9, 0.2, 1.3, 4.6 }) + { + 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}"); + + // On [-3, -1], left of every pole here, F(-1) - F(-3) is the integral. + double At(Entity e, double point) => (double)e.Substitute("x", point).EvalNumerical().RealPart; + const int steps = 200; + var (from, to) = (-3.0, -1.0); + var width = (to - from) / steps; + var simpson = At(original, from) + At(original, to); + for (var step = 1; step < steps; step++) + simpson += (step % 2 == 1 ? 4 : 2) * At(original, from + step * width); + simpson *= width / 3; + var difference = antiderivative.Substitute("x", to).EvalNumerical() - antiderivative.Substitute("x", from).EvalNumerical(); + Assert.True(Math.Abs((double)difference.ImaginaryPart) < 1e-9, $"F(-1) - F(-3) is {difference} for {integrand}"); + Assert.True(Math.Abs((double)difference.RealPart - simpson) < 1e-6 * Math.Max(1, Math.Abs(simpson)), + $"F(-1) - F(-3) is {difference} for {integrand}, where Simpson's rule gives {simpson}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/PowersOfExponentialsOfOneRateIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/PowersOfExponentialsOfOneRateIntegralTest.cs new file mode 100644 index 000000000..776709de0 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/PowersOfExponentialsOfOneRateIntegralTest.cs @@ -0,0 +1,61 @@ +// +// 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 +{ + /// + /// Two powers of exponentials written apart whose logarithmic derivatives are one constant, + /// each read as a constant multiple of the one exponential: Rubi's 2.3 + /// (a + b (F^(e (c + d x)))^n)^p (G^(h (f + g x)))^(d e n ln(F)/(g h ln(G))), which the + /// derivative-of-a-product search ran past its budget on, is + /// (G^(h (f + g x)))^(...) (a + b (F^(e (c + d x)))^n)^(p + 1)/(b d e n (p + 1) ln(F) (F^(e (c + d x)))^n). + /// #233 + /// + /// + /// Pinned to moderate values: with every symbol near 5 the powers pass the range the + /// evaluation keeps, and both sides read as infinite. + /// + [Trait("Area", "Calculus")] + public sealed class PowersOfExponentialsOfOneRateIntegralTest + { + [Theory] + [InlineData("(a + b*(F^(k*(c + d*x)))^n)^p*(G^(h*(f + g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))", false)] + [InlineData("(G^(h*(f + g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))/(a + b*(F^(k*(c + d*x)))^n)^2", false)] + [InlineData("(G^(h*(f + g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))*sqrt(a + b*(F^(k*(c + d*x)))^n)", false)] + [InlineData("(a + b*(F^(k*(c + d*x)))^n)^p*(F^(h*(f + g*x)))^(d*k*n/(g*h))", false)] + [InlineData("(a + b*(F^(k*(c + d*x)))^n)^p*(G^(h*(f + g*x)))^(d*k*n*ln(F)/(g*h*ln(G)))", true)] + public void IsIntegrated(string integrand, bool negativeBases) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 1500, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("F", negativeBases ? -1.9 : 1.9) + .Substitute("k", 0.4).Substitute("c", 0.3).Substitute("d", 0.6).Substitute("n", 1.7).Substitute("p", 0.45) + .Substitute("G", negativeBases ? -2.6 : 2.6).Substitute("h", 0.35).Substitute("f", -0.2).Substitute("g", 0.8); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -2.3, -1.1, -0.4, 0.3, 0.9, 1.7 }) + { + 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}"); + } + } +}