Skip to content

Commit 84dc0e7

Browse files
Reduce a polynomial times a high half-odd power of a quadratic with numbers for coefficients (#1913)
* Reduce a polynomial times a high half-odd power of a quadratic with numbers for coefficients The reduction to R Q^(k + 1/2) + K/sqrt(Q) took a numerator of degree 8 and a power of 9/2 at most, and solved its system as expressions. With numbers for every coefficient the system is now solved in exact rationals, kept in lowest terms, to degree 24 and a power of 25/2, and a coefficient written as an expression that evaluates to a rational is that rational rather than a name. 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 were refused at the bound, and solved as expressions their values came to half a million nodes. Part of #718 Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the high half-odd powers in the breaking changes Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --------- Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent 1e79b0e commit 84dc0e7

3 files changed

Lines changed: 138 additions & 2 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,19 @@ A not-equal written `<>` is read as `not (a = b)`, and printed that way.
4242

4343
## Unreleased — since 2.5.0
4444

45+
### A high half-odd power of a quadratic beside a power of another is reduced in exact rationals
46+
47+
**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
48+
were declined. The pieces over the first quadratic are a polynomial times its half-odd power, reduced to
49+
`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
50+
numbers for every coefficient the solve is now done in exact rationals, to degree 24 and a power of 25/2
51+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
52+
53+
| Input | Was (2.5.0) | Now |
54+
|---|---|---|
55+
| `"1/((3-2*x+x^2)^(11/2)*(1+x+2*x^2)^5)".ToEntity().Integrate("x")` | `integral(...)` | 1,853 characters |
56+
| `"1/((3-2*x+x^2)^(21/2)*(1+x+2*x^2)^10)".ToEntity().Integrate("x")` | `integral(...)` | 3,475 characters |
57+
4558
### `|` is divisibility, and was disjunction
4659

4760
**Read this one first.** It is the only entry here that changes what existing input *means* rather

‎Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs‎

Lines changed: 69 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17910,6 +17910,60 @@ static Dictionary<EInteger, Entity> Scaled(Dictionary<EInteger, Entity> poly, En
1791017910
return answer;
1791117911
}
1791217912

17913+
/// <summary>
17914+
/// The solution of a linear system whose entries are all rational numbers, by Gauss-Jordan
17915+
/// elimination in exact rationals kept in lowest terms -- PeterO's rationals are not reduced
17916+
/// on their own, and unreduced the solve for Welz's piece took two minutes -- an unknown no
17917+
/// row pivots on set to zero; null where an
17918+
/// entry is not a rational number or the system is inconsistent.
17919+
/// </summary>
17920+
private static Entity[]? TrySolveInRationals(Entity[][] matrix, Entity[] rhs)
17921+
{
17922+
static ERational? Read(Entity entry) => entry is Number.Rational r ? r.ERational : entry.Evaled is Number.Rational e ? e.ERational : null;
17923+
var rows = matrix.Length;
17924+
var columns = rows == 0 ? 0 : matrix[0].Length;
17925+
var m = new ERational[rows][];
17926+
for (var i = 0; i < rows; i++)
17927+
{
17928+
m[i] = new ERational[columns + 1];
17929+
for (var j = 0; j < columns; j++)
17930+
if (Read(matrix[i][j]) is { } value) m[i][j] = value; else return null;
17931+
if (Read(rhs[i]) is { } right) m[i][columns] = right; else return null;
17932+
}
17933+
var pivotOf = new int[columns];
17934+
for (var j = 0; j < columns; j++) pivotOf[j] = -1;
17935+
var row = 0;
17936+
for (var column = 0; column < columns && row < rows; column++)
17937+
{
17938+
var pivot = -1;
17939+
for (var i = row; i < rows; i++)
17940+
if (!m[i][column].IsZero) { pivot = i; break; }
17941+
if (pivot < 0)
17942+
continue;
17943+
(m[row], m[pivot]) = (m[pivot], m[row]);
17944+
var lead = m[row][column];
17945+
for (var j = column; j <= columns; j++)
17946+
m[row][j] = m[row][j].Divide(lead).ToLowestTerms();
17947+
for (var i = 0; i < rows; i++)
17948+
{
17949+
if (i == row || m[i][column].IsZero)
17950+
continue;
17951+
var factor = m[i][column];
17952+
for (var j = column; j <= columns; j++)
17953+
m[i][j] = m[i][j].Subtract(factor.Multiply(m[row][j])).ToLowestTerms();
17954+
}
17955+
pivotOf[column] = row;
17956+
row++;
17957+
}
17958+
for (var i = row; i < rows; i++)
17959+
if (!m[i][columns].IsZero)
17960+
return null;
17961+
var values = new Entity[columns];
17962+
for (var j = 0; j < columns; j++)
17963+
values[j] = pivotOf[j] < 0 ? Number.Integer.Zero : Number.Rational.Create(m[pivotOf[j]][columns]);
17964+
return values;
17965+
}
17966+
1791317967
/// <summary>
1791417968
/// A polynomial times an odd half power of a quadratic, <c>P Q^(m/2)</c>, reduced to
1791517969
/// <c>R Q^(k + 1/2) + K/sqrt(Q)</c> by one linear solve, with <c>K/sqrt(Q)</c> the
@@ -17968,7 +18022,13 @@ static Dictionary<EInteger, Entity> Scaled(Dictionary<EInteger, Entity> poly, En
1796818022
|| pRead.Values.Any(coefficient => coefficient.ContainsNode(x)))
1796918023
return null;
1797018024
var degreeOfP = pRead.Count == 0 ? 0 : pRead.Keys.Max()!.ToInt32Unchecked();
17971-
if (degreeOfP > 8 || m > 9 || m < -9)
18025+
// The solve is one linear system of about degreeOfP + |m| unknowns; with numbers for every
18026+
// coefficient it stays small to far higher powers, and Welz's
18027+
// `1/((3 - 2 x + x^2)^(11/2) (1 + x + 2 x^2)^5)` left a piece of degree 9 over the
18028+
// eleventh half power, and its `^(21/2)`, `^10` twin one of degree 19.
18029+
var numeric = !radicand.Vars.Any(v => v != x) && !polynomial.Vars.Any(v => v != x);
18030+
var (degreeCap, powerCap) = numeric ? (24, 25) : (8, 9);
18031+
if (degreeOfP > degreeCap || m > powerCap || m < -powerCap)
1797218032
return null;
1797318033

1797418034
// Each coefficient that is neither a number nor a symbol is named for the solve and
@@ -17981,6 +18041,10 @@ Entity Named(Entity value)
1798118041
{
1798218042
if (value is Number or Variable)
1798318043
return value;
18044+
// A number written as an expression, `-2` as `-1 * 2`, is that number: named, it kept
18045+
// the system out of the exact rational solve below.
18046+
if (!value.Vars.Any() && value.Evaled is Number.Rational number)
18047+
return number;
1798418048
foreach (var (name, named) in names)
1798518049
if (named == value)
1798618050
return name;
@@ -18034,7 +18098,10 @@ Dictionary<EInteger, Entity> PowerOfQ(int j)
1803418098
matrix[row][column] = columns[column].TryGetValue(power, out var entry) ? entry : Number.Integer.Zero;
1803518099
rhs[row] = target.TryGetValue(power, out var wanted) ? wanted : Number.Integer.Zero;
1803618100
}
18037-
if (!Functions.PartialFractions.TrySolveLinear(matrix, rhs, out var values) || values is null)
18101+
// With numbers for every entry the system is solved in exact rationals: eliminated as
18102+
// expressions, the values of Welz's degree-9 piece over the eleventh half power came to
18103+
// half a million nodes of unfolded fractions and were refused at the bound.
18104+
if ((TrySolveInRationals(matrix, rhs) ?? (Functions.PartialFractions.TrySolveLinear(matrix, rhs, out var solved) ? solved : null)) is not { } values)
1803818105
return null;
1803918106
Entity r = Number.Integer.Zero;
1804018107
for (var k = 0; k <= degreeOfR; k++)
Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
//
2+
// Copyright (c) 2019-2026 Angouri.
3+
// AngouriMath is licensed under MIT.
4+
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
5+
// Website: https://am.angouri.org.
6+
//
7+
8+
using System;
9+
using AngouriMath.Extensions;
10+
using Xunit;
11+
12+
namespace AngouriMath.Tests.Calculus
13+
{
14+
/// <summary>
15+
/// A high half-odd power of a quadratic beside a power of another: the pieces over the first are a
16+
/// polynomial times its half-odd power, reduced by one linear solve, which with numbers for
17+
/// coefficients is solved in exact rationals. Welz's <c>1/((3 - 2 x + x^2)^(11/2) (1 + x + 2 x^2)^5)</c>.
18+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
19+
/// </summary>
20+
[Trait("Area", "Calculus")]
21+
public sealed class AHighHalfPowerOfAQuadraticBesideAPowerOfAnotherIntegralTest
22+
{
23+
private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
24+
{
25+
var integral = integrand.ToEntity().Integrate("x");
26+
Assert.DoesNotContain("integral(", integral.Stringize());
27+
28+
var pinned = integral.Substitute("C", 0);
29+
Entity original = integrand.ToEntity();
30+
foreach (var (name, value) in pins)
31+
{
32+
pinned = pinned.Substitute(name, value);
33+
original = original.Substitute(name, value);
34+
}
35+
var derivative = pinned.Differentiate("x");
36+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
37+
{
38+
var got = derivative.Substitute("x", at).EvalNumerical();
39+
var want = original.Substitute("x", at).EvalNumerical();
40+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
41+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
42+
Assert.True(difference / scale < 1e-8,
43+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
44+
}
45+
}
46+
47+
// d with the sign that keeps the argument real at every sampled x, and |d| below 2/2.3^2
48+
// so that the inverse sine and cosine stay inside [-1, 1]; off it the exact checks hold too.
49+
[Theory]
50+
[InlineData("1/((3-2*x+x^2)^(11/2)*(1+x+2*x^2)^5)")]
51+
[InlineData("1/((3-2*x+x^2)^(21/2)*(1+x+2*x^2)^10)")]
52+
[InlineData("(1+x^9)*(3-2*x+x^2)^(-11/2)")]
53+
public void IsAnsweredAndDifferentiatesBack(string integrand)
54+
=> DifferentiatesBack(integrand);
55+
}
56+
}

0 commit comments

Comments
 (0)