Skip to content

A root of a quotient is a quotient of roots where that is exact, and a root gathers with a whole power of its own base #1300

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’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Sep 12, 2026
Merged

A root of a quotient is a quotient of roots where that is exact, and a root gathers with a whole power of its own base #1300

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
145 changes: 120 additions & 25 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 @@ -5051,9 +5051,12 @@ private static bool TryReadAPowerTimesARadicalQuadratic(Entity expr, Entity.Vari
// Two square roots holding the variable, or there is nothing to combine; counted
// before any of the rewriting below is paid for, since this runs on every
// sub-integrand of the chain.
if (expr.Nodes.Count(node => node is Powf(var @base, Number.Rational half)
// A root of a quotient counts as the two roots it splits into.
if (expr.Nodes.Sum(node => node is Powf(var @base, Number.Rational half)
&& half is not Number.Integer && half.ERational.Denominator.Equals(EInteger.FromInt32(2))
&& @base.ContainsNode(x)) < 2)
&& @base.ContainsNode(x)
? (@base is Divf || Functions.SingleQuotient.Of(@base).Denominator != Number.Integer.One ? 2 : 1)
: 0) < 2)
return null;
var rewritten = CombineRadicalsIn(expr, x);
if (rewritten == expr)
Expand All @@ -5069,7 +5072,48 @@ private static bool TryReadAPowerTimesARadicalQuadratic(Entity expr, Entity.Vari
/// reads as it stands.
/// </summary>
private static Entity CombineRadicalsIn(Entity expr, Entity.Variable x)
=> expr.Replace(node => node is Mulf or Divf ? CombineRadicalsInAQuotient(node, x) ?? node : node);
=> SplitRootsOfQuotientsIn(expr, x).Replace(node => node is Mulf or Divf ? CombineRadicalsInAQuotient(node, x) ?? node : node);

/// <summary>
/// Every root of a quotient of polynomials in <paramref name="expr"/> written as a
/// quotient of roots, where that is exact: <c>sqrt(P/Q) = sqrt(P)/sqrt(Q)</c> on the
/// principal branch unless <c>Q &lt; 0 &lt; P</c>, where the left is <c>i sqrt(P/|Q|)</c>
/// and the right <c>-i sqrt(P/|Q|)</c> -- with both negative the two <c>i</c>s cancel and
/// with <c>P</c> negative alone they agree. So the rule asks whether there is a real
/// <c>x</c> with <c>Q(x) &lt; 0 &lt; P(x)</c>, on one point of each interval between the
/// real roots, and splits only where there is not.
/// </summary>
/// <remarks>
/// Charlwood's <c>arcsin(x/sqrt(1 - x^2))</c> by parts against one leaves
/// <c>x (1 - x^2)^(-3/2)/sqrt((1 - 2x^2)/(1 - x^2))</c>, a root of a quotient that no
/// rule reads; split, it is <c>x/((1 - x^2) sqrt(1 - 2x^2))</c>, which Euler answers.
/// <c>1 - x^2 &lt; 0</c> only past <c>1</c>, where <c>1 - 2x^2</c> is negative too.
/// </remarks>
private static Entity SplitRootsOfQuotientsIn(Entity expr, Entity.Variable x)
=> expr.Replace(node =>
{
if (node is not Powf(var @base, Number.Rational power) || power is Number.Integer
|| !power.ERational.Denominator.Equals(EInteger.FromInt32(2)) || !@base.ContainsNode(x))
return node;
var (above, below) = Functions.SingleQuotient.Of(@base);
if (below == Number.Integer.One || !below.ContainsNode(x))
return node;
if (!TreeAnalyzer.TryGetPolynomial(above, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(below, x, out var belowRead)
|| above.Vars.Any(v => v != x) || below.Vars.Any(v => v != x))
return node;
// Quadratics at most on either side: those are what the rules behind this read,
// and a root of a sextic set free is a search that ends nowhere --
// `sin(x)/sqrt(1 - sin(x)^6)` under the tangent is one, and went from a
// twenty-millisecond decline to a twenty-second one.
static int Degree(Dictionary<EInteger, Entity> read) => read.Count == 0 ? 0 : read.Keys.Max()!.ToInt32Checked();
if (Degree(aboveRead) > 2 || Degree(belowRead) > 2)
return node;
// Kept as written, so that the base of the root below matches the same base
// written beside it as a whole power, for the gathering that follows.
if (!OnEveryRealInterval(new List<Entity> { above, below }, x, signs => !(signs[1] < 0 && signs[0] > 0)))
return node;
return MathS.Pow(above, power) / MathS.Pow(below, power);
});

/// <summary>
/// A factor of a product that is a small power of a sum holding a square root of a
Expand Down Expand Up @@ -5124,41 +5168,84 @@ Entity WrittenOut(Entity side)
var cancelled = CancelEqualSumFactors(numerator, denominator, out numerator, out denominator) | writtenOut;

// Every factor with an odd half-power of a polynomial base is a whole power of the
// base times one square root of it; the square roots are what combine.
// base times one square root of it; the square roots are what combine. A whole
// power of the same base, written beside it, is gathered in first: `P * P^(-1/2)`
// is `P^(1/2)` exactly, since a whole power combines with any principal power of
// the same base -- `(1 - x^2) sqrt(1 - 2x^2)/sqrt(1 - x^2)` is
// `sqrt(1 - x^2) sqrt(1 - 2x^2)`, and what a root of a quotient splits into is
// often beside its own base that way.
var bases = new List<Entity>();
Entity above = Number.Integer.One;
Entity below = Number.Integer.One;
foreach (var factor in Mulf.LinearChildren(numerator))
var halves = new Dictionary<Entity, int>(); // base -> exponent in halves
var others = new List<(Entity Factor, bool Below)>();
var gathered = false;
foreach (var (side, isBelow) in new[] { (numerator, false), (denominator, true) })
foreach (var factor in Mulf.LinearChildren(side))
{
var sign = isBelow ? -1 : 1;
if (TryReadASquareRootOfAPolynomial(factor, x, out var @base, out var wholePower))
halves[@base] = halves.TryGetValue(@base, out var so) ? so + sign * (2 * wholePower + 1) : sign * (2 * wholePower + 1);
else if (factor is Powf(var b, Number.Integer k) && k.EInteger.CanFitInInt32() && halves.ContainsKey(b))
{
halves[b] += sign * 2 * k.EInteger.ToInt32Unchecked();
gathered = true;
}
else if (factor.ContainsNode(x) && halves.ContainsKey(factor))
{
halves[factor] += sign * 2;
gathered = true;
}
else
others.Add((factor, isBelow));
}
// A whole power read before its root came: gathered on a second pass.
foreach (var (factor, isBelow) in others.ToList())
{
if (TryReadASquareRootOfAPolynomial(factor, x, out var @base, out var wholePower))
var sign = isBelow ? -1 : 1;
if (factor is Powf(var b, Number.Integer k) && k.EInteger.CanFitInInt32() && halves.ContainsKey(b))
{
bases.Add(@base);
above = above * MathS.Pow(@base, wholePower);
halves[b] += sign * 2 * k.EInteger.ToInt32Unchecked();
others.Remove((factor, isBelow));
gathered = true;
}
else if (factor.ContainsNode(x) && halves.ContainsKey(factor))
{
halves[factor] += sign * 2;
others.Remove((factor, isBelow));
gathered = true;
}
else
above = above * factor;
}
foreach (var factor in Mulf.LinearChildren(denominator))
foreach (var (factor, isBelow) in others)
if (isBelow) below = below * factor;
else above = above * factor;
foreach (var pair in halves)
{
if (TryReadASquareRootOfAPolynomial(factor, x, out var @base, out var wholePower))
var n = pair.Value;
if (n % 2 == 0)
{
// 1/(B^q sqrt(B)) = sqrt(B)/B^(q+1)
bases.Add(@base);
below = below * MathS.Pow(@base, wholePower + 1);
// A whole power after all: no root left of this base.
var whole = n / 2;
if (whole > 0) above = above * MathS.Pow(pair.Key, whole);
else if (whole < 0) below = below * MathS.Pow(pair.Key, -whole);
continue;
}
else
below = below * factor;
// n = 2q + 1 for every odd n, negative ones included: -1 = 2(-1) + 1.
var q = (n - 1) / 2;
bases.Add(pair.Key);
if (q > 0) above = above * MathS.Pow(pair.Key, q);
else if (q < 0) below = below * MathS.Pow(pair.Key, -q);
}
if (bases.Count < 2)
if (bases.Count < 2 || !AtMostOneIsNegativeOnTheReals(bases, x))
{
if (!cancelled)
// Nothing to combine, or roots that must not be: what was cancelled or
// gathered is still worth handing on, each root on its own.
if (!cancelled && !gathered)
return null;
foreach (var @base in bases)
above = above * MathS.Sqrt(@base);
bases.Clear();
}
else if (!AtMostOneIsNegativeOnTheReals(bases, x))
return null;

Entity radical = Number.Integer.One;
if (bases.Count > 0)
Expand Down Expand Up @@ -5231,6 +5318,15 @@ private static bool TryReadASquareRootOfAPolynomial(Entity factor, Entity.Variab
/// <see langword="false"/> where the roots cannot be had.
/// </summary>
private static bool AtMostOneIsNegativeOnTheReals(List<Entity> bases, Entity.Variable x)
=> OnEveryRealInterval(bases, x, signs => signs.Count(sign => sign < 0) <= 1);

/// <summary>
/// Whether <paramref name="holds"/> is true of the signs of the polynomials
/// <paramref name="bases"/> on every interval between their real roots -- one point of
/// each interval decides it, since no base changes sign inside one. <see langword="false"/>
/// where the roots cannot be had.
/// </summary>
private static bool OnEveryRealInterval(List<Entity> bases, Entity.Variable x, System.Func<List<int>, bool> holds)
{
var roots = new List<double>();
foreach (var @base in bases)
Expand Down Expand Up @@ -5259,15 +5355,14 @@ private static bool AtMostOneIsNegativeOnTheReals(List<Entity> bases, Entity.Var
}
foreach (var at in samples)
{
var negative = 0;
var signs = new List<int>();
foreach (var @base in bases)
{
if (@base.Substitute(x, at).Evaled is not Number.Real value || !value.IsFinite)
return false;
if (value < 0)
negative++;
signs.Add(value < 0 ? -1 : value > 0 ? 1 : 0);
}
if (negative > 1)
if (!holds(signs))
return false;
}
return true;
Expand Down
34 changes: 34 additions & 0 deletions Sources/Tests/UnitTests/Calculus/CombinedRadicalsTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,40 @@ private static void DifferentiatesBack(string integrand, double[] points)
[InlineData("ln(x + sqrt(x^2 - 1))/(1 + x^2)^(3/2)")]
public void OneBaseMayBeNegative(string integrand) => DifferentiatesBack(integrand, new[] { 1.2, 1.6, 2.3, 3.1 });

/// <summary>
/// A root of a quotient of polynomials written as a quotient of roots, where that is
/// exact: <c>sqrt(P/Q) = sqrt(P)/sqrt(Q)</c> unless <c>Q &lt; 0 &lt; P</c> somewhere on the
/// reals. Charlwood's <c>arcsin(x/sqrt(1 - x^2))</c> by parts against one leaves
/// <c>x (1 - x^2)^(-3/2)/sqrt((1 - 2x^2)/(1 - x^2))</c>, which no rule reads and which
/// split is <c>x/((1 - x^2) sqrt(1 - 2x^2))</c>. Checked where the arcsine is real,
/// <c>|x| &lt; 1/sqrt(2)</c>, on both sides of zero.
/// </summary>
[Theory]
[InlineData("arcsin(x/sqrt(1 - x^2))")]
[InlineData("x*(1 - x^2)^(-3/2)/sqrt((1 - 2*x^2)/(1 - x^2))")]
public void ARootOfAQuotientIsSplit(string integrand) => DifferentiatesBack(integrand, new[] { -0.6, -0.3, 0.2, 0.45, 0.65 });

/// <summary>
/// <c>sqrt((1 + x)/(3 + 2x))</c> is not answered: <c>3 + 2x &lt; 0 &lt; 1 + x</c> nowhere, so
/// the split is exact, but the two roots it leaves are both negative below
/// <c>-3/2</c> and no rule takes them further; and <c>x/((1 - x^2) sqrt((1 - 2x^2)/(1 - x^2)))</c>
/// splits and gathers to <c>x/(sqrt(1 - x^2) sqrt(1 - 2x^2))</c>, two roots both
/// negative past <c>1</c>, where the integrand is real, so those are not combined
/// either. Either verdict but a wrong answer.
/// </summary>
[Theory]
[InlineData("sqrt((1 + x)/(3 + 2*x))")]
[InlineData("sqrt((1 + x^2)/(1 - x^2))")]
[InlineData("x/((1 - x^2)*sqrt((1 - 2*x^2)/(1 - x^2)))")]
public void ARootOfAQuotientTheRulesBehindCannotUse(string integrand)
{
var integral = integrand.ToEntity().Integrate("x");
Assert.DoesNotContain("NaN", integral.Stringize());
if (integral.Stringize().Contains("integral("))
return;
DifferentiatesBack(integrand, new[] { -0.7, -0.3, 0.2, 0.55, 0.85 });
}

/// <summary>
/// A small power of a sum of radicals, as a factor, is written out first:
/// Bondarenko's <c>1/(sqrt(1 - x) + sqrt(1 + x))^2</c> is <c>1/(2 + 2 sqrt(1 - x^2))</c>,
Expand Down
Loading