diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 4620ca12e..afed7e112 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1535,6 +1535,31 @@ by the linear's coefficients, `(b d + 2c d x)^9`, which ran past a minute | `"(b + 2*c*x)^9/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | past a minute; 76,153 characters on the unreleased master | 607 characters | | `"(b*d + 2*c*d*x)^9/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | past a minute | 611 characters | +### A polynomial over a power of a symbolic quadratic beside the root of another is reduced in one solve, and split over the quadratic's roots + +**Answers where there were none, and shorter ones.** `(a + b x + c x^2)^(3/2)/(d + k x + f x^2)^3` was declined, +and so were the third and fourth powers of `d + b x + c x^2` below `sqrt(a + b x + c x^2)`, Rubi's 1.2.1.5. The +reduction took the power of the quadratic `A` below the bar down one at a time and carried each step's quotients into +the next as written, past its bound on their size at the second step; what it answered was large on the unreleased +master, `(a + b x + c x^2)^(3/2)/(d + k x + f x^2)^2` in 140,668 characters. Up to the fourth power, with symbols in +the coefficients, it is one linear solve over the polynomials in them now, each coefficient over one denominator in +lowest terms, that denominator written as powers of the discriminant of `A` and of the resultant of the two +quadratics. And `(g + h x)/(A sqrt(B))`, where the discriminant of `A` has a symbol in it, is split over the two roots +of `A`, each piece the logarithm Rubi writes as an artanh of a linear over the root of `B`, where it was two +arctangents with the root of the resultant in every coefficient. Of Rubi's 1.2.1.5, 3 answers are new and 11 shorter; +of its 1.2.1.6, 80 are shorter; every other answer of the two files, and of 1.2.1.4, is unchanged +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a + b*x + c*x^2)^(3/2)/(d + k*x + f*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 2,160 characters | +| `"1/((d + b*x + c*x^2)^3*sqrt(a + b*x + c*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 1,379 characters | +| `"1/((d + b*x + c*x^2)^4*sqrt(a + b*x + c*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 2,734 characters | +| `"(a + b*x + c*x^2)^(3/2)/(d + k*x + f*x^2)^2".ToEntity().Integrate("x")` | `integral(...)`; 140,668 characters on the unreleased master | 5,686 characters | +| `"1/((a + b*x + c*x^2)^(1/2)*(d + k*x + f*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)`; 38,378 characters on the unreleased master | 2,158 characters | +| `"1/((a + b*x + c*x^2)^(1/2)*(d + k*x + f*x^2))".ToEntity().Integrate("x")` | `integral(...)`; 4,982 characters on the unreleased master | 863 characters | +| `"(A + B*x)/((a*k + b*k*x + b*f*x^2)^2*sqrt(d + k*x + f*x^2))".ToEntity().Integrate("x")` | `integral(...)`; 87,464 characters on the unreleased master | 2,641 characters | + ### A quadratic's half-odd power beside a power of another quadratic is integrated **Answers where there were none.** `sqrt(3 - x + 2 x^2)/(2 + 3 x + 5 x^2)^2` was declined, and so were @@ -1551,7 +1576,7 @@ quadratic with numbers in it, although beside the root each piece is a shape the | `"1/((3 - x + 2*x^2)^(5/2)*(2 + 3*x + 5*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 1,380 characters | | `"1/((3 - x + 2*x^2)^(5/2)*(2 + 3*x + 5*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,629 characters | | `"(2 + 5*x + x^2)*sqrt(3 + 2*x + 5*x^2)/(1 + 4*x - 7*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 5,242 characters | -| `"(a + b*x + c*x^2)^(1/2)/(d + e*x + f*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 50,758 characters, six symbols | +| `"(a + b*x + c*x^2)^(1/2)/(d + e*x + f*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 1,218 characters, six symbols | ### A power of a perfect square with a symbolic leading coefficient is a power of its root diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index d45e1812e..105e047ad 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -1906,7 +1906,7 @@ private static bool TrySolveOverRationals(Entity[][] matrix, Entity[] rhs, [NotN /// d the last pivot -- one fraction per unknown, no back-substitution to /// compound them. /// - private static bool TrySolveOverPolynomials(Entity[][] matrix, Entity[] rhs, [NotNullWhen(true)] out Entity[]? values) + internal static bool TrySolveOverPolynomials(Entity[][] matrix, Entity[] rhs, [NotNullWhen(true)] out Entity[]? values) { values = null; var symbols = matrix.SelectMany(row => row).Concat(rhs).SelectMany(entry => entry.Vars).Distinct() diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 8dec90348..bcb89f856 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -20811,6 +20811,16 @@ Dictionary PowerOfQ(int j) var qSquared = (p * p - e * f).Expand().InnerSimplified; if (VanishesIdentically(qSquared)) return null; + // Over the roots of A where they are symbols, as Rubi writes it: a few hundred + // characters where the form below, with q in every coefficient, is thousands. + if ((!VanishesIdentically(e) || !VanishesIdentically(f)) + && OverTheRootsOfTheQuadraticBelow(g, h, a2, a1, a0, b2, b1, b0, p, e, radicand, x, expr) is { } overTheRoots) + { + var divided = (overTheRoots / constant).InnerSimplified; + if (!divided.Nodes.Any(node => node is Number.Complex { IsNaN: true }) + && Functions.PartialFractions.DerivativeHoldsAtSampledPoints(divided, expr, x)) + return divided; + } var q = MathS.Sqrt(qSquared).InnerSimplified; var roots = !VanishesIdentically(e) ? new[] { (Lambda: p + q, Mu: e), (Lambda: p - q, Mu: e) } : !VanishesIdentically(f) ? new[] { (Lambda: f, Mu: p + q), (Lambda: f, Mu: p - q) } @@ -20873,6 +20883,92 @@ Dictionary PowerOfQ(int j) return answer; } + /// + /// (g + h x)/(A sqrt(B)) over the two roots of A, where its discriminant is + /// not a number: with s = sqrt(A1^2 - 4 A0 A2) and r = (-A1 +- s)/(2 A2), + /// (g + h x)/A is ((g + h r1)/(x - r1) - (g + h r2)/(x - r2))/s, and each + /// 1/((x - r) sqrt(B)) is -ln((N + 2 sqrt(B(r)) sqrt(B))/(x - r))/sqrt(B(r)) + /// with N = 2 B0 + B1 r + (B1 + 2 B2 r) x, Rubi's artanh of N/(2 sqrt(B(r)) sqrt(B)) + /// written with N once. Null where the discriminant is a number. + /// + /// + /// + /// With y = x - r, N is 2 B(r) + B'(r) y and y N' - N is + /// -2 B(r) - sqrt(B(r)) N/sqrt(B), so the derivative of the logarithm is + /// -sqrt(B(r))/(y sqrt(B)) on either branch of sqrt(B(r)), which appears + /// twice: the piece holds where B(r) is negative or off the real line, and where + /// r is. Where the roots are real, the argument of the logarithm is real and of + /// one sign, or off the real line, on the whole of each interval between the roots of + /// A and B, since (N + 2 sqrt(B(r)) sqrt(B)) (N - 2 sqrt(B(r)) sqrt(B)) + /// is (B1^2 - 4 B0 B2) y^2, zero only at the pole: continuous on each. Where they + /// are a conjugate pair that is not shown; at four pinnings with B positive + /// everywhere, F(b) - F(a) agreed with Simpson's rule on intervals up to + /// [-5, 5]. B(r) is written as the + /// remainder of B by A at r, (P + E r)/A2 with + /// P = A2 B0 - A0 B2 and E = A2 B1 - A1 B2, as Rubi writes it. + /// + /// + /// The roots, s, B(r) and a compound g or h are symbols of + /// their own while the answer is built, written back after. Over the linears for which + /// A is a sum of multiples of B and their squares, Rubi's 1.2.1.5 + /// 1/(sqrt(a + b x + c x^2) (d + e x + f x^2)) was 4,982 characters, the root of the + /// resultant of the two in every coefficient; it is 863. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + private static Entity? OverTheRootsOfTheQuadraticBelow(Entity g, Entity h, Entity a2, Entity a1, Entity a0, + Entity b2, Entity b1, Entity b0, Entity p, Entity e, Entity radicand, Entity.Variable x, Entity context) + { + var discriminant = (a1 * a1 - 4 * a0 * a2).Expand().InnerSimplified; + if (discriminant.Evaled is Number.Complex || !discriminant.Vars.Any()) + return null; + // Expanded, so that `c^2 b - b c c` is the zero it is. + p = p.Expand().InnerSimplified; + e = VanishesIdentically(e) ? Number.Integer.Zero : e.Expand().InnerSimplified; + var named = context; + Variable Fresh(string prefix) + { + var fresh = Variable.CreateUnique(named, prefix); + named = named + fresh; + return fresh; + } + var s = Fresh("s_roots"); + var written = new List<(Variable Name, Entity Value)>(); + Entity Named(Entity value, string prefix) + { + if (value is Number or Variable) + return value; + var name = Fresh(prefix); + written.Add((name, value)); + return name; + } + var (top, slope) = (Named(g, "g_roots"), Named(h, "h_roots")); + var root = MathS.Sqrt(radicand); + Entity? answer = null; + var roots = new List<(Variable Name, Entity Value)>(); + var atTheRoots = new List<(Variable Name, Entity Value)>(); + foreach (var (sign, side) in new[] { (-1, 1), (1, -1) }) + { + var r = Fresh("r_roots"); + var w = Fresh("w_roots"); + roots.Add((r, (side * s - a1) / (2 * a2))); + atTheRoots.Add((w, (p + e * r) / a2)); + // -ln((N + 2 sqrt(w) sqrt(B))/(x - r))/sqrt(w), over s with g + h r in front. + var linear = 2 * b0 + b1 * r + (b1 + 2 * b2 * r) * x; + var piece = sign * (top + slope * r) * MathS.Ln((linear + 2 * MathS.Sqrt(w) * root) / (x - r)) / (s * MathS.Sqrt(w)); + answer = answer is null ? piece : answer + piece; + } + answer = answer!.InnerSimplified; + foreach (var (name, value) in atTheRoots) + answer = answer.Substitute(name, value); + foreach (var (name, value) in roots) + answer = answer.Substitute(name, value); + answer = answer.Substitute(s, MathS.Sqrt(discriminant)); + foreach (var (name, value) in written) + answer = answer.Substitute(name, value); + return answer; + } + /// /// The largest the linear p + q x of one step of /// is let grow, @@ -20977,6 +21073,18 @@ static Entity WithoutConditions(Entity coefficient) for (var i = 0; i < coefficients.Length; i++) coefficients[i] = monomials.TryGetValue(EInteger.FromInt32(i), out var monomial) ? WithoutConditions(monomial) : Number.Integer.Zero; + // With symbols in it, the whole reduction as one solve over the polynomials in them, + // every coefficient in lowest terms; step by step, each step's quotients were carried + // through the next as written. + if (expr.Vars.Any(symbol => symbol != x) + && ReducedOverThePolynomialsInItsSymbols(coefficients, power, a2, a1, a0, b2, b1, b0, quadratic, radicand, constant, x, expr) is { } reduced) + { + var divided = reduced.InnerSimplified; + if (!divided.Nodes.Any(node => node is Number.Complex { IsNaN: true }) + && Functions.PartialFractions.DerivativeHoldsAtSampledPoints(divided, expr, x)) + return divided; + } + var root = MathS.Sqrt(radicand); Entity answer = Number.Integer.Zero; for (var j = power; j >= 1; j--) @@ -21069,6 +21177,282 @@ static Entity WithoutConditions(Entity coefficient) return answer; } + /// + /// The reduction of + /// as one linear system over the polynomials in the symbols: P/(A^n sqrt(B)) is + /// (Q sqrt(B)/A^(n - 1))' + (g + h x)/(A sqrt(B)) + sigma/sqrt(B) with Q of + /// degree 2n - 3, which is + /// 2 P = 2 Q' A B + Q B' A - 2 (n - 1) Q A' B + 2 (g + h x) A^(n - 1) + 2 sigma A^n, + /// one equation for each power of x up to 2n in as many unknowns. Null where a + /// coefficient is not a polynomial in at most seven symbols, or the solve declines. + /// + /// + /// + /// The decomposition is the one the reduction a power at a time makes, since it is + /// unique where A and B have no root in common; written so, each unknown + /// is one quotient of polynomials, where a step at a time carried each step's quotients + /// through the next as written. Solved by fraction-free elimination, the values are put + /// over their least common denominator in lowest terms, and that is written as the + /// powers of the discriminant of A, the resultant of A and B and the + /// symbols it is, Rubi's e^2 - 4 d f and + /// (c d - a f)^2 - (b d - a e) (c e - b f). A's content comes out first: over + /// a monic quadratic the integrand arrives over d f + e f x + f^2 x^2. The two + /// last pieces are asked with g, h and sigma as symbols of their own. + /// + /// + /// A step at a time, Rubi's 1.2.1.5 (a + b x + c x^2)^(3/2)/(d + e x + f x^2)^3 was + /// declined past , and + /// (a + b x + c x^2)^(3/2)/(d + e x + f x^2)^2 was 140,668 characters; they are 2,160 + /// and 5,686. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// + /// The highest power of the quadratic below the bar that + /// takes: with six symbols the fifth + /// ran past twenty seconds in the elimination, where the reduction a power at a time + /// declines it in half a second. + /// + private const int LargestPowerSolvedAtOnce = 4; + + private static Entity? ReducedOverThePolynomialsInItsSymbols(Entity[] coefficients, int power, + Entity a2, Entity a1, Entity a0, Entity b2, Entity b1, Entity b0, Entity quadratic, Entity radicand, Entity constant, Entity.Variable x, Entity expr) + { + // A named constant, e or pi, is a symbol of its own here and put back at the end: the + // polynomials are read in the symbols, and `d + e x + f x^2` was not one. + var standIns = new List<(Variable StandIn, Variable Constant)>(); + var spelled = expr; + foreach (var named in new[] { a2, a1, a0, b2, b1, b0, constant }.Concat(coefficients) + .SelectMany(coefficient => coefficient.VarsAndConsts).Where(symbol => symbol.IsConstant).Distinct().ToList()) + { + var standIn = Variable.CreateUnique(spelled, "c_named"); + spelled = spelled + standIn; + standIns.Add((standIn, named)); + } + Entity WithStandIns(Entity coefficient) + { + foreach (var (standIn, named) in standIns) + coefficient = coefficient.Substitute(named, standIn); + return coefficient; + } + if (standIns.Count > 0) + { + (a2, a1, a0, b2, b1, b0, constant) = (WithStandIns(a2), WithStandIns(a1), WithStandIns(a0), WithStandIns(b2), WithStandIns(b1), WithStandIns(b0), WithStandIns(constant)); + coefficients = coefficients.Select(WithStandIns).ToArray(); + } + var symbols = new[] { a2, a1, a0, b2, b1, b0 }.Concat(coefficients).SelectMany(coefficient => coefficient.Vars).Distinct() + .OrderBy(symbol => symbol.Name, System.StringComparer.Ordinal).ToList(); + if (symbols.Count == 0 || symbols.Contains(x) || symbols.Count >= Functions.MultivariatePolynomial.MaxVariables + || power > LargestPowerSolvedAtOnce) + return null; + var indices = new Dictionary(symbols.Count + 1); + for (var i = 0; i < symbols.Count; i++) + indices[symbols[i]] = i; + indices[x] = symbols.Count; + var everything = symbols.Append(x).ToList(); + var order = Enumerable.Range(0, everything.Count).ToList(); + Functions.MultivariatePolynomial? Read(Entity polynomial) => Functions.MultivariatePolynomial.TryParse(polynomial, indices); + Functions.MultivariatePolynomial? Gcd(Functions.MultivariatePolynomial? left, Functions.MultivariatePolynomial? right) + => left is null || right is null ? null : Functions.PolynomialGcd.Gcd(left, right, order, 0); + if (Read(a2) is not { } p2 || Read(a1) is not { } p1 || Read(a0) is not { } p0 + || Read(b2) is not { } q2 || Read(b1) is not { } q1 || Read(b0) is not { } q0) + return null; + // A's content out of it: over a monic quadratic the integrand comes as one over + // `d f + k f x + f^2 x^2`, whose roots would be written with sqrt(f^2 (k^2 - 4 d f)). + Entity content = Number.Integer.One; + if (Gcd(Gcd(p2, p1), p0) is { IsConstant: false } common + && p2.DivideExact(common) is { } reduced2 && p1.DivideExact(common) is { } reduced1 && p0.DivideExact(common) is { } reduced0) + { + (p2, p1, p0) = (reduced2, reduced1, reduced0); + (a2, a1, a0) = (p2.ToEntity(everything), p1.ToEntity(everything), p0.ToEntity(everything)); + quadratic = a2 * MathS.Sqr(x) + a1 * x + a0; + content = common.ToEntity(everything); + } + // A root of A that is one of B's, or a root of A twice, is not the generic case the + // solve is for: the system is singular, and a solution of it is not in lowest terms. + if (p1.Multiply(p1) is not { } squared || p0.Multiply(p2) is not { } product + || squared.Subtract(product.ScaleBy(ERational.FromInt32(4))) is not { IsZero: false } discriminant + || p2.Multiply(q0) is not { } r1 || p0.Multiply(q2) is not { } r2 || p2.Multiply(q1) is not { } r3 + || p1.Multiply(q2) is not { } r4 || p1.Multiply(q0) is not { } r5 || p0.Multiply(q1) is not { } r6 + || r1.Subtract(r2) is var first && first.Multiply(first) is not { } firstSquared + || r3.Subtract(r4).Multiply(r5.Subtract(r6)) is not { } cross + || firstSquared.Subtract(cross) is not { IsZero: false } resultant) + return null; + var a = a2 * MathS.Sqr(x) + a1 * x + a0; + var b = b2 * MathS.Sqr(x) + b1 * x + b0; + var aSlope = 2 * a2 * x + a1; + var bSlope = 2 * b2 * x + b1; + var size = 2 * power + 1; + var below = power - 1; + // The columns: Q's powers of x, then g, h and sigma. + var columns = new List(); + for (var i = 0; i < 2 * power - 2; i++) + { + Entity column = MathS.Pow(x, i) * (bSlope * a - 2 * below * aSlope * b); + if (i > 0) + column = column + 2 * i * MathS.Pow(x, i - 1) * a * b; + columns.Add(column); + } + columns.Add(2 * MathS.Pow(a, below)); + columns.Add(2 * x * MathS.Pow(a, below)); + columns.Add(2 * MathS.Pow(a, power)); + Entity[]? ByPowers(Entity polynomial) + { + if (Read(polynomial) is not { } parsed) + return null; + var byPowers = new Entity[size]; + for (var i = 0; i < size; i++) + byPowers[i] = Number.Integer.Zero; + foreach (var pair in parsed.CoefficientsIn(symbols.Count)) + { + if (pair.Key >= size) + return null; + byPowers[pair.Key] = pair.Value.ToEntity(everything); + } + return byPowers; + } + var matrix = new Entity[size][]; + for (var row = 0; row < size; row++) + matrix[row] = new Entity[size]; + for (var column = 0; column < size; column++) + { + if (ByPowers(columns[column]) is not { } entries) + return null; + for (var row = 0; row < size; row++) + matrix[row][column] = entries[row]; + } + Entity numerator = Number.Integer.Zero; + for (var i = 0; i < coefficients.Length; i++) + numerator = numerator + 2 * coefficients[i] * MathS.Pow(x, i); + if (ByPowers(numerator) is not { } rhs + || !Functions.PartialFractions.TrySolveOverPolynomials(matrix, rhs, out var values)) + return null; + + // Every value over one denominator, the least common multiple of theirs, each in + // lowest terms first: the solve cancels a small quotient only, and its last pivot + // carries powers of the factors below that the values do not. + var tops = new Functions.MultivariatePolynomial[values.Length]; + var bottom = Functions.MultivariatePolynomial.One(everything.Count); + var bottoms = new Functions.MultivariatePolynomial[values.Length]; + for (var i = 0; i < values.Length; i++) + { + var (written, writtenBelow) = Functions.SingleQuotient.Of(values[i]); + if (Read(written) is not { } top || Read(writtenBelow) is not { IsZero: false } over) + return null; + if (!top.IsZero && Gcd(top, over) is { IsConstant: false } shared + && top.DivideExact(shared) is { } lowerTop && over.DivideExact(shared) is { } lowerOver) + (top, over) = (lowerTop, lowerOver); + (tops[i], bottoms[i]) = (top, over); + if (top.IsZero) + continue; + if (Gcd(bottom, over) is not { } both || over.DivideExact(both) is not { } extra || bottom.Multiply(extra) is not { } wider) + return null; + bottom = wider; + } + var numerators = new Functions.MultivariatePolynomial[values.Length]; + for (var i = 0; i < values.Length; i++) + if (tops[i].IsZero) + numerators[i] = tops[i]; + else if (bottom.DivideExact(bottoms[i]) is not { } scale || tops[i].Multiply(scale) is not { } raised) + return null; + else + numerators[i] = raised; + // The constant below the bar and A's content to its power go below too, and then out + // goes what every numerator shares with it: a power of A's leading coefficient, which + // the integrand over a monic quadratic has above as well. + if (Read(constant) is not { IsZero: false } constantRead || constantRead.Multiply(bottom) is not { } withTheConstant) + return null; + bottom = withTheConstant; + if (content != Number.Integer.One) + { + if (Read(content) is not { } contentRead || contentRead.Power(power) is not { } contentPower || bottom.Multiply(contentPower) is not { } withTheContent) + return null; + bottom = withTheContent; + } + var sharedByAll = bottom; + foreach (var value in numerators) + if (!value.IsZero && Gcd(sharedByAll, value) is { } narrowed) + sharedByAll = narrowed; + if (!sharedByAll.IsConstant && bottom.DivideExact(sharedByAll) is { } lowerBottom) + { + var lowered = new Functions.MultivariatePolynomial[numerators.Length]; + for (var i = 0; i < numerators.Length; i++) + if (numerators[i].IsZero) + lowered[i] = numerators[i]; + else if (numerators[i].DivideExact(sharedByAll) is { } lowerValue) + lowered[i] = lowerValue; + else + return null; + (numerators, bottom) = (lowered, lowerBottom); + } + + // The denominator as a product of what the reduction divides by: the discriminant + // of A, the resultant of A and B, and the symbols, Rubi's `(e^2 - 4 d f)` and + // `(c d - a f)^2 - (b d - a e) (c e - b f)`. + var candidates = new List { discriminant, resultant }; + for (var i = 0; i < symbols.Count; i++) + candidates.Add(Functions.MultivariatePolynomial.Monomial(everything.Count, i)); + Entity denominator = Number.Integer.One; + var rest = bottom; + foreach (var candidate in candidates) + { + if (candidate.IsConstant) + continue; + var normalized = candidate.Normalized(); + var multiplicity = 0; + while (rest.DivideExact(normalized) is { } divided) + { + rest = divided; + multiplicity++; + } + if (multiplicity > 0) + denominator = denominator * (multiplicity == 1 ? normalized.ToEntity(everything) : MathS.Pow(normalized.ToEntity(everything), multiplicity)); + } + denominator = rest.ToEntity(everything) * denominator; + + static bool IsZero(Functions.MultivariatePolynomial value) => value.IsZero; + var root = MathS.Sqrt(radicand); + Entity? polynomial = null; + for (var i = 0; i < 2 * power - 2; i++) + if (!IsZero(numerators[i])) + { + var coefficient = numerators[i].ToEntity(everything); + var term = i == 0 ? coefficient : coefficient * MathS.Pow(x, i); + polynomial = polynomial is null ? term : polynomial + term; + } + var (g, h, sigma) = (numerators[2 * power - 2], numerators[2 * power - 1], numerators[2 * power]); + Entity answer = polynomial is null ? Number.Integer.Zero : polynomial * root / MathS.Pow(quadratic, below); + Variable Fresh(string prefix) + { + var fresh = Variable.CreateUnique(spelled, prefix); + spelled = spelled + fresh; + return fresh; + } + if (!IsZero(g) || !IsZero(h)) + { + var (top, slope) = (IsZero(g) ? (Entity)Number.Integer.Zero : Fresh("g_reduced"), IsZero(h) ? (Entity)Number.Integer.Zero : Fresh("h_reduced")); + if (SolveALinearOverAQuadraticBesideTheRootOfAnother((top + slope * x) / (quadratic * root), x) is not { } closing) + return null; + if (top is Variable topName) + closing = closing.Substitute(topName, g.ToEntity(everything)); + if (slope is Variable slopeName) + closing = closing.Substitute(slopeName, h.ToEntity(everything)); + answer = answer + closing; + } + if (!IsZero(sigma)) + { + var over = Fresh("k_reduced"); + if (Integration.ComputeIndefiniteIntegral(over / root, x, integrateByParts: false) is not { } overTheRoot + || overTheRoot.Nodes.Any(node => node is Integralf)) + return null; + answer = answer + overTheRoot.Substitute(over, sigma.ToEntity(everything)); + } + answer = answer / denominator; + foreach (var (standIn, named) in standIns) + answer = answer.Substitute(standIn, named); + return answer; + } + /// /// A polynomial over a power of a linear beside the square root of a quadratic, /// P/((x - p)^k sqrt(Q)) with P of degree below k, by the reciprocal diff --git a/Sources/Tests/UnitTests/Calculus/AQuadraticOverItsRootsBesideTheRootOfAnotherIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AQuadraticOverItsRootsBesideTheRootOfAnotherIntegralTest.cs new file mode 100644 index 000000000..096947673 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AQuadraticOverItsRootsBesideTheRootOfAnotherIntegralTest.cs @@ -0,0 +1,91 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A polynomial over a power of a symbolic quadratic beside the root of another, reduced by + /// one solve over the polynomials in the symbols, and what is left over the quadratic split + /// over its two roots, Rubi's writing. Rubi's 1.2.1.5: + /// (a + b x + c x^2)^(3/2)/(d + e x + f x^2)^3 was declined, and + /// (a + b x + c x^2)^(3/2)/(d + e x + f x^2)^2 was 140,668 characters. + /// #233 + /// + /// + /// Differentiated back at points on both sides of the quadratic's roots where they are real, + /// with a root where the radicand is negative too, and where they are a conjugate pair: + /// the pieces over the roots hold on either branch of the root of the radicand there. + /// + [Trait("Area", "Calculus")] + public sealed class AQuadraticOverItsRootsBesideTheRootOfAnotherIntegralTest + { + // d + k x + f x^2 and d + b x + c x^2 with the roots -0.22 and -2.28, beside a radicand + // positive everywhere. + private static readonly (string, double)[] RealRoots = + { ("a", 4), ("b", 5), ("c", 2), ("d", 1), ("k", 5), ("f", 2), ("g", 0.7), ("h", 1.3) }; + private static readonly double[] AroundTheRealRoots = { -3, -2.5, -1.5, -1, -0.1, 0.5, 1.5 }; + + // The same roots of d + k x + f x^2, the radicand negative at -0.22 and positive at -2.28. + private static readonly (string, double)[] ARootWhereTheRadicandIsNegative = + { ("a", -1), ("b", 0.5), ("c", 1), ("d", 1), ("k", 5), ("f", 2), ("g", 0.7), ("h", 1.3) }; + private static readonly double[] WhereTheRadicandIsPositive = { -3, -2.5, -2, -1.5, 1, 1.5, 2.5 }; + + // Both quadratics below the bar with a conjugate pair of roots. + private static readonly (string, double)[] ComplexRoots = + { ("a", 1), ("b", 1), ("c", 1), ("d", 2), ("k", 1), ("f", 1.5), ("g", 0.7), ("h", 1.3) }; + private static readonly double[] AlongTheLine = { -2, -1, 0, 0.5, 1.5 }; + + private static void DifferentiatesBack(string integrand, Entity integral, double[] points, (string, double)[] pins) + { + var derivative = integral.Substitute("C", 0).Differentiate("x"); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + derivative = derivative.Substitute(name, value); + original = original.Substitute(name, value); + } + var compared = 0; + foreach (var at in points) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)); + Assert.True(difference / scale < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 4, $"only {compared} points could be compared for {integrand}"); + } + + [Theory] + [InlineData("1/((a + b*x + c*x^2)^(1/2)*(d + k*x + f*x^2))", 1000)] + [InlineData("(g + h*x)/((d + k*x + f*x^2)*sqrt(a + b*x + c*x^2))", 2500)] + [InlineData("1/((d + b*x + c*x^2)*sqrt(a + b*x + c*x^2))", 1000)] + [InlineData("1/((a + b*x + c*x^2)^(1/2)*(d + k*x + f*x^2)^2)", 2500)] + [InlineData("1/((d + b*x + c*x^2)^3*sqrt(a + b*x + c*x^2))", 2000)] + [InlineData("(a + b*x + c*x^2)^(1/2)/(d + k*x + f*x^2)^2", 1500)] + [InlineData("(a + b*x + c*x^2)^(3/2)/(d + k*x + f*x^2)^2", 7000)] + [InlineData("(a + b*x + c*x^2)^(3/2)/(d + k*x + f*x^2)^3", 2500)] + public void OverTheRootsOfTheQuadratic(string integrand, int longest) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < longest, $"{text.Length} characters of answer for {integrand}"); + DifferentiatesBack(integrand, integral, AroundTheRealRoots, RealRoots); + DifferentiatesBack(integrand, integral, WhereTheRadicandIsPositive, ARootWhereTheRadicandIsNegative); + DifferentiatesBack(integrand, integral, AlongTheLine, ComplexRoots); + } + } +}