Skip to content

Commit bcfefc2

Browse files
Solve a polynomial through its factors rather than by its degree (#272) (#742)
* Solve a polynomial through its factors rather than by its degree (#272) The solver read a polynomial equation only as a polynomial: it gathered the monomials, took the degree, and applied the formula for that degree -- Cardano at three, Ferrari at four, nothing at all above four. A factorization was never looked for, and one handed over outright was expanded away before the degree was taken. So (x - 1)(x^2 - 3) = 0, whose roots are written on its face, came back as two nested cube roots of 26 + 18i, and a quintic that splits into a cubic's worth of whole roots and a quadratic fell through to the numeric solver. Two things are added. A product is now solved factor by factor, since a product is zero exactly where one of its factors is; and where the equation is written out as a sum, its rational roots are divided out first, which leaves a factor of lower degree for the formulas to answer exactly. Zeroing one factor only makes the product zero where the others are defined, so each root is checked against the whole product before it is kept -- x * (1/x) is undefined at 0, not zero there, and the verification the solvers share cannot catch that, since it drops a root only on positive evidence and a NaN residual is not evidence. Two-term polynomials are left alone. a*x^n + b is answered whole by inverting x^n = -b/a, which gives the roots as a + bi; splitting it off at its rational roots leaves the rest to be dug out of a quotient, and x^3 - 8 would read (-2 -+ sqrt(-12))/2 rather than (-1/2 +- i*sqrt(3)/2)*2. A test pins that. lim of the measured effect, against 40ce109: x^3 - x^2 - 3x + 3 = 0 -(-1 + (26 + 18i)^(1/3) + 10/(26 + 18i)^(1/3))/3, ... -> { 1, sqrt(3), -sqrt(3) } 160 ms -> 9 ms x^5 - x^4 - 7x^3 + 7x^2 + 12x - 12 = 0 { -2, -1.73205080756887741... - 2.3e-15i, 1, ... } -> { 1, 2, -2, sqrt(3), -sqrt(3) } 153 ms -> 12 ms (x - 1)(x^2 - 3) = 0 -> { 1, sqrt(3), -sqrt(3) } 138 ms -> 8 ms x^3 - 6x^2 + 11x - 6 = 0-> { 1, 2, 3 } 160 ms -> 1 ms 4929 unit tests and 130 F# tests pass; corpus 112/117 with 0 wrong, 0 error and 0 timeout, unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Pin the roots the split turns out to recover A sweep over 572 polynomials built from factors whose roots are known in advance -- every product and every expansion of two or three factors drawn from a pool of linear and irreducible quadratic ones -- says the fix does more than tidy the output. Before it, ten of those came back **short of roots**, and after it none do; no case gained a root that is not one, in either direction. Every one of the ten had the same shape, a rational root times a quadratic times x^2 + x + 1: x^5 - 2x^3 - x^2 + 2 = 0 { 1, sqrt(2), -sqrt(2) } -> all five, the pair (-1 +- i*sqrt(3))/2 included Above degree four the equation went to the numeric solver whole, and that searches the real axis, so the complex conjugate pair was simply absent -- and a set two roots short reads exactly like a complete one. Splitting off the rational root leaves a quartic, which is answered in one piece. This is the stronger half of #272 and it was not visible from the issue, which reports the answers as ugly rather than as incomplete. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
1 parent 61af68a commit bcfefc2

3 files changed

Lines changed: 274 additions & 0 deletions

File tree

‎Sources/AngouriMath/Functions/Continuous/Solvers/EquationSolver/AnalyticalEquationSolver.cs‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,10 +140,53 @@ static Entity TryDowncast(Entity equation, Variable x, Entity root)
140140
return Solve(expression, x, compensateSolving).Filter(predicate, x);
141141
case Piecewise p:
142142
return EquationSolver.SolvePiecewise(p, x, (e, x) => Solve(e, x, compensateSolving));
143+
// A product is zero exactly where one of its factors is. Without this the
144+
// factors are expanded back into a polynomial and the equation is answered
145+
// by the general formula for whatever degree that turns out to be, so
146+
// (x - 1)(x^2 - 3) = 0 -- handed over already factored -- came back as two
147+
// nested cube roots of 26 + 18i rather than as 1 and +-sqrt(3).
148+
// https://github.com/asc-community/AngouriMath/issues/272
149+
//
150+
// A factor free of x contributes no roots. It could still be zero itself,
151+
// which would make every x a solution, but that is a reading the solver
152+
// does not take anywhere else either: a * (x^2 - 3) = 0 already answered
153+
// +-sqrt(3) by dividing a out, so nothing changes here.
154+
case Mulf(var left, var right) when left.ContainsNode(x) || right.ContainsNode(x):
155+
{
156+
var factors = new List<Entity>(2);
157+
if (left.ContainsNode(x)) factors.Add(left);
158+
if (right.ContainsNode(x)) factors.Add(right);
159+
var product = expr;
160+
// Zeroing one factor makes the product zero only where the others are
161+
// defined: x * (1/x) is undefined at 0, not zero there. This has to be
162+
// checked here rather than left to the verification every solver shares,
163+
// because that deliberately keeps a root whose residual does not
164+
// evaluate to a number at all -- a root is only dropped on positive
165+
// evidence, and NaN is not evidence (see SolveStatement.IsSpurious).
166+
bool ProductIsDefinedAt(Entity root)
167+
=> root.Vars.Any()
168+
|| product.Substitute(x, root).Evaled is not Complex value
169+
|| value.IsFinite;
170+
var fromFactors = factors.Select(factor => Solve(factor, x, compensateSolving)).Unite();
171+
// Collapsed first: uniting sets builds a Unionf, and the roots cannot be
172+
// looked at one at a time until it is a finite set again.
173+
return fromFactors.InnerSimplified is FiniteSet found
174+
? found.Where(ProductIsDefinedAt).ToSet()
175+
: fromFactors;
176+
}
143177
default:
144178
break;
145179
}
146180

181+
// The same equation written out as a sum. Where a rational root can be divided
182+
// out, what is left is a polynomial of lower degree, and the factors are then
183+
// solved one at a time by the case above -- which is how x^3 - x^2 - 3x + 3 is
184+
// answered exactly, and how a quintic is answered at all.
185+
// https://github.com/asc-community/AngouriMath/issues/272
186+
if (Functions.PolynomialFactoring.TrySplitOffRationalRoots(expr, x, out var split)
187+
&& split is Mulf)
188+
return Solve(split, x, compensateSolving);
189+
147190
if (PolynomialSolver.SolveAsPolynomial(expr, x, out var isIdentity) is { } poly)
148191
return poly.Select(e => TryDowncast(expr, x, e.InnerSimplified)).ToSet();
149192
if (isIdentity)

‎Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,62 @@ internal static bool TryFactor(Entity expr, Variable x, [NotNullWhen(true)] out
105105
return degree > 1;
106106
}
107107

108+
/// <summary>
109+
/// <paramref name="expr"/> written as the linear factors of its rational roots
110+
/// times whatever will not divide, or <see langword="false"/> where it has no
111+
/// rational root at all.
112+
/// </summary>
113+
/// <remarks>
114+
/// <para>
115+
/// <see cref="TryFactor"/> answers a different question and answers it more
116+
/// strictly: it offers a factorization only when the polynomial splits completely
117+
/// into whole linear factors, because a partly factored sum is not obviously a
118+
/// better way of *writing* the same thing. Solving is not a question of how to
119+
/// write it. <c>(x - 1)(x^2 - 3)</c> is no nicer to read than
120+
/// <c>x^3 - x^2 - 3x + 3</c>, but it is the difference between answering with
121+
/// <c>1</c> and <c>+-sqrt(3)</c> and answering with two nested cube roots of
122+
/// <c>26 + 18i</c> -- and above degree four, between answering at all and handing
123+
/// the equation to the numeric solver.
124+
/// https://github.com/asc-community/AngouriMath/issues/272
125+
/// </para>
126+
/// <para>
127+
/// Fractional roots are kept here where <see cref="TryFactor"/> declines them, for
128+
/// the same reason: a root of <c>1/2</c> is as exact an answer as any, whatever
129+
/// the factored form looks like on the page.
130+
/// </para>
131+
/// </remarks>
132+
internal static bool TrySplitOffRationalRoots(
133+
Entity expr, Variable x, [NotNullWhen(true)] out Entity? factored)
134+
{
135+
factored = null;
136+
if (!TryGetRationalCoefficients(expr, x, out var coefficients))
137+
return false;
138+
139+
// Two terms is a x^n + b, which is answered whole by inverting x^n = -b/a, and
140+
// that gives the n roots in the form a + bi. Splitting it would divide out the
141+
// rational ones and leave the rest to be dug out of a quotient: x^3 - 8 reads
142+
// as 2 and (-1/2 +- i*sqrt(3)/2)*2 that way, and as 2 and (-2 -+ sqrt(-12))/2
143+
// this way. Nothing here improves on that, so it is left alone.
144+
if (coefficients.Count(coefficient => !coefficient.IsZero) <= 2)
145+
return false;
146+
147+
var factors = new List<Entity>();
148+
foreach (var candidate in RootCandidates(coefficients))
149+
while (coefficients.Length > 1 && Evaluate(coefficients, candidate).IsZero)
150+
{
151+
factors.Add(LinearFactor(candidate, x));
152+
coefficients = DivideByLinear(coefficients, candidate);
153+
}
154+
if (factors.Count == 0)
155+
return false;
156+
157+
var remainder = BuildPolynomial(coefficients, x);
158+
if (remainder != Integer.One)
159+
factors.Add(remainder);
160+
factored = factors.Aggregate((left, right) => left * right);
161+
return true;
162+
}
163+
108164
/// <summary>
109165
/// One step of a partial fraction decomposition, at a rational root of the
110166
/// denominator: <c>N/D</c> becomes <c>A/(x - r) + R/Q</c>, where <c>D = (x - r)Q</c>.
Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
//
2+
// Copyright (c) 2019-2022 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.Linq;
9+
using AngouriMath;
10+
using AngouriMath.Extensions;
11+
using Xunit;
12+
13+
namespace AngouriMath.Tests.Algebra
14+
{
15+
/// <summary>
16+
/// A polynomial equation that has a factor the solver can answer exactly must be
17+
/// answered through that factor, rather than by handing the whole polynomial to the
18+
/// general formula for its degree -- or, above degree four, to the numeric solver.
19+
/// https://github.com/asc-community/AngouriMath/issues/272
20+
/// </summary>
21+
public sealed class FactoredPolynomialSolveTest
22+
{
23+
private static Entity.Set.FiniteSet Roots(string equation) =>
24+
(Entity.Set.FiniteSet)equation.ToEntity().Solve("x");
25+
26+
private static void AssertRootsAre(string equation, params string[] expected)
27+
{
28+
var roots = Roots(equation);
29+
foreach (var root in expected)
30+
Assert.Contains(root.ToEntity().InnerSimplified, roots);
31+
Assert.Equal(expected.Length, roots.Count);
32+
}
33+
34+
/// <summary>
35+
/// The clearest statement of the defect, and the one that does not depend on how a
36+
/// radical is spelled: an equation whose roots are all real must not be answered
37+
/// with expressions built out of complex numbers. Cardano's formula reaches the
38+
/// real roots of x^3 - x^2 - 3x + 3 through (26 + 18i)^(1/3).
39+
/// </summary>
40+
private static void AssertNoImaginaryPartsInTheWorking(string equation)
41+
{
42+
foreach (var root in Roots(equation))
43+
foreach (var node in root.Nodes)
44+
if (node is Entity.Number.Complex and not Entity.Number.Real)
45+
Assert.True(false,
46+
$"{equation} has only real roots, but {root.Stringize()} reaches one " +
47+
$"through {node.Stringize()}");
48+
}
49+
50+
/// <summary>A root that is a decimal is a root the solver gave up on.</summary>
51+
private static void AssertEveryRootIsExact(string equation)
52+
{
53+
foreach (var root in Roots(equation))
54+
foreach (var node in root.Nodes)
55+
if (node is Entity.Number.Real and not Entity.Number.Rational)
56+
Assert.True(false,
57+
$"{equation} was answered with the decimal {node.Stringize()} " +
58+
$"inside {root.Stringize()}");
59+
}
60+
61+
// Handed the factors outright, the solver still expanded them into a cubic and ran
62+
// Cardano over it, so the answer to an equation whose roots are 1 and +-sqrt(3) came
63+
// back as two nested cube roots of 26 + 18i.
64+
[Fact]
65+
public void AProductIsSolvedThroughItsFactors() =>
66+
AssertRootsAre("(x - 1) * (x ^ 2 - 3) = 0", "1", "sqrt(3)", "-sqrt(3)");
67+
68+
[Fact]
69+
public void AProductOfSeveralFactorsIsSolvedThroughThemToo() =>
70+
AssertRootsAre("(x - 1) * (x - 2) * (x + 2) * (x ^ 2 - 3) = 0",
71+
"1", "2", "-2", "sqrt(3)", "-sqrt(3)");
72+
73+
// The same equation written out. A rational root is there to be found, and dividing
74+
// it out leaves a quadratic that is answered exactly.
75+
[Fact]
76+
public void AnExpandedCubicWithARationalRootIsSplitAtIt() =>
77+
AssertRootsAre("x ^ 3 - x ^ 2 - 3 * x + 3 = 0", "1", "sqrt(3)", "-sqrt(3)");
78+
79+
// Above degree four there is no general formula, so the whole equation used to fall
80+
// through to the numeric solver -- which returned decimals, one of them carrying a
81+
// spurious -2.3e-15i.
82+
[Fact]
83+
public void AQuinticThatSplitsIntoSolvableFactorsIsSolvedExactly() =>
84+
AssertRootsAre("x ^ 5 - x ^ 4 - 7 * x ^ 3 + 7 * x ^ 2 + 12 * x - 12 = 0",
85+
"1", "2", "-2", "sqrt(3)", "-sqrt(3)");
86+
87+
[Theory]
88+
[InlineData("(x - 1) * (x ^ 2 - 3) = 0")]
89+
[InlineData("x ^ 3 - x ^ 2 - 3 * x + 3 = 0")]
90+
[InlineData("x ^ 3 - 6 * x ^ 2 + 11 * x - 6 = 0")]
91+
public void RealRootsAreNotReachedThroughComplexNumbers(string equation) =>
92+
AssertNoImaginaryPartsInTheWorking(equation);
93+
94+
[Theory]
95+
[InlineData("x ^ 3 - x ^ 2 - 3 * x + 3 = 0")]
96+
[InlineData("x ^ 5 - x ^ 4 - 7 * x ^ 3 + 7 * x ^ 2 + 12 * x - 12 = 0")]
97+
[InlineData("(x - 1) * (x - 2) * (x + 2) * (x ^ 2 - 3) = 0")]
98+
public void RootsThatCanBeWrittenExactlyAre(string equation) =>
99+
AssertEveryRootIsExact(equation);
100+
101+
private static bool IsNear(Entity root, double re, double im)
102+
{
103+
if (root.Vars.Any())
104+
return false;
105+
var value = root.EvalNumerical();
106+
return System.Math.Abs(value.RealPart.EDecimal.ToDouble() - re) < 1e-9
107+
&& System.Math.Abs(value.ImaginaryPart.EDecimal.ToDouble() - im) < 1e-9;
108+
}
109+
110+
/// <summary>
111+
/// Worse than an unreadable answer: an incomplete one. Above degree four the whole
112+
/// equation went to the numeric solver, which searches for real roots, so
113+
/// (x - 1)(x^2 - 2)(x^2 + x + 1) = 0 came back as { 1, sqrt(2), -sqrt(2) } -- the
114+
/// two roots of its x^2 + x + 1 factor were missing, and a set that is short of two
115+
/// roots does not look any different from one that is not.
116+
/// Splitting off the rational root leaves a quartic, which is solved whole.
117+
/// </summary>
118+
[Theory]
119+
[InlineData("x ^ 5 - 2 * x ^ 3 - x ^ 2 + 2 = 0")]
120+
[InlineData("(x - 1) * (x ^ 2 - 2) * (x ^ 2 + x + 1) = 0")]
121+
public void AQuinticKeepsTheRootsOfItsIrreducibleFactor(string equation)
122+
{
123+
var roots = Roots(equation);
124+
var sqrt2 = System.Math.Sqrt(2);
125+
var halfSqrt3 = System.Math.Sqrt(3) / 2;
126+
foreach (var (re, im) in new[] { (1.0, 0.0), (sqrt2, 0.0), (-sqrt2, 0.0),
127+
(-0.5, halfSqrt3), (-0.5, -halfSqrt3) })
128+
Assert.Contains(roots, root => IsNear(root, re, im));
129+
Assert.Equal(5, roots.Count);
130+
}
131+
132+
// What must not change. A polynomial with no rational root has nothing to split off,
133+
// and its roots are what they always were.
134+
[Fact]
135+
public void APolynomialWithoutARationalRootIsUnaffected() =>
136+
AssertRootsAre("x ^ 2 - 2 = 0", "sqrt(2)", "-sqrt(2)");
137+
138+
// A two-term polynomial is answered whole by inverting x^n = c, which gives its
139+
// roots as a + bi. Splitting off the rational ones instead leaves the others to be
140+
// dug out of the quotient, and x^3 - 8 comes back as (-2 -+ sqrt(-12))/2 rather
141+
// than as (-1/2 +- i*sqrt(3)/2)*2. Nothing is gained, so those are left alone.
142+
[Theory]
143+
[InlineData("x ^ 3 - 8 = 0")]
144+
[InlineData("x ^ 6 - 1 = 0")]
145+
public void ABinomialKeepsTheAnswerInversionGivesIt(string equation)
146+
{
147+
foreach (var root in Roots(equation))
148+
foreach (var node in root.Nodes)
149+
if (node is Entity.Powf(Entity.Number.Real { IsNegative: true } radicand, Entity.Number.Rational))
150+
Assert.True(false,
151+
$"{equation} was answered through a root of the negative " +
152+
$"{radicand.Stringize()}, in {root.Stringize()}");
153+
}
154+
155+
// Repeated factors are one root, not several.
156+
[Fact]
157+
public void ARepeatedFactorContributesItsRootOnce() =>
158+
AssertRootsAre("(x - 1) ^ 2 * (x + 3) = 0", "1", "-3");
159+
160+
// A factor free of the variable has no roots to contribute, and must not make the
161+
// equation look unsolvable either.
162+
[Fact]
163+
public void AConstantFactorContributesNothing() =>
164+
AssertRootsAre("5 * (x ^ 2 - 3) = 0", "sqrt(3)", "-sqrt(3)");
165+
166+
// Solving factor by factor must not smuggle in a value that is not a root of the
167+
// whole equation. Here x = 0 zeroes the first factor but makes the second undefined.
168+
[Fact]
169+
public void AFactorsRootThatTheWholeEquationRejectsIsDropped()
170+
{
171+
var roots = Roots("x * (1 / x) = 0");
172+
Assert.DoesNotContain(Entity.Number.Integer.Create(0), roots);
173+
}
174+
}
175+
}

0 commit comments

Comments
 (0)