New in Maple 2021 is the PolynomialTools:-Approximate subpackage. It includes three main commands for computing approximate operations on polynomials in many variables: Factor, GCD, and Divide, as well as commands for computing the matrices used for these operations.

These new commands are useful for cases where no exact answer is possible, and a solution is instead found for a nearby problem which does have a solution. Typically, this assumes that some amount of error or noise was introduced into the coefficients of the input polynomial(s) that destroyed algebraic structure, and which these commands attempt to recover.

The SNAP package has provided the gcd and division with remainder operations for univariate polynomials but those commands typically require the user to provide a numerical tolerance.  These new commands do not take a tolerance as input, but instead always return an approximate solution, even if it has very large backward error following the model in the papers of Gao et al. and Kaltofen et al in the references below.

This new package supports univariate computations of quotients and gcd under that model (univariate inputs to Factor are just sent to the top-level command factor).

factor( x^4-3.0*x^2-4.0 );

$\left(x+2.\right)\left(x-2.\right)\left(x^2+1.\right)$

factor( x^4-3.0*x^2-4.0*I );

$\left(x+1.85321006568802+0.291797565955056\mathrm{I}\right)\left(x+0.627404403732890-0.861903714978381\mathrm{I}\right)\left(x-0.627404403732890+0.861903714978381\mathrm{I}\right)\left(x-1.85321006568802-0.291797565955056\mathrm{I}\right)$

SNAP:-QRGCD( x^4-3.0*x^2-4.0, x^2-4.0, x, 10e-10);

$0,0.2425356250,-0.970142500145332+0.242535625036333x^2$

SNAP:-Quotient( x^4-3.0*x^2-4.0, x^2-4.0, x );

$1.00000000000000x^2+1.00000000000000$

New functionality is described below.

with(PolynomialTools):

with(Approximate);

$\left[\mathrm{ConvolutionMatrix}\,\mathrm{Divide}\,\mathrm{Factor}\,\mathrm{GCD}\,\mathrm{RuppertMatrix}\,\mathrm{SylvesterMatrix}\right]$

The commands CoefficientVector and FromCoefficientVector in the main PolynomialTools package were both enhanced to work with multivariate polynomials.

v := CoefficientVector(y^3-y^2*x+(y-2)*x^2+x^3, [x,y]);

$v≔\left[\begin{array}{c}0\\ 0\\ 0\\ \mathrm{-2}\\ 0\\ 0\\ 1\\ 1\\ \mathrm{-1}\\ 1\end{array}\right]$

FromCoefficientVector( v, [a, b]);

$a^3+a^{2}b-a{b}^2+b^3-2a^2$

Approximate division using a LeastSquares computation on the matrix representing polynomial division.

f := x^2+y^2-1;

$f≔x^2+y^2-1$

ConvolutionMatrix(f, [x,y], 2);

infolevel[PolynomialTools] := 1;

${\mathrm{infolevel}}_{\mathrm{PolynomialTools}}≔1$

d := sqrt(2.0)*(x^2+y^2-1);

$d≔1.414213562x^2+1.414213562y^2-1.414213562$

G := sort( (x^2-y^3+1), [x,y]); F := sort( expand( d*G ), [x,y]);

$G≔-y^3+x^2+1$

$F≔-1.414213562x^2{y}^3-1.414213562y^5+1.414213562x^4+1.414213562x^2{y}^2+1.414213562y^3+1.414213562y^2-1.414213562$

ad_8 := Divide( expand( F + 10^(-8)*x*y ), expand(G+10^(-8)*x^2*y), [x,y] );

Divide:   computed approximate quotient has backward error 2.291288e-08
Divide:   computed approximate quotient has backward error 2.291288e-08

$\mathrm{ad\_8}≔-1.41421356376777+1.72890601399040\times {10}^{\mathrm{-16}}x+4.71404493218586\times {10}^{\mathrm{-9}}y+1.41421356730330x^2+3.33333340088403\times {10}^{\mathrm{-9}}xy+1.41421356200000y^2$

ad_8 := sort( tcoeff(d)*fnormal(expand(ad_8)/(tcoeff(ad_8))), [x,y]);

$\mathrm{ad\_8}≔1.414213566x^2+1.414213560y^2-1.414213562$

ilog10(norm(expand(d-ad_8),2) / norm(d,2));

$\mathrm{-9}$

Approximate GCD computation using a low-rank approximation of a generalized Sylvester matrix.

d := sqrt(2.0)*(x^2+y^2-1);

$d≔1.414213562x^2+1.414213562y^2-1.414213562$

F := sort( expand( d*(x^3-y^3+1) ), [x,y]);

$F≔1.414213562x^5+1.414213562x^3{y}^2-1.414213562x^2{y}^3-1.414213562y^5-1.414213562x^3+1.414213562y^3+1.414213562x^2+1.414213562y^2-1.414213562$

G := sort( expand( d*(x^2-y^3+1) ), [x,y]);

$G≔-1.414213562x^2{y}^3-1.414213562y^5+1.414213562x^4+1.414213562x^2{y}^2+1.414213562y^3+1.414213562y^2-1.414213562$

S1 := SylvesterMatrix(F, G, [x,y], 3);

Maximal rank, means degree( gcd(F,G) ) < 3

min(upperbound(S1)) - LinearAlgebra:-Rank(S1);

$0$

ad_8 := GCD( expand( F + 10^(-8)*x*y ), expand(G+10^(-8)*x^2*y), [x,y] );

GCD:   approximate GCD inputs of total degrees 5 and 5 in [x, y]
GCD:   approximate GCD determined with an error of 7.484964e-10

$\mathrm{ad\_8}≔0.569494797046587+3.40008238745826\times {10}^{\mathrm{-10}}x-1.02787949443617\times {10}^{\mathrm{-9}}y-0.569494797402477x^2-1.49456632732006\times {10}^{\mathrm{-10}}xy-0.569494797883116y^2$

ad_8 := sort( tcoeff(d)*fnormal(expand(ad_8)/(tcoeff(ad_8))), [x,y]);

$\mathrm{ad\_8}≔1.414213563x^2+1.414213563y^2-1.414213562$

ilog10(norm(expand(d-ad_8),2) / norm(d,2));

$\mathrm{-10}$

Approximate factorization computation using a low-rank approximation of a Ruppert matrix.

F := sort( expand( (x^2+y^2-1)*(x^3-y^3+1) ), [x,y]);

$F≔x^5+x^3{y}^2-x^2{y}^3-y^5-x^3+y^3+x^2+y^2-1$

R1 := RuppertMatrix( F, [x,y]);

The rank deficiency of the Ruppert matrix is equal to the number of factors.

min(upperbound(R1)) - LinearAlgebra:-Rank(R1);

$2$

aF_8 := Factor( expand( F + 10^(-8)*x*y ), [x,y] );

Factor:   approximate factorization input of degree [5, 5] in [x, y]
Factor:   factorization square-free input of degree [5, 5] in [x, y]
Factor:   polynomial has 2 approximate factors at numerical tolerance 10e-9

$\mathrm{aF\_8}≔-3.51077852332865\left(-0.553434072558778+1.77801493208625\times {10}^{\mathrm{-9}}x+7.90820151946415\times {10}^{\mathrm{-10}}y+0.553434076452427x^2+3.77682411351434\times {10}^{\mathrm{-9}}xy+0.553434072372129y^2\right)\left(-0.514672153977448-2.34835368967946\times {10}^{\mathrm{-9}}x-5.64898113877529\times {10}^{\mathrm{-9}}y-1.17024429658888\times {10}^{\mathrm{-9}}{x}^2-1.49078210800354\times {10}^{\mathrm{-9}}xy+2.96989995359614\times {10}^{\mathrm{-9}}{y}^2-0.514672156404694x^3-5.26729024754203\times {10}^{\mathrm{-10}}{x}^{2}y-2.25599642616380\times {10}^{\mathrm{-9}}{y}^{2}x+0.514672152857353y^3\right)$

sort( fnormal( expand(aF_8) ), [x,y] );

$1.000000007x^5+1.000000004x^3{y}^2-0.9999999990x^2{y}^3-0.9999999926y^5-0.9999999953x^3+1.000000004y^3+0.9999999999x^2+1.000000001y^2-0.9999999951$

ilog10( norm(expand(F-aF_8),2) / norm(F,2) );

$\mathrm{-9}$

Gao, S.; Kaltofen, E.; May, J.; Yang, Z.; and Zhi, L. "Approximate factorization of multivariate polynomials via differential equations." Proceedings of the 2004 International Symposium on Symbolic and Algebraic Computation (ISSAC 2004),  pp. 167-174. Ed. J. Gutierrez. ACM Press, 2004.

Kaltofen, E.; May, J.; Yang, Z.; and Zhi, L. "Approximate factorization of multivariate polynomials using singular value decomposition." Journal of Symbolic Computation Vol. 43(5), (2008): 359-376.