# improve_factorization.mpl
#
# By: John May
#
# Created: March 15th, 2005
# Last Revision: July 22nd, 2005

# provides backward_error:
#read("backwards_error.mpl");


# This routine converts an input bivariate polynomial into a vector
 #  Written by John P May
 
 tv2veco := proc(f,ord)
        local i, j, k, tmpcoeff, td, fvec,x,y,z;
        x:=ord[1];
        y:=ord[2];
        z:=ord[3];        
    #    tmpcoeff:=map(PolynomialTools[CoefficientList], map(PolynomialTools[CoefficientList], PolynomialTools[CoefficientList](f,x) , y), z);
        td := degree(f, {x,y,z});
        fvec := [];
        for i from td to 0 by -1 do
            for j from i to 0 by -1 do
                for k from 0 to i - j do 
                # Insert coeff of x^j*z^k*y^(i-j-k)
                    fvec := [op(fvec), coeff(coeff(coeff(f,x,j),z,k),y,i-j-k)];
                end do;    
           end do;
       end do;
       Vector(fvec);
end proc;


# This routine converts the output of bv2vec back into a bivariate polynomial
#  Written by John P May

vec2tvo := proc(v,ord)
        local i, j, k, l, deg, f,x,y,z;
        x:=ord[1];
        y:=ord[2];
        z:=ord[3];
        
        l := LinearAlgebra[Dimension](v);
        f := 0;
        i:=-1;
        for deg from 0 to l do
           for j from 0 to deg do
              for k from 0 to deg-j do            
                        try
                                f := f + v[i]*x^(j)*y^(k)*z^(deg-j-k);
                                i := i - 1;
                        catch:
                                break;
                        end try;
               end do;         
            end do;
        end do;
        f;
end proc;



multiplication_map_jacobian := proc(fs::list(polynom),  # List of polynomials
                vars::list(name)    # the variables of fs as a list
                )

 local 	Fs,			# Symbolic version of fs
 		ns, 		# degrees of fs
 		F, 			# The coeff vector of the symbolic product of the Fs
 		C, 			# The Jacobian of F
 		eqnvars, i,
 		numfacs;

	# This method uses total degree support of the polynomials

	# Computes the jacobian of the map (fi)->prod(f_i) evaluated at the given (f_i)
  
  
	numfacs := nops(fs);#print('numfacs'=numfacs);
	ns := [seq(degree(fs[i],vars),i=1..numfacs)];#print('degrees'=ns);

	# Build symbolic polynomials, so we can take the derivative.
	# This line could be modified to support > 3 vars
	
	Fs := [seq(add(add(add(a[k,j,l,i-j-l]*vars[1]^j*vars[3]^(i-j-l)*vars[2]^l,l=0..i-j),j=0..i),i=0..ns[k]), k=1..numfacs)];


	# Build the list of variables in the right order
	# Need a version of bv2veco for 3 variables!
	eqnvars := convert(Vector([seq(tv2veco(Fs[i],vars),i=1..numfacs)]), list);

	# Build multiplication map in terms of symbolic coefficients

	F := Vector([tv2veco(expand(product(Fs[i],i=1..numfacs)),vars)]);

	C := VectorCalculus[Jacobian](F, eqnvars); 

	# Specialize the Jacobian to the inputs
	# This line could be modifed to support >2 vars
	subs({ seq(seq(seq(seq(a[k,j,l,i-j-l] = coeff(coeff(coeff(fs[k],vars[1],j),vars[2],l),vars[3],i-j-l),l=0..i-j),j=0..i), i=0..ns[k]),k=1..numfacs)},C);

end proc: # multiplication_map_jacobian


improve_factorization := proc(
			f::polynom, 			# Bivariate polynomial being factored
			facts::list(polynom), 	# a list of approximate factors of f
			vars::list(name)		# the variables of f and facs
			)

	local c, n, beo, l, facs, ii, zn, kk, jj;
	
	# Backward error of the input factorization: 
	c,beo:=backward_error(f, convert(facts,`*`),vars);
	userinfo(1, {appfac,imfac}, print('org_error'=beo));
	userinfo(3, {appfac,imfac}, print('c'=c));
	
	n := nops(facts);
	
	# compute the lengths of the coeff. vectors of the factors
	# this is a stupid way to do this - there is a formula to compute this
	l:=[seq(op(1,tv2veco(facts[i], vars)),i=1..n)];
	
	userinfo(4, {appfac,imfac}, print('l' = l));
	
	Digits := floor(evalhf(Digits));
	
	# Process the factors: remove extraneous imaginary parts
	facs:=map(z1->vec2tvo(map(z2->`if`(Im(z2)=0,Re(z2),z2), tv2veco(z1,vars)), vars),facts):
		
	# multiply each factor by part of the optimal scaling computed by backward_error()
	# possibly better ways to do this.
	
	if c < 0 then c := -c; facs[1]:=-1.*facs[1] end if;
	
	facs:=map(z->evalf(c^(1/n)) * z, facs):

	# Gauss-Newton Iteration to refine the factors
	for ii from 1 to 15 do

		userinfo(4, {appfac,imfac}, print(ii,'bwerror'=norm(expand(f-convert(facs,`*`)),2))) ;
		
		# Using using LeastSquares should be more stable, but doesn't seem to work at all.
		zn := LinearAlgebra:-MatrixInverse( multiplication_map_jacobian(facs,vars)) . tv2veco(expand(convert(facs,`*`) - f),vars);
		#zn := LinearAlgebra:-LeastSquares( multiplication_map_jacobian(facs,vars) , tv2veco(expand(convert(facs,`*`) - f),vars));		

		userinfo(4, {appfac,imfac}, print(ii,'zn'=LinearAlgebra:-Norm(zn,2)));

  		kk:=0;  
  		for jj from 1 to n do
			facs[jj]:=facs[jj]-vec2tvo(zn[ (kk+1)..(kk+l[jj]) ], vars);
			kk:=kk+l[jj];
		end do;

		if(LinearAlgebra:-Norm(zn,2) < 5*10^(-12) ) 
			then break end if; 

	end do:
	
	userinfo(1, {appfac,imfac}, print('iters'=ii));

	userinfo(1, {appfac,imfac}, print('new_error'= backward_error(f, convert(facs,`*`), vars)[2], 'improvement' = beo  / backward_error(f, convert(facs,`*`), vars)[2]));

	return facs;

end proc: # improve_factorization