diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a1f9d34ed..61d0ed849 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -42,6 +42,19 @@ A not-equal written `<>` is read as `not (a = b)`, and printed that way. ## Unreleased — since 2.5.0 +### A high half-odd power of a quadratic beside a power of another is reduced in exact rationals + +**Answers where there were none.** Welz's `1/((3 - 2 x + x^2)^(11/2) (1 + x + 2 x^2)^5)` and its `^(21/2)`, `^10` twin +were declined. The pieces over the first quadratic are a polynomial times its half-odd power, reduced to +`R Q^(k + 1/2) + K/sqrt(Q)` by one linear solve, which took a numerator of degree 8 and a power of 9/2 at most. With +numbers for every coefficient the solve is now done in exact rationals, to degree 24 and a power of 25/2 +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/((3-2*x+x^2)^(11/2)*(1+x+2*x^2)^5)".ToEntity().Integrate("x")` | `integral(...)` | 1,853 characters | +| `"1/((3-2*x+x^2)^(21/2)*(1+x+2*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 3,475 characters | + ### `|` is divisibility, and was disjunction **Read this one first.** It is the only entry here that changes what existing input *means* rather diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 176887f8a..ded68afa2 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -17910,6 +17910,60 @@ static Dictionary Scaled(Dictionary poly, En return answer; } + /// + /// The solution of a linear system whose entries are all rational numbers, by Gauss-Jordan + /// elimination in exact rationals kept in lowest terms -- PeterO's rationals are not reduced + /// on their own, and unreduced the solve for Welz's piece took two minutes -- an unknown no + /// row pivots on set to zero; null where an + /// entry is not a rational number or the system is inconsistent. + /// + private static Entity[]? TrySolveInRationals(Entity[][] matrix, Entity[] rhs) + { + static ERational? Read(Entity entry) => entry is Number.Rational r ? r.ERational : entry.Evaled is Number.Rational e ? e.ERational : null; + var rows = matrix.Length; + var columns = rows == 0 ? 0 : matrix[0].Length; + var m = new ERational[rows][]; + for (var i = 0; i < rows; i++) + { + m[i] = new ERational[columns + 1]; + for (var j = 0; j < columns; j++) + if (Read(matrix[i][j]) is { } value) m[i][j] = value; else return null; + if (Read(rhs[i]) is { } right) m[i][columns] = right; else return null; + } + var pivotOf = new int[columns]; + for (var j = 0; j < columns; j++) pivotOf[j] = -1; + var row = 0; + for (var column = 0; column < columns && row < rows; column++) + { + var pivot = -1; + for (var i = row; i < rows; i++) + if (!m[i][column].IsZero) { pivot = i; break; } + if (pivot < 0) + continue; + (m[row], m[pivot]) = (m[pivot], m[row]); + var lead = m[row][column]; + for (var j = column; j <= columns; j++) + m[row][j] = m[row][j].Divide(lead).ToLowestTerms(); + for (var i = 0; i < rows; i++) + { + if (i == row || m[i][column].IsZero) + continue; + var factor = m[i][column]; + for (var j = column; j <= columns; j++) + m[i][j] = m[i][j].Subtract(factor.Multiply(m[row][j])).ToLowestTerms(); + } + pivotOf[column] = row; + row++; + } + for (var i = row; i < rows; i++) + if (!m[i][columns].IsZero) + return null; + var values = new Entity[columns]; + for (var j = 0; j < columns; j++) + values[j] = pivotOf[j] < 0 ? Number.Integer.Zero : Number.Rational.Create(m[pivotOf[j]][columns]); + return values; + } + /// /// A polynomial times an odd half power of a quadratic, P Q^(m/2), reduced to /// R Q^(k + 1/2) + K/sqrt(Q) by one linear solve, with K/sqrt(Q) the @@ -17968,7 +18022,13 @@ static Dictionary Scaled(Dictionary poly, En || pRead.Values.Any(coefficient => coefficient.ContainsNode(x))) return null; var degreeOfP = pRead.Count == 0 ? 0 : pRead.Keys.Max()!.ToInt32Unchecked(); - if (degreeOfP > 8 || m > 9 || m < -9) + // The solve is one linear system of about degreeOfP + |m| unknowns; with numbers for every + // coefficient it stays small to far higher powers, and Welz's + // `1/((3 - 2 x + x^2)^(11/2) (1 + x + 2 x^2)^5)` left a piece of degree 9 over the + // eleventh half power, and its `^(21/2)`, `^10` twin one of degree 19. + var numeric = !radicand.Vars.Any(v => v != x) && !polynomial.Vars.Any(v => v != x); + var (degreeCap, powerCap) = numeric ? (24, 25) : (8, 9); + if (degreeOfP > degreeCap || m > powerCap || m < -powerCap) return null; // Each coefficient that is neither a number nor a symbol is named for the solve and @@ -17981,6 +18041,10 @@ Entity Named(Entity value) { if (value is Number or Variable) return value; + // A number written as an expression, `-2` as `-1 * 2`, is that number: named, it kept + // the system out of the exact rational solve below. + if (!value.Vars.Any() && value.Evaled is Number.Rational number) + return number; foreach (var (name, named) in names) if (named == value) return name; @@ -18034,7 +18098,10 @@ Dictionary PowerOfQ(int j) matrix[row][column] = columns[column].TryGetValue(power, out var entry) ? entry : Number.Integer.Zero; rhs[row] = target.TryGetValue(power, out var wanted) ? wanted : Number.Integer.Zero; } - if (!Functions.PartialFractions.TrySolveLinear(matrix, rhs, out var values) || values is null) + // With numbers for every entry the system is solved in exact rationals: eliminated as + // expressions, the values of Welz's degree-9 piece over the eleventh half power came to + // half a million nodes of unfolded fractions and were refused at the bound. + if ((TrySolveInRationals(matrix, rhs) ?? (Functions.PartialFractions.TrySolveLinear(matrix, rhs, out var solved) ? solved : null)) is not { } values) return null; Entity r = Number.Integer.Zero; for (var k = 0; k <= degreeOfR; k++) diff --git a/Sources/Tests/UnitTests/Calculus/AHighHalfPowerOfAQuadraticBesideAPowerOfAnotherIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AHighHalfPowerOfAQuadraticBesideAPowerOfAnotherIntegralTest.cs new file mode 100644 index 000000000..812f12522 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AHighHalfPowerOfAQuadraticBesideAPowerOfAnotherIntegralTest.cs @@ -0,0 +1,56 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A high half-odd power of a quadratic beside a power of another: the pieces over the first are a + /// polynomial times its half-odd power, reduced by one linear solve, which with numbers for + /// coefficients is solved in exact rationals. Welz's 1/((3 - 2 x + x^2)^(11/2) (1 + x + 2 x^2)^5). + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class AHighHalfPowerOfAQuadraticBesideAPowerOfAnotherIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + + var pinned = integral.Substitute("C", 0); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + pinned = pinned.Substitute(name, value); + original = original.Substitute(name, value); + } + var derivative = pinned.Differentiate("x"); + foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 }) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + 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-8, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + // d with the sign that keeps the argument real at every sampled x, and |d| below 2/2.3^2 + // so that the inverse sine and cosine stay inside [-1, 1]; off it the exact checks hold too. + [Theory] + [InlineData("1/((3-2*x+x^2)^(11/2)*(1+x+2*x^2)^5)")] + [InlineData("1/((3-2*x+x^2)^(21/2)*(1+x+2*x^2)^10)")] + [InlineData("(1+x^9)*(3-2*x+x^2)^(-11/2)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + => DifferentiatesBack(integrand); + } +}