Skip to content

Commit cba4232

Browse files
Take a rational function past the twelfth degree by Hermite's reduction (#1893)
* Take a rational function past the twelfth degree by Hermite's reduction A denominator up to the twenty-fourth degree is read, the squarefree part Hermite leaves still held to the twelfth for the resultant; and such a rational function with numbers for coefficients is asked before the search, its terms over the least common multiple of their denominators. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Hold the rational integrator's old degree below a substitution, and ask the early route only at the question asked Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the answers Hermite's reduction past the twelfth degree gives, and build on netstandard2.0 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 34780ec commit cba4232

5 files changed

Lines changed: 131 additions & 3 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -447,6 +447,18 @@ case is most of its length ([#718](https://github.com/asc-community/AngouriMath/
447447
| `"x^4/((a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 28,551 characters |
448448
| `"1/(x^2*(a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 10,194 characters |
449449

450+
### A rational function past the twelfth degree is reduced by Hermite's method
451+
452+
**Answers where there were none.** The integrator read a denominator of degree twelve at most, so a derivative over
453+
powers of `x^4 + x + 3` up to the fourth was declined, after seven to fifteen seconds. A rational function with numbers
454+
for coefficients over a denominator past the twelfth degree now takes Hermite's reduction, as long as what is left for the
455+
logarithms stays within twelve ([#718](https://github.com/asc-community/AngouriMath/issues/718)).
456+
457+
| Input | Was (2.5.0) | Now |
458+
|---|---|---|
459+
| `"3*(-47+228*x+120*x^2+19*x^3)/(3+x+x^4)^4+(42-320*x-75*x^2-8*x^3)/(3+x+x^4)^3+30*x/(3+x+x^4)^2".ToEntity().Integrate("x")` | `integral(...)` | 142 characters |
460+
| `"(-3+10*x+4*x^3-30*x^5)/(3+x+x^4)^3-3*(1+4*x^3)*(2-3*x+5*x^2+x^4-5*x^6)/(3+x+x^4)^4".ToEntity().Integrate("x")` | `integral(...)` | 142 characters |
461+
450462
### A power of a linear over a power of a quadratic is reduced a power at a time
451463

452464
**Answers where there were none.** `sqrt(d + k x)/(a + b x + c x^2)^3` took thirty seconds for 26,437 characters on the

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

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2621,6 +2621,50 @@ node is Powf(var radicand, Number.Rational r) && r is not Number.Integer && radi
26212621
/// <see cref="RothsteinTrager"/>. After <see cref="SolveByPartialFractions"/>, so that
26222622
/// everything that answers keeps the form it gives.
26232623
/// </summary>
2624+
internal static Entity? SolveAHighDegreeRationalFunctionByHermite(Entity expr, Entity.Variable x)
2625+
{
2626+
// A sum of quotients over powers of one denominator, Rubi's spelling of a derivative,
2627+
// is one quotient for this. At the question asked only, which is where nothing else
2628+
// answered one.
2629+
if (!Integration.AnsweringTheQuestionAsked || expr.Vars.Any(symbol => symbol != x))
2630+
return null;
2631+
// Over the least common multiple of the terms' denominators, each base to its highest
2632+
// power, rather than their product: three terms over the fourth, third and second
2633+
// powers of one quartic are one quotient over its fourth power.
2634+
var terms = Sumf.LinearChildren(expr).Select(term => Functions.SingleQuotient.Of(term)).ToList();
2635+
var highest = new Dictionary<Entity, int>();
2636+
var powersOf = new List<Dictionary<Entity, int>>();
2637+
foreach (var (_, below) in terms)
2638+
{
2639+
var powers = new Dictionary<Entity, int>();
2640+
foreach (var factor in Mulf.LinearChildren(below))
2641+
{
2642+
var (@base, whole) = factor is Powf(var b0, Number.Integer { EInteger: var w }) && w.Sign > 0 && w.CanFitInInt32()
2643+
? (b0, w.ToInt32Unchecked()) : (factor, 1);
2644+
powers[@base] = powers.TryGetValue(@base, out var already) ? already + whole : whole;
2645+
}
2646+
foreach (var pair in powers)
2647+
highest[pair.Key] = System.Math.Max(highest.TryGetValue(pair.Key, out var already) ? already : 0, pair.Value);
2648+
powersOf.Add(powers);
2649+
}
2650+
Entity denominator = Number.Integer.One;
2651+
foreach (var pair in highest)
2652+
denominator = denominator * MathS.Pow(pair.Key, Number.Integer.Create(pair.Value));
2653+
Entity numerator = Number.Integer.Zero;
2654+
for (var i = 0; i < terms.Count; i++)
2655+
{
2656+
Entity cofactor = Number.Integer.One;
2657+
foreach (var pair in highest)
2658+
{
2659+
var missing = pair.Value - (powersOf[i].TryGetValue(pair.Key, out var have) ? have : 0);
2660+
if (missing > 0)
2661+
cofactor = cofactor * MathS.Pow(pair.Key, Number.Integer.Create(missing));
2662+
}
2663+
numerator = numerator + terms[i].Numerator * cofactor;
2664+
}
2665+
return RothsteinTrager.IntegrateAHighDegree(numerator, denominator, x);
2666+
}
2667+
26242668
internal static Entity? SolveByRothsteinTrager(Entity expr, Entity.Variable x)
26252669
{
26262670
if (!TryReadAsQuotient(expr, out var numerator, out var denominator))

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -666,6 +666,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
666666
// the question asked only, and the polynomial term below writes `1/u` as `u^(-1)` and
667667
// asks that a level down: `1/(e^atanh(a x) (c - c/(a x))^(3/2))` reached it nowhere.
668668
if ((answer = IndefiniteIntegralSolver.SolveByWritingEachLinearOnce(expr, x, integrateByParts)) is { }) return answer;
669+
// A rational function with numbers for coefficients over a denominator past the twelfth
670+
// degree, by Hermite's reduction before anything searches: nothing below answered one.
671+
if ((answer = IndefiniteIntegralSolver.SolveAHighDegreeRationalFunctionByHermite(expr, x)) is { }) return answer;
669672
// A power of a linear beside the derivative of a quadratic over a whole power of it, by
670673
// parts once, before anything takes the symbolic quadratic apart.
671674
if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearBesideTheDerivativeOfAQuadraticOverItsPower(expr, x, integrateByParts)) is { }) return answer;

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

Lines changed: 31 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,16 +63,44 @@ internal static class RothsteinTrager
6363
/// </summary>
6464
private const int MaxDegree = 12;
6565

66+
/// <summary>
67+
/// The largest denominator degree read at all: Hermite's reduction is a few Diophantine
68+
/// solves over the rationals, and what it leaves for the resultant is the squarefree part,
69+
/// held to <see cref="MaxDegree"/>. Rubi's 1.3.1 <c>3 (19 x^3 + 120 x^2 + 228 x - 47)/(x^4 + x + 3)^4 + ...</c>
70+
/// is a derivative over a sixteenth-degree denominator with nothing left for the logarithms,
71+
/// and was declined after a minute and a half of searching elsewhere.
72+
/// </summary>
73+
private const int MaxReducedDegree = 24;
74+
6675
/// <summary>
6776
/// The integral of <paramref name="numerator"/> over <paramref name="denominator"/>,
6877
/// both polynomials in <paramref name="x"/> with rational coefficients and the fraction
6978
/// proper; <see langword="null"/> where they are not, or where the answer does not
7079
/// differentiate back to the integrand at the sampled points.
7180
/// </summary>
81+
/// <summary>
82+
/// <see cref="Integrate"/>, for a denominator past <see cref="MaxDegree"/> only: what the
83+
/// rest of the integrator never answered, asked before it searches.
84+
/// </summary>
85+
internal static Entity? IntegrateAHighDegree(Entity numerator, Entity denominator, Variable x)
86+
=> PolynomialFactoring.TryGetRationalCoefficients(denominator, x, 1, 1, MaxReducedDegree, out var belowCoefficients)
87+
&& RationalPolynomial.Create(belowCoefficients).Degree > MaxDegree
88+
? IntegrateUpToDegree(numerator, denominator, x, MaxReducedDegree)
89+
: null;
90+
91+
/// <summary>
92+
/// The integral of <paramref name="numerator"/> over <paramref name="denominator"/>,
93+
/// held to <see cref="MaxDegree"/> unless asked otherwise: below a substitution a
94+
/// denominator past it -- <c>1/(cos(x) + cos(3x))^5</c> under the tangent's half angle --
95+
/// was reduced past the budget where it used to be declined at once.
96+
/// </summary>
7297
internal static Entity? Integrate(Entity numerator, Entity denominator, Variable x)
98+
=> IntegrateUpToDegree(numerator, denominator, x, MaxDegree);
99+
100+
private static Entity? IntegrateUpToDegree(Entity numerator, Entity denominator, Variable x, int degree)
73101
{
74-
if (!PolynomialFactoring.TryGetRationalCoefficients(numerator, x, 1, 0, MaxDegree, out var aboveCoefficients)
75-
|| !PolynomialFactoring.TryGetRationalCoefficients(denominator, x, 1, 1, MaxDegree, out var belowCoefficients))
102+
if (!PolynomialFactoring.TryGetRationalCoefficients(numerator, x, 1, 0, degree, out var aboveCoefficients)
103+
|| !PolynomialFactoring.TryGetRationalCoefficients(denominator, x, 1, 1, degree, out var belowCoefficients))
76104
return null;
77105
var above = RationalPolynomial.Create(aboveCoefficients);
78106
var below = RationalPolynomial.Create(belowCoefficients);
@@ -87,7 +115,7 @@ internal static class RothsteinTrager
87115
var answer = rationalPart;
88116
if (!logAbove.IsZero)
89117
{
90-
if (LogarithmicPart(logAbove, logBelow, x) is not { } logarithmic)
118+
if (logBelow.Degree > MaxDegree || LogarithmicPart(logAbove, logBelow, x) is not { } logarithmic)
91119
return null;
92120
answer = answer + logarithmic;
93121
}
Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
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 rational function with numbers for coefficients over a denominator past the twelfth
16+
/// degree, by Hermite's reduction: Rubi's 1.3.1 derivatives over powers of
17+
/// <c>x^4 + x + 3</c>, written as sums of quotients over different powers of it.
18+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
19+
/// </summary>
20+
[Trait("Area", "Calculus")]
21+
public sealed class AHighDegreeRationalFunctionIntegralTest
22+
{
23+
[Theory]
24+
[InlineData("3*(-47+228*x+120*x^2+19*x^3)/(3+x+x^4)^4+(42-320*x-75*x^2-8*x^3)/(3+x+x^4)^3+30*x/(3+x+x^4)^2")]
25+
[InlineData("(-3+10*x+4*x^3-30*x^5)/(3+x+x^4)^3-3*(1+4*x^3)*(2-3*x+5*x^2+x^4-5*x^6)/(3+x+x^4)^4")]
26+
public void IsAnsweredAndDifferentiatesBack(string integrand)
27+
{
28+
var integral = integrand.ToEntity().Integrate("x");
29+
Assert.DoesNotContain("integral(", integral.Stringize());
30+
var derivative = integral.Substitute("C", 0).Differentiate("x");
31+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
32+
{
33+
var got = derivative.Substitute("x", at).EvalNumerical();
34+
var want = integrand.ToEntity().Substitute("x", at).EvalNumerical();
35+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
36+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart));
37+
Assert.True(difference / scale < 1e-8, $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
38+
}
39+
}
40+
}
41+
}

0 commit comments

Comments
 (0)