diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 391fb969b..b96e19a3f 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -2923,6 +2923,55 @@ biquadratic that the rule above makes, and the terms in `1/(x^2 - r)^j` go by th | `"sqrt(c + d*tan(x))/(a + b*tan(x))^2".ToEntity().Integrate("x")` | `integral(...)`; 40,540 characters on the unreleased master | 1,465 characters | | `"(c + d*tan(x))^(3/2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)`; 71,293 characters on the unreleased master | 2,786 characters | +### A power of x that is not whole beside binomials in a power of it, and thirteen more algebraic shapes, are integrated + +**Answers where there were none.** Rubi's algebraic integrands of these shapes were declined or ran out of time +([#233](https://github.com/asc-community/AngouriMath/issues/233)): +- `x^q R(x^n)`, `q` not whole and `R` over whole powers of two binomials in `x^n`, split into partial fractions in + `w = x^n` with `x^r` kept beside each term: `x^(9/2)/((a + b x^2)(c + d x^2)^2)` is `sqrt(x)` times + `w^2/((a + b w)(c + d w)^2)` written apart, three binomial differentials. +- `x^q P(x^n) (a + b x^n)^p`, and the same over a trinomial in `x^n`, for a whole `p ≤ -2`, by the reduction that + raises `p` to `-1`: `1/(x^(3/2) (a + c x^4)^2)`, `1/(x^(3/2) (a + b x^2 + c x^4)^2)`. +- `R(Q) Q^(k/2)` for a quadratic `Q`, under `t = Q'/sqrt(Q)`, where `Q = D/(4 c - t^2)` and `dx/sqrt(Q) = 2 dt/(4 c - t^2)` + with `D = 4 a c - b^2`, past the fourth power the reduction over a power of the quadratic answers: + `1/((d + b x + c x^2)^5 sqrt(a + b x + c x^2))`. +- A linear below the cube roots of two linears whose product is a quadratic, read as the cube root of the quadratic: + `1/((a + b x)(c + d x)^(1/3)(b c + a d + 2 b d x)^(1/3))`. Two equal half-odd powers of linears beside a power of a + quadratic, read as one power of their product: `1/((a + b x + c x^2)^2 sqrt(1 - d x) sqrt(1 + d x))`. +- A linear over a linear beside the reciprocal root of an even quartic, at the coefficients where the integral is an + inverse hyperbolic tangent or an arctangent of `N^2/(q sqrt(P))`: + `(1 + x - sqrt(3))/((1 + x + sqrt(3)) sqrt(-4 + x^4 + 4 sqrt(3) x^2))`. +- A rational function beside the reciprocal root of a palindromic quartic, in `u = x + 1/x`: + `(k f - k f x^2)/((a d + b d x + a d x^2) sqrt(a + b x + c x^2 + b x^3 + a x^4))`. +- `P(x) (a + b x)^m (c + d x)^n` with `m + n = -1`, by the reduction that leaves the two powers alone: + `(k + f x)^3/((a + b x)^(1/3)(c + d x)^(2/3))`. +- A rational function with symbols in it beside a half-odd power of a quadratic, split first over the factors below + the bar as written: `(a + b x + c x^2)^(3/2)/((d + k x)(f + g x)^2)`. +- A scaling of the variable that leaves numbers only, `x = (c/d) y`: `1/((c + d x)(2 c^3 + d^3 x^3)^(1/3))`. +- A power below `-1` of an odd root of a binomial in `x^n`, beside another binomial, divided by the root first: + `1/(x^9 (a + b x^3)^(4/3)(c + d x^3))`. +- A product under a nested root written out: `sqrt(x (a x + b sqrt(-a/b^2 + a^2 x^2/b^2)))/(x sqrt(-a/b^2 + a^2 x^2/b^2))`. +- A sum whose terms share a power of one polynomial that is not whole, written as that power times the rest: + `(1 - x^6)^(2/3) + (1 - x^6)^(2/3)/x^6` is `-(1 - x^6)^(5/3)/(5 x^5)` differentiated. + +The depth of nested questions is sixteen where it was twelve. Twenty-three of 1.1.2.4's +`x^(k/2)/((a + b x^2)^j (c + d x^2)^l)`, which 2.5.0 declined and the unreleased master answered, change form: +`x^(7/2)/((a + b x^2)(c + d x^2)^2)` is 799 characters there and 1,246 here, and +`1/(x^(3/2) (a + b x^2)(c + d x^2)^2)` is 1,577 there and 1,367 here. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^(9/2)/((a + b*x^2)*(c + d*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,313 characters | +| `"1/(x^(3/2)*(a + c*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,237 characters | +| `"1/((d + b*x + c*x^2)^5*sqrt(a + b*x + c*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 10,616 characters | +| `"1/((a + b*x)*(c + d*x)^(1/3)*(b*c + a*d + 2*b*d*x)^(1/3))".ToEntity().Integrate("x")` | `integral(...)` | 1,120 characters | +| `"(1 + x - sqrt(3))/((1 + x + sqrt(3))*sqrt(-4 + x^4 + 4*x^2*sqrt(3)))".ToEntity().Integrate("x")` | `integral(...)` | 241 characters | +| `"(k*f - k*f*x^2)/((a*d + b*d*x + a*d*x^2)*sqrt(a + b*x + c*x^2 + b*x^3 + a*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 1,460 characters | +| `"(k + f*x)^3/((a + b*x)^(1/3)*(c + d*x)^(2/3))".ToEntity().Integrate("x")` | `integral(...)` | 1,413 characters | +| `"1/((c + d*x)*(2*c^3 + d^3*x^3)^(1/3))".ToEntity().Integrate("x")` | `integral(...)` | 754 characters | +| `"1/(x^9*(a + b*x^3)^(4/3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 3,205 characters | +| `"x^(7/2)/((a+b*x^2)*(c+d*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,246 characters, where the unreleased master gives 799 | + ### A quotient in `x^2` over a power of `x` and a biquadratic is split in `x^2` **Shorter answers, sooner.** `cot(x)^(13/2) (a + b tan(x))^(5/2) (A + B tan(x))` was declined in 2.5.0 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 602536f0d..3106762f3 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -7686,6 +7686,60 @@ node is Powf(var xAgain, Number.Integer power) && xAgain == x && power.EInteger. return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// asked at the question asked before the polynomial + /// term, for a power below -1 of a cube root or a higher odd one of a binomial in x^n + /// beside a power of x of x^(-n) or below and a whole power of a second binomial in + /// x^n below the bar: Rubi's 1/(x^9 (a + b x^3)^(4/3) (c + d x^3)). + /// + /// + /// Asked where it stands, a level down from the polynomial term and the cancellation with + /// functions as indeterminates, its rational function in u = x/(a + b x^3)^(1/3) ran past + /// the budget there, where at the question asked it is answered in two seconds. Only this + /// shape, where the division answered nothing in time; with square roots, and with the + /// power of x above x^(-n), the integrand is answered as it was. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveByDividingByAHighPowerOfARootBesideAnotherBinomialFirst(Entity expr, Entity.Variable x) + { + if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x)) + return null; + var power = ERational.Zero; + var n = 0; + ERational? root = null; + var others = 0; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + continue; + if (factor == x || factor is Powf(var variable, Number.Rational) && variable == x) + { + var exponent = factor is Powf(_, Number.Rational { ERational: var e }) ? e : ERational.One; + power = underneath ? power.Subtract(exponent) : power.Add(exponent); + continue; + } + var (@base, written) = factor is Powf(var raised, Number.Rational { ERational: var p }) ? (raised, p) : (factor, ERational.One); + if (underneath) + written = written.Negate(); + if (!TryReadAsABinomialIn(@base, x, out _, out _, out var degree) || n != 0 && degree != n) + return null; + n = degree; + if (!written.Denominator.Equals(EInteger.One)) + { + if (root is not null) + return null; + root = written; + } + else if (written.Sign < 0) + others++; + else + return null; + } + return root is { } r && n >= 3 && n % 2 == 1 && r.Denominator.Equals(EInteger.FromInt32(n)) && r.CompareTo(ERational.FromInt32(-1)) < 0 + && others > 0 && power.CompareTo(ERational.FromInt32(-n)) <= 0 + ? SolveByDividingByTheRoot(expr, x) : null; + } + /// /// with each base that is a sum written as its polynomial in /// , and each factor then written alike above and below the bar @@ -7908,6 +7962,87 @@ Entity Symbolic(Entity side) => Mulf.LinearChildren(side).Where(factor => factor return number.Sign > 0 ? (MathS.Sqrt(Number.Rational.Create(number)) * root).InnerSimplified : null; } + /// + /// A linear over a linear beside the reciprocal of the square root of an even quartic, + /// N/(M sqrt(P)) with P = a + b x^2 + c x^4, at the coefficients where its + /// integral is K artanh(N^2/(q sqrt(P))) or K atan(N^2/(q sqrt(P))): Rubi's 1.2.2.3 + /// (1 + x - sqrt(3))/((1 + x + sqrt(3)) sqrt(-4 + x^4 + 4 sqrt(3) x^2)). + /// + /// + /// + /// The derivative of artanh(N^2/(q sqrt(P))) is + /// q N (2 N' P - N P'/2)/(sqrt(P) (q^2 P - N^4)), and that of the arctangent is the same + /// with -q^2 for q^2 and -N^4 for N^4. So the integrand is K times + /// one of them where sigma T = rho P - N^4 as polynomials, for T = (2 N' P - N P'/2) M, + /// with sigma = K q and rho = q^2, or sigma = -K q and rho = -q^2. + /// The top terms of 2 N' P - N P'/2 cancel, so both sides have degree four, and the + /// five coefficients are a linear system in sigma and rho: solved from the first and the last, + /// it holds on the rest exactly where Rubi's three conditions on the coefficients do, and + /// otherwise the integral is elliptic. The answer is differentiated back at sampled points, + /// which is where the rest of the system is checked. These were declined: no rule read the + /// square of the numerator over the root. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquare(Entity expr, Entity.Variable x) + { + Entity constant = Number.Integer.One; + Entity? numerator = null, denominator = null, quartic = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + constant = underneath ? constant / factor : constant * factor; + else if (factor is Powf(var @base, Number.Rational { ERational: var power }) && quartic is null + && power.CompareTo(underneath ? ERational.Create(1, 2) : ERational.Create(-1, 2)) == 0) + quartic = @base; + else if (!underneath && numerator is null && factor is not Powf) + numerator = factor; + else if (underneath && denominator is null && factor is not Powf) + denominator = factor; + else + return null; + } + if (numerator is null || denominator is null || quartic is null + || !TreeAnalyzer.TryGetPolyLinear(numerator, x, out var n1, out var n0) || !TreeAnalyzer.TryGetPolyLinear(denominator, x, out var m1, out var m0) + || !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) + || new[] { n0, n1, m0, m1 }.Any(coefficient => coefficient.ContainsNode(x))) + 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) || Zero(n1) || Zero(m1) || Zero(n0 * m1 + n1 * m0)) + return null; + // From the constant terms and the fourth: sigma 2 n1 a m0 - rho a = -n0^4 and + // -sigma 2 n0 c m1 - rho c = -n1^4. + var determinant = -2 * a * c * (n1 * m0 + n0 * m1); + var sigma = ((c * MathS.Pow(n0, 4) - a * MathS.Pow(n1, 4)) / determinant).InnerSimplified; + var rho = ((-2 * a * m0 * MathS.Pow(n1, 5) - 2 * c * m1 * MathS.Pow(n0, 5)) / determinant).InnerSimplified; + // Numbers written as they simplify: `2 - sqrt(3)` for `((1 - sqrt(3))^4 + 4)/16 ...`. + if (!sigma.Vars.Any()) + sigma = sigma.Simplify(); + if (!rho.Vars.Any()) + rho = rho.Simplify(); + if (Zero(rho)) + return null; + var square = MathS.Sqr(numerator) / MathS.Sqrt(quartic); + Entity answer; + if (!rho.Vars.Any() && rho.Evaled is Number.Real { IsNegative: true }) + { + var q = MathS.Sqrt((-rho).InnerSimplified); + answer = constant * (-sigma / q) * MathS.Arctan(square / q); + } + else + { + var q = MathS.Sqrt(rho); + answer = constant * (sigma / q) * MathS.Hyperbolic.Artanh(square / q); + } + answer = answer.InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + /// /// 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 @@ -14963,9 +15098,11 @@ void Gather(Entity node, int times) /// /// The most negative whole power /// takes down - /// the recurrence: each step is one term more in the answer. + /// the recurrence: each step is one term more in the answer. Sixteen, since Rubi's + /// 1/(x^4 (b x^(2/3) + a x)^(3/2)) is 3/(u^13 (b + a u)^(3/2)) under u = x^(1/3), and + /// past the bound it ran past the budget in the cancellation with functions as indeterminates. /// - private const int MaximumRecurrenceDepth = 12; + private const int MaximumRecurrenceDepth = 16; /// /// An integrand holding a fractional power of something linear in the variable, @@ -25625,7 +25762,10 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres return null; var three = Number.Integer.Create(3); var root3 = MathS.Sqrt(three); - if (VanishesIdentically((c * dShifted * dShifted + 9 * aShifted * e * e).Expand().InnerSimplified.Expand())) + // Each condition decided as a value: simplified, a true one was not told from zero where + // the quadratic is a product of two linears with four symbols in them. + static bool Holds(Entity condition) => VanishesIdentically(condition.Expand().InnerSimplified.Expand()) || Functions.PartialFractions.IsZeroAsAValue(condition); + if (Holds(c * dShifted * dShifted + 9 * aShifted * e * e)) { var k = MathS.Pow(aShifted, Number.Rational.Create(1, 3)); var cubeRoot = MathS.Pow(quadratic, Number.Rational.Create(1, 3)); @@ -25638,7 +25778,7 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres closed = closed.InnerSimplified; return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(closed, expr, x) ? closed : null; } - if (!VanishesIdentically((c * dShifted * dShifted - 3 * aShifted * e * e).Expand().InnerSimplified.Expand())) + if (!Holds(c * dShifted * dShifted - 3 * aShifted * e * e)) return null; var q = MathS.Pow(6 * c * c * e * e / (dShifted * dShifted), Number.Rational.Create(1, 3)); var r = MathS.Pow(quadratic, Number.Rational.Create(1, 3)); @@ -25652,6 +25792,297 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; } + /// + /// A linear below the bar beside the cube roots of two other linears whose product is a + /// quadratic, K/(L (A^(1/3) B^(1/3))), as the linear below the cube root of + /// A B that answers, + /// with A^(1/3) B^(1/3) written back for that root in the answer: Rubi's + /// 1/((a + b x) (c + d x)^(1/3) (b c + a d + 2 b d x)^(1/3)) and + /// 1/((d - 3 e x)^(1/3) (d + e x) (d + 3 e x)^(1/3)). + /// + /// + /// The closed form is F(x, R) with logarithms and an arctangent of rational functions + /// of x and R, and its derivative is the integrand wherever R^3 = A B and + /// R' = (A B)' R/(3 A B); A^(1/3) B^(1/3) is such an R on every branch, + /// where (A B)^(1/3) is it only where A and B are not both negative. Both were declined: + /// the closed form reads one cube root, and , which + /// joins two, does it for roots it can sign. The answer is differentiated back. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadratic(Entity expr, Entity.Variable x) + { + Entity constant = Number.Integer.One; + Entity? linear = null; + var roots = new List(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (!underneath && factor is Powf(var inverse, Number.Rational { ERational: var minus }) && minus.CompareTo(ERational.Create(-1, 3)) == 0) + roots.Add(inverse); + else if (underneath && factor is Powf(var @base, Number.Rational { ERational: var third }) && third.CompareTo(ERational.Create(1, 3)) == 0) + roots.Add(@base); + else if (underneath && linear is null && !(factor is Powf)) + linear = factor; + else + return null; + } + if (linear is null || roots.Count != 2 || roots[0] == roots[1] + || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out var intercept) || slope.ContainsNode(x) || intercept.ContainsNode(x) + || roots.Any(root => !TreeAnalyzer.TryGetPolyLinear(root, x, out var s, out var i) || s.ContainsNode(x) || i.ContainsNode(x) + || Functions.PartialFractions.IsZeroAsAValue(s))) + return null; + var quadratic = Functions.PartialFractions.Bare((roots[0] * roots[1]).Expand().InnerSimplified); + if (SolveALinearBelowACubeRootOfAQuadraticInClosedForm(constant / (linear * MathS.Pow(quadratic, Number.Rational.Create(1, 3))), x) is not { } joined) + return null; + var asTwo = MathS.Pow(roots[0], Number.Rational.Create(1, 3)) * MathS.Pow(roots[1], Number.Rational.Create(1, 3)); + var answer = WithARadicandWrittenBack(joined, quadratic, x, 3, power => MathS.Pow(asTwo, Number.Integer.Create(power))).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// + /// with every power P^(k/m) of a polynomial that is + /// as a value, coefficient by coefficient, written as + /// of k, the power of the m-th root it is. + /// + private static Entity WithARadicandWrittenBack(Entity answer, Entity radicand, Entity.Variable x, int order, System.Func writtenFor) + { + if (!TreeAnalyzer.TryGetPolynomial(radicand, x, out var reference)) + return answer; + return answer.Replace(node => + { + if (node is not Powf(var @base, Number.Rational { ERational: var power }) || !@base.ContainsNode(x)) + return node; + var scaled = power.Multiply(ERational.FromInt32(order)).ToLowestTerms(); + if (!scaled.Denominator.Equals(EInteger.One) || !scaled.Numerator.CanFitInInt32() + || !TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Count != reference.Count + || read.Any(pair => !reference.TryGetValue(pair.Key, out var theirs) || !Functions.PartialFractions.IsZeroAsAValue(pair.Value - theirs))) + return node; + return writtenFor(scaled.Numerator.ToInt32Checked()); + }); + } + + /// + /// An algebraic function whose symbols a scaling of the variable takes out, x = s y + /// with s = c/d for a linear c + d x in it, or s = (a/b)^(1/n) for an + /// odd n and a binomial a + b x^n: where every polynomial in it is then a constant + /// times one with rational coefficients, the integral is a constant times one with numbers, + /// which the logarithmic ansatz for a cube root answers. Rubi's + /// 1/((c + d x) (2 c^3 + d^3 x^3)^(1/3)) is 1/(c d) int dy/((1 + y) (2 + y^3)^(1/3)), + /// and (a + b x^3)^(1/3)/(a - b x^3) is a^(-2/3) (a/b)^(1/3) int (1 + y^3)^(1/3)/(1 - y^3) dy. + /// + /// + /// + /// With P(s y) = mu P~(y) for each polynomial P, the integrand is + /// s prod mu^e g(y), where a root P^(k/m) stands as R^k for + /// R = P^(1/m)/mu^(1/m), whose m-th power is P~(y) on every branch. The answer + /// in y is a function of y and the roots of the P~ whose derivative holds wherever their + /// m-th powers are the P~, so it is written back with each root of a P~ as that + /// R and y as x/s, rather than as a root of P(x)/mu, which is it only + /// on the principal branch where mu is positive. A polynomial above the bar that does + /// not scale to numbers is written out in y, each term asked on its own. + /// + /// + /// The ansatz reads one cube root with numbers for coefficients, and these were declined + /// with symbols in them; with numbers, 1/((1 + 2^(1/3) x) (1 + x^3)^(2/3)) ran past the + /// budget in it, where in y = 2^(1/3) x it is 1/((1 + y) (2 + y^3)^(2/3)) and + /// is answered in a third of a second. Late, before the ansatz, and at the question asked; + /// the answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveByScalingTheVariableToNumbers(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + // Algebraic in x, with a root of a polynomial in it. + foreach (var node in expr.Nodes) + if (node.ContainsNode(x) && node is not (Variable or Sumf or Minusf or Mulf or Divf or Powf) + || node is Powf(var @base, var exponent) && @base.ContainsNode(x) && (exponent.ContainsNode(x) || exponent is not Number.Rational)) + return null; + Entity constant = Number.Integer.One; + var factors = new List<(Entity Base, ERational Power, Dictionary Read)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, power) = factor is Powf(var raised, Number.Rational { ERational: var p }) ? (raised, p) : (factor, ERational.One); + if (!TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Keys.Any(degree => degree.Sign < 0 || !degree.CanFitInInt32()) + || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + factors.Add((@base, underneath ? power.Negate() : power, read)); + } + if (!factors.Any(factor => !factor.Power.Denominator.Equals(EInteger.One))) + return null; + var symbolic = expr.Vars.Any(symbol => symbol != x); + // The scalings to try: one from each linear, and one from each binomial in an odd power. + var scalings = new List<(Entity Scale, int Order)>(); + foreach (var (_, _, read) in factors) + { + if (read.Count != 2 || !read.TryGetValue(EInteger.Zero, out var low) || Functions.PartialFractions.IsZeroAsAValue(low)) + continue; + var (degree, high) = read.Select(pair => (pair.Key, pair.Value)).First(pair => !pair.Key.IsZero); + if (Functions.PartialFractions.IsZeroAsAValue(high)) + continue; + if (degree.Equals(EInteger.One)) + { + // With numbers, only a linear whose ratio is irrational: the ansatz reads a + // rational one as it stands. + if (symbolic || (low / high).InnerSimplified is not Number.Rational) + scalings.Add((Functions.PartialFractions.Bare((low / high).InnerSimplified), 1)); + } + else if (symbolic && degree.CompareTo(EInteger.FromInt32(MostOrderScaledToNumbers)) <= 0 && !degree.IsEven) + scalings.Add((Functions.PartialFractions.Bare((low / high).InnerSimplified), degree.ToInt32Checked())); + } + var y = Entity.Variable.CreateUnique(expr, "y_scaled"); + foreach (var (ratio, order) in scalings.Distinct()) + if (Scaled(ratio, order) is { } answer) + return answer; + return null; + + Entity? Scaled(Entity ratio, int order) + { + // s^k for a whole k: a power of the ratio where k is a multiple of the order, and + // otherwise a power of its root. + Entity ToThe(EInteger k) + => k.IsZero ? Number.Integer.One + : k.Remainder(EInteger.FromInt32(order)).IsZero + ? MathS.Pow(ratio, Number.Integer.Create(k.Divide(EInteger.FromInt32(order)))) + : MathS.Pow(ratio, Number.Rational.Create(k, EInteger.FromInt32(order))); + var scale = order == 1 ? ratio : MathS.Pow(ratio, Number.Rational.Create(1, order)); + Entity inFront = scale * constant; + Entity inY = Number.Integer.One; + var aboveInY = new List(); + var written = new List<(Entity Base, Entity Mu, Dictionary Read, int Order)>(); + foreach (var (@base, power, read) in factors) + { + var bare = read.Count == 1 && read.TryGetValue(EInteger.One, out var unit) && Functions.PartialFractions.IsZeroAsAValue(unit - 1); + if (bare) + { + // x^p is s^p y^p. + inFront *= MathS.Pow(scale, Number.Rational.Create(power)); + inY *= MathS.Pow(y, Number.Rational.Create(power)); + continue; + } + var lowest = read.Keys.Min()!; + var mu = read[lowest] * ToThe(lowest); + var inNumbers = new Dictionary(); + var allRational = true; + foreach (var (degree, coefficient) in (read).Select(pair => (pair.Key, pair.Value))) + { + if (order != 1 && !degree.Remainder(EInteger.FromInt32(order)).IsZero) + { + allRational = false; + break; + } + if (AsARationalNumber(coefficient * ToThe(degree) / mu) is not { } rational) + { + allRational = false; + break; + } + inNumbers[degree] = rational; + } + if (!allRational) + { + // A polynomial above the bar, written out in y. + if (power.Denominator.Equals(EInteger.One) && power.Sign > 0) + { + aboveInY.Add(MathS.Pow(@base.Substitute(x, scale * y), Number.Integer.Create(power.Numerator))); + continue; + } + return null; + } + Entity tilde = Number.Integer.Zero; + foreach (var (degree, rational) in (inNumbers.OrderBy(pair => pair.Key)).Select(pair => (pair.Key, pair.Value))) + tilde += degree.IsZero ? rational : rational * MathS.Pow(y, Number.Integer.Create(degree)); + tilde = tilde.InnerSimplified; + inFront *= MathS.Pow(mu, Number.Rational.Create(power)); + inY *= MathS.Pow(tilde, Number.Rational.Create(power)); + if (!power.Denominator.Equals(EInteger.One)) + { + if (!power.Denominator.CanFitInInt32()) + return null; + written.Add((@base, mu, inNumbers, power.Denominator.ToInt32Checked())); + } + } + if (written.Count == 0) + return null; + // Two roots of one P~ are two roots that need not be the same branch: the integrand in y + // multiplies them into a power of P~ itself, and the answer has no root of it left to + // write back. `1/(sqrt(-2 - b x) sqrt(2 + b x))` is `1/(sqrt(-2) sqrt(2) (1 + y))` only + // where `1 + y > 0`, and its derivative had the other sign on the rest of the line. + for (var i = 0; i < written.Count; i++) + for (var j = i + 1; j < written.Count; j++) + if (written[i].Read.Count == written[j].Read.Count + && written[i].Read.All(pair => written[j].Read.TryGetValue(pair.Key, out var other) && other == pair.Value)) + return null; + // The polynomial above, its terms in y with their coefficients in front. + var terms = new List<(Entity Coefficient, Entity Question)>(); + if (aboveInY.Count == 0) + terms.Add((Number.Integer.One, inY)); + else + { + var product = aboveInY.Skip(1).Aggregate(aboveInY[0], (left, right) => left * right); + if (!TreeAnalyzer.TryGetPolynomial(product, y, out var above) || above.Count > MostTermsScaledToNumbers + || above.Values.Any(coefficient => coefficient.ContainsNode(y))) + return null; + foreach (var (degree, coefficient) in (above).Select(pair => (pair.Key, pair.Value))) + if (!Functions.PartialFractions.IsZeroAsAValue(coefficient)) + terms.Add((LowestOverTheSymbols(coefficient), degree.IsZero ? inY : MathS.Pow(y, Number.Integer.Create(degree)) * inY)); + } + if (terms.Any(term => term.Question.Vars.Any(symbol => symbol != y))) + return null; + Entity sum = Number.Integer.Zero; + foreach (var (coefficient, question) in terms) + { + if (Integration.ComputeAsAQuestionOfItsOwn(question.InnerSimplified, y, integrateByParts) is not { } integral + || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += coefficient * integral; + } + // Each root of a P~ as P^(1/m)/mu^(1/m), and y as x/s. + foreach (var (@base, mu, inNumbers, m) in written) + { + Entity tilde = Number.Integer.Zero; + foreach (var (degree, rational) in (inNumbers).Select(pair => (pair.Key, pair.Value))) + tilde += degree.IsZero ? rational : rational * MathS.Pow(y, Number.Integer.Create(degree)); + var root = MathS.Pow(@base, Number.Rational.Create(1, m)) / MathS.Pow(mu, Number.Rational.Create(1, m)); + sum = WithARadicandWrittenBack(sum, tilde, y, m, k => MathS.Pow(root, Number.Integer.Create(k))); + } + var answer = (inFront * sum.Substitute(y, x / scale)).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || answer.ContainsNode(y) + || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + } + + /// + /// as the rational number it is, where it is one: without symbols, + /// as it simplifies, and with them, its value at pinned symbols where the difference from + /// that is zero as a value. + /// + private static Number.Rational? AsARationalNumber(Entity value) + { + var simplified = Functions.PartialFractions.Bare(value.InnerSimplified); + if (!simplified.Vars.Any()) + return simplified.Evaled as Number.Rational; + var pinned = simplified; + var index = 0; + foreach (var symbol in simplified.Vars.ToList()) + pinned = pinned.Substitute(symbol, Number.Rational.Create(EInteger.FromInt32(137 + 100 * index++), EInteger.FromInt32(100))); + return pinned.Evaled is Number.Rational rational && Functions.PartialFractions.IsZeroAsAValue(simplified - rational) ? rational : null; + } + + /// The largest odd degree of a binomial scales by, and the most terms it writes out. + private const int MostOrderScaledToNumbers = 9, MostTermsScaledToNumbers = 6; + /// /// Powers that are not whole of polynomials that are constant multiples of powers of one /// linear, gathered into one power of it: B^p = (B^p/L^(k p)) L^(k p) for @@ -28030,67 +28461,33 @@ Entity Without(Entity side) } /// - /// A half-odd power of a constant over a polynomial whose sign changes, among the factors - /// of the integrand, written apart with the factor that keeps it right on both sides of the - /// polynomial's zeros: for a real Q other than zero, (A/Q)^r is - /// K Q^(-r) with K = (A/Q)^r Q^r, whose derivative is zero wherever it has one, - /// and for a positive number A that is A^r times the sign of Q, the - /// root of a negative being i times the root of its modulus on both sides of the bar. - /// - /// - /// sqrt(1/(1 - x^2)) is sgn(1 - x^2)/sqrt(1 - x^2), the arcsine's derivative - /// where it is real, and was declined, with every such power over a polynomial that is not - /// positive: the reading that writes a power of a quotient apart does so only where it is. - /// Late in the chain, after the substitutions, which answer some of these without the - /// factor -- x sqrt(1/(4 - x^2)) is -(1/(4 - x^2))^(-1/2) -- and asked as the - /// same question, as the sign taken out above is. - /// https://github.com/asc-community/AngouriMath/issues/718 - /// - internal static Entity? SolveByWritingAHalfOddPowerOfAReciprocalWithItsSign(Entity expr, Entity.Variable x, bool integrateByParts) - { - foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) - { - if (underneath || factor is not Powf(var @base, Number.Rational exponent) - || !exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) || !@base.ContainsNode(x)) - continue; - var (above, below) = Functions.SingleQuotient.Of(@base); - if (below == Number.Integer.One || above.ContainsNode(x) || !below.ContainsNode(x) - || !TreeAnalyzer.TryGetPolynomial(below, x, out var read) || read.Count == 0 - || read.Keys.Max()!.CompareTo(EInteger.FromInt32(4)) > 0 || IsPositiveForReal(below, x)) - continue; - var positive = above.Evaled is Number.Real { IsPositive: true }; - var inFront = positive ? MathS.Signum(below) : factor * MathS.Pow(below, exponent); - var replacement = positive ? MathS.Pow(above, exponent) * MathS.Pow(below, -exponent) : MathS.Pow(below, -exponent); - var written = expr.Replace(node => node == factor ? replacement : node); - return Integration.ComputeAsTheSameQuestion(written, x, integrateByParts) is { } answer ? inFront * answer : null; - } - return null; - } - - /// - /// x^(n - 1) g(x^n) with a symbolic n, which is g(u)/n under - /// u = x^n: the power in front is the derivative of the power inside up to the - /// constant n, and the rest mentions x only as x^n. Rubi's 6.5.2 and - /// 6.6.2 -- (e x)^(n - 1) (a + b sech(c + d x^n))^p -- write the power in front as - /// (e x)^(n - 1), which is e^(n - 1) x^(n - 1) for a positive e, and - /// the answer says so. A concrete n is the general substitution's; a symbolic - /// one was nobody's, and x^(n - 1) e^(x^n) was declined. - /// https://github.com/asc-community/AngouriMath/issues/718 + /// A power of x that is not whole beside a polynomial in x^n over a whole power + /// below -1 of a binomial in x^n with a symbol in it, x^q P(x^n) (a + b x^n)^p, + /// each term of P by Rubi's reduction that raises p by one, down to + /// int x^(q + k n)/(a + b x^n), which is asked on its own: + /// + /// int x^m B^p = -x^(m + 1) B^(p + 1)/(a n (p + 1)) + (m + 1 + n (p + 1))/(a n (p + 1)) int x^m B^(p + 1) + /// + /// Rubi's 1/(x^(3/2) (a + c x^4)^2) and (A + B x^3)/(x^(3/2) (a + b x^3)^3). /// /// - /// x^(k n - 1) with k half an odd number is the same substitution by half - /// the power: under u = x^(n/2), x^(k n - 1) dx is (2/n) u^(2k - 1) du - /// and g(x^n) is g(u^2). Beside an exponential of x^n that is a moment - /// of the Gaussian: Rubi's 2.3, f^(a + b x^n) x^(-1 + 5n/2), and 6.1.3's - /// x^(-1 + n/2) sinh(a + b x^n). - /// https://github.com/asc-community/AngouriMath/issues/1501 + /// The derivative of x^(m + 1) B^(p + 1) is x^m B^p ((m + 1 + n (p + 1)) B - a n (p + 1)) + /// for any m, which is the step. These were declined: the substitution under the root makes + /// 2/(u^2 (a + c u^8)^2) of the first, and the partial fractions over the eighth roots of + /// a symbolic binomial gave out past the first power, which is answered, as is + /// x^(-3/2)/(a + c x^4) itself in a tenth of a second. 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? SolveByAPowerOfTheVariableTimesAFunctionOfItsPower(Entity expr, Entity.Variable x, bool integrateByParts) + internal static Entity? SolveAFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOne(Entity expr, Entity.Variable x, bool integrateByParts) { - Entity? inFront = null; - Entity? scale = null; - Entity? rest = null; + if (!expr.Vars.Any(symbol => symbol != x)) + return null; Entity constant = Number.Integer.One; + var q = ERational.Zero; + Entity? binomial = null; + var p = 0; + Entity above = Number.Integer.One; foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) { if (!factor.ContainsNode(x)) @@ -28098,22 +28495,606 @@ Entity Without(Entity side) constant = underneath ? constant / factor : constant * factor; continue; } - // The power in front: x^m, or (e x)^m with e free of x, to a symbolic m. - if (inFront is null && !underneath && factor is Powf(var @base, var exponent) && !exponent.ContainsNode(x) - && exponent.Evaled is not Number && exponent.Vars.Any() - && (@base == x || @base is Mulf(var l, var r) && (l == x && !r.ContainsNode(x) || r == x && !l.ContainsNode(x)))) + if (factor == x || factor is Powf(var variable, Number.Rational) && variable == x) { - inFront = exponent; - scale = @base == x ? null : @base is Mulf(var l2, var r2) ? (l2 == x ? r2 : l2) : null; + var exponent = factor is Powf(_, Number.Rational { ERational: var e }) ? e : ERational.One; + q = underneath ? q.Subtract(exponent) : q.Add(exponent); continue; } - var piece = underneath ? Number.Integer.One / factor : factor; - rest = rest is null ? piece : rest * piece; - } - if (inFront is null || rest is null) - return null; - // The rest mentions x only as x^n, one n throughout. - Entity? n = null; + var (raised, whole) = factor is Powf(var @base, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (@base, w.ToInt32Checked()) : (factor, 1); + if (underneath) + whole = -whole; + if (whole < 0) + { + if (binomial is not null && binomial != raised) + return null; + binomial = raised; + p += whole; + } + else + above *= factor; + } + q = q.ToLowestTerms(); + if (binomial is null || q.Denominator.Equals(EInteger.One) || p > -2 || p < -MostBinomialPowersRaisedBesideAFractionalPower + || !TryReadAsABinomialIn(binomial, x, out var a, out var b, out var n) + || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b) + || !TreeAnalyzer.TryGetPolynomial(above, x, out var top) + || top.Any(pair => pair.Value.ContainsNode(x) || pair.Key.Sign < 0 || !pair.Key.Remainder(EInteger.FromInt32(n)).IsZero) + || top.Count > MostBinomialPowersRaisedBesideAFractionalPower) + return null; + Entity sum = Number.Integer.Zero; + foreach (var (degree, coefficient) in (top).Select(pair => (pair.Key, pair.Value))) + { + if (Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + var m = q.Add(ERational.FromEInteger(degree)); + var raisedX = MathS.Pow(x, Number.Rational.Create(m.Add(ERational.One))); + Entity factor = coefficient; + for (var k = p; k <= -2; k++) + { + var below = a * n * (k + 1); + sum += LowestOverTheSymbols(-factor / below) * raisedX * MathS.Pow(binomial, Number.Integer.Create(k + 1)); + factor = LowestOverTheSymbols(factor * (Number.Rational.Create(m.Add(ERational.One)) + n * (k + 1)) / below); + } + if (Functions.PartialFractions.IsZeroAsAValue(factor)) + continue; + if (Integration.ComputeAsAQuestionOfItsOwn(MathS.Pow(x, Number.Rational.Create(m)) / binomial, x, integrateByParts) is not { } integral + || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += factor * 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, and the most terms above it takes. + private const int MostBinomialPowersRaisedBesideAFractionalPower = 12; + + /// + /// A power of x that is not whole beside a whole power below -1 of a trinomial + /// T = a + b x^n + c x^(2 n) with a symbol in it, x^q T^p, by Rubi's reduction + /// that raises p by one, down to int x^(q + k n)/T, each asked on its own: + /// + /// int x^m T^p = -x^(m + 1) (b^2 - 2 a c + b c x^n) T^(p + 1)/(a n (p + 1) (b^2 - 4 a c)) + /// + int x^m T^(p + 1) (b^2 (m + n (p + 1) + 1) - 2 a c (m + 2 n (p + 1) + 1) + b c (m + n (2 p + 3) + 1) x^n)/(a n (p + 1) (b^2 - 4 a c)) + /// + /// Rubi's 1/(x^(3/2) (a + b x^2 + c x^4)^2). + /// + /// + /// The step is the derivative of its first term, for any m, as for the binomial in + /// . It was + /// declined, where x^(-3/2)/(a + b x^2 + c x^4) and x^(1/2)/(a + b x^2 + c x^4) are + /// answered. Late, beside the binomial's; the answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveAFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOne(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + var q = ERational.Zero; + Entity? trinomial = null; + var p = 0; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor == x || factor is Powf(var variable, Number.Rational) && variable == x) + { + var exponent = factor is Powf(_, Number.Rational { ERational: var e }) ? e : ERational.One; + q = underneath ? q.Subtract(exponent) : q.Add(exponent); + continue; + } + var (raised, whole) = factor is Powf(var @base, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (@base, w.ToInt32Checked()) : (factor, 1); + if (underneath) + whole = -whole; + if (whole >= 0 || trinomial is not null && trinomial != raised) + return null; + trinomial = raised; + p += whole; + } + q = q.ToLowestTerms(); + if (trinomial is null || q.Denominator.Equals(EInteger.One) || p > -2 || p < -MostTrinomialPowersRaisedBesideAFractionalPower + || !TreeAnalyzer.TryGetPolynomial(trinomial, x, out var read) || read.Count != 3 || read.Values.Any(value => value.ContainsNode(x)) + || !read.TryGetValue(EInteger.Zero, out var a) || read.Keys.Max() is not { } top || !top.CanFitInInt32() || top.IsZero || !top.IsEven) + return null; + var n = top.ToInt32Checked() / 2; + if (!read.TryGetValue(EInteger.FromInt32(n), out var b) || !read.TryGetValue(top, out var c)) + return null; + var discriminant = (b * b - 4 * a * c).InnerSimplified; + if (new[] { a, b, c, discriminant }.Any(Functions.PartialFractions.IsZeroAsAValue)) + return null; + var xToTheN = MathS.Pow(x, Number.Integer.Create(n)); + // The integrals still to reduce at the current power, by the power of x beside them. + var pending = new Dictionary { [q] = Number.Integer.One }; + Entity sum = Number.Integer.Zero; + for (var k = p; k <= -2; k++) + { + var below = a * n * (k + 1) * discriminant; + var next = new Dictionary(); + foreach (var (m, coefficient) in (pending).Select(pair => (pair.Key, pair.Value))) + { + var mPlusOne = Number.Rational.Create(m.Add(ERational.One)); + sum += LowestOverTheSymbols(-coefficient / below) * MathS.Pow(x, mPlusOne) * (b * b - 2 * a * c + b * c * xToTheN) + * MathS.Pow(trinomial, Number.Integer.Create(k + 1)); + var first = b * b * (mPlusOne + n * (k + 1)) - 2 * a * c * (mPlusOne + 2 * n * (k + 1)); + var second = b * c * (mPlusOne + n * (2 * k + 3)); + foreach (var (at, part) in new[] { (m, first), (m.Add(ERational.FromInt32(n)), second) }) + { + var added = LowestOverTheSymbols(coefficient * part / below); + var key = next.Keys.FirstOrDefault(existing => existing.CompareTo(at) == 0) ?? at; + next[key] = next.TryGetValue(key, out var already) ? LowestOverTheSymbols(already + added) : added; + } + } + pending = next; + } + foreach (var (m, coefficient) in (pending).Select(pair => (pair.Key, pair.Value))) + { + if (Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + if (Integration.ComputeAsAQuestionOfItsOwn(MathS.Pow(x, Number.Rational.Create(m)) / trinomial, x, integrateByParts) is not { } integral + || integral.Nodes.Any(node => node is Integralf)) + 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 MostTrinomialPowersRaisedBesideAFractionalPower = 4; + + /// + /// The same half-odd power of two linears beside a power of a quadratic below the bar, + /// A^r B^r R(x), as one power of their product, the integral of (A B)^r R(x) asked, + /// and sqrt(A) sqrt(B) written back in it for the root of A B: Rubi's + /// 1/((a + b x + c x^2)^2 sqrt(1 - d x) sqrt(1 + d x)). + /// + /// + /// The answer is a function of x and the root whose derivative holds for any root whose + /// square is A B, and sqrt(A) sqrt(B) is one on every branch; so it holds as + /// written back, where (A B)^r and A^r B^r differ in sign. Asked of the two roots, + /// the polynomial term and the substitutions after it ran past the budget, where of the one + /// it is answered in two seconds. Only beside a quadratic to the second power below the bar + /// or past it, where it was declined, and at the question asked; the answer is differentiated + /// back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveTwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadratic(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var roots = new List<(Entity Base, ERational Power)>(); + Entity? above = null, below = null; + var aSquaredQuadratic = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (factor is Powf(var @base, Number.Rational { ERational: var power }) && power.Denominator.Equals(EInteger.FromInt32(2)) && @base.ContainsNode(x)) + { + if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out var intercept) || slope.ContainsNode(x) || intercept.ContainsNode(x)) + return null; + roots.Add((@base, underneath ? power.Negate() : power)); + continue; + } + if (factor.ContainsNode(x)) + { + var (raised, whole) = factor is Powf(var inner, Number.Integer { EInteger: var w }) && w.CanFitInInt32() ? (inner, w.ToInt32Checked()) : (factor, 1); + if (!TreeAnalyzer.TryGetPolynomial(raised, x, out var read) || read.Values.Any(coefficient => coefficient.ContainsNode(x)) || read.Keys.Any(degree => degree.Sign < 0)) + return null; + if (underneath && whole >= 2 && read.Keys.Max()!.Equals(EInteger.FromInt32(2))) + aSquaredQuadratic = true; + } + if (underneath) + below = below is null ? factor : below * factor; + else + above = above is null ? factor : above * factor; + } + if (roots.Count != 2 || roots[0].Power.CompareTo(roots[1].Power) != 0 || roots[0].Base == roots[1].Base || !aSquaredQuadratic) + return null; + var product = Functions.PartialFractions.Bare((roots[0].Base * roots[1].Base).Expand().InnerSimplified); + if (!TreeAnalyzer.TryGetPolynomial(product, x, out var readProduct) || readProduct.Keys.Max()!.CompareTo(EInteger.FromInt32(2)) != 0 + || Functions.PartialFractions.IsZeroAsAValue(readProduct[EInteger.FromInt32(2)])) + return null; + // As one quotient, the root below the bar with the rest of what is there: written as a + // product of reciprocals it went to Euler's substitution, sixteen seconds past the budget. + var shared = roots[0].Power; + var root = MathS.Pow(product, Number.Rational.Create(shared.Abs())); + var joinedIntegrand = shared.Sign < 0 + ? (above ?? Number.Integer.One) / (below is null ? root : below * root) + : (above is null ? root : above * root) / (below ?? Number.Integer.One); + if (Integration.ComputeAsAQuestionOfItsOwn(joinedIntegrand, x, integrateByParts) is not { } joined + || joined.Nodes.Any(node => node is Integralf)) + return null; + var asTwo = MathS.Sqrt(roots[0].Base) * MathS.Sqrt(roots[1].Base); + var answer = WithARadicandWrittenBack(joined, product, x, 2, k => MathS.Pow(asTwo, Number.Integer.Create(k))).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// + /// A rational function of a quadratic Q beside a half-odd power of it, + /// R(Q) Q^(k/2) with R's factors Q + delta, under t = Q'/sqrt(Q): with + /// D = 4 a c - b^2, Q = D/(4 c - t^2) and dx/sqrt(Q) = 2 dt/(4 c - t^2), so it + /// is rational in t. Rubi's 1/((d + b x + c x^2)^3 sqrt(a + b x + c x^2)) is + /// 2 (4 c - t^2)^2/(D + 4 c (d - a) - (d - a) t^2)^3. + /// + /// + /// + /// t^2 = (b + 2 c x)^2/Q = 4 c - D/Q, which gives Q, and dt/dx = D/(2 Q^(3/2)); each + /// holds wherever Q is not zero, on either sign of it. The quadratics below the bar differ from + /// the radicand by a constant, so they are linear in 1/(4 c - t^2), and the integral in t is + /// a polynomial over a power of a binomial in t^2, which is answered at once. After the + /// reduction over a power of the quadratic beside the root of another, which answers the third + /// and fourth powers in a third of the characters and declines the fifth; ahead of it, this + /// wrote 4620 characters where that writes 1439. Only the third power below the bar or past it, + /// the lower ones being answered as they are; the answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRoot(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity constant = Number.Integer.One; + Entity? radicand = null; + var half = ERational.Zero; + var quadratics = new List<(Entity Base, int Power)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor is Powf(var root, Number.Rational { ERational: var power }) && power.Denominator.Equals(EInteger.FromInt32(2))) + { + if (radicand is not null) + return null; + radicand = root; + half = underneath ? power.Negate() : power; + continue; + } + var (raised, whole) = factor is Powf(var @base, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (@base, w.ToInt32Checked()) : (factor, 1); + quadratics.Add((raised, underneath ? -whole : whole)); + } + if (radicand is null || quadratics.Count == 0 || !quadratics.Any(pair => pair.Power <= -FewestPowersUnderTheDerivativeOverTheRoot) + || !TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var c, out var b, out var a) || c.ContainsNode(x) || b.ContainsNode(x) || a.ContainsNode(x) + || Functions.PartialFractions.IsZeroAsAValue(c)) + return null; + var discriminant = (4 * a * c - b * b).InnerSimplified; + if (Functions.PartialFractions.IsZeroAsAValue(discriminant)) + return null; + var t = Entity.Variable.CreateUnique(expr, "t_root"); + var below = 4 * c - MathS.Sqr(t); + // Q^(half) dx = Q^(half + 1/2) dx/sqrt(Q), and each power of Q is D/(4 c - t^2). + var whole2 = half.Add(ERational.Create(1, 2)).ToLowestTerms(); + if (!whole2.Denominator.Equals(EInteger.One) || !whole2.Numerator.CanFitInInt32()) + return null; + var exponentOfW = -1 - whole2.Numerator.ToInt32Checked(); + var inFront = constant * 2 * MathS.Pow(discriminant, Number.Integer.Create(whole2.Numerator)); + // Each compound coefficient in t named by a symbol of its own, and written back: the + // integral in t with `4 a c - b^2 + 4 c (d - a)` in it was declined, and with a symbol + // there it is answered at once. + var named = new List<(Entity.Variable Name, Entity Value)>(); + Entity Named(Entity value) + { + value = Functions.PartialFractions.Bare(LowestOverTheSymbols(value.Expand())); + if (value is Entity.Variable or Number) + return value; + var name = Entity.Variable.CreateUnique(expr + t + named.Aggregate((Entity)Number.Integer.Zero, (sum, pair) => sum + pair.Name), "k_root"); + named.Add((name, value)); + return name; + } + Entity inT = Number.Integer.One; + foreach (var (@base, power) in quadratics) + { + if (!TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var c2, out var b2, out var a2) + || !Functions.PartialFractions.IsZeroAsAValue(c2 - c) || !Functions.PartialFractions.IsZeroAsAValue(b2 - b)) + return null; + var delta = Functions.PartialFractions.Bare((a2 - a).InnerSimplified); + // Q + delta = (D + delta (4 c - t^2))/(4 c - t^2). + var written = Functions.PartialFractions.IsZeroAsAValue(delta) + ? Named(discriminant) : Named(discriminant + 4 * c * delta) - Named(delta) * MathS.Sqr(t); + inT *= MathS.Pow(written, Number.Integer.Create(power)); + exponentOfW -= power; + } + inT = Functions.PartialFractions.Bare((inT * MathS.Pow(Named(4 * c) - MathS.Sqr(t), Number.Integer.Create(exponentOfW))).InnerSimplified); + if (Integration.ComputeAsAQuestionOfItsOwn(inT, t, integrateByParts) is not { } inTermsOfT || inTermsOfT.Nodes.Any(node => node is Integralf)) + return null; + foreach (var (name, value) in named) + inTermsOfT = inTermsOfT.Substitute(name, value); + var answer = (inFront * inTermsOfT.Substitute(t, (b + 2 * c * x) / MathS.Sqrt(radicand))).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// The lowest power of a quadratic below the bar takes. + private const int FewestPowersUnderTheDerivativeOverTheRoot = 3; + + /// + /// A sum whose terms share a power that is not whole of one polynomial, up to whole powers + /// of it, written as that power times the sum of the rest: Rubi's + /// (1 - x^6)^(2/3) + (1 - x^6)^(2/3)/x^6 is (1 - x^6)^(2/3) (1 + x^6)/x^6, the + /// derivative of -(1 - x^6)^(5/3)/(5 x^5), where neither term apart is elementary. + /// + /// + /// Linearity takes the sum apart first, and a sum of two binomial differentials outside + /// Chebyshev's cases was declined term by term, while the product is answered by the rules + /// for a polynomial beside a power of a binomial in a thirtieth of a second. At the question + /// asked and late, where every rule reading the sum as written has declined it; the product + /// is asked on its own. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveASumOfTermsSharingAPowerOfAPolynomialAsOneProduct(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || expr is not (Sumf or Minusf)) + return null; + var terms = Sumf.LinearChildren(expr); + if (terms.Count < 2 || terms.Count > MostTermsSharingAPower) + return null; + // The powers that are not whole in the first term, each a candidate. + foreach (var (candidate, _) in FactorsOfTheIntegrand(terms[0])) + { + if (candidate is not Powf(var radicand, Number.Rational { ERational: var first }) || first.Denominator.Equals(EInteger.One) + || !radicand.ContainsNode(x)) + continue; + var rests = new List<(Entity Remaining, ERational Power)>(); + foreach (var term in terms) + { + var power = ERational.Zero; + Entity rest = Number.Integer.One; + var found = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(term)) + if (factor is Powf(var @base, Number.Rational { ERational: var exponent }) && @base == radicand) + { + power = underneath ? power.Subtract(exponent) : power.Add(exponent); + found = true; + } + else + rest = underneath ? rest / factor : rest * factor; + power = power.ToLowestTerms(); + if (!found || !power.Subtract(first).ToLowestTerms().Denominator.Equals(EInteger.One)) + break; + rests.Add((rest, power)); + } + if (rests.Count != terms.Count) + continue; + var lowest = rests.Select(pair => pair.Power).Aggregate((least, power) => power.CompareTo(least) < 0 ? power : least); + Entity sum = Number.Integer.Zero; + foreach (var (rest, power) in rests) + { + var whole = power.Subtract(lowest).ToLowestTerms().Numerator; + sum += whole.IsZero ? rest : rest * MathS.Pow(radicand, Number.Integer.Create(whole)); + } + var product = MathS.Pow(radicand, Number.Rational.Create(lowest)) * Functions.SingleQuotient.Combine(sum); + return Integration.ComputeAsAQuestionOfItsOwn(product, x, integrateByParts); + } + return null; + } + + /// The most terms takes. + private const int MostTermsSharingAPower = 8; + + /// + /// A rational function beside the reciprocal of the square root of a palindromic quartic + /// P = p0 + p1 x + p2 x^2 + p1 x^3 + p0 x^4, where it is a function of + /// u = x + 1/x: P/x^2 = p0 (u^2 - 2) + p1 u + p2 = Q(u), and with + /// S = sqrt(P)/x, a root of Q, R(x) dx/sqrt(P) = h(u) du/S for + /// h = R(x) x/(x^2 - 1). Rubi's + /// (e f - e f x^2)/((a d + b d x + a d x^2) sqrt(a + b x + c x^2 + b x^3 + a x^4)) is + /// -(e f/d) int du/((a u + b) sqrt(a u^2 + b u + c - 2 a)). + /// + /// + /// + /// h is a function of u where it is the same at x and at 1/x, the two roots of + /// X^2 - u X + 1; that is read by writing its numerator and denominator modulo that + /// quadratic, each as beta x + alpha with alpha and beta polynomials in u, the quotient + /// free of x exactly where beta_N alpha_D = alpha_N beta_D, decided as a value. + /// + /// + /// The integral in u is written back with each power of the root of Q as the same power of + /// S = sqrt(P)/x rather than of sqrt(Q(x + 1/x)): S is that root for a positive x + /// and its negative for a negative one, and the answer holds for any root whose square is Q, + /// so with S it holds on both sides of zero, where |x| would have to be taken out + /// otherwise. These were declined: no rule read the quartic's symmetry. The answer is + /// differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocal(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity constant = Number.Integer.One; + Entity? quartic = null; + Entity rational = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor is Powf(var @base, Number.Rational { ERational: var power }) && power.Denominator.Equals(EInteger.FromInt32(2))) + { + if (quartic is not null || power.CompareTo(underneath ? ERational.Create(1, 2) : ERational.Create(-1, 2)) != 0) + return null; + quartic = @base; + continue; + } + var raised = factor is Powf(var inner, Number.Integer) ? inner : factor; + if (!TreeAnalyzer.TryGetPolynomial(raised, x, out var read) || read.Keys.Any(degree => degree.Sign < 0) + || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + rational = underneath ? rational / factor : rational * factor; + } + if (quartic is null || rational == Number.Integer.One || !TreeAnalyzer.TryGetPolynomial(quartic, x, out var p) + || p.Keys.Any(degree => degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(4)) > 0) || p.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + Entity P(int k) => p.TryGetValue(EInteger.FromInt32(k), out var value) ? value : Number.Integer.Zero; + static bool Zero(Entity value) => Functions.PartialFractions.IsZeroAsAValue(value); + if (Zero(P(0)) || !Zero(P(0) - P(4)) || !Zero(P(1) - P(3))) + return null; + // h = R x/(x^2 - 1) over one bar, each side modulo X^2 - u X + 1. + var (top, bottom) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(rational * x / (MathS.Sqr(x) - 1))); + if (!TreeAnalyzer.TryGetPolynomial(top, x, out var above) || !TreeAnalyzer.TryGetPolynomial(bottom, x, out var below) + || above.Concat(below).Any(pair => pair.Key.Sign < 0 || !pair.Key.CanFitInInt32() || pair.Value.ContainsNode(x)) + || above.Keys.Concat(below.Keys).Max()!.CompareTo(EInteger.FromInt32(MostPalindromicDegreeRead)) > 0) + return null; + var u = Entity.Variable.CreateUnique(expr, "u_palindromic"); + var (betaN, alphaN) = Reduced(above); + var (betaD, alphaD) = Reduced(below); + if (!TreeAnalyzer.TryGetPolynomial((betaN * alphaD - alphaN * betaD).Expand(), u, out var cross) || !cross.Values.All(Zero)) + return null; + var (overU, underU) = TreeAnalyzer.TryGetPolynomial(betaD.Expand(), u, out var readBetaD) && readBetaD.Values.All(Zero) ? (alphaN, alphaD) : (betaN, betaD); + var h = Functions.PartialFractions.Bare(LowestOverTheSymbols(overU / underU)); + if (h.ContainsNode(x)) + return null; + var q = (P(0) * (MathS.Sqr(u) - 2) + P(1) * u + P(2)).Expand().InnerSimplified; + if (Integration.ComputeAsAQuestionOfItsOwn(constant * h / MathS.Sqrt(q), u, integrateByParts) is not { } inU || inU.Nodes.Any(node => node is Integralf)) + return null; + var s = MathS.Sqrt(quartic) / x; + var answer = WithARadicandWrittenBack(inU, q, u, 2, k => MathS.Pow(s, Number.Integer.Create(k))).Substitute(u, x + 1 / x).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) || answer.ContainsNode(u) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + + // A polynomial in x modulo X^2 - u X + 1, by Horner's rule from the top: `(b x + a) x + k` + // is `(b u + a) x + (k - b)` there. As (beta, alpha) for `beta x + alpha`. + (Entity, Entity) Reduced(Dictionary polynomial) + { + Entity beta = Number.Integer.Zero, alpha = Number.Integer.Zero; + for (var k = polynomial.Keys.Max()!.ToInt32Checked(); k >= 0; k--) + { + var coefficient = polynomial.TryGetValue(EInteger.FromInt32(k), out var value) ? value : Number.Integer.Zero; + (beta, alpha) = (beta * u + alpha, coefficient - beta); + } + return (beta.Expand(), alpha.Expand()); + } + } + + /// The largest degree reads. + private const int MostPalindromicDegreeRead = 12; + + /// + /// A root of a product of a polynomial and a sum with a root in it, the product written out + /// under the root: Rubi's sqrt(x (a x + b sqrt(-a/b^2 + a^2 x^2/b^2)))/(x sqrt(-a/b^2 + a^2 x^2/b^2)) + /// is sqrt(a x^2 + b x sqrt(...)) over the same, which is answered. + /// + /// + /// Written as a product, the rule that takes a power of x out of a root took x out and left + /// sqrt(x) sqrt(a x + b sqrt(...)), which nothing read. The two radicands are one + /// polynomial in x and the inner root, so the identity is exact. At the question asked and + /// late, after every rule reading the product as written; the rewritten integrand is asked on + /// its own. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveByWritingOutAProductUnderANestedRoot(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var changed = false; + var rewritten = expr.Replace(node => + { + if (node is not Powf(Mulf product, Number.Rational { ERational: var power }) || power.Denominator.Equals(EInteger.One) + || !Mulf.LinearChildren(product).Any(factor => factor is Sumf or Minusf + && factor.Nodes.Any(inner => inner is Powf(var innerBase, Number.Rational { ERational: var innerPower }) && !innerPower.Denominator.Equals(EInteger.One) && innerBase.ContainsNode(x))) + || Mulf.LinearChildren(product).Count(factor => factor.ContainsNode(x)) < 2) + return node; + changed = true; + return MathS.Pow(product.Expand(), Number.Rational.Create(power)); + }); + return changed && rewritten != expr ? Integration.ComputeAsAQuestionOfItsOwn(rewritten, x, integrateByParts) : null; + } + + /// + /// A half-odd power of a constant over a polynomial whose sign changes, among the factors + /// of the integrand, written apart with the factor that keeps it right on both sides of the + /// polynomial's zeros: for a real Q other than zero, (A/Q)^r is + /// K Q^(-r) with K = (A/Q)^r Q^r, whose derivative is zero wherever it has one, + /// and for a positive number A that is A^r times the sign of Q, the + /// root of a negative being i times the root of its modulus on both sides of the bar. + /// + /// + /// sqrt(1/(1 - x^2)) is sgn(1 - x^2)/sqrt(1 - x^2), the arcsine's derivative + /// where it is real, and was declined, with every such power over a polynomial that is not + /// positive: the reading that writes a power of a quotient apart does so only where it is. + /// Late in the chain, after the substitutions, which answer some of these without the + /// factor -- x sqrt(1/(4 - x^2)) is -(1/(4 - x^2))^(-1/2) -- and asked as the + /// same question, as the sign taken out above is. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByWritingAHalfOddPowerOfAReciprocalWithItsSign(Entity expr, Entity.Variable x, bool integrateByParts) + { + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (underneath || factor is not Powf(var @base, Number.Rational exponent) + || !exponent.ERational.Denominator.Equals(EInteger.FromInt32(2)) || !@base.ContainsNode(x)) + continue; + var (above, below) = Functions.SingleQuotient.Of(@base); + if (below == Number.Integer.One || above.ContainsNode(x) || !below.ContainsNode(x) + || !TreeAnalyzer.TryGetPolynomial(below, x, out var read) || read.Count == 0 + || read.Keys.Max()!.CompareTo(EInteger.FromInt32(4)) > 0 || IsPositiveForReal(below, x)) + continue; + var positive = above.Evaled is Number.Real { IsPositive: true }; + var inFront = positive ? MathS.Signum(below) : factor * MathS.Pow(below, exponent); + var replacement = positive ? MathS.Pow(above, exponent) * MathS.Pow(below, -exponent) : MathS.Pow(below, -exponent); + var written = expr.Replace(node => node == factor ? replacement : node); + return Integration.ComputeAsTheSameQuestion(written, x, integrateByParts) is { } answer ? inFront * answer : null; + } + return null; + } + + /// + /// x^(n - 1) g(x^n) with a symbolic n, which is g(u)/n under + /// u = x^n: the power in front is the derivative of the power inside up to the + /// constant n, and the rest mentions x only as x^n. Rubi's 6.5.2 and + /// 6.6.2 -- (e x)^(n - 1) (a + b sech(c + d x^n))^p -- write the power in front as + /// (e x)^(n - 1), which is e^(n - 1) x^(n - 1) for a positive e, and + /// the answer says so. A concrete n is the general substitution's; a symbolic + /// one was nobody's, and x^(n - 1) e^(x^n) was declined. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + /// x^(k n - 1) with k half an odd number is the same substitution by half + /// the power: under u = x^(n/2), x^(k n - 1) dx is (2/n) u^(2k - 1) du + /// and g(x^n) is g(u^2). Beside an exponential of x^n that is a moment + /// of the Gaussian: Rubi's 2.3, f^(a + b x^n) x^(-1 + 5n/2), and 6.1.3's + /// x^(-1 + n/2) sinh(a + b x^n). + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + internal static Entity? SolveByAPowerOfTheVariableTimesAFunctionOfItsPower(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity? inFront = null; + Entity? scale = null; + Entity? rest = null; + Entity constant = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + // The power in front: x^m, or (e x)^m with e free of x, to a symbolic m. + if (inFront is null && !underneath && factor is Powf(var @base, var exponent) && !exponent.ContainsNode(x) + && exponent.Evaled is not Number && exponent.Vars.Any() + && (@base == x || @base is Mulf(var l, var r) && (l == x && !r.ContainsNode(x) || r == x && !l.ContainsNode(x)))) + { + inFront = exponent; + scale = @base == x ? null : @base is Mulf(var l2, var r2) ? (l2 == x ? r2 : l2) : null; + continue; + } + var piece = underneath ? Number.Integer.One / factor : factor; + rest = rest is null ? piece : rest * piece; + } + if (inFront is null || rest is null) + return null; + // The rest mentions x only as x^n, one n throughout. + Entity? n = null; foreach (var node in rest.Nodes) if (node is Powf(var b, var e) && b == x && !e.ContainsNode(x)) { @@ -32442,6 +33423,362 @@ static bool IsAPowerOf(EInteger value, int order) } } + /// + /// A power of x that is not whole beside whole powers of two binomials in x^n or + /// more, x^q R(x^n) with R over (a + b x^n)^j (c + d x^n)^k, split into + /// partial fractions in w = x^n, with x^r for q = r + n m, 0 < r < n, + /// kept beside each term: Rubi's x^(9/2)/((a + b x^2) (c + d x^2)^2) is + /// sqrt(x) times w^2/((a + b w) (c + d w)^2) written apart, three binomial + /// differentials. + /// + /// + /// + /// The split is over the binomials as linears in w, so its coefficients are rational in the + /// symbols, and it is checked at sampled points before it is used; each term is a constant + /// times x^(r + n i)/(a + b x^n)^l, which the rules for one binomial answer in a tenth + /// of a second. These were declined past the budget: the polynomial term makes each + /// binomial monic, and the substitution u = sqrt(x) and the partial fractions after + /// it went on to take u^10/((a + b u^4) (c + d u^4)^2) apart over the fourth roots of + /// both. With a symbol only, at the question asked, before the polynomial term, and with one of + /// the binomials below the bar past its first power: over two first powers the routes as + /// written answer. A whole power of x is a rational function, and the rational integrator's. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAPowerOfXBesideBinomialsInAPowerOfItSplitInThatPower(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + var power = ERational.Zero; + var n = 0; + var binomials = new List<(Entity Base, Entity C, Entity D, int Power)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor == x || factor is Powf(var variable, Number.Rational) && variable == x) + { + var exponent = factor is Powf(_, Number.Rational { ERational: var e }) ? e : ERational.One; + power = underneath ? power.Subtract(exponent) : power.Add(exponent); + continue; + } + var (raised, whole) = factor is Powf(var @base, Number.Integer { EInteger: var w }) && w.CanFitInInt32() + ? (@base, w.ToInt32Checked()) : (factor, 1); + if (underneath) + whole = -whole; + if (!TryReadAsABinomialIn(raised, x, out var c, out var d, out var degree) || degree < 2 || n != 0 && degree != n + || Functions.PartialFractions.IsZeroAsAValue(c) || Functions.PartialFractions.IsZeroAsAValue(d)) + return null; + n = degree; + var index = binomials.FindIndex(binomial => binomial.Base == raised); + if (index >= 0) + binomials[index] = (raised, c, d, binomials[index].Power + whole); + else + binomials.Add((raised, c, d, whole)); + } + power = power.ToLowestTerms(); + // Two below the bar at least, one of them past its first power: over two first powers the + // routes as written answer in a second. + if (n == 0 || power.Denominator.Equals(EInteger.One) || binomials.Count(binomial => binomial.Power < 0) < 2 + || binomials.Where(binomial => binomial.Power < 0).Sum(binomial => -binomial.Power) < 3 + || binomials.Any(binomial => System.Math.Abs(binomial.Power) > MostBinomialPowersSplitInTheirPower)) + return null; + // q = r + n m, with m the floor of q/n. + var over = power.Denominator.Multiply(EInteger.FromInt32(n)); + var m = power.Numerator.Divide(over); + if (power.Numerator.Sign < 0 && !m.Multiply(over).Equals(power.Numerator)) + m -= EInteger.One; + if (!m.CanFitInInt32() || System.Math.Abs(m.ToInt32Checked()) > MostBinomialPowersSplitInTheirPower) + return null; + var r = power.Subtract(ERational.FromEInteger(m.Multiply(EInteger.FromInt32(n)))).ToLowestTerms(); + var u = Entity.Variable.CreateUnique(expr, "w_power"); + Entity above = m.Sign > 0 ? MathS.Pow(u, Number.Integer.Create(m)) : Number.Integer.One; + Entity? below = m.Sign < 0 ? MathS.Pow(u, Number.Integer.Create(m.Negate())) : null; + foreach (var (_, c, d, whole) in binomials) + { + var inU = c + d * u; + if (whole > 0) + above *= whole == 1 ? inU : MathS.Pow(inU, Number.Integer.Create(whole)); + else if (whole < 0) + { + var written = whole == -1 ? inU : MathS.Pow(inU, Number.Integer.Create(-whole)); + below = below is null ? written : below * written; + } + } + if (below is null || !Functions.PartialFractions.TrySplitOverWrittenFactors(above.Expand(), below, u, out var decomposition)) + return null; + if (decomposition is Divf(var split, var under) && !under.ContainsNode(u)) + (decomposition, constant) = (split, constant / under); + // Each term over its binomial, its numerator a polynomial in w written out term by term: + // `x^r w^i/(c + d w)^l` at `w = x^n`. + var root = MathS.Pow(x, Number.Rational.Create(r)); + var questions = new List(); + foreach (var term in Sumf.LinearChildren(decomposition)) + { + if (TreeAnalyzer.IsZero(term)) + continue; + var (top, bottom) = Functions.SingleQuotient.Of(term); + if (!TreeAnalyzer.TryGetPolynomial(top, u, out var read) || read.Keys.Any(degree => degree.Sign < 0 || !degree.CanFitInInt32()) + || read.Values.Any(coefficient => coefficient.ContainsNode(u))) + return null; + var inX = bottom.Substitute(u, MathS.Pow(x, Number.Integer.Create(n))); + foreach (var (degree, coefficient) in (read).Select(pair => (pair.Key, pair.Value))) + { + if (Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + var exponent = r.Add(ERational.FromEInteger(degree.Multiply(EInteger.FromInt32(n)))); + questions.Add(coefficient * MathS.Pow(x, Number.Rational.Create(exponent)) / inX); + } + } + if (questions.Count < 2) + return null; + Entity sum = Number.Integer.Zero; + foreach (var question in questions) + { + if (Integration.ComputeAsAQuestionOfItsOwn(question, x, integrateByParts) is not { } integral + || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += integral; + } + var answer = constant * sum; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + + /// The largest power of a binomial, and of x over the degree, takes. + private const int MostBinomialPowersSplitInTheirPower = 8; + + /// + /// A polynomial of degree two or more beside two powers of linears whose exponents add up + /// to -1, P(x) (a + b x)^m (c + d x)^n with m + n = -1, by the reduction that + /// leaves a constant times the two powers alone: + /// int P L1^m L2^n dx = R L1^(m + 1) L2^(n + 1) + C int L1^m L2^n dx, R a polynomial + /// of degree one less than P. Rubi's (e + f x)^3/((a + b x)^(1/3) (c + d x)^(2/3)). + /// + /// + /// + /// The derivative of R L1^(m + 1) L2^(n + 1) is + /// L1^m L2^n (R' L1 L2 + R ((m + 1) b L2 + (n + 1) d L1)), which takes x^j to + /// j a c x^(j - 1) + (j (a d + b c) + (m + 1) b c + (n + 1) a d) x^j + (j + 1) b d x^(j + 1); + /// the last coefficient is never zero, so R is read off P's coefficients from the top down, + /// and what is left of P is the constant C. + /// + /// + /// Written out term by term, each power of x beside the two roots was answered by the + /// substitution for a quotient of two linear radicals in two to six seconds, and the cube + /// past the budget, where the two powers alone take a tenth of one. Only for a polynomial of + /// degree two or more, which is where that ran out of time, and with a symbol in it; the + /// answer is differentiated back at sampled points. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAPolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOne(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + Entity polynomial = Number.Integer.One; + var powers = new List<(Entity Linear, ERational Exponent)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor is Powf(var @base, Number.Rational exponent) && exponent is not Number.Integer + && TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out var intercept) + && !slope.ContainsNode(x) && !intercept.ContainsNode(x)) + { + var written = underneath ? exponent.ERational.Negate() : exponent.ERational; + var index = powers.FindIndex(pair => pair.Linear == @base); + if (index >= 0) + powers[index] = (@base, powers[index].Exponent.Add(written)); + else + powers.Add((@base, written)); + if (powers.Count > 2) + return null; + continue; + } + var raised = factor is Powf(var inner, Number.Integer { EInteger.Sign: > 0 }) ? inner : factor; + if (underneath || !TreeAnalyzer.TryGetPolynomial(raised, x, out var read) + || read.Keys.Any(degree => degree.Sign < 0) || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + polynomial *= factor; + } + if (powers.Count != 2 || powers[0].Exponent.Add(powers[1].Exponent).CompareTo(ERational.FromInt32(-1)) != 0 + || !TreeAnalyzer.TryGetPolynomial(polynomial, x, out var p) || p.Values.Any(coefficient => coefficient.ContainsNode(x)) + || p.Keys.Any(degree => degree.Sign < 0 || !degree.CanFitInInt32()) || p.Keys.Max() is not { } top + || top.CompareTo(EInteger.FromInt32(2)) < 0 || top.CompareTo(EInteger.FromInt32(MostPolynomialDegreeReducedBesideTwoLinearPowers)) > 0 + || !TreeAnalyzer.TryGetPolyLinear(powers[0].Linear, x, out var b, out var a) + || !TreeAnalyzer.TryGetPolyLinear(powers[1].Linear, x, out var d, out var c)) + return null; + var ((first, m), (second, n)) = (powers[0], powers[1]); + if (Functions.PartialFractions.IsZeroAsAValue(b) || Functions.PartialFractions.IsZeroAsAValue(d) + || Functions.PartialFractions.IsZeroAsAValue(a * d - b * c)) + return null; + var k = top.ToInt32Checked(); + Entity Coefficient(int i) => p.TryGetValue(EInteger.FromInt32(i), out var value) ? value : Number.Integer.Zero; + var mPlusOne = Number.Rational.Create(m.Add(ERational.One)); + var nPlusOne = Number.Rational.Create(n.Add(ERational.One)); + var alpha = mPlusOne * b * c + nPlusOne * a * d; + // R's coefficients r_0 .. r_(k - 1), from the top: the coefficient of x^i in the + // operator applied to R is i b d r_(i - 1) + (i (a d + b c) + alpha) r_i + (i + 1) a c r_(i + 1). + var rs = new Entity[k + 1]; + for (var i = 0; i <= k; i++) + rs[i] = Number.Integer.Zero; + for (var i = k; i >= 1; i--) + { + var left = Coefficient(i) - (i * (a * d + b * c) + alpha) * rs[i] - (i + 1) * a * c * (i + 1 <= k ? rs[i + 1] : Number.Integer.Zero); + rs[i - 1] = LowestOverTheSymbols(left / (i * b * d)); + } + var remainder = LowestOverTheSymbols(Coefficient(0) - alpha * rs[0] - a * c * rs[1]); + Entity reduced = Number.Integer.Zero; + for (var j = 0; j < k; j++) + if (!Functions.PartialFractions.IsZeroAsAValue(rs[j])) + reduced += rs[j] * (j == 0 ? Number.Integer.One : j == 1 ? x : MathS.Pow(x, Number.Integer.Create(j))); + Entity sum = reduced * MathS.Pow(first, mPlusOne) * MathS.Pow(second, nPlusOne); + if (!Functions.PartialFractions.IsZeroAsAValue(remainder)) + { + if (Integration.ComputeAsAQuestionOfItsOwn(MathS.Pow(first, Number.Rational.Create(m)) * MathS.Pow(second, Number.Rational.Create(n)), x, integrateByParts) is not { } integral + || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += remainder * integral; + } + var answer = constant * sum; + return answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? null : answer; + } + + /// The largest degree of the polynomial takes. + private const int MostPolynomialDegreeReducedBesideTwoLinearPowers = 12; + + /// + /// A rational function with symbols in it beside a half-odd power of a quadratic, + /// R(x) Q^(k/2), its polynomial part divided out and the rest split into partial + /// fractions over the factors below the bar as written, each term beside the power the + /// question of one factor: Rubi's (a + b x + c x^2)^(3/2)/((d + e x) (f + g x)^2) and + /// x^4/((a + b x + c x^2)^(3/2) (d - f x^2)). + /// + /// + /// The terms are x^i Q^(k/2) and x^i Q^(k/2)/F^j for a linear or a quadratic F, each + /// answered by the rules for one factor beside the root, where the whole was declined: the + /// rules for a rational function beside the root of a quadratic read one factor below the + /// bar. The split is checked at sampled points before it is used. Late, at the question + /// asked, and with a symbol, as + /// is for the root of a linear. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirst(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x)) + return null; + Entity constant = Number.Integer.One; + Entity? radicand = null; + var half = ERational.Zero; + Entity above = Number.Integer.One; + Entity? below = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor is Powf(var root, Number.Rational { ERational: var power }) && power.Denominator.Equals(EInteger.FromInt32(2))) + { + if (radicand is not null) + return null; + radicand = root; + half = underneath ? power.Negate() : power; + continue; + } + var raised = factor is Powf(var inner, Number.Integer { EInteger.Sign: > 0 }) ? inner : factor; + if (!TreeAnalyzer.TryGetPolynomial(raised, x, out var read) || read.Keys.Any(degree => degree.Sign < 0) + || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + if (underneath) + below = below is null ? factor : below * factor; + else + above *= factor; + } + if (radicand is null || below is null || !TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var leading, out _, out _) + || leading.ContainsNode(x) || Functions.PartialFractions.IsZeroAsAValue(leading) + || Mulf.LinearChildren(below).Any(factor => (factor is Powf(var b0, Number.Integer) ? b0 : factor) == radicand) + || !TreeAnalyzer.TryGetPolynomial(above, x, out var top) || !TreeAnalyzer.TryGetPolynomial(below, x, out var bottom) + || top.Concat(bottom).Any(pair => !pair.Key.CanFitInInt32() || pair.Value.ContainsNode(x))) + return null; + // The polynomial part, by long division with the coefficients in lowest terms. + var degreeBelow = bottom.Keys.Max()!.ToInt32Checked(); + var lead = bottom[EInteger.FromInt32(degreeBelow)]; + var remainder = top.ToDictionary(pair => pair.Key.ToInt32Checked(), pair => pair.Value); + var quotient = new Dictionary(); + for (var degree = remainder.Keys.DefaultIfEmpty(-1).Max(); degree >= degreeBelow; degree--) + { + if (!remainder.TryGetValue(degree, out var coefficient) || Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + var step = LowestOverTheSymbols(coefficient / lead); + quotient[degree - degreeBelow] = step; + foreach (var (k, value) in (bottom).Select(pair => (pair.Key, pair.Value))) + { + var at = k.ToInt32Checked() + degree - degreeBelow; + remainder[at] = LowestOverTheSymbols((remainder.TryGetValue(at, out var already) ? already : Number.Integer.Zero) - step * value); + } + remainder.Remove(degree); + } + var questions = new List(); + var halfOdd = MathS.Pow(radicand, Number.Rational.Create(half)); + foreach (var (k, coefficient) in (quotient).Select(pair => (pair.Key, pair.Value))) + if (!Functions.PartialFractions.IsZeroAsAValue(coefficient)) + questions.Add(coefficient * (k == 0 ? halfOdd : MathS.Pow(x, Number.Integer.Create(k)) * halfOdd)); + Entity left = Number.Integer.Zero; + foreach (var (k, coefficient) in (remainder).Select(pair => (pair.Key, pair.Value))) + if (!Functions.PartialFractions.IsZeroAsAValue(coefficient)) + left += k == 0 ? coefficient : coefficient * MathS.Pow(x, Number.Integer.Create(k)); + if (left != Number.Integer.Zero) + { + IReadOnlyList terms; + if (Mulf.LinearChildren(below).Count(factor => factor.ContainsNode(x)) >= 2) + { + if (!Functions.PartialFractions.TrySplitOverWrittenFactors(left, below, x, out var decomposition)) + return null; + if (decomposition is Divf(var split, var under) && !under.ContainsNode(x)) + (decomposition, constant) = (split, constant / under); + terms = Sumf.LinearChildren(decomposition); + } + else + terms = new[] { left / below }; + foreach (var term in terms) + { + if (TreeAnalyzer.IsZero(term)) + continue; + var (termAbove, termBelow) = Functions.SingleQuotient.Of(term); + if (!TreeAnalyzer.TryGetPolynomial(termAbove, x, out var read) || read.Values.Any(value => value.ContainsNode(x))) + return null; + foreach (var (k, coefficient) in (read).Select(pair => (pair.Key, pair.Value))) + if (!Functions.PartialFractions.IsZeroAsAValue(coefficient)) + questions.Add(coefficient * (k.IsZero ? halfOdd : MathS.Pow(x, Number.Integer.Create(k)) * halfOdd) / termBelow); + } + } + if (questions.Count < 2 || questions.Count > MostQuestionsBesideAHalfOddPowerOfAQuadratic) + return null; + Entity sum = Number.Integer.Zero; + foreach (var question in questions) + { + if (Integration.ComputeAsAQuestionOfItsOwn(question, x, integrateByParts) is not { } integral || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += integral; + } + var answer = constant * sum; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + + /// The most terms asks. + private const int MostQuestionsBesideAHalfOddPowerOfAQuadratic = 8; + /// /// A rational function over one block quadratic in a power of x, a + b x^n + c x^(2n) /// with a symbol in it, split at the block's roots in u = x^n: the single block diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 80893e834..1ff48ce49 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -741,6 +741,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // And over the cube root of a quadratic beside another about the same axis, in the // variable that takes the linear terms away: Rubi's 1.2.1.6 pseudo-elliptic 232 is 231 there. if ((answer = IndefiniteIntegralSolver.SolveOverTheCubeRootOfAQuadraticBesideAnotherAboutOneAxis(expr, x, integrateByParts)) is { }) return answer; + // And beside the cube roots of two linears whose product is a quadratic, with the two roots + // written back for the root of the product. + if ((answer = IndefiniteIntegralSolver.SolveALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadratic(expr, x)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveTwoLinearsToOnePowerAsOneQuadratic(expr, x, integrateByParts)) is { }) return answer; // A linear and a quadratic to one power as one binomial, `(1 + x) (1 - x + x^2)` as // `1 + x^3`, which the rules for a power of a binomial answer. @@ -810,6 +813,16 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // term hands the reciprocal of the product down a level, and the sign of the complement // wrote it there with sgn(sqrt(cos(y)) - 1), owing cos(y) >= 0. if ((answer = IndefiniteIntegralSolver.SolveAHalfOddPowerOfACosineBesideTheRootOfAMinusACosine(expr, x)) is { }) return answer; + // The same half-odd power of two linears beside a squared quadratic, as one power of their + // product: `1/((a + b x + c x^2)^2 sqrt(1 - d x) sqrt(1 + d x))`, before the polynomial term. + if ((answer = IndefiniteIntegralSolver.SolveTwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadratic(expr, x, integrateByParts)) is { }) return answer; + // A power below -1 of a cube root of a binomial beside a high power of x and another binomial, + // divided by the root here, where it is answered: a level down it ran past the budget. + if ((answer = IndefiniteIntegralSolver.SolveByDividingByAHighPowerOfARootBesideAnotherBinomialFirst(expr, x)) is { }) return answer; + // A power of x that is not whole beside binomials in x^n, split over them in w = x^n first, + // before the polynomial term makes them monic and the substitution under the root takes + // `u^10/((a + b u^4) (c + d u^4)^2)` apart over fourth roots: `x^(9/2)/((a + b x^2) (c + d x^2)^2)`. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXBesideBinomialsInAPowerOfItSplitInThatPower(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveAsPolynomialTerm(expr, x, integrateByParts)) is { }) return answer; // A product of powers of the variable with a power of a constant multiple of it // among them, `(c x)^m x^n`, by the power rule with the written power kept as @@ -1055,6 +1068,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // `(1 + x + sqrt(3))/((1 + x - sqrt(3)) sqrt(1 + x^3))`. Before the reduction over the // linear, which writes it as `1` and `2 sqrt(3)` over the linear, neither elementary. if ((answer = IndefiniteIntegralSolver.SolveAPseudoEllipticQuotientOverTheRootOfACubicBinomial(expr, x)) is { }) return answer; + // And a linear over a linear beside the root of an even quartic, where it is a function of + // the numerator's square over the root: Rubi's `(1 + x - sqrt(3))/((1 + x + sqrt(3)) sqrt(-4 + x^4 + 4 sqrt(3) x^2))`. + if ((answer = IndefiniteIntegralSolver.SolveALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquare(expr, x)) is { }) return answer; // x over a cubic binomial beside the root of another, at the two ratios where that is // elementary: Welz's `x/((4 - x^3) sqrt(1 - x^3))`. if ((answer = IndefiniteIntegralSolver.SolveXOverACubicBinomialBesideTheRootOfAnother(expr, x, integrateByParts)) is { }) return answer; @@ -1368,6 +1384,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // times the length of the one answer. if (IndefiniteIntegralSolver.HoldsAPowerOfASumWithARoot(expr, x) && (answer = IndefiniteIntegralSolver.SolveByEulerSubstitution(expr, x, onlyBySumItWrites: true)) is { }) return answer; + // A polynomial beside two powers of linears whose exponents add up to -1, reduced to the two + // powers alone, before the split writes the polynomial out: each term of + // `(e + f x)^3/((a + b x)^(1/3) (c + d x)^(2/3))` took seconds by the quotient of the roots. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOne(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveBySplittingSum(expr, x, integrateByParts)) is { }) return answer; // A sum whose terms all carry the same power of one linear form, as that power times // the sum of the rest: `x cosh(x)^(3/2) - x sqrt(cosh(x))/3` is `x` times a pair @@ -1456,6 +1476,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) if ((answer = IndefiniteIntegralSolver.SolveALinearOverAQuadraticBesideTheRootOfAnother(expr, x)) is { }) return answer; // And over a power of the quadratic, a power at a time down to that one. if ((answer = IndefiniteIntegralSolver.SolveAPolynomialOverAPowerOfAQuadraticBesideTheRootOfAnother(expr, x)) is { }) return answer; + // A rational function of a quadratic beside a half-odd power of it, under t = Q'/sqrt(Q): + // after the reduction over a power of the quadratic, which answers the third and fourth + // powers of `1/((d + b x + c x^2)^k sqrt(a + b x + c x^2))` in a third of the characters. + if ((answer = IndefiniteIntegralSolver.SolveAFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRoot(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionBesideTheRootOfAQuadratic(expr, x)) is { }) return answer; // A linear below the bar that divides the radicand, written over it: after the two // rules above, which read a linear beside the root as a pole and find nothing to take @@ -1473,6 +1497,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // beside the one with numbers: the rules since then answer what they share with it // more shortly, `1/(a - b x^4)^(1/4)` by two arctangents where this gives four terms. if ((answer = IndefiniteIntegralSolver.SolveABinomialDifferentialWithSymbols(expr, x)) is { }) return answer; + // An algebraic function whose symbols a scaling of x takes out, as a constant times one with + // numbers: `1/((c + d x) (2 c^3 + d^3 x^3)^(1/3))` in y = d x/c, which the ansatz below answers; + // and before the ansatz, which ran past the budget on `1/((1 + 2^(1/3) x) (1 + x^3)^(2/3))`. + if ((answer = IndefiniteIntegralSolver.SolveByScalingTheVariableToNumbers(expr, x, integrateByParts)) is { }) return answer; // A rational function of x and one cube root of a polynomial: no substitution // rationalises it, and the elementary ones are logarithms of `L - y` for linear L // whose cube agrees with the polynomial at the poles, found by an ansatz. @@ -1571,6 +1599,22 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // One block quadratic in a power of x, split at its roots, which are complex where its // discriminant is negative: as late, so that a rule answering the block whole goes first. if ((answer = IndefiniteIntegralSolver.SolveByOneBlockInAPowerOfXAtItsRoots(expr, x, integrateByParts)) is { }) return answer; + // A rational function with symbols in it beside a half-odd power of a quadratic, its polynomial + // part divided out and the rest split over the factors below the bar first. + if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirst(expr, x, integrateByParts)) is { }) return answer; + // A power of x that is not whole beside a polynomial in x^n over a whole power below -1 of a + // binomial in x^n, each term reduced to the first power: `1/(x^(3/2) (a + c x^4)^2)`. + if ((answer = IndefiniteIntegralSolver.SolveAFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOne(expr, x, integrateByParts)) is { }) return answer; + // And beside a trinomial in x^n: `1/(x^(3/2) (a + b x^2 + c x^4)^2)`. + if ((answer = IndefiniteIntegralSolver.SolveAFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOne(expr, x, integrateByParts)) is { }) return answer; + // A sum whose terms share a power that is not whole of one polynomial, as that power times + // the rest: `(1 - x^6)^(2/3) + (1 - x^6)^(2/3)/x^6`, whose terms apart are not elementary. + if ((answer = IndefiniteIntegralSolver.SolveASumOfTermsSharingAPowerOfAPolynomialAsOneProduct(expr, x, integrateByParts)) is { }) return answer; + // A rational function over the root of a palindromic quartic, under u = x + 1/x, the root of + // the quartic over x standing for the root of the quadratic in u. + if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocal(expr, x, integrateByParts)) is { }) return answer; + // A root of a product with a root in one of its factors, written out under the root. + if ((answer = IndefiniteIntegralSolver.SolveByWritingOutAProductUnderANestedRoot(expr, x, integrateByParts)) is { }) return answer; // A rational function with symbols in it beside a root of a linear, split into partial // fractions first, each term the question of one factor beside the root. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionBesideARootOfALinearSplitFirst(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOneIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOneIntegralTest.cs new file mode 100644 index 000000000..4ed48f67e --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOneIntegralTest.cs @@ -0,0 +1,50 @@ +// +// 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 that is not whole beside a polynomial in x^n over a whole power below -1 of a + /// binomial in x^n, by the reduction that raises the power of the binomial. Rubi's 1.1.3.2 and 1.1.3.4. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AFractionalPowerOfXBesideABinomialToAWholePowerBelowMinusOneIntegralTest + { + [Theory] + [InlineData("1/(x^(3/2)*(a + c*x^4)^2)")] + [InlineData("(A + B*x^3)/(x^(3/2)*(a + b*x^3)^3)")] + [InlineData("(A + B*x^3)/(x^(5/2)*(a + b*x^3)^3)")] + public void TheBinomialsPowerRaisedToMinusOne(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 6000, $"{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("A", 0.9).Substitute("B", 1.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOneIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOneIntegralTest.cs new file mode 100644 index 000000000..6707809d7 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOneIntegralTest.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 that is not whole beside a whole power below -1 of a trinomial in x^n, by the + /// reduction that raises the power of the trinomial. Rubi's 1.2.2.2 1/(x^(3/2) (a + b x^2 + c x^4)^2). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AFractionalPowerOfXBesideATrinomialToAWholePowerBelowMinusOneIntegralTest + { + [Theory] + [InlineData("1/(x^(3/2)*(a + b*x^2 + c*x^4)^2)")] + [InlineData("sqrt(x)/(a + b*x^2 + c*x^4)^2")] + public void TheTrinomialsPowerRaisedToMinusOne(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); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/AFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRootIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRootIntegralTest.cs new file mode 100644 index 000000000..c51b2a111 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRootIntegralTest.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 rational function of a quadratic beside a half-odd power of it, under t = Q'/sqrt(Q), where + /// it is rational. Rubi's 1.2.1.5 1/((d + b x + c x^2)^3 sqrt(a + b x + c x^2)). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AFunctionOfAQuadraticBesideItsRootUnderItsDerivativeOverTheRootIntegralTest + { + [Theory] + [InlineData("1/((d + b*x + c*x^2)^3*sqrt(a + b*x + c*x^2))")] + [InlineData("1/((d + b*x + c*x^2)^4*sqrt(a + b*x + c*x^2))")] + public void RationalUnderTheDerivativeOverTheRoot(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 16000, $"{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.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/ALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadraticIntegralTest.cs new file mode 100644 index 000000000..3d0bd08d2 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadraticIntegralTest.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 linear below the bar beside the cube roots of two linears whose product is a quadratic, + /// answered as the cube root of the quadratic with the two roots written back. Rubi's 1.1.1.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ALinearBesideCubeRootsOfTwoLinearsAsOneOfAQuadraticIntegralTest + { + [Theory] + [InlineData("1/((a + b*x)*(c + d*x)^(1/3)*(b*c + a*d + 2*b*d*x)^(1/3))")] + [InlineData("1/((d - 3*k*x)^(1/3)*(d + k*x)*(d + 3*k*x)^(1/3))")] + public void TheTwoRootsWrittenBack(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[] { -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/ALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquareIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquareIntegralTest.cs new file mode 100644 index 000000000..b299d7776 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquareIntegralTest.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 linear over a linear beside the root of an even quartic where its integral is a constant + /// times the inverse tangent, or its hyperbolic one, of the numerator's square over the root. Rubi's 1.2.2.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ALinearOverALinearBesideTheRootOfAnEvenQuarticAsAFunctionOfTheSquareIntegralTest + { + [Theory] + [InlineData("(1 + x - sqrt(3))/((1 + x + sqrt(3))*sqrt(-4 + x^4 + 4*x^2*sqrt(3)))")] + [InlineData("(1 + x + sqrt(3))/((1 + x - sqrt(3))*sqrt(-4 + x^4 - 4*x^2*sqrt(3)))")] + [InlineData("(1 + 2*x - sqrt(3))/((1 + 2*x + sqrt(3))*sqrt(-1 + 4*x^4 + 4*x^2*sqrt(3)))")] + [InlineData("(1 + 2*x + sqrt(3))/((1 + 2*x - sqrt(3))*sqrt(-1 + 4*x^4 - 4*x^2*sqrt(3)))")] + public void AnInverseTangentOfTheSquareOverTheRoot(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + 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[] { -2.1, -1.4, 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/APolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOneIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOneIntegralTest.cs new file mode 100644 index 000000000..a73933906 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOneIntegralTest.cs @@ -0,0 +1,50 @@ +// +// 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 beside two powers of linears whose exponents add up to -1, reduced to a constant + /// times the two powers alone. Rubi's 1.1.1.3 (e + f x)^3/((a + b x)^(1/3) (c + d x)^(2/3)). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APolynomialBesideTwoLinearPowersWhoseExponentsAddUpToMinusOneIntegralTest + { + [Theory] + [InlineData("(k + f*x)^3/((a + b*x)^(1/3)*(c + d*x)^(2/3))")] + [InlineData("(k + f*x)^2/((a + b*x)^(1/3)*(c + d*x)^(2/3))")] + [InlineData("x^3/((a + b*x)^(3/4)*(c + d*x)^(1/4))")] + public void ReducedToTheTwoPowersAlone(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).Substitute("f", 1.7); + 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/APowerOfXBesideBinomialsInAPowerOfItSplitInThatPowerIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfXBesideBinomialsInAPowerOfItSplitInThatPowerIntegralTest.cs new file mode 100644 index 000000000..171c199f0 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfXBesideBinomialsInAPowerOfItSplitInThatPowerIntegralTest.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 power of x that is not whole beside binomials in x^n, split into partial fractions in + /// w = x^n with the power of x kept beside each term. Rubi's 1.1.2.4 ran past the budget. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfXBesideBinomialsInAPowerOfItSplitInThatPowerIntegralTest + { + [Theory] + [InlineData("x^(9/2)/((a + b*x^2)*(c + d*x^2)^2)")] + [InlineData("x^(5/2)/((a + b*x^2)^2*(c + d*x^2)^2)")] + [InlineData("x^(7/2)/((a + b*x^2)^2*(c + d*x^2)^3)")] + [InlineData("sqrt(x)/((a + b*x^2)^2*(c + d*x^2)^3)")] + [InlineData("1/(x^(3/2)*(a + b*x^2)^2*(c + d*x^2)^2)")] + public void SplitOverTheBinomialsInXToTheN(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 8000, $"{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[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirstIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirstIntegralTest.cs new file mode 100644 index 000000000..c39d6f1fa --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirstIntegralTest.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 rational function beside a half-odd power of a quadratic, its polynomial part divided out + /// and the rest split over the factors below the bar. Rubi's 1.2.1.4 and 1.2.1.6. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ARationalFunctionBesideAHalfOddPowerOfAQuadraticSplitFirstIntegralTest + { + [Theory] + [InlineData("(a + b*x + c*x^2)^(3/2)/((d + k*x)*(f + g*x)^2)")] + [InlineData("x^4/((a + b*x + c*x^2)^(3/2)*(d - f*x^2))")] + public void SplitBesideTheRoot(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 35000, $"{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).Substitute("f", 1.7).Substitute("g", 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/ARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocalIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocalIntegralTest.cs new file mode 100644 index 000000000..869a4566f --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocalIntegralTest.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 rational function beside the root of a palindromic quartic, under u = x + 1/x, with the + /// root of the quartic over x for the root of the quadratic in u. Rubi's 1.3.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ARationalFunctionOverTheRootOfAPalindromicQuarticUnderXPlusItsReciprocalIntegralTest + { + [Theory] + [InlineData("(k*f - k*f*x^2)/((a*d + b*d*x + a*d*x^2)*sqrt(a + b*x + c*x^2 + b*x^3 + a*x^4))")] + [InlineData("(k*f - k*f*x^2)/((-a*d + b*d*x - a*d*x^2)*sqrt(-a + b*x + c*x^2 + b*x^3 - a*x^4))")] + public void UnderXPlusItsReciprocal(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("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.3).Substitute("k", 0.9).Substitute("f", 1.7); + 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/ASumOfTermsSharingAPowerOfAPolynomialAsOneProductIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ASumOfTermsSharingAPowerOfAPolynomialAsOneProductIntegralTest.cs new file mode 100644 index 000000000..506ccefdc --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumOfTermsSharingAPowerOfAPolynomialAsOneProductIntegralTest.cs @@ -0,0 +1,48 @@ +// +// 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 sum whose terms share a power that is not whole of one polynomial, written as that power + /// times the rest, where the terms apart are not elementary. Rubi's 1.3.2 (1 - x^6)^(2/3) + (1 - x^6)^(2/3)/x^6. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ASumOfTermsSharingAPowerOfAPolynomialAsOneProductIntegralTest + { + [Theory] + [InlineData("(1 - x^6)^(2/3) + (1 - x^6)^(2/3)/x^6")] + public void TheSharedPowerTakenOut(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 500, $"{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.6, -0.3, 0.3, 0.6, 0.8 }) + { + 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/ByDividingByAHighPowerOfARootBesideAnotherBinomialFirstIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ByDividingByAHighPowerOfARootBesideAnotherBinomialFirstIntegralTest.cs new file mode 100644 index 000000000..b61370ca2 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ByDividingByAHighPowerOfARootBesideAnotherBinomialFirstIntegralTest.cs @@ -0,0 +1,50 @@ +// +// 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 below -1 of a cube root of a binomial beside a high power of x and a second binomial, + /// divided by the root at the question asked. Rubi's 1.1.3.4 1/(x^9 (a + b x^3)^(4/3) (c + d x^3)). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ByDividingByAHighPowerOfARootBesideAnotherBinomialFirstIntegralTest + { + [Theory] + [InlineData("1/(x^9*(a + b*x^3)^(4/3)*(c + d*x^3))")] + [InlineData("1/(x^6*(a + b*x^3)^(4/3)*(c + d*x^3))")] + public void DividedByTheRootAtTheQuestionAsked(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 8000, $"{text.Length} characters of answer for {integrand}"); + // Exact pins: with decimals the answer's coefficients, up to c^45, came out NaN in the check. + Entity Pinned(Entity e) => e.Substitute("a", "13/10".ToEntity()).Substitute("b", "7/10".ToEntity()).Substitute("c", "2/5".ToEntity()).Substitute("d", "3/10".ToEntity()); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/ByScalingTheVariableToNumbersIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ByScalingTheVariableToNumbersIntegralTest.cs new file mode 100644 index 000000000..87588d544 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ByScalingTheVariableToNumbersIntegralTest.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 +{ + /// + /// An algebraic function whose symbols a scaling of the variable takes out, integrated as a + /// constant times one with numbers for coefficients. Rubi's 1.3.2 and 1.1.3.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ByScalingTheVariableToNumbersIntegralTest + { + [Theory] + [InlineData("1/((c + d*x)*(-c^3 + d^3*x^3)^(1/3))")] + [InlineData("1/((c + d*x)*(2*c^3 + d^3*x^3)^(1/3))")] + [InlineData("1/((c + d*x)*(2*c^3 + d^3*x^3)^(2/3))")] + [InlineData("(k + f*x)/((c + d*x)*(-c^3 + d^3*x^3)^(1/3))")] + [InlineData("(a + b*x^3)^(1/3)/(a - b*x^3)")] + [InlineData("1/((1 + 2^(1/3)*x)*(1 + x^3)^(2/3))")] + public void AConstantTimesAnIntegralWithNumbers(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).Substitute("f", 1.7); + 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/ByWritingOutAProductUnderANestedRootIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ByWritingOutAProductUnderANestedRootIntegralTest.cs new file mode 100644 index 000000000..1ee94b8d7 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ByWritingOutAProductUnderANestedRootIntegralTest.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 root of a product of a polynomial and a sum with a root in it, the product written out under + /// the root. Rubi's 1.3.2 sqrt(x (a x + b sqrt(-a/b^2 + a^2 x^2/b^2)))/(x sqrt(-a/b^2 + a^2 x^2/b^2)). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ByWritingOutAProductUnderANestedRootIntegralTest + { + [Theory] + [InlineData("sqrt(x*(a*x + b*sqrt(-a/b^2 + a^2*x^2/b^2)))/(x*sqrt(-a/b^2 + a^2*x^2/b^2))")] + [InlineData("sqrt(x*(-a*x + b*sqrt(a/b^2 + a^2*x^2/b^2)))/(x*sqrt(a/b^2 + a^2*x^2/b^2))")] + public void TheProductWrittenOut(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs index 3a331febc..fe15f953d 100644 --- a/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs @@ -42,6 +42,17 @@ public void ByTheRecurrence(string integrand) public void BesideAMultipleOfTheSameLinear(string integrand) => DifferentiatesBack(integrand, new[] { -1.2, -0.6, 0.3, 0.8, 1.5, 2.4 }); + /// + /// Past the twelfth power, which the recurrence used to stop at: Rubi's 1.1.4.2 + /// 1/(x^4 (b x^(2/3) + a x)^(3/2)) is 3/(u^13 (b + a u)^(3/2)) under u = x^(1/3). + /// #233 + /// + [Theory] + [InlineData("1/(x^13*(c + d*x)^(3/2))")] + [InlineData("1/(x^4*(c*x^(2/3) + d*x)^(3/2))")] + public void PastTheTwelfthPower(string integrand) + => DifferentiatesBack(integrand, new[] { 0.3, 0.8, 1.3, 1.9, 2.6, 3.1 }); + private static void DifferentiatesBack(string integrand, double[] points) { var integral = integrand.ToEntity().Integrate("x"); diff --git a/Sources/Tests/UnitTests/Calculus/TwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadraticIntegralTest.cs new file mode 100644 index 000000000..ea4d25638 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadraticIntegralTest.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 +{ + /// + /// The same half-odd power of two linears beside a squared quadratic, as one power of their + /// product with the two roots written back. Rubi's 1.2.1.4 1/((a + b x + c x^2)^2 sqrt(1 - d x) sqrt(1 + d x)). + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class TwoEqualHalfOddPowersOfLinearsAsOneBesideAPowerOfAQuadraticIntegralTest + { + [Theory] + [InlineData("1/((a + b*x + c*x^2)^2*sqrt(1 - d*x)*sqrt(1 + d*x))")] + public void AsOnePowerOfTheProduct(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 60000, $"{text.Length} characters of answer for {integrand}"); + // By central differences of the answer: its derivative, taken symbolically with the + // numbers pinned, carries a condition that is false at these pins, though the answer and + // its derivative agree with the integrand to every digit there. + Entity Pinned(Entity e) => e.Substitute("a", "13/10".ToEntity()).Substitute("b", "7/10".ToEntity()).Substitute("c", "2/5".ToEntity()).Substitute("d", "3/10".ToEntity()); + var antiderivative = Pinned(integral.Substitute("C", 0)); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + const double h = 1e-6; + foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + var got = (antiderivative.Substitute("x", at + h).EvalNumerical() - antiderivative.Substitute("x", at - h).EvalNumerical()) / (2 * h); + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-6 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"the antiderivative of {integrand} changes at the rate {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}