Factor a polynomial in two variables, by Kronecker's substitution

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

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
 

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;
+            }
+
+
+            /// <summary>Whether a constant polynomial is 1, so that it need not be printed.</summary>
+            private static bool SameAsOne(MultivariatePolynomial poly)
+                => poly.IsConstant && poly.CoefficientOf(0).CompareTo(ERational.One) == 0;
+
+            /// <summary>
+            /// The factorisation of a polynomial in exactly two variables, as an expression.
+            /// </summary>
+            /// <remarks>
+            /// Kronecker's substitution: see <see cref="BivariateFactorization"/> for what it
+            /// does, what it refuses, and why a wrong answer is not among the things it can do.
+            /// </remarks>
+            private static Entity? Bivariate(
+                MultivariatePolynomial poly, IReadOnlyList<Variable> variables,
+                IReadOnlyDictionary<Variable, int> 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<Entity>();
+                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;
             }
 
             /// <summary>

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
+{
+    /// <summary>
+    /// Factorisation of a polynomial in two variables, by reducing it to one variable and
+    /// putting the answer back.
+    /// </summary>
+    /// <remarks>
+    /// <para>
+    /// <b>Kronecker's substitution.</b> A factor of a polynomial of degree <c>d</c> in <c>x</c>
+    /// has degree at most <c>d</c> in <c>x</c>, so with <c>s = d + 1</c> the map
+    /// <c>x^i y^j -> t^(i + s*j)</c> is injective on every monomial that can appear in the
+    /// polynomial or in any of its factors: <c>i</c> is the remainder and <c>j</c> the quotient of
+    /// the exponent by <c>s</c>, 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.
+    /// </para>
+    /// <para>
+    /// <b>A candidate is a guess and is checked by division.</b> 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.
+    /// </para>
+    /// <para>
+    /// <b>What it will not do.</b> The image has degree <c>d + s * e</c> for a polynomial of
+    /// degree <c>e</c> in <c>y</c>, and the one-variable factoriser stops at
+    /// <see cref="IntegerPolynomial.MaxDegree"/> — 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 (<a href="https://github.com/asc-community/AngouriMath/issues/746">#746</a> item 43).
+    /// </para>
+    /// </remarks>
+    internal static class BivariateFactorization
+    {
+        /// <summary>
+        /// 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.
+        /// </summary>
+        private const int MaxImageFactors = 12;
+
+        /// <summary>
+        /// The factors of <paramref name="poly"/> in <paramref name="x"/> and <paramref name="y"/>,
+        /// each to the first power and each of positive degree in <paramref name="x"/>, or
+        /// <see langword="null"/> where nothing could be settled. A single factor means the
+        /// polynomial did not factor, which is an answer.
+        /// </summary>
+        internal static IReadOnlyList<MultivariatePolynomial>? 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<IntegerPolynomial>();
+            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);
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        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);
+
+        /// <summary>
+        /// Two variables again, reading each exponent as its remainder and quotient by
+        /// <paramref name="stride"/> — or <see langword="null"/> where a power is past what a
+        /// packed monomial holds.
+        /// </summary>
+        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;
+        }
+
+        /// <summary>
+        /// Every subset of the image's irreducible factors, smallest first, kept where its
+        /// product divides what is left of the polynomial.
+        /// </summary>
+        /// <remarks>
+        /// 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.
+        /// </remarks>
+        private static IReadOnlyList<MultivariatePolynomial>? Recombine(
+            MultivariatePolynomial poly, List<IntegerPolynomial> irreducibles, int x, int y, int stride)
+        {
+            var found = new List<MultivariatePolynomial>();
+            var remaining = poly;
+            var available = new List<IntegerPolynomial>(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;
+        }
+
+        /// <summary>
+        /// Whether the factors multiply back to the polynomial up to a rational constant, which
+        /// is as far as a factorisation is ever fixed.
+        /// </summary>
+        private static bool DividesBackExactly(MultivariatePolynomial product, MultivariatePolynomial poly)
+            => !product.IsZero
+               && poly.DivideExact(product) is { } quotient
+               && quotient.IsConstant;
+
+        /// <summary>The index subsets of a given size, in a fixed order.</summary>
+        private static IEnumerable<List<int>> Subsets(int count, int size)
+        {
+            var chosen = new List<int>(size);
+            return Walk(0);
+
+            IEnumerable<List<int>> Walk(int from)
+            {
+                if (chosen.Count == size)
+                {
+                    yield return new List<int>(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);
+                }
+            }
+        }
+    }
+}