diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index b815876e1..022dc7d2d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -47,6 +47,7 @@ read first. | | `"x + 1 // done".ToEntity()`, and any input whose last line ends in a `//` comment | `UnhandledParseException: extraneous input '/'` | `x + 1` — the comment is skipped, as the block form already was | | | `MathS.Polynomials.Factor("x * y + y", "x")`, and any polynomial whose coefficients in the named variable share a common divisor | `null` — a refusal | `y * (x + 1)` | | | `MathS.Polynomials.SquareFreePart("(x - y) ^ 2 * (x + y)", "x")`, and any polynomial in more than one variable | `null` — a refusal | `x ^ 2 - y ^ 2` | +| | `MathS.Polynomials.Factor("x ^ 2 - y ^ 2", "x")`, and polynomials in two variables of small enough bidegree | `null` — a refusal | `(x + y) * (x - y)` | ### An equation nothing settled is no longer answered with the empty set @@ -304,6 +305,40 @@ bracketing for. The change only ever adds `\left(`/`\right)` groups, which CShar parses, so nothing downstream needs a matching change ([#822](https://github.com/asc-community/AngouriMath/issues/822)). +### `Factor` factors a polynomial in two variables + +After the content is taken out, what remains may still have polynomial coefficients — and where it +is in two variables it can be factored anyway, by **Kronecker's substitution**. A factor of a +polynomial of degree `d` in `x` has degree at most `d` in `x`, so with `s = d + 1` the map +`x^i y^j → t^(i + s*j)` is injective on every monomial that can appear in the polynomial or in any +of its factors. The one-variable image is factored by the existing factoriser, and each subset of +its irreducible factors names a candidate. + +| | 2.3.0 | now | +|---|---|---| +| `Factor("x ^ 2 - y ^ 2", "x")` | `null` | `(x + y) * (x - y)` | +| `Factor("x ^ 2 + 2 * x * y + y ^ 2", "x")` | `null` | `(x + y) ^ 2` | +| `Factor("x ^ 3 - y ^ 3", "x")` | `null` | `(x - y) * (x ^ 2 + x * y + y ^ 2)` | +| `Factor("x ^ 4 - y ^ 4", "x")` | `null` | `(x + y) * (x ^ 2 + y ^ 2) * (x - y)` | +| `Factor("x ^ 2 * y ^ 2 - 1", "x")` | `null` | `(x * y + 1) * (x * y - 1)` | +| `Factor("x ^ 2 - y ^ 2 + 2 * x + 1", "x")` | `null` | `(x + y + 1) * (x - y + 1)` | +| `Factor("x ^ 2 + y ^ 2", "x")` | `null` | `null` — irreducible over ℚ | +| `Factor("x * y + z", "x")` | `null` | `null` — three variables | + +**It cannot answer wrongly.** The substitution is injective on monomials but not on factorisations, +so the image may factor further than the polynomial does and a candidate is a guess. Every one is +tested by exact division before it is kept, and the assembled factors are divided back into the +input, so the failure mode is a refusal. + +**What it refuses.** The image has degree `d + s*e` for degree `e` in the second variable, and the +one-variable factoriser stops at 32 — so this reaches bidegrees like (2, 10), (3, 7) and (5, 4) and +refuses past them. The recombination is over subsets, so the image's factor count is capped too. +Lifting that ceiling is Hensel lifting with an evaluation homomorphism, which is a different piece +of work. + +`MathS.Polynomials.Factor` has no caller inside the library, so no simplification, solution or +integral changes with it. + ### The square-free part is taken where the coefficients are polynomials `MathS.Polynomials.SquareFreePart` refused every polynomial in more than one variable, for the same diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 73fa5efed..10473aee5 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -27,6 +27,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Core/Transformations/*.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[AngouriMath/Functions/Algebra/Polynomials/BivariateFactorization.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [AngouriMath/Functions/Algebra/Groebner/*.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Convenience/MathS.Polynomials.cs b/Sources/AngouriMath/Convenience/MathS.Polynomials.cs index 4bdc514d2..daa0af086 100644 --- a/Sources/AngouriMath/Convenience/MathS.Polynomials.cs +++ b/Sources/AngouriMath/Convenience/MathS.Polynomials.cs @@ -152,15 +152,72 @@ public static class Polynomials for (var i = 0; i < variables.Count; i++) if (i != main) others.Add(i); - if (PolynomialGcd.ContentIn(poly, main, others, 0) is not { } content - || content.IsConstant) + if (PolynomialGcd.ContentIn(poly, main, others, 0) is not { } content) return null; if (poly.DivideExact(content) is not { } primitive) return null; + + // What is left may still have polynomial coefficients, and where it is in two + // variables it can be factored anyway -- see BivariateFactorization. var rest = Assemble( PolynomialFactorization.FactorComplete(primitive.ToEntity(variables), variable), - variable); - return rest is null ? null : content.ToEntity(variables) * rest; + variable) + ?? Bivariate(primitive, variables, index, variable); + if (rest is null) + return null; + return content.IsConstant && content.DivideExact(content) is not null + && SameAsOne(content) + ? rest + : content.ToEntity(variables) * rest; + } + + + /// Whether a constant polynomial is 1, so that it need not be printed. + private static bool SameAsOne(MultivariatePolynomial poly) + => poly.IsConstant && poly.CoefficientOf(0).CompareTo(ERational.One) == 0; + + /// + /// The factorisation of a polynomial in exactly two variables, as an expression. + /// + /// + /// Kronecker's substitution: see for what it + /// does, what it refuses, and why a wrong answer is not among the things it can do. + /// + private static Entity? Bivariate( + MultivariatePolynomial poly, IReadOnlyList variables, + IReadOnlyDictionary index, Variable variable) + { + if (variables.Count != 2) + return null; + var main = index[variable]; + var other = main == 0 ? 1 : 0; + if (BivariateFactorization.Factor(poly, main, other) is not { } factors + || factors.Count < 2) + return null; + // Repeated factors are collected into a power, as the one-variable path does: + // the recombination finds a square as the same factor twice, and printing it + // twice would be a different answer to the same question depending on which + // path answered it. + Entity? product = null; + var pieces = new List(); + foreach (var factor in factors) + pieces.Add(factor.ToEntity(variables)); + var taken = new bool[pieces.Count]; + for (var i = 0; i < pieces.Count; i++) + { + if (taken[i]) + continue; + var multiplicity = 1; + for (var j = i + 1; j < pieces.Count; j++) + if (!taken[j] && pieces[i] == pieces[j]) + { + taken[j] = true; + multiplicity++; + } + var piece = multiplicity > 1 ? pieces[i].Pow(multiplicity) : pieces[i]; + product = product is null ? piece : product * piece; + } + return product; } /// diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/BivariateFactorization.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/BivariateFactorization.cs new file mode 100644 index 000000000..dc8b8bd48 --- /dev/null +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/BivariateFactorization.cs @@ -0,0 +1,245 @@ +// +// 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 PeterO.Numbers; +using System.Collections.Generic; + +namespace AngouriMath.Functions +{ + /// + /// Factorisation of a polynomial in two variables, by reducing it to one variable and + /// putting the answer back. + /// + /// + /// + /// Kronecker's substitution. A factor of a polynomial of degree d in x + /// has degree at most d in x, so with s = d + 1 the map + /// x^i y^j -> t^(i + s*j) is injective on every monomial that can appear in the + /// polynomial or in any of its factors: i is the remainder and j the quotient of + /// the exponent by s, and neither can be confused with another pair. So a factorisation + /// of the one-variable image can be read back, and each subset of its irreducible factors + /// names a candidate. + /// + /// + /// A candidate is a guess and is checked by division. The image can factor further than + /// the polynomial does — the substitution is injective on monomials, not on factorisations — + /// so a subset whose product reads back as a polynomial need not divide the original. Every + /// one is tested with exact division before it is kept, which is why this cannot answer + /// wrongly: the worst it does is fail to find a factorisation that exists and say so. + /// + /// + /// What it will not do. The image has degree d + s * e for a polynomial of + /// degree e in y, and the one-variable factoriser stops at + /// — so this reaches bidegrees like (2, 10), + /// (3, 7) and (5, 4) and refuses beyond them. The recombination is over subsets, so the + /// number of irreducible factors of the image is capped as well. Both limits are refusals, + /// never wrong answers, and neither is the algorithm one would write to lift this ceiling: + /// that is Hensel lifting with an evaluation homomorphism, and it is a different piece of + /// work (#746 item 43). + /// + /// + internal static class BivariateFactorization + { + /// + /// Beyond this the subset search is refused rather than paid for: the recombination is + /// over every subset of the image's irreducible factors, and a polynomial this reducible + /// is not what the substitution is for. + /// + private const int MaxImageFactors = 12; + + /// + /// The factors of in and , + /// each to the first power and each of positive degree in , or + /// where nothing could be settled. A single factor means the + /// polynomial did not factor, which is an answer. + /// + internal static IReadOnlyList? Factor( + MultivariatePolynomial poly, int x, int y) + { + var degreeInX = poly.DegreeIn(x); + var degreeInY = poly.DegreeIn(y); + if (poly.IsZero || degreeInX < 1 || degreeInY < 1) + return null; + + var stride = degreeInX + 1; + var imageDegree = degreeInX + stride * degreeInY; + if (imageDegree > IntegerPolynomial.MaxDegree) + return null; + + if (ToImage(poly, x, y, stride, imageDegree) is not { } image) + return null; + if (PolynomialFactorization.FactorPrimitive(image.PrimitivePart()) is not { } parts) + return null; + + // The multiplicities are flattened: a square in the image may or may not be a square + // in two variables, and the recombination settles that by division rather than by + // carrying the exponent across the substitution. + var irreducibles = new List(); + foreach (var part in parts) + for (var i = 0; i < part.Multiplicity; i++) + { + if (irreducibles.Count == MaxImageFactors) + return null; + irreducibles.Add(part.Factor); + } + if (irreducibles.Count < 2) + return new[] { poly }; + + return Recombine(poly, irreducibles, x, y, stride); + } + + /// + /// The one-variable image, with the denominators cleared — the constant they came to is + /// not wanted, since a rational multiple of a factor is the same factor. + /// + private static IntegerPolynomial? ToImage( + MultivariatePolynomial poly, int x, int y, int stride, int imageDegree) + { + var coefficients = new ERational[imageDegree + 1]; + for (var i = 0; i < coefficients.Length; i++) + coefficients[i] = ERational.Zero; + + foreach (var byX in poly.CoefficientsIn(x)) + foreach (var byY in byX.Value.CoefficientsIn(y)) + { + // Anything left is a third variable, and this is the two-variable case. + if (!byY.Value.IsConstant) + return null; + var at = byX.Key + stride * byY.Key; + if (at > imageDegree) + return null; + coefficients[at] = coefficients[at].Add(byY.Value.CoefficientOf(0)); + } + + var denominator = EInteger.One; + foreach (var coefficient in coefficients) + denominator = Lcm(denominator, coefficient.Denominator); + var whole = new EInteger[coefficients.Length]; + for (var i = 0; i < whole.Length; i++) + whole[i] = coefficients[i].Numerator + .Multiply(denominator.Divide(coefficients[i].Denominator)); + return IntegerPolynomial.Create(whole); + } + + private static EInteger Lcm(EInteger left, EInteger right) + => left.Divide(left.Gcd(right)).Multiply(right); + + /// + /// Two variables again, reading each exponent as its remainder and quotient by + /// — or where a power is past what a + /// packed monomial holds. + /// + private static MultivariatePolynomial? FromImage( + IntegerPolynomial image, int variableCount, int x, int y, int stride) + { + var result = MultivariatePolynomial.Zero(variableCount); + for (var power = 0; power <= image.Degree; power++) + { + if (image[power].IsZero) + continue; + var inX = power % stride; + var inY = power / stride; + if (inX > MultivariatePolynomial.MaxDegree || inY > MultivariatePolynomial.MaxDegree) + return null; + if (MultivariatePolynomial.Monomial(variableCount, x).Power(inX) is not { } partX + || MultivariatePolynomial.Monomial(variableCount, y).Power(inY) is not { } partY + || partX.Multiply(partY) is not { } monomial) + return null; + result = result.Add(monomial.ScaleBy(ERational.Create(image[power], EInteger.One))); + } + return result.IsZero ? null : result; + } + + /// + /// Every subset of the image's irreducible factors, smallest first, kept where its + /// product divides what is left of the polynomial. + /// + /// + /// Smallest first so that what is taken out is irreducible: a subset that divides and + /// whose proper subsets do not is a factor with nothing inside it. The loop restarts + /// after each success because the remaining polynomial has changed. + /// + private static IReadOnlyList? Recombine( + MultivariatePolynomial poly, List irreducibles, int x, int y, int stride) + { + var found = new List(); + var remaining = poly; + var available = new List(irreducibles); + + var progress = true; + while (progress && available.Count > 0) + { + progress = false; + for (var size = 1; size <= available.Count / 2 && !progress; size++) + foreach (var subset in Subsets(available.Count, size)) + { + var product = IntegerPolynomial.One; + foreach (var index in subset) + if (product.Multiply(available[index]) is { } multiplied) + product = multiplied; + else + return null; + if (FromImage(product, poly.VariableCount, x, y, stride) is not { } candidate + || candidate.DegreeIn(x) < 1 + || remaining.DivideExact(candidate) is not { } quotient) + continue; + found.Add(candidate.Normalized()); + remaining = quotient; + for (var i = subset.Count - 1; i >= 0; i--) + available.RemoveAt(subset[i]); + progress = true; + break; + } + } + + if (found.Count == 0) + return new[] { poly }; + if (!remaining.IsConstant) + found.Add(remaining.Normalized()); + // Nothing is returned that does not multiply back to what was asked about. + var check = MultivariatePolynomial.One(poly.VariableCount); + foreach (var factor in found) + if (check.Multiply(factor) is { } multiplied) + check = multiplied; + else + return null; + return DividesBackExactly(check, poly) ? found : null; + } + + /// + /// Whether the factors multiply back to the polynomial up to a rational constant, which + /// is as far as a factorisation is ever fixed. + /// + private static bool DividesBackExactly(MultivariatePolynomial product, MultivariatePolynomial poly) + => !product.IsZero + && poly.DivideExact(product) is { } quotient + && quotient.IsConstant; + + /// The index subsets of a given size, in a fixed order. + private static IEnumerable> Subsets(int count, int size) + { + var chosen = new List(size); + return Walk(0); + + IEnumerable> Walk(int from) + { + if (chosen.Count == size) + { + yield return new List(chosen); + yield break; + } + for (var index = from; index < count; index++) + { + chosen.Add(index); + foreach (var subset in Walk(index + 1)) + yield return subset; + chosen.RemoveAt(chosen.Count - 1); + } + } + } + } +} diff --git a/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialSurfaceTest.cs b/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialSurfaceTest.cs index eada3b0e2..ac7d2f853 100644 --- a/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialSurfaceTest.cs +++ b/Sources/Tests/UnitTests/Algebra/Polynomials/PolynomialSurfaceTest.cs @@ -156,6 +156,57 @@ public void ASquareFreePartHasNoRepeatedRoot(string input) Assert.Equal(Integer.Create(0), part!.Substitute("x", root).Simplify()); } + /// + /// A polynomial in two variables is factored by Kronecker's substitution: with + /// s one more than its degree in x, the map x^i y^j -> t^(i + s*j) + /// is injective on every monomial that can appear in it or in any of its factors, so a + /// factorisation of the image reads back and each subset of its irreducible factors + /// names a candidate. Every candidate is checked by exact division, so the failure mode + /// is a refusal and not a wrong answer. + /// + /// + /// x ^ 2 - y ^ 2 is the case #746 item 43 names, and the one this test exists for. + /// Compared numerically, for the reason the content test above gives. + /// + [Theory] + [InlineData("x ^ 2 - y ^ 2", 2)] + [InlineData("x ^ 2 + 2 * x * y + y ^ 2", 1)] + [InlineData("x ^ 3 - y ^ 3", 2)] + [InlineData("x ^ 2 * y ^ 2 - 1", 2)] + [InlineData("x ^ 4 - y ^ 4", 3)] + [InlineData("x ^ 2 - y ^ 2 + 2 * x + 1", 2)] + public void APolynomialInTwoVariablesIsFactored(string input, int distinctFactors) + { + var expr = input.ToEntity(); + var factored = MathS.Polynomials.Factor(expr, "x"); + Assert.NotNull(factored); + Assert.NotEqual(expr, factored); + + // As many distinct factors as the mathematics has, so a partial factorisation + // reported as a whole one fails rather than passing quietly. + Assert.Equal(distinctFactors, CountFactors(factored!)); + + var variables = expr.Vars.Concat(factored!.Vars).Distinct().ToArray(); + var random = new Random(20260825); + for (var trial = 0; trial < 20; trial++) + { + Entity before = expr, after = factored; + foreach (var variable in variables) + { + Entity value = Math.Round(random.NextDouble() * 6 - 3, 4); + before = before.Substitute(variable, value); + after = after.Substitute(variable, value); + } + Assert.Equal( + before.EvalNumerical().RealPart.EDecimal.ToDouble(), + after.EvalNumerical().RealPart.EDecimal.ToDouble(), + 9); + } + } + + private static int CountFactors(Entity product) + => product is Mulf(var left, var right) ? CountFactors(left) + CountFactors(right) : 1; + /// /// The square-free part is p / gcd(p, dp/dx) whatever ring the coefficients live /// in, so it is not univariate for any reason but the representation it used to be @@ -216,10 +267,9 @@ public void ASquareFreePartOutsideTheLayerIsRefused(string input) /// not exist, which is a wrong answer and not a graceful failure. /// [Theory] - [InlineData("x ^ 2 * y ^ 2 - 1")] // multivariate, and the content is 1 - [InlineData("x ^ 2 - y ^ 2")] // needs factorisation over Q(y), which this is not + [InlineData("x ^ 2 + y ^ 2")] // irreducible over Q in both variables [InlineData("x ^ 2 - a")] // a symbolic coefficient is not rational - [InlineData("x * y + z")] // the coefficients are coprime + [InlineData("x * y + z")] // three variables: the substitution is over two [InlineData("sin(x) + 1")] // not a polynomial [InlineData("x ^ 33 - 1")] // past the degree bound of the factoriser public void FactorisationRefusesRatherThanReturningTheInput(string input)