From 6d75aa6d0e4eb828a93c10d65c3d8172ac161e70 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 21 Oct 2026 16:04:55 +0000 Subject: [PATCH] Integrate a rational function of x/sqrt(P) for an even quartic P, high powers of two linears, and nine more algebraic shapes Eleven shapes of Rubi's family 1, each answered where the rules declined or ran out of time: - a rational function beside a half-odd power of P = a + b x^2 + c x^4, where it is one of t = x/sqrt(P), since dt/dx = (a - c x^4)/P^(3/2); - (c x)^t, t not whole, as x^t beside (c x)^t x^(-t), constant wherever defined; - two roots whose radicands' quotient is a constant times an even power of x, as one root up to a factor whose square is constant; - a polynomial over high whole powers of two linears, under t = L1/L2; - an algebraic function of 1/x alone, under x = 1/u; - x^m (a + b x^n)^p, m and n whole and p <= -2, by the reduction to p = -1; - a power of a sum of powers of x beside a symbolic power of a binomial, multiplied out into binomial differentials; - a difference of radicands that is one power of x, written as one; - a rational function beside an odd root of a linear past the square, split first where a linear below meets the root at an irrational value; - an odd integrand beside a high power of a quadratic left to u = x^2, and a quartic in one shifted linear to u = d + k x. Chain against master f7f2a7a8: pocket of 1,079 rows of family 1, 1,013 -> 1,038 solved, 0 wrong; samples unchanged but family 1's 40-a-file, 1,375 -> 1,379; the 25 rows that moved, alone, 0 -> 25 solved, each in 18-772 ms; no row's allocation grew past 1.3x. x^2/((1 - a x)^16 (1 + a x)^11), which master answers in 20 s and 8.5 GB as 16,292 characters (exactly right; the probe's double evaluation calls it wrong), is 1,380 in 0.12 s. Part of #233. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 40 ++ .../Integration/IndefiniteIntegralSolver.cs | 669 +++++++++++++++++- .../Integration/Integration.Definition.cs | 25 + ...lToAWholePowerBelowMinusOneIntegralTest.cs | 51 ++ ...erOfAMultipleOfXAsAPowerOfXIntegralTest.cs | 53 ++ ...tOfSumsBesideASymbolicPowerIntegralTest.cs | 52 ++ .../AQuarticInOneShiftedLinearIntegralTest.cs | 49 ++ ...onBesideAnIrrationalOddRootIntegralTest.cs | 53 ++ ...ootsWithAMonomialDifferenceIntegralTest.cs | 49 ++ ...cFunctionOfTheReciprocalOfXIntegralTest.cs | 49 ++ ...onalFunctionBesideAPowerOfXIntegralTest.cs | 51 ++ ...HighWholePowersOfTwoLinearsIntegralTest.cs | 51 ++ ...eOverTheRootOfAnEvenQuarticIntegralTest.cs | 62 ++ ...otsAsOneUpToAConstantFactorIntegralTest.cs | 53 ++ 14 files changed, 1305 insertions(+), 2 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/ABinomialToAWholePowerBelowMinusOneIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/APowerOfAMultipleOfXAsAPowerOfXIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/AProductOfSumsBesideASymbolicPowerIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/AQuarticInOneShiftedLinearIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAnIrrationalOddRootIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/ASumOfRootsWithAMonomialDifferenceIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/AnAlgebraicFunctionOfTheReciprocalOfXIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/AnOddRationalFunctionBesideAPowerOfXIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/HighWholePowersOfTwoLinearsIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/TheVariableOverTheRootOfAnEvenQuarticIntegralTest.cs create mode 100644 Sources/Tests/UnitTests/Calculus/TwoRootsAsOneUpToAConstantFactorIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 0909ba42e..b59eb4d7d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -2963,6 +2963,46 @@ without a common factor | `"(c + d*tan(x))^(3/2)/(a + i*a*tan(x))^2".ToEntity().Integrate("x")` | `integral(...)`; 470,119 characters on the unreleased master | 77,377 characters | | `"(c + d*tan(x))^(1/2)/(a + i*a*tan(x))^2".ToEntity().Integrate("x")` | `integral(...)`; 258,943 characters on the unreleased master | 71,266 characters | +### A rational function of `x/sqrt(P)` for an even quartic `P`, high powers of two linears, and nine more algebraic shapes are integrated + +**Answers where there were none, and shorter ones.** Rubi's algebraic integrands of these shapes were declined or ran +out of time ([#233](https://github.com/asc-community/AngouriMath/issues/233)): +- A rational function beside a half-odd power of `P = a + b x^2 + c x^4`, where it is a rational function of + `t = x/sqrt(P)`: `dt/dx = (a - c x^4)/P^(3/2)`, so `(a - c x^4)/((a d + a e x^2 + c d x^4) sqrt(P))` is + `dt/(d + (a e - b d) t^2)`. Which integrands are is read in `y = x^2` modulo `y^2 - w y + a/c`, `w = y + a/(c y)`; + an elliptic integral does not pass. +- `(c x)^t`, `t` not whole, written `x^t` beside `(c x)^t x^(-t)`, whose derivative is zero wherever it is defined: + `(c x)^(-1 - 2 n - n p) (a + b x^n)^p`. +- Two square roots whose radicands have a constant times an even power of x for their quotient are one root, up to a + factor whose square is constant: `(g + h x) sqrt(a + b x + c x^2)/(a d + b d x + c d x^2)^(3/2)`. +- A polynomial over high whole powers of two linears, under `t = L1/L2`: `x^2/((1 - a x)^16 (1 + a x)^11)`. +- An algebraic function of `1/x` alone, under `x = 1/u`: `1/((a + b/x)^(5/2) (c + d/x)^2)`. +- `x^m (a + b x^n)^p` for whole `m` and `n` and `p ≤ -2`, by the reduction that raises `p` to `-1`: + `1/(x^2 (a + b x^3)^5)`. +- A power of a sum of powers of x beside a symbolic power of a binomial, multiplied out into binomial differentials: + `(a + b x^n)^3 (c + d x^n)^(-4 - 1/n)`. +- A difference of radicands that is one power of x, written as one: `x^4/(sqrt(a + b x) + sqrt(a + c x))^3`. +- A rational function beside an odd root of a linear past the square, split first where a linear below meets the root + at an irrational value: `1/(x^3 (1 - x^2)^(1/3) (3 + x^2))`. +- An odd integrand beside a high power of a quadratic is left to `u = x^2`, and a quartic in one shifted linear to + `u = d + k x`, where the partial fractions spent the budget or wrote pages. + +Answers the unreleased master gave change length: `x^2/((1 - a x)^16 (1 + a x)^11)` was 16,292 characters there, +after twenty seconds, and is 1,380 in a tenth of a second; `1/((a + b/x)^(5/2) (c + d/x)^2)` was 16,675 and is 4,667. +One is longer: `1/(x (1 - x^2)^(1/3) (3 + x^2))`, 342 characters there and 495 now. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a - c*x^4)/((a*d + a*k*x^2 + c*d*x^4)*sqrt(a + b*x^2 + c*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 1,949 characters | +| `"(g + h*x)*sqrt(a + b*x + c*x^2)/(a*d + b*d*x + c*d*x^2)^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | 4,376 characters | +| `"x^2/((1 - a*x)^16*(1 + a*x)^11)".ToEntity().Integrate("x")` | `integral(...)` | 1,380 characters | +| `"1/((a + b/x)^(5/2)*(c + d/x)^2)".ToEntity().Integrate("x")` | `integral(...)` | 4,667 characters | +| `"(c*x)^(-1 - 2*n - n*p)*(a + b*x^n)^p".ToEntity().Integrate("x")` | `integral(...)` | 255 characters | +| `"1/(x^2*(a + b*x^3)^5)".ToEntity().Integrate("x")` | `integral(...)` | 741 characters | +| `"(d + k*x)^4/(a + b*(d + k*x)^2 + c*(d + k*x)^4)".ToEntity().Integrate("x")` | no answer in twenty seconds | 452 characters | +| `"1/(x^5*(a + b*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 652 characters, where the unreleased master gave 67,828 | +| `"1/(x*(1 - x^2)^(1/3)*(3 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | 503 characters, where the unreleased master gave 350 | + ### A polynomial with symbols in it that is a binomial or an even quartic in a shifted variable is integrated in it **Answers where there were none.** `1/(c^2 x^3 + 3 b c x^2 + 3 b^2 x + 3 a b)` was left unevaluated, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 31dc3156e..6f991d533 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -7909,6 +7909,315 @@ Entity Terms(bool positive, Entity s) ? (constant, cubic) : null; } + /// + /// A rational function beside a half-odd power of an even quartic, R(x) P^(m/2) with + /// P = a + b x^2 + c x^4, as a rational function of t = x/sqrt(P) where it is + /// one: Rubi's (a - c x^4)/((a d + a e x^2 + c d x^4) sqrt(a + b x^2 + c x^4)) is + /// int dt/(d + (a e - b d) t^2), and sqrt(a + b x^2 + c x^4)/(a d - c d x^4) is + /// int dt/(d ((1 - b t^2)^2 - 4 a c t^4)). + /// + /// + /// + /// dt/dx = (a - c x^4)/P^(3/2), so the integrand is H dt for + /// H = R P^((m + 3)/2)/(a - c x^4), a rational function of x. With y = x^2, + /// t^2 = y/(a + b y + c y^2) is 1/(c w + b) for w = y + a/(c y); so + /// H is a rational function of t where it is even, H = h(x^2), and h is a + /// rational function of w -- the same at y and at a/(c y), the two roots of + /// Y^2 - w Y + a/c. That is read by writing the numerator and the denominator of + /// h modulo that quadratic, each as u_1 y + u_0 with the u_i polynomials + /// in w: the quotient is free of y exactly where n_1 d_0 = n_0 d_1, decided as a value, + /// and is then n_1/d_1. An elliptic integral does not pass: 1/sqrt(P) has + /// h(a/(c y)) = -h(y). + /// + /// + /// Each step is an identity for the principal root wherever P is not zero, of either sign: + /// P^(m/2) P^(3/2) = P^((m + 3)/2) for an odd m, t^2 = x^2/P, and the derivative + /// of P^(-1/2) on an interval where P is negative. These were declined, as no rule + /// writes the variable over the root; with numbers for coefficients too, + /// x^2/((1 + x^4) sqrt(x^4 - 1)) is -int t^2/(1 + 4 t^4) dt. Asked late and at + /// the question asked, so that every rule reading the integrand as written goes first, and + /// the answer is differentiated back at sampled points before it is returned. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveByTheVariableOverTheRootOfAnEvenQuartic(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + // The one radicand, its total power, and the rest, which is to be rational in x. + Entity? quartic = null; + var power = ERational.Zero; + Entity rest = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (factor is Powf(var @base, Number.Rational exponent) && exponent is not Number.Integer && @base.ContainsNode(x)) + { + if (!exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) || quartic is not null && quartic != @base) + return null; + quartic = @base; + power = (underneath ? power.Subtract(exponent.ERational) : power.Add(exponent.ERational)).ToLowestTerms(); + continue; + } + rest = underneath ? rest / factor : rest * factor; + } + if (quartic is null || !power.Denominator.Equals(EInteger.FromInt32(2)) || !power.Numerator.CanFitInInt32() + || rest.Nodes.Any(node => node.ContainsNode(x) && node is not (Sumf or Minusf or Mulf or Divf or Variable) + && !(node is Powf(_, Number.Integer))) + || !TreeAnalyzer.TryGetPolynomial(quartic, x, out var read) + || read.Keys.Any(degree => !(degree.IsZero || degree.Equals(EInteger.FromInt32(2)) || degree.Equals(EInteger.FromInt32(4)))) + || read.Values.Any(coefficient => coefficient.ContainsNode(x)) + || !read.TryGetValue(EInteger.Zero, out var a) || !read.TryGetValue(EInteger.FromInt32(4), out var c)) + return null; + var b = read.TryGetValue(EInteger.FromInt32(2), out var middle) ? middle : Number.Integer.Zero; + static bool Zero(Entity value) => Functions.PartialFractions.IsZeroAsAValue(value); + if (Zero(a) || Zero(c)) + return null; + // H, over one bar, its two sides even in x once a common power of x is taken off. + var raisedBy = Number.Integer.Create(power.Numerator.Add(EInteger.FromInt32(3)).ToInt32Checked() / 2); + var h = Functions.SingleQuotient.Combine(rest * MathS.Pow(quartic, raisedBy) / (a - c * MathS.Pow(x, 4))); + var (top, bottom) = Functions.SingleQuotient.Of(h); + if (!TreeAnalyzer.TryGetPolynomial(top.Expand(), x, out var above) || !TreeAnalyzer.TryGetPolynomial(bottom.Expand(), x, out var below) + || above.Count == 0 || below.Count == 0 + || above.Concat(below).Any(pair => pair.Key.Sign < 0 || !pair.Key.CanFitInInt32() || pair.Value.ContainsNode(x))) + return null; + var shift = EInteger.Min(above.Keys.Min()!, below.Keys.Min()!); + if (above.Keys.Concat(below.Keys).Any(degree => !(degree - shift).IsEven)) + return null; + // The two sides in y = x^2 modulo `y^2 - w y + a/c`, each as `u_1 y + u_0`, the u_i + // polynomials in w by their coefficients. + var ratio = Functions.PartialFractions.Bare((a / c).InnerSimplified); + var (n1, n0) = Reduced(above); + var (d1, d0) = Reduced(below); + if (!PolynomialDifference(PolynomialProduct(n1, d0), PolynomialProduct(n0, d1)).Values.All(Zero)) + return null; + var (overW, underW) = d1.Values.All(Zero) ? (n0, d0) : (n1, d1); + if (underW.Values.All(Zero)) + return null; + // At `w = (1 - b t^2)/(c t^2)`, each side times `(c t^2)` to its degree in w. + var t = Variable.CreateUnique(expr, "t_quartic"); + var (overT, overDegree) = InT(overW); + var (underT, underDegree) = InT(underW); + var inT = Functions.PartialFractions.Bare(Functions.SingleQuotient.Combine( + overT / underT * MathS.Pow(c * MathS.Sqr(t), Number.Integer.Create(underDegree - overDegree))).InnerSimplified); + if (inT.ContainsNode(x)) + return null; + if (Integration.ComputeAsAQuestionOfItsOwn(inT, t, integrateByParts) is not { } inTermsOfT) + return null; + var answer = inTermsOfT.Substitute(t, x / MathS.Sqrt(quartic)).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + + // A polynomial in w at `w = (1 - b t^2)/(c t^2)` times `(c t^2)^K`, K its degree. + (Entity, int) InT(Dictionary inW) + { + var degree = inW.Where(pair => !Zero(pair.Value)).Select(pair => pair.Key.ToInt32Checked()).DefaultIfEmpty(0).Max(); + Entity sum = Number.Integer.Zero; + foreach (var pair in inW) + { + if (Zero(pair.Value)) + continue; + var k = pair.Key.ToInt32Checked(); + sum += pair.Value * MathS.Pow(1 - b * MathS.Sqr(t), Number.Integer.Create(k)) * MathS.Pow(c * MathS.Sqr(t), Number.Integer.Create(degree - k)); + } + return (sum.Expand().InnerSimplified, degree); + } + + // The polynomial in y = x^2, read off the even powers past the shift, modulo + // `y^2 - w y + a/c`, by Horner's rule from the top: `(p y + q) y + k` is + // `(p w + q) y - p a/c + k` there. As `(p, q)` for `p y + q`. + (Dictionary, Dictionary) Reduced(Dictionary polynomial) + { + var p = new Dictionary(); + var q = new Dictionary(); + var byW = new Dictionary { [EInteger.One] = Number.Integer.One }; + var byRatio = new Dictionary { [EInteger.Zero] = (-ratio).InnerSimplified }; + var highest = (polynomial.Keys.Max()! - shift).ToInt32Checked() / 2; + for (var k = highest; k >= 0; k--) + { + var term = new Dictionary(); + if (polynomial.TryGetValue(shift + EInteger.FromInt32(2 * k), out var coefficient)) + term[EInteger.Zero] = coefficient; + (p, q) = (Bared(PolynomialSum(PolynomialProduct(p, byW), q)), Bared(PolynomialSum(term, PolynomialProduct(p, byRatio)))); + } + return (p, q); + } + + // Each coefficient without the condition the simplification of a quotient attaches. + static Dictionary Bared(Dictionary polynomial) + => polynomial.ToDictionary(pair => pair.Key, pair => Functions.PartialFractions.Bare(pair.Value)); + } + + /// + /// A factor (c x)^t of the integrand, t not whole, written as x^t with + /// (c x)^t x^(-t) in front: Rubi's (c x)^(-1 - 2 n - n p) (a + b x^n)^p and + /// (d x)^(-1 - 2 n (1 + p)) (a^2 + 2 a b x^n + b^2 x^(2 n))^p, whose powers of x alone are + /// answered. + /// + /// + /// The derivative of (c x)^t x^(-t) is t (c x)^t x^(-t) (c/(c x) - 1/x), zero + /// wherever it is defined, since (c x)^(t - 1) is (c x)^t/(c x) for the principal + /// power: it is a constant on each interval of the domain, for either sign of c and of x, and + /// goes in front of the integral of the rest. + /// writes c^t x^t for a positive c only, and reads (c x^n)^t for an n that is + /// not whole. Late and at the question asked, so that a rule reading (c x)^t as + /// written goes first; the answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveByWritingAPowerOfAMultipleOfXAsAPowerOfX(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + Entity inFront = Number.Integer.One; + Entity written = Number.Integer.One; + var changed = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + var rewritten = factor; + if (factor is Powf(var multiple, var t) && multiple != x && !t.ContainsNode(x) && t.Evaled is not Number.Integer + && TreeAnalyzer.TryGetPolyLinear(multiple, x, out var slope, out var intercept) + && !slope.ContainsNode(x) && !intercept.ContainsNode(x) + && Functions.PartialFractions.IsZeroAsAValue(intercept) && !Functions.PartialFractions.IsZeroAsAValue(slope)) + { + rewritten = MathS.Pow(x, t); + var ratio = MathS.Pow(multiple, t) * MathS.Pow(x, -t); + inFront = underneath ? inFront / ratio : inFront * ratio; + changed = true; + } + written = underneath ? written / rewritten : written * rewritten; + } + if (!changed || Integration.ComputeAsAQuestionOfItsOwn(written, x, integrateByParts) is not { } integral) + return null; + var answer = (inFront * integral).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// + /// Two square roots written as one up to a factor whose square is a constant: where the + /// quotient of the radicands is a constant times an even power of x, sqrt(A)/sqrt(B) + /// is K x^k, and where their product is free of a root of x that each holds, + /// sqrt(A) sqrt(B) is K sqrt(A B). Rubi's + /// (g + h x) sqrt(a + b x + c x^2)/(a d + b d x + c d x^2)^(3/2), + /// sqrt(b - a/x^2)/sqrt(a - b x^2) and + /// 1/(sqrt(1 + x) sqrt(-1 - sqrt(x)) sqrt(-1 + sqrt(x))). + /// + /// + /// + /// K is the quotient of the two sides, sqrt(A)/(x^k sqrt(B)) or + /// sqrt(A) sqrt(B)/sqrt(A B), so the rewriting is exact; and its square is + /// A/(x^(2k) B) or 1, free of x, so K is one of two values and constant on each + /// interval where it is continuous: it goes in front of the integral of the rest, as Rubi + /// writes it. For powers A^(p/2) B^(q/2) with odd p and q, the quotient is + /// B^((p + q)/2) (sqrt(A)/sqrt(B))^p, and the product needs p = q. Each was + /// declined: combines two roots of polynomials it + /// can sign, and reads no root inside a radicand, and + /// wants a logarithm to name + /// the radicand. Late and at the question asked, so that every rule reading the roots as + /// written goes first. The rewriting is checked at sampled points rather than the answer + /// differentiated back: the decimal evaluation of the derivative of K, zero, disagreed where + /// a root of x is imaginary, and the answer there was right. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveByWritingTwoRootsAsOneUpToAConstantFactor(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var factors = FactorsOfTheIntegrand(expr); + var roots = new List<(int At, Entity Base, int Odd)>(); + for (var i = 0; i < factors.Count; i++) + if (factors[i].Factor is Powf(var @base, Number.Rational exponent) && exponent is not Number.Integer && @base.ContainsNode(x) + && exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) && exponent.ERational.Numerator.CanFitInInt32()) + { + var odd = exponent.ERational.Numerator.ToInt32Checked(); + roots.Add((i, @base, factors[i].Underneath ? -odd : odd)); + } + if (roots.Count < 2 || roots.Count > 4) + return null; + for (var i = 0; i < roots.Count; i++) + for (var j = 0; j < roots.Count; j++) + { + if (i == j || roots[i].Base == roots[j].Base) + continue; + var (a, p) = (roots[i].Base, roots[i].Odd); + var (b, q) = (roots[j].Base, roots[j].Odd); + Entity written, factor; + if (QuotientAsAConstantTimesAnEvenPowerOfX(a, b) is { } k) + { + // A^(p/2) B^(q/2) = B^((p + q)/2) (K x^k)^p. + written = Raised(b, (p + q) / 2) * Raised(x, k * p); + factor = Raised(MathS.Sqrt(a) / (Raised(x, k) * MathS.Sqrt(b)), p); + } + else if (i < j && p == q && HoldsARootOfX(a) && HoldsARootOfX(b) + && Functions.PartialFractions.Bare((a * b).Expand().InnerSimplified) is var product && !HoldsARootOfX(product) + && product.ContainsNode(x)) + { + // A^(p/2) B^(p/2) = (K sqrt(A B))^p. + written = MathS.Pow(product, Number.Rational.Create(p, 2)); + factor = Raised(MathS.Sqrt(a) * MathS.Sqrt(b) / MathS.Sqrt(product), p); + } + else + continue; + Entity rest = Number.Integer.One; + for (var m = 0; m < factors.Count; m++) + if (m != roots[i].At && m != roots[j].At) + rest = factors[m].Underneath ? rest / factors[m].Factor : rest * factors[m].Factor; + var rewritten = Functions.PartialFractions.Bare((rest * written).InnerSimplified); + // The rewriting checked, rather than the answer differentiated back: the decimal + // evaluation of the derivative of K, zero, disagreed at points where a root of + // x is imaginary, and the answer there is right. + bool rewrittenHolds; + using (MathS.Settings.DowncastingEnabled.Set(false)) + rewrittenHolds = Functions.PartialFractions.HoldsAtSampledPoints(factor * rewritten, expr, x); + if (!rewrittenHolds || Integration.ComputeAsAQuestionOfItsOwn(rewritten, x, integrateByParts) is not { } integral) + return null; + var answer = (factor * integral).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + return null; + + static Entity Raised(Entity e, int power) => power switch + { + 0 => Number.Integer.One, + 1 => e, + _ => MathS.Pow(e, Number.Integer.Create(power)), + }; + + bool HoldsARootOfX(Entity e) => e.Nodes.Any(node => node is Powf(var @base, Number.Rational r) && r is not Number.Integer && @base.ContainsNode(x)); + + // k where A/B is a constant times x^(2k): each side a polynomial over a power of x, the + // polynomials a constant multiple of one another. + int? QuotientAsAConstantTimesAnEvenPowerOfX(Entity first, Entity second) + { + if (HoldsARootOfX(first) || HoldsARootOfX(second) + || OverAPowerOfX(first) is not var (firstAbove, firstPower) || OverAPowerOfX(second) is not var (secondAbove, secondPower) + || (secondPower - firstPower) % 2 != 0 + || TryReadAsAConstantMultiple(firstAbove, secondAbove, x, allowOne: true) is null) + return null; + return (secondPower - firstPower) / 2; + } + + // A polynomial over a whole power of x, as the polynomial and the power. + (Entity, int)? OverAPowerOfX(Entity e) + { + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(e)); + var power = 0; + foreach (var part in Mulf.LinearChildren(below)) + { + if (!part.ContainsNode(x)) + above /= part; + else if (part == x) + power++; + else if (part is Powf(var @base, Number.Integer whole) && @base == x && whole.EInteger.CanFitInInt32()) + power += whole.EInteger.ToInt32Checked(); + else + return null; + } + return TreeAnalyzer.TryGetPolynomial(above, x, out var read) && read.Keys.All(degree => degree.Sign >= 0) && read.Values.All(value => !value.ContainsNode(x)) + ? (above, power) : null; + } + } + /// /// A root of a quadratic binomial beside another, 1/((A + B x^2)^(1/3) (C + D x^2)) /// with B C + 3 A D = 0 or B C - 9 A D = 0, and @@ -18345,6 +18654,82 @@ private static bool IsAFunctionToRaise(Entity factor, Entity.Variable x) return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; } + /// + /// A power of x beside a whole power below -1 of a binomial with a symbol in it, + /// x^m (a + b x^n)^p for whole m, n and p <= -2, by Rubi's reduction that + /// raises p by one, down to int x^m/(a + b x^n), which is asked on its own: + /// + /// J(p) = -x^(m + 1) (a + b x^n)^(p + 1)/(a n (p + 1)) + (m + 1 + n (p + 1))/(a n (p + 1)) J(p + 1) + /// + /// + /// + /// The derivative of x^(m + 1) B^(p + 1) is x^m B^p ((m + 1 + n (p + 1)) B - a n (p + 1)), + /// which is the step, exact for whole powers. Rubi's 1/(x^2 (a + b x^3)^5), which the + /// root of its perfect square leaves from 1/(x^2 (a^2 + 2 a b x^3 + b^2 x^6)^(5/2)), was + /// declined: the partial fractions over the cube roots of a symbolic binomial gave out past the + /// first power, which is answered. Late, so that every rule answering the power as written goes + /// first; the answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveABinomialToAWholePowerBelowMinusOneByTheRecurrence(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + var m = 0; + Entity? binomial = null; + var p = 0; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, whole) = factor is Powf(var raised, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (raised, w.ToInt32Checked()) : (factor, 1); + if (underneath) + whole = -whole; + if (@base == x) + { + m += whole; + continue; + } + if (binomial is not null && binomial != @base || @base is not (Sumf or Minusf)) + return null; + binomial = @base; + p += whole; + } + if (binomial is null || p > -2 || p < -MostBinomialPowersReduced + || !TreeAnalyzer.TryGetPolynomial(binomial, x, out var read) || read.Count != 2 + || !read.TryGetValue(EInteger.Zero, out var a) || read.Keys.Max() is not { } degree || !degree.CanFitInInt32() + || degree.Sign <= 0 || read.Values.Any(value => value.ContainsNode(x)) || Functions.PartialFractions.IsZeroAsAValue(a)) + return null; + var n = degree.ToInt32Checked(); + Entity sum = Number.Integer.Zero; + Entity coefficient = Number.Integer.One; + var raisedX = m + 1 == 0 ? Number.Integer.One : m + 1 == 1 ? (Entity)x : MathS.Pow(x, Number.Integer.Create(m + 1)); + for (var q = p; q <= -2; q++) + { + var below = a * n * (q + 1); + sum += LowestOverTheSymbols(-coefficient / below) * raisedX * MathS.Pow(binomial, Number.Integer.Create(q + 1)); + coefficient = LowestOverTheSymbols(coefficient * (m + 1 + n * (q + 1)) / below); + } + if (!Functions.PartialFractions.IsZeroAsAValue(coefficient)) + { + var last = (m == 0 ? Number.Integer.One : MathS.Pow(x, Number.Integer.Create(m))) / binomial; + if (Integration.ComputeAsAQuestionOfItsOwn(last, x, integrateByParts) is not { } integral) + return null; + sum += coefficient * integral; + } + var answer = (constant * sum).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// The most powers raises. + private const int MostBinomialPowersReduced = 16; + /// /// A binomial differential x^m (a + b x^n)^p with a symbol in its exponents, in /// Chebyshev's third case at a whole (m + 1)/n + p + 1 = -k of at most 0: Rubi's @@ -18575,6 +18960,127 @@ static Entity Simplified(Entity exponent) return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// A polynomial over high whole powers of two linears, P(x)/(L1^m L2^n) with + /// m + n = k past and P of degree + /// at most k - 2, under t = L1/L2: a polynomial in t over a power of it. + /// Rubi's x^2/((1 - a x)^16 (1 + a x)^11). + /// + /// + /// + /// With D = a1 b2 - a2 b1, x = (a1 - a2 t)/(b2 t - b1), L2 = D/(b2 t - b1) + /// and dx = -D dt/(b2 t - b1)^2, the integrand is -D^(1 - k) t^(-m) S(t) dt for + /// S(t) = sum_j p_j (a1 - a2 t)^j (b2 t - b1)^(k - 2 - j), as + /// has it for exponents + /// that are not whole. Here they are, so each t^(i - m) is L1^(i - m) L2^(m - i) + /// exactly, and the one at i = m - 1 integrates to ln(L1) - ln(L2). + /// + /// + /// The partial fractions take the two linears apart with a symbol in them over k terms, and + /// ran past the budget at k = 27; below the bound they answer in their own form, and + /// this leaves them to it, as it leaves numbers to Hermite's reduction. The answer is + /// differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAPolynomialOverHighWholePowersOfTwoLinears(Entity expr, Entity.Variable x) + { + // With numbers only, Hermite's reduction above has it. + if (!expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + Entity above = Number.Integer.One; + var linears = new List<(Entity Linear, int Power)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, whole) = factor is Powf(var raised, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (raised, w.ToInt32Checked()) : (factor, 1); + if (underneath) + whole = -whole; + if (whole < 0 && TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out var intercept) + && !slope.ContainsNode(x) && !intercept.ContainsNode(x)) + { + var index = linears.FindIndex(pair => pair.Linear == @base); + if (index >= 0) + linears[index] = (@base, linears[index].Power - whole); + else if (linears.Count == 2) + return null; + else + linears.Add((@base, -whole)); + continue; + } + if (whole <= 0) + return null; + above *= factor; + } + if (linears.Count != 2) + return null; + var ((first, m), (second, n)) = (linears[0], linears[1]); + var k = m + n; + if (k < FewestLinearPowersWrittenInTheirQuotient || k > MostLinearPowersWrittenInTheirQuotient + || !TreeAnalyzer.TryGetPolynomial(above.Expand(), x, out var coefficients) + || coefficients.Any(pair => pair.Key.Sign < 0 || pair.Key.CompareTo(EInteger.FromInt32(k - 2)) > 0 || pair.Value.ContainsNode(x)) + || !TreeAnalyzer.TryGetPolyLinear(first, x, out var b1, out var a1) || !TreeAnalyzer.TryGetPolyLinear(second, x, out var b2, out var a2)) + return null; + // A power of x beside one linear is a binomial differential, which the rules for one + // answer at once. + var d = (a1 * b2 - a2 * b1).InnerSimplified; + if (Functions.PartialFractions.IsZeroAsAValue(d) || Functions.PartialFractions.IsZeroAsAValue(b1) || Functions.PartialFractions.IsZeroAsAValue(b2) + || Functions.PartialFractions.IsZeroAsAValue(a1) || Functions.PartialFractions.IsZeroAsAValue(a2)) + return null; + // S by its coefficients, from the binomial theorem: `(a1 - a2 t)^j (b2 t - b1)^(k - 2 - j)` + // has at t^i the sum over r + s = i of C(j, r) C(k - 2 - j, s) a1^(j - r) (-a2)^r b2^s (-b1)^(k - 2 - j - s). + // Expanded as a tree, a nineteenth power of `d t - b` came back unread. + var terms = new Dictionary>(); + foreach (var pair in coefficients) + { + var j = pair.Key.ToInt32Checked(); + var rest = k - 2 - j; + for (var r = 0; r <= j; r++) + for (var q = 0; q <= rest; q++) + { + var term = Number.Integer.Create(Binomial(j, r) * Binomial(rest, q)) * pair.Value + * Raised(a1, j - r) * Raised(-a2, r) * Raised(b2, q) * Raised(-b1, rest - q); + if (!terms.TryGetValue(r + q, out var list)) + terms[r + q] = list = new List(); + list.Add(term); + } + } + var ofT = terms.ToDictionary(pair => EInteger.FromInt32(pair.Key), pair => pair.Value.Skip(1).Aggregate(pair.Value[0], (sum, term) => sum + term)); + Entity sum = Number.Integer.Zero; + foreach (var pair in ofT.OrderBy(pair => pair.Key)) + { + var i = pair.Key.ToInt32Checked(); + var s = LowestOverTheSymbols(pair.Value); + if (s.Evaled is Number.Complex { IsZero: true }) + continue; + sum += i == m - 1 + ? s * (MathS.Ln(first) - MathS.Ln(second)) + : LowestOverTheSymbols(s / (i + 1 - m)) * MathS.Pow(first, Number.Integer.Create(i + 1 - m)) * MathS.Pow(second, Number.Integer.Create(m - 1 - i)); + } + var answer = (LowestOverTheSymbols(-constant * MathS.Pow(d, 1 - k)) * sum).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + + static Entity Raised(Entity e, int power) => power == 0 ? Number.Integer.One : power == 1 ? e : MathS.Pow(e, Number.Integer.Create(power)); + + static EInteger Binomial(int n, int r) + { + var value = EInteger.One; + for (var i = 1; i <= r; i++) + value = value * EInteger.FromInt32(n - r + i) / EInteger.FromInt32(i); + return value; + } + } + + /// The fewest and the most powers takes. + private const int FewestLinearPowersWrittenInTheirQuotient = 20, MostLinearPowersWrittenInTheirQuotient = 60; + /// /// The most negative sum of the two exponents /// takes: the polynomial @@ -23809,6 +24315,15 @@ Entity TakeOut(Entity node) // a second and ten levels of this ran past the budget. if (Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic) || linear is null && VanishesIdentically(a1)) return null; + // Nor beside a power of x where the integrand is odd and the power of the quadratic high: + // it is then x times a function of x^2, the substitution's under u = x^2, and peeling + // fifteen powers of x off `(b + 2 c x^2)/(x^15 (b + c x^2)^8)` here ran past the budget, + // where in u it is two linears. An even one, and a low power, are left here, where they + // were answered. + if (linear is not null && power >= 6 && VanishesIdentically(a1) && TreeAnalyzer.TryGetPolyLinear(linear, x, out _, out var interceptOfTheLinear) + && VanishesIdentically(interceptOfTheLinear) + && top.Keys.All(degree => degree.IsEven != (linearPower % 2 == 0))) + return null; // A quadratic that factors over linears is the partial fractions' to take apart, and // they do it cheaply: `(d + e x)^10/((d + e x) (a e + c d x))^4`, written out, went from // two seconds past the budget here. @@ -26283,6 +26798,65 @@ private static (Entity Constant, Entity Base, Entity ExponentInU, Variable U, En return (constant, exponential.Base, exponentInU, u, linear, slope, power, beside ? (above.InnerSimplified, below.InnerSimplified, slopePower) : null); } + /// + /// An algebraic function of 1/x alone with a root in it, under x = 1/u: + /// f(1/u) (-1/u^2) du. Rubi's 1/((a + b/x)^(5/2) (c + d/x)^2) is + /// -int du/(u^2 (a + b u)^(5/2) (c + d u)^2), two linears in u beside a root of one. + /// + /// + /// x = 1/u is one to one on each side of zero, and d/dx F(1/x) = -F'(1/x)/x^2, + /// so the answer in u written back holds wherever the integrand is defined. Only where x + /// stands nowhere but below a bar, as b/x^k in two sums or more, each below the bar + /// itself, and some power of such a sum is not whole: a rational function of x is the + /// partial fractions', one sum beside a power of x is read by the binomial rules as it + /// stands, and with a sum above the bar the routes this pre-empts answer, where in u + /// sqrt(a + b u)/(u^2 (c + d u)^3) ran past the budget. wrote the + /// root over its linears, and the partial fractions it reached ran past the budget. At the + /// question asked, before that rule. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveAnAlgebraicFunctionOfTheReciprocalOfX(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + // Algebraic in x, with a power that is not whole of something holding it. + var aRoot = false; + foreach (var node in expr.Nodes) + { + if (node.ContainsNode(x) && node is not (Variable or Sumf or Minusf or Mulf or Divf or Powf)) + return null; + if (node is Powf(var @base, var exponent) && @base.ContainsNode(x) && exponent.ContainsNode(x)) + return null; + } + // A root to a power past the first below the bar: over a square root alone the routes in + // x answered, two of them faster than in u. + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + if (underneath && factor is Powf(var @base, Number.Rational { ERational: var power } exponent) && exponent is not Number.Integer + && @base.ContainsNode(x) && power.CompareTo(ERational.Create(3, 2)) >= 0) + aRoot = true; + // Two sums of 1/x or more, every factor with x below the bar, the shape the linears + // written once took apart: one beside powers of x is answered as it stands, and with a sum + // above the bar the routes it pre-empts answer, where in u two ran past the budget. + if (!aRoot || expr.Nodes.Where(node => node is Sumf or Minusf && node.ContainsNode(x)).Distinct().Count() < 2 + || FactorsOfTheIntegrand(expr).Any(pair => !pair.Underneath && pair.Factor.ContainsNode(x))) + return null; + var u = Variable.CreateUnique(expr, "u_reciprocal"); + var inU = expr.Replace(node => node switch + { + Divf(var above, Powf(var below, Number.Integer k)) when below == x && k.EInteger.Sign > 0 && !above.ContainsNode(x) => above * MathS.Pow(u, k), + Divf(var above, var below) when below == x && !above.ContainsNode(x) => above * u, + Powf(var below, Number.Integer { EInteger.Sign: < 0 } k) when below == x => MathS.Pow(u, (-k).InnerSimplified), + _ => node + }); + if (inU.ContainsNode(x)) + return null; + var integrand = Functions.PartialFractions.Bare((-inU / MathS.Sqr(u)).InnerSimplified); + if (Integration.ComputeAsAQuestionOfItsOwn(integrand, u, integrateByParts) is not { } inTermsOfU) + return null; + var answer = inTermsOfU.Substitute(u, 1 / x).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + /// /// Sines, cosines and the rest of a function of the reciprocal of a linear, beside a whole /// power of the linear: L^m f(a + b/L^k) with L = c + d x. Under u = 1/L, @@ -30770,9 +31344,30 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons /// in it: with numbers the substitution's rational function is split over the rationals. /// internal static Entity? SolveARationalFunctionBesideARootOfALinearSplitFirst(Entity expr, Entity.Variable x, bool integrateByParts) + => expr.Vars.Any(symbol => symbol != x) ? SplitFirstBesideARootOfALinear(expr, x, integrateByParts, withNumbers: false) : null; + + /// + /// The same with numbers for coefficients, beside a cube root or a higher one, where a linear + /// below the bar has the radicand at its root no rational power of that order: Rubi's + /// 1/(u^2 (1 - u)^(1/3) (3 + u)), which x^3 (1 - x^2)^(1/3) (3 + x^2) is under + /// u = x^2. In v = (1 - u)^(1/3) the linear 3 + u is 4 - v^3, and + /// the partial fractions over the cube root of 4 beside the rest ran past the budget, where + /// split over the rationals first each term beside the root is answered in a tenth of a + /// second. Before the cancellation with functions as indeterminates, which spent five seconds + /// declining it, and the substitution, which this pre-empts only there. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveARationalFunctionBesideAnIrrationalOddRootOfALinearSplitFirst(Entity expr, Entity.Variable x, bool integrateByParts) { - if (!expr.Vars.Any(symbol => symbol != x)) + if (expr.Vars.Any(symbol => symbol != x)) return null; + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + return SplitFirstBesideARootOfALinear(numerator / denominator, x, integrateByParts, withNumbers: true); + } + + private static Entity? SplitFirstBesideARootOfALinear(Entity expr, Entity.Variable x, bool integrateByParts, bool withNumbers) + { + var linearsBelow = new List(); Entity? radicand = null; var exponent = ERational.Zero; Entity constant = Number.Integer.One; @@ -30803,12 +31398,21 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons { below = below == Number.Integer.One ? written : below * written; belowFactors++; + if (read.Keys.Max()!.Equals(EInteger.One)) + linearsBelow.Add(@base); } else above = above == Number.Integer.One ? written : above * written; } if (radicand is null) return null; + // With numbers, only an odd root past the square of a radicand that is no rational power + // of its order at the root of some linear below: there the substitution's rational function + // is over an irrational root, and elsewhere it is answered whole. + if (withNumbers && (exponent.Denominator.CompareTo(EInteger.FromInt32(3)) < 0 || exponent.Denominator.IsEven + || !exponent.Denominator.CanFitInInt32() + || !linearsBelow.Any(linear => RadicandIsNoRationalPowerAtTheRootOf(linear, exponent.Denominator.ToInt32Checked())))) + return null; // The whole part of the power goes into the rational function, and what is left of it // is a root below the bar, `L^q = L^k / L^(k - q)` with `k` the least whole number past // `q`: `1/(x (1 + x^2) (a + b x)^(3/2))` is split with `a + b x` among the factors and @@ -30844,6 +31448,25 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons sum += integral; } return constant * sum; + + // Whether the radicand, at the root of the linear, is a rational number with no rational + // root of the order: 4 at the root of `3 + u` beside `(1 - u)^(1/3)`. + bool RadicandIsNoRationalPowerAtTheRootOf(Entity linear, int order) + { + if (!TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var intercept) + || radicand!.Substitute(x, -intercept / slope).Evaled is not Number.Rational { ERational: var value } || value.IsZero) + return false; + return !IsAPowerOf(value.Numerator.Abs(), order) || !IsAPowerOf(value.Denominator, order); + } + + static bool IsAPowerOf(EInteger value, int order) + { + var root = EInteger.FromInt64((long)System.Math.Round(System.Math.Pow(EDecimal.FromEInteger(value).ToDouble(), 1.0 / order))); + for (var step = -1; step <= 1; step++) + if (root + EInteger.FromInt32(step) is var candidate && candidate.Sign >= 0 && candidate.Pow(order).Equals(value)) + return true; + return false; + } } /// @@ -30883,16 +31506,31 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons /// x^2 Shi(a + b x) sinh(a + b x) leaves a run of them beside powers of x, and /// written out each term is a search of its own. /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// And beside one power of a sum of powers of x with a symbol in its exponent and in the + /// powers, kept beside each term of eight at most: Rubi's + /// (a + b x^n)^3 (c + d x^n)^(-4 - 1/n) is four binomial differentials in Chebyshev's + /// third case, each answered, and as a product it was declined. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// /// internal static Entity? SolveByExpandingAProductOfSums(Entity expr, Entity.Variable x, bool integrateByParts) { if (expr is not Mulf || !Integration.AnsweringTheQuestionAsked) return null; var sums = 0; + // And one power of a sum of powers of x with a symbol in its exponent, kept beside each + // term: `(a + b x^n)^3 (c + d x^n)^(-4 - 1/n)` is four terms `x^(k n) (c + d x^n)^(-4 - 1/n)`, + // each a binomial differential in Chebyshev's third case, where the product is not. + Entity? kept = null; + Entity expanded = Number.Integer.One; foreach (var factor in Mulf.LinearChildren(expr)) { if (!factor.ContainsNode(x)) + { + expanded *= factor; continue; + } switch (factor) { case Sumf or Minusf when IsASumOfPowersOfX(factor, x): @@ -30902,15 +31540,27 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons case Entity.Variable when factor == x: case Powf(var @base, var exponent) when @base == x && !exponent.ContainsNode(x): break; + case Powf((Sumf or Minusf) and var sum, var exponent) when kept is null && exponent.Vars.Any() && !exponent.ContainsNode(x) + && sum.Nodes.Any(node => node is Powf(var @base, var power) && @base == x && power.Vars.Any()) && IsASumOfPowersOfX(sum, x): + kept = factor; + continue; default: return null; } + expanded *= factor; } if (sums == 0) return null; - var written = expr.Expand(); + var written = (kept is null ? expr : expanded).Expand(); if (written is not Sumf and not Minusf) return null; + if (kept is not null) + { + var terms = Sumf.LinearChildren(written).ToList(); + if (terms.Count > 8) + return null; + written = terms.Skip(1).Aggregate(terms[0] * kept, (sum, term) => sum + term * kept); + } return Integration.ComputeIndefiniteIntegral(written, x, integrateByParts); // Each term a constant times a power of x, the power's exponent free of x. @@ -33035,6 +33685,12 @@ private static (Entity Numerator, Entity Denominator)? WithASharedPowerCancelled // or a multiple of one root, and neither is this rule's. if (product.Evaled is Number.Complex { IsZero: true } || product.Simplify().Evaled is Number.Complex { IsZero: true }) return null; + // A difference that is one power of x, written as one: `a + b x - (a + c x)` is + // `(b - c) x`, and left as the inner simplification writes it, the cube below + // `x^4 (sqrt(a + b x) - sqrt(a + c x))^3` ran a substitution past the budget. + if (TreeAnalyzer.TryGetPolynomial(product.Expand(), x, out var difference) && difference.Count == 1 + && difference.Single() is { Key: var degree, Value: var coefficient } && degree.Sign > 0 && !coefficient.ContainsNode(x)) + product = Functions.PartialFractions.Bare(coefficient.InnerSimplified) * (degree.Equals(EInteger.One) ? x : MathS.Pow(x, Number.Integer.Create(degree))); // An even power of the conjugate has the root of the product of the radicands // in it, `2 sqrt(A B)` in the square, which over a constant difference leaves a // root of a quadratic where the power as written was answered: the squares over @@ -36331,6 +36987,15 @@ bool HoldsItsValue(Entity factor) if (Mulf.LinearChildren(denominator).Count(factor => factor.ContainsNode(x)) != 1 || Mulf.LinearChildren(denominator).First(factor => factor.ContainsNode(x)) is not (Sumf or Minusf)) return null; + // Nor one written in powers of a single shifted linear, `(d + e x)^4/(a + b (d + e x)^2 + c (d + e x)^4)`: + // in u = d + e x it is the even quartic the rules for a trinomial answer, and the shifted + // linear rule writes it so; split here first, the shift it asked of the split form ran on + // past the budget. + if (expr.Nodes.OfType().Where(power => power.Base.ContainsNode(x)).Select(power => power.Base).Distinct().ToList() is { Count: 1 } bases + && TreeAnalyzer.TryGetPolyLinear(bases[0], x, out _, out var offsetOfTheLinear) && !offsetOfTheLinear.ContainsNode(x) + && !VanishesIdentically(offsetOfTheLinear) + && !expr.Replace(node => node == bases[0] ? Variable.CreateUnique(expr, "u_shifted") : node).ContainsNode(x)) + return null; if (!TreeAnalyzer.TryGetPolynomial(denominator, x, out var read) || read.Keys.Any(key => key.Sign < 0 || key.CompareTo(EInteger.FromInt32(4)) > 0) || !read.ContainsKey(EInteger.FromInt32(4)) || read.Values.Any(value => value.ContainsNode(x)) || !TreeAnalyzer.TryGetPolynomial(numerator, x, out var readAbove) || readAbove.Values.Any(value => value.ContainsNode(x))) diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index b771c2129..5d00475e6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -726,6 +726,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A whole power of a linear in the inverse sine, cosine or hyperbolic cosine of // `c + d x^2` with `c^2 = 1`, by its recurrence, which closes: by parts it went round. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(expr, x)) is { }) return answer; + // An algebraic function of 1/x alone with a root in it, under x = 1/u, before the linears + // are written once: `1/((a + b/x)^(5/2) (c + d/x)^2)` ran past the budget that way. + if ((answer = IndefiniteIntegralSolver.SolveAnAlgebraicFunctionOfTheReciprocalOfX(expr, x, integrateByParts)) is { }) return answer; // A product of powers of quotients of linears with one linear in two of them, written // over each linear once with the constant that takes in front. First, since it answers // the question asked only, and the polynomial term below writes `1/u` as `u^(-1)` and @@ -757,6 +760,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A rational function with numbers for coefficients over a denominator past the twelfth // degree, by Hermite's reduction before anything searches: nothing below answered one. if ((answer = IndefiniteIntegralSolver.SolveAHighDegreeRationalFunctionByHermite(expr, x)) is { }) return answer; + // And a polynomial over high whole powers of two linears with a symbol in them, under + // t = L1/L2, where the partial fractions ran past the budget: `x^2/((1 - a x)^16 (1 + a x)^11)`. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialOverHighWholePowersOfTwoLinears(expr, x)) is { }) return answer; // A whole power of the square of a polynomial written out, written over its root: the rules // for a power of a binomial read `(a + b x^2)^4` where `(a^2 + 2 a b x^2 + b^2 x^4)^2` ran // past the budget. @@ -1022,6 +1028,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // split, since Hearn's `e^(x^2)/x + 2x e^(x^2) ln(x) + ...` is a sum whose terms are // not elementary apart; at the top only, and only for a tower no closed rule reads. if ((answer = IndefiniteIntegralSolver.SolveByRischNormanAnsatz(expr, x)) is { }) return answer; + // A rational function with numbers beside an odd root of a linear past the square, split over + // the rationals first where a linear below meets the root at an irrational value: + // `1/(u^2 (1 - u)^(1/3) (3 + u))`, whose substitution is over the cube root of 4, and which the + // cancellation below spent five seconds declining. + if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionBesideAnIrrationalOddRootOfALinearSplitFirst(expr, x, integrateByParts)) is { }) return answer; // A quotient of polynomials in x and the functions of x in it, cancelled with those // functions as indeterminates and its denominator factored the same way: // `(x^2 + 2x sin(x) + sin(x)^2)/(x + sin(x))^2` is `1`, and nothing else reads it. @@ -1537,6 +1548,20 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A half-odd power of a constant over a polynomial of either sign, written apart with // its sign in front. After the substitutions, which answer some of these without it. if ((answer = IndefiniteIntegralSolver.SolveByWritingAHalfOddPowerOfAReciprocalWithItsSign(expr, x, integrateByParts)) is { }) return answer; + // A power of x beside a whole power below -1 of a binomial with a symbol in it, by the + // reduction that raises the power to -1: `1/(x^2 (a + b x^3)^5)`, which the partial + // fractions declined past the first power. + if ((answer = IndefiniteIntegralSolver.SolveABinomialToAWholePowerBelowMinusOneByTheRecurrence(expr, x, integrateByParts)) is { }) return answer; + // A rational function beside a half-odd power of an even quartic P, as a rational + // function of x/sqrt(P) where it is one: Rubi's pseudo-elliptic + // `(a - c x^4)/((a d + a e x^2 + c d x^4) sqrt(a + b x^2 + c x^4))`. Late, since it fires + // on nothing a rule reading the root as written answers. + if ((answer = IndefiniteIntegralSolver.SolveByTheVariableOverTheRootOfAnEvenQuartic(expr, x, integrateByParts)) is { }) return answer; + // A factor `(c x)^t` as `x^t`, `(c x)^t x^(-t)` in front, constant on each interval. + if ((answer = IndefiniteIntegralSolver.SolveByWritingAPowerOfAMultipleOfXAsAPowerOfX(expr, x, integrateByParts)) is { }) return answer; + // Two square roots as one, up to a factor whose square is a constant, in front: + // `sqrt(b - a/x^2)/sqrt(a - b x^2)` is that factor over x. + if ((answer = IndefiniteIntegralSolver.SolveByWritingTwoRootsAsOneUpToAConstantFactor(expr, x, integrateByParts)) is { }) return answer; // The sign of a real-valued factor is constant between its zeros, and goes in // front of the antiderivative of the rest. After every rule that reads the sign // where it stands: `cos(x) sgn(sin(x))` is `|sin(x)|` by the substitution, and diff --git a/Sources/Tests/UnitTests/Calculus/ABinomialToAWholePowerBelowMinusOneIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ABinomialToAWholePowerBelowMinusOneIntegralTest.cs new file mode 100644 index 000000000..aa88687c8 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ABinomialToAWholePowerBelowMinusOneIntegralTest.cs @@ -0,0 +1,51 @@ +// +// 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 beside a whole power below -1 of a binomial with a symbol in it, by the reduction + /// that raises the power to -1, and the root of a perfect square that leaves one: Rubi's 1.2.3.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ABinomialToAWholePowerBelowMinusOneIntegralTest + { + [Theory] + [InlineData("1/(x^2*(a + b*x^3)^5)")] + [InlineData("1/(x^3*(a + b*x^3)^5)")] + [InlineData("1/(x^2*(a^2 + 2*a*b*x^3 + b^2*x^6)^(5/2))")] + [InlineData("1/(x^3*(a^2 + 2*a*b*x^3 + b^2*x^6)^(5/2))")] + public void RaisedToTheFirstPower(string integrand) + { + 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).Substitute("c", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/APowerOfAMultipleOfXAsAPowerOfXIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAMultipleOfXAsAPowerOfXIntegralTest.cs new file mode 100644 index 000000000..7fc06d3e8 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfAMultipleOfXAsAPowerOfXIntegralTest.cs @@ -0,0 +1,53 @@ +// +// 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 factor (c x)^t written as x^t, with (c x)^t x^(-t) in front, whose + /// derivative is zero: Rubi's 1.1.3.2 and 1.2.3.2 with a symbol for the multiple of x. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfAMultipleOfXAsAPowerOfXIntegralTest + { + [Theory] + [InlineData("(c*x)^(-1 - 2*n - n*p)*(a + b*x^n)^p")] + [InlineData("(c*x)^(-1 - 3*n - n*p)*(a + b*x^n)^p")] + [InlineData("(c*x)^(-1 - 4*n - n*p)*(a + b*x^n)^p")] + [InlineData("(d*x)^(-1 - 2*n*(1 + p))*(a^2 + 2*a*b*x^n + b^2*x^(2*n))^p")] + public void WrittenAsAPowerOfX(string integrand) + { + 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}"); + // n past the wrap, so that x^n at a negative x is not confined to a half-plane. + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", -0.4).Substitute("d", 0.3) + .Substitute("n", 2.3).Substitute("p", 0.6); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/AProductOfSumsBesideASymbolicPowerIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AProductOfSumsBesideASymbolicPowerIntegralTest.cs new file mode 100644 index 000000000..6e30ebdb6 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AProductOfSumsBesideASymbolicPowerIntegralTest.cs @@ -0,0 +1,52 @@ +// +// 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 whole power of a binomial in x^n beside a power of another with a symbol in its + /// exponent, written out with the power kept beside each term: each term is a binomial + /// differential in Chebyshev's third case. Rubi's 1.1.3.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AProductOfSumsBesideASymbolicPowerIntegralTest + { + [Theory] + [InlineData("(a + b*x^n)^3*(c + d*x^n)^(-4 - 1/n)")] + [InlineData("(a + b*x^n)^2*(c + d*x^n)^(-4 - 1/n)")] + public void WrittenOutBesideThePower(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 10000, $"{text.Length} characters of answer for {integrand}"); + // n past the wrap, so that x^n at a negative x is not confined to a half-plane. + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.3) + .Substitute("n", 2.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/AQuarticInOneShiftedLinearIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AQuarticInOneShiftedLinearIntegralTest.cs new file mode 100644 index 000000000..de298ee17 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AQuarticInOneShiftedLinearIntegralTest.cs @@ -0,0 +1,49 @@ +// +// 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 a linear over a trinomial in its square, written in the linear rather than split + /// over the quartic's quadratics first: Rubi's 1.2.3.2, which ran past the budget split. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AQuarticInOneShiftedLinearIntegralTest + { + [Theory] + [InlineData("(d + k*x)^4/(a + b*(d + k*x)^2 + c*(d + k*x)^4)")] + [InlineData("(d + k*x)^2/(a + b*(d + k*x)^2 + c*(d + k*x)^4)")] + public void WrittenInTheShiftedLinear(string integrand) + { + 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).Substitute("c", 0.4).Substitute("d", 0.9).Substitute("k", 1.1); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/ARationalFunctionBesideAnIrrationalOddRootIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAnIrrationalOddRootIntegralTest.cs new file mode 100644 index 000000000..afd0d1a68 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAnIrrationalOddRootIntegralTest.cs @@ -0,0 +1,53 @@ +// +// 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 rational function with numbers beside a cube root of a linear, split over the rationals + /// first where a linear below meets the root at an irrational value: Rubi's 1.1.2.4 and 1.1.3.4 + /// under the substitution of a power of x, which ran past the budget over the cube root of 4. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ARationalFunctionBesideAnIrrationalOddRootIntegralTest + { + [Theory] + [InlineData("1/(x^3*(1 - x^2)^(1/3)*(3 + x^2))")] + [InlineData("1/(x^5*(1 - x^2)^(1/3)*(3 + x^2))")] + [InlineData("1/(x^3*(1 - x^2)^(1/3)*(3 + x^2)^2)")] + [InlineData("1/(x^5*(1 - x^2)^(1/3)*(3 + x^2)^2)")] + [InlineData("1/(x^4*(1 - x^3)^(4/3)*(1 + x^3))")] + public void SplitOverTheRationalsFirst(string integrand) + { + 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; + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -0.9, -0.4, 0.3, 0.6, 0.8, 1.7, 2.9 }) + { + 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/ASumOfRootsWithAMonomialDifferenceIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ASumOfRootsWithAMonomialDifferenceIntegralTest.cs new file mode 100644 index 000000000..ebf0a3a0a --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumOfRootsWithAMonomialDifferenceIntegralTest.cs @@ -0,0 +1,49 @@ +// +// 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 over a power of a sum of two roots whose radicands differ by a multiple of x, + /// rationalised with the difference written as that multiple: Rubi's 1.3.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ASumOfRootsWithAMonomialDifferenceIntegralTest + { + [Theory] + [InlineData("x^4/(sqrt(a + b*x) + sqrt(a + c*x))^3")] + [InlineData("x^2/(sqrt(a + b*x) + sqrt(a + c*x))^2")] + public void RationalisedOverAPowerOfX(string integrand) + { + 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).Substitute("c", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/AnAlgebraicFunctionOfTheReciprocalOfXIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AnAlgebraicFunctionOfTheReciprocalOfXIntegralTest.cs new file mode 100644 index 000000000..0bb84eed2 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AnAlgebraicFunctionOfTheReciprocalOfXIntegralTest.cs @@ -0,0 +1,49 @@ +// +// 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 +{ + /// + /// An algebraic function of 1/x alone with a root in it, under x = 1/u: Rubi's 1.1.3.3, which ran + /// past the budget with its linears written once. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AnAlgebraicFunctionOfTheReciprocalOfXIntegralTest + { + [Theory] + [InlineData("1/((a + b/x)^(5/2)*(c + d/x)^2)")] + [InlineData("1/((a + b/x)^(5/2)*(c + d/x)^3)")] + public void UnderTheReciprocal(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 10000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -2.9, -1.2, 0.3, 0.8, 1.3, 2.9 }) + { + 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/AnOddRationalFunctionBesideAPowerOfXIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AnOddRationalFunctionBesideAPowerOfXIntegralTest.cs new file mode 100644 index 000000000..899ad4fb3 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AnOddRationalFunctionBesideAPowerOfXIntegralTest.cs @@ -0,0 +1,51 @@ +// +// 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 +{ + /// + /// An odd rational function over a power of x and a whole power of a binomial in x^2, left to + /// the substitution under u = x^2 rather than lowered a power of x at a time: Rubi's 1.1.2.2 and + /// 1.1.3.4, which ran past the budget. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AnOddRationalFunctionBesideAPowerOfXIntegralTest + { + [Theory] + [InlineData("1/(x^7*(a + b*x^2)^10)")] + [InlineData("(b + 2*c*x^2)/(x^15*(b + c*x^2)^8)")] + [InlineData("1/(x^5*(a + b*x^2)^10)")] + public void AnsweredUnderTheSquare(string integrand) + { + 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).Substitute("c", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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/HighWholePowersOfTwoLinearsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/HighWholePowersOfTwoLinearsIntegralTest.cs new file mode 100644 index 000000000..5aa106fca --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/HighWholePowersOfTwoLinearsIntegralTest.cs @@ -0,0 +1,51 @@ +// +// 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 high whole powers of two linears with a symbol in them, under the + /// quotient of the linears: Rubi's 1.1.1.3, where the partial fractions ran past the budget. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class HighWholePowersOfTwoLinearsIntegralTest + { + [Theory] + [InlineData("x^2/((1 - a*x)^16*(1 + a*x)^11)")] + [InlineData("x^2/((1 - a*x)^14*(1 + a*x)^10)")] + [InlineData("(c + x)^3/((a + b*x)^10*(c + d*x)^11)")] + [InlineData("1/((a + b*x)^10*(c + d*x)^11)")] + public void WrittenInTheQuotientOfTheLinears(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 10000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", -0.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.6, 1.3, 2.9 }) + { + 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/TheVariableOverTheRootOfAnEvenQuarticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TheVariableOverTheRootOfAnEvenQuarticIntegralTest.cs new file mode 100644 index 000000000..f3aa02c86 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TheVariableOverTheRootOfAnEvenQuarticIntegralTest.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 rational function beside a half-odd power of an even quartic P, as a rational function + /// of t = x/sqrt(P): the pseudo-elliptic integrals of Rubi's 1.2.2.7 and 1.3.2, which + /// no rule wrote in t. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class TheVariableOverTheRootOfAnEvenQuarticIntegralTest + { + [Theory] + [InlineData("(a - c*x^4)/((a*d + a*k*x^2 + c*d*x^4)*sqrt(a + b*x^2 + c*x^4))")] + [InlineData("(a - c*x^4)/((a*d + a*k*x^2 + c*d*x^4)*sqrt(a - b*x^2 + c*x^4))")] + [InlineData("(a - c*x^4)/((a*k + c*d*x^2)*(d + k*x^2)*sqrt(a + b*x^2 + c*x^4))")] + [InlineData("sqrt(a + b*x^2 + c*x^4)/(a*d - c*d*x^4)")] + [InlineData("sqrt(a + b*x^2 - c*x^4)/(a*d + c*d*x^4)")] + [InlineData("sqrt(a + b*x^4)/(a*c - b*c*x^4)")] + [InlineData("sqrt(a - b*x^4)/(a*c + b*c*x^4)")] + [InlineData("x^2/((1 + x^4)*sqrt(-1 + x^4))")] + public void WrittenInTheVariableOverTheRoot(string integrand) + { + 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).Substitute("c", 0.4).Substitute("d", 0.3).Substitute("k", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -2.9, -1.2, -0.7, 0.3, 0.8, 1.7, 2.9 }) + { + 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}"); + } + + [Theory] + [InlineData("1/sqrt(a + b*x^2 + c*x^4)")] + [InlineData("x^2/sqrt(1 + x^4)")] + public void AnEllipticIntegralIsLeftAlone(string integrand) + => Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize()); + } +} diff --git a/Sources/Tests/UnitTests/Calculus/TwoRootsAsOneUpToAConstantFactorIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TwoRootsAsOneUpToAConstantFactorIntegralTest.cs new file mode 100644 index 000000000..478c886c9 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TwoRootsAsOneUpToAConstantFactorIntegralTest.cs @@ -0,0 +1,53 @@ +// +// 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 square roots written as one up to a factor whose square is a constant, the factor in + /// front: where their radicands' quotient is a constant times an even power of x, or their + /// product is free of a root of x each holds. Rubi's 1.1.3.3, 1.2.1.6 and 1.3.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class TwoRootsAsOneUpToAConstantFactorIntegralTest + { + [Theory] + [InlineData("(g + h*x)*sqrt(a + b*x + c*x^2)/(a*d + b*d*x + c*d*x^2)^(3/2)")] + [InlineData("sqrt(b - a/x^2)/sqrt(a - b*x^2)")] + [InlineData("1/(sqrt(1 + x)*sqrt(-1 - sqrt(x))*sqrt(-1 + sqrt(x)))")] + [InlineData("1/(sqrt(a^2 + b^2*x)*sqrt(a - b*sqrt(x))*sqrt(a + b*sqrt(x)))")] + public void WrittenAsOneRoot(string integrand) + { + 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).Substitute("c", 0.4).Substitute("d", -0.3) + .Substitute("g", 0.9).Substitute("h", 1.1); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + 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}"); + } + } +}