In[]:=
f=2*w^3*x*z+3*x^2*y*z^2-2*w*x*z+6*y*z-8*w+10;g=-w^4*x^3*y^2*z+8*x*y^2*z-6*w*x+3*x-1h=f^2+g^2;k=Expand[h]Factor[k,ExtensionI](*PolynomialSumOfSquaresList[k,{x,y}]//EchoTiming*)(*Total[%^2//Factor]*)
Out[]=
-1+3x-6wx+8xz-z
2
y
4
w
3
x
2
y
Out[]=
101-160w+64-6x+12wx+9-36w+36-40wxz+32xz+40xz-32xz+120yz-96wyz-16xz+48z-96wz+2z-6z+12z+4-8+4-24wxy+24xy+60y-48wy+36+64-16+-12wy+12y+36+9
2
w
2
x
2
x
2
w
2
x
2
w
3
w
4
w
2
y
2
x
2
y
2
x
2
y
4
w
3
x
2
y
4
w
4
x
2
y
5
w
4
x
2
y
2
w
2
x
2
z
4
w
2
x
2
z
6
w
2
x
2
z
2
z
3
w
2
z
2
x
2
z
2
x
2
z
2
y
2
z
2
x
4
y
2
z
4
w
4
x
4
y
2
z
8
w
6
x
4
y
2
z
3
x
3
z
3
w
3
x
3
z
2
x
2
y
3
z
4
x
2
y
4
z
Out[]=
((1-10)+8w-3x+6wx+2wxz-2xz-6yz-8xz+z-3y)((1+10)-8w-3x+6wx-2wxz+2xz+6yz-8xz+z+3y)
3
w
2
y
4
w
3
x
2
y
2
x
2
z
3
w
2
y
4
w
3
x
2
y
2
x
2
z
In[]:=