read NearestNpoly;

f1 := 92*x^3+44*x^2*y+40*x*y^2-67*y^3+8*x^2+66*x*y+68*y^2-95*x-62*y-18: 
f2 := 68*x^3+39*x^2*y+20*x*y^2+45*y^3-65*x^2-67*x*y+93*y^2+43*x+8*y+6: 
f := expand(f1*f2)/100;
F := evalf(f,3);

# The distance from F to the reducible polynomial f that we started with:
Distance(F, f);

# Note that starting from the approximate polynomial F, the polynomial f
# is unlikely to be the nearest reducible polynomial.
#
# The tangent space V of reducible polynomials at f is spanned by:
#   x^i*y^j*f1, x^i*y^j*f2,  0 <= i+j <= 3.
# and has dimension 10 + 10 - 1   (the -1 is because f was counted twice).
#
# So reducible polynomials that are close (distance O(epsilon)) to f
# are very close (distance O(epsilon^2)) to the affine space f+V.
#
# So if the distance from F to f is O(epsilon) then the distance from F
# to its nearest reducible polynomial is very close (within O(epsilon^2))
# to the distance from F to f+V.
#
# This way we can (if the distance from F to f is small) give a very
# good estimate for the distance from F to its nearest reducible
# polynomial, we expect this estimate to be off by only O(epsilon^2).

V := {	seq(seq(x^i * y^j * f1, i=0..3-j),j=0..3),
	seq(seq(x^i * y^j * f2, i=0..3-j),j=0..2)}:
V := {seq(CoeffVector(evalf(expand(i)), x,y,6), i=V)}:
v := CoeffVector(F, x,y,6):
v := v-ProjectOnV(v, V):
estimate := sqrt(LinearAlgebra:-DotProduct(v,v));

NF := NearestNpoly(F, 9);
Distance(F, NF);

# The polynomial NF we found is reducible (at least, up to the accuracy
# determined by the setting of Digits) and its distance from F is very
# close to the estimated distance from F to its nearest reducible polynomial.
#
# This means that the code does indeed find "the nearest reducible" polynomial
# in R[x,y],  where "the nearest" should be OK within O(epsilon^2)-tolerance,
# and "reducible" should be OK within  10^(-Digits)-tolerance, where Digits
# can be set arbitrarily by the user, and epsilon is the distance from F
# to its nearest reducible polynomial.