Skip to content

A quotient in x^2 over a power of x and a biquadratic is split in x^2 #1806

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Oct 6, 2026
Merged

A quotient in x^2 over a power of x and a biquadratic is split in x^2 #1806

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 18 additions & 0 deletions BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,6 +336,24 @@ quotient of two such linears the sum is `(b - d t^2)^2 + (c t^2 - a)^2`. Rubi's
| `"(c + d*tan(x))^(3/2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)` | the same |
| `"1/((a + b*tan(x))^(3/2)*(c + d*tan(x))^(3/2))".ToEntity().Integrate("x")` | `integral(...)`; past a minute on the unreleased master | in the root of the quotient of the two, in a second |

### 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
and answered in a million characters after eighteen seconds since; it is answered in two thousand,
in under a second. Under `u = tan(x)` and the root of the quotient `t = sqrt(u/(a + b u))`, Rubi's
half-odd powers of the cotangent beside one of `a + b tan(x)` leave a polynomial in `t^2` over
`t^(2k) Q(t^2)`, `Q` a biquadratic, which was split in `t`: the repeated `t^(2k)` beside the quartic
went to the conjugates of a sum of two squares, at length, or, with the power of a product
distributed, to the Hermite reduction, which declined it after seconds. With a power of `x` on both
sides cancelled, such a quotient is split in `w = x^2` now: the terms in `w^(-j)` are the expansion
of the remainder over `Q` at `w = 0`, and what is left is `(d + e w)/Q`, which the biquadratic's rule
answers by arctangents ([#718](https://github.com/asc-community/AngouriMath/issues/718)).

| Input | Was (2.5.0) | Now |
|---|---|---|
| `"cot(x)^(13/2)*(a + b*tan(x))^(5/2)*(A + B*tan(x))".ToEntity().Integrate("x")` | `integral(...)`; 1,018,377 characters after 18 seconds on the unreleased master | 2,305 characters, in under a second |
| `"cot(x)^(5/2)/(a + b*tan(x))^(3/2)".ToEntity().Integrate("x")` | `integral(...)`; 4,779 characters after 6 seconds on the unreleased master | 955 characters |

### A function of `x^n` for a symbolic `n` beside a power of `x` is integrated in a power of `x`

**Answers where there were none.** `x^(-1 + 4n)/(a + b x^n + c x^(2n))` and the rest of Rubi's
Expand Down
92 changes: 69 additions & 23 deletions Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2422,7 +2422,15 @@ over is Entity.Powf(var @base, var power) ?
/// <c>(d + e r_1)/(q (x^2 - r_1)) - (d + e r_2)/(q (x^2 - r_2))</c>, and
/// <c>1/(x^2 - r)</c> is <c>atan(x/sqrt(-r))/sqrt(-r)</c> for any complex <c>r</c>. A
/// higher even degree is divided down first, and an odd numerator is left to the
/// substitution <c>u = x^2</c>.
/// substitution <c>u = x^2</c>. And over <c>x^(2k)</c> times the biquadratic, a power of
/// <c>x</c> on both sides cancelled first: in <c>w = x^2</c> the remainder <c>R</c> over
/// <c>w^k Q</c> is the expansion of <c>R/Q</c> at <c>w = 0</c> up to <c>w^(k-1)</c>, over
/// <c>w^k</c>, plus <c>(d + e w)/Q</c>, so the terms in <c>x^(-2j)</c> go by the power
/// rule and the rest as above. Split in <c>x</c>, the repeated factor <c>x^(2k)</c> beside
/// the quartic went to the Hermite reduction or to the conjugates of a sum of two squares,
/// and <c>2 a (1 - b t^2)^4/(t^4 (1 - 2 b t^2 + (a^2 + b^2) t^4))</c>, which the root of the
/// quotient makes of Rubi's <c>cot(c + d x)^(5/2)/(a + b tan(c + d x))^(3/2)</c>, was
/// declined after seconds in one spelling and answered at length in the other.
/// </summary>
/// <remarks>
/// <para>
Expand All @@ -2441,46 +2449,76 @@ over is Entity.Powf(var @base, var power) ?
{
if (!TryReadAsQuotient(expr, out var numerator, out var denominator))
return null;
if (!TreeAnalyzer.TryGetPolynomial(denominator, x, out var below) || below.Count == 0
|| !below.ContainsKey(EInteger.FromInt32(4)) || below.Keys.Any(power => !(power.IsZero || power.Equals(EInteger.FromInt32(2)) || power.Equals(EInteger.FromInt32(4))))
|| below.Values.Any(coefficient => coefficient.ContainsNode(x))
if (!TreeAnalyzer.TryGetPolynomial(denominator, x, out var belowAsWritten) || belowAsWritten.Count == 0
|| belowAsWritten.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32())
|| belowAsWritten.Values.Any(coefficient => coefficient.ContainsNode(x))
|| !denominator.Vars.Any(v => v != x))
return null;
var a = below.TryGetValue(EInteger.Zero, out var a0) ? a0 : Number.Integer.Zero;
var b = below.TryGetValue(EInteger.FromInt32(2), out var b0) ? b0 : Number.Integer.Zero;
var c = below[EInteger.FromInt32(4)];
if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var aboveAsWritten) || aboveAsWritten.Count == 0
|| aboveAsWritten.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32())
|| aboveAsWritten.Values.Any(coefficient => coefficient.ContainsNode(x)))
return null;
// A power of x on both sides is cancelled first, and what is left below is x^(2k) times the
// biquadratic: `2 a (1 - b t^2)^4 t/(t^5 ((1 - b t^2)^2 + a^2 t^4))` is how the root of
// a quotient writes Rubi's `cot(c + d x)^(5/2)/(a + b tan(c + d x))^(3/2)`.
var common = EInteger.Min(aboveAsWritten.Keys.Min()!, belowAsWritten.Keys.Min()!).ToInt32Unchecked();
var below = belowAsWritten.ToDictionary(pair => pair.Key.ToInt32Unchecked() - common, pair => pair.Value);
var above = aboveAsWritten.ToDictionary(pair => pair.Key.ToInt32Unchecked() - common, pair => pair.Value);
var lowest = below.Keys.Min();
if (lowest % 2 != 0 || !below.ContainsKey(lowest + 4)
|| below.Keys.Any(power => power != lowest && power != lowest + 2 && power != lowest + 4)
|| above.Keys.Any(power => power % 2 != 0))
return null;
var k = lowest / 2;
var a = below[lowest];
var b = below.TryGetValue(lowest + 2, out var b0) ? b0 : Number.Integer.Zero;
var c = below[lowest + 4];
// Zero as a value, and not only as written: the coefficients are read off the
// denominator as it is written, and `((c - u^2)/c - 1)(u^2 - c)` is `-u^2 (u^2 - c)/c`,
// with `-(c/c - 1) c` for its constant term. Read as a biquadratic, one of its roots
// in `u^2` was that zero, and `atan(u/sqrt(0))/sqrt(0)` made the answer NaN.
// https://github.com/asc-community/AngouriMath/issues/1665
if (Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(c))
return null;
if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var above) || above.Count == 0
|| above.Keys.Any(power => power.Sign < 0 || !power.IsEven) || above.Values.Any(coefficient => coefficient.ContainsNode(x)))
return null;

var discriminant = Functions.PartialFractions.Bare((b * b - 4 * a * c).Simplify());
if (Functions.PartialFractions.IsZeroAsAValue(discriminant))
return null;

// The numerator in w = x^2, divided down by c w^2 + b w + a to a remainder d + e w.
var degree = above.Keys.Max()!.ToInt32Checked() / 2;
var inW = new Entity[degree + 1];
for (var k = 0; k <= degree; k++)
inW[k] = above.TryGetValue(EInteger.FromInt32(2 * k), out var at) ? at : Number.Integer.Zero;
// The numerator in w = x^2, divided down by w^k (c w^2 + b w + a) to a remainder of a
// degree below k + 2.
var degree = above.Keys.Max() / 2;
var inW = new Entity[System.Math.Max(degree, k + 1) + 1];
for (var j = 0; j < inW.Length; j++)
inW[j] = above.TryGetValue(2 * j, out var at) ? at : Number.Integer.Zero;
Entity polynomialPart = Number.Integer.Zero;
for (var k = degree; k >= 2; k--)
for (var j = degree; j >= k + 2; j--)
{
var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k] / c);
var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j] / c);
if (lead == Number.Integer.Zero)
continue;
polynomialPart += lead * MathS.Pow(x, 2 * (k - 2));
inW[k - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k - 1] - lead * b);
inW[k - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k - 2] - lead * a);
}
var d = inW[0];
var e = degree >= 1 ? inW[1] : Number.Integer.Zero;
polynomialPart += lead * MathS.Pow(x, 2 * (j - k - 2));
inW[j - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 1] - lead * b);
inW[j - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 2] - lead * a);
}
// The remainder R over w^k Q is sum r_m w^(m - k) over m < k, the expansion of R/Q at
// w = 0, plus (d + e w)/Q: r_m = (R_m - b r_(m-1) - c r_(m-2))/a, and d + e w is what
// (R - Q sum r_m w^m)/w^k leaves, every lower coefficient cancelling by the recurrence.
var pole = new Entity[k];
for (var m = 0; m < k; m++)
{
var term = inW[m];
if (m >= 1) term -= b * pole[m - 1];
if (m >= 2) term -= c * pole[m - 2];
pole[m] = Functions.PartialFractions.InLowestTermsOverTheSymbols(term / a);
}
var dTerm = inW[k];
if (k >= 1) dTerm -= b * pole[k - 1];
if (k >= 2) dTerm -= c * pole[k - 2];
var eTerm = inW[k + 1];
if (k >= 1) eTerm -= c * pole[k - 1];
var d = k == 0 ? inW[0] : Functions.PartialFractions.InLowestTermsOverTheSymbols(dTerm);
var e = k == 0 ? (degree >= 1 ? inW[1] : Number.Integer.Zero) : Functions.PartialFractions.InLowestTermsOverTheSymbols(eTerm);

var q = MathS.Sqrt(discriminant);
var firstRoot = (-b + q) / (2 * c);
Expand All @@ -2492,6 +2530,14 @@ over is Entity.Powf(var @base, var power) ?
return null;
total += whole;
}
// r_m x^(2m - 2k), whose exponent is odd and so never minus one.
for (var m = 0; m < k; m++)
{
if (pole[m] == Number.Integer.Zero)
continue;
var exponent = 2 * (m - k) + 1;
total += pole[m] * MathS.Pow(x, exponent) / exponent;
}
// `1/(x^2 - r)` is `atan(x/s)/s` with `s = sqrt(-r)` for every complex `r` but zero:
// `d/dx atan(x/s)/s` is `1/(s^2 + x^2)` whatever `s` is, and for a positive `r` the
// arctangent of an imaginary argument is the hyperbolic one, `-atanh(x/sqrt(r))/sqrt(r)`,
Expand Down
59 changes: 59 additions & 0 deletions Sources/Tests/UnitTests/Calculus/EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
//
// 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
{
/// <summary>
/// A polynomial in <c>x^2</c> over <c>x^(2k)</c> times a biquadratic with a symbol in it, a power
/// of <c>x</c> on both sides cancelled first, split in <c>w = x^2</c>: the terms in
/// <c>w^(-j)</c> are the expansion of the remainder over the biquadratic at <c>w = 0</c>, and
/// what is left is the biquadratic's own. A half-odd power of the cotangent beside one of
/// <c>a + b tan(x)</c> comes to one under <c>u = tan(x)</c> and the root of the quotient
/// <c>t = sqrt(u/(a + b u))</c>; split in <c>t</c>, the answer to the last row ran to a million
/// characters. Rubi's 4.3.2.1 and 4.3.3.1.
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
/// </summary>
/// <remarks>
/// The cotangent's rows are real where <c>tan(x)</c> and <c>a + b tan(x)</c> are positive,
/// which the points are. The pinned answer is differentiated as it is, since simplifying
/// answers of this shape costs minutes.
/// </remarks>
[Trait("Area", "Calculus")]
public sealed class EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest
{
[Theory]
[InlineData("1/(x^2*(a + b*x^2 + c*x^4))")]
[InlineData("(1 + x^10)/(x^2*(a + b*x^2 + c*x^4))")]
[InlineData("2*a*(1 - b*x^2)^4*x/(x^5*(1 - 2*b*x^2 + (a^2 + b^2)*x^4))")]
[InlineData("cot(x)^(5/2)/(a + b*tan(x))^(3/2)")]
[InlineData("cot(x)^(9/2)*sqrt(a + b*tan(x))*(A + B*tan(x))")]
[InlineData("cot(x)^(13/2)*(a + b*tan(x))^(5/2)*(A + B*tan(x))")]
public void SplitInTheSquare(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", 2.9)
.Substitute("A", 0.4).Substitute("B", 1.1);
var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
var original = Pinned(integrand.ToEntity());
foreach (var at in new[] { 0.2, 0.5, 0.9, 1.2 })
{
var want = original.Substitute("x", at).EvalNumerical();
var got = derivative.Substitute("x", at).EvalNumerical();
Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
< 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)),
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
}
}
}
}
Loading