#------------------------------------------------------------------------------------------------------------------------- # module : multifac.mpl # author : Lihong Zhi Zhengfeng Yang # date : April 6, 2004. #-------------------------------------------------------------------------------------------------------------------------- # This routine compute the factorization of squarefree or nonsquarefree multivariate polynomials multifac:=proc(f,ord) local fx,e,ff,FF,cofa,e0,ddd,RES,vars,tt; # vars:=mainvariable(f,ord); vars:=ord; fx:=diff(f,vars[1]); if degree(f)<= 1 then RETURN([[f,1]]);fi; if nargs=3 then e:=args[3]; tt:=testsquarefree(evalf(expand(f)),ord,e); if tt=0 then lprint(`the polynomial is squarefree`); FF:=appfac(f,vars); FF:=makelist(FF); RETURN(FF); fi; ff:=multigcd(f,fx,vars,e); userinfo(1,multifac,print('ff[1]'=ff[1])); if ff[1]<>0 then cofa:=mcofactor(f,ff[2],vars,e); FF:=appfac(cofa,vars); RES:=multiplicityoffactor(ff[2],FF,vars,e); RES:=orderlist(RES); RETURN(RES); else FF:=appfac(f,vars); FF:=makelist(FF); RETURN(FF); fi; else tt:=testsquarefree(evalf(expand(f)),ord); if tt=0 then lprint(`the polynomial is squarefree`); FF:=appfac(f,vars); FF:=makelist(FF); RETURN(FF); fi; (e,ddd):=epsilonsqrf(f,vars); #print(`e=`,e); userinfo(1,multifac,print('e'=e)); ff:=multigcd(f,fx,vars,e); userinfo(1,multifac,print('ff[1]'=ff[1])); if ff[1]<>0 then cofa:=mcofactor(f,ff[2],vars,e); FF:=appfac(cofa,vars); RES:=multiplicityoffactor(ff[2],FF,vars,e); RES:=orderlist(RES); RETURN(RES); else FF:=appfac(f,vars); FF:=makelist(FF); RETURN(FF); fi; fi; end proc; # This routine compute the factorization of nonsquarefree polynomials. # Wenshin's interpolation routine was used to compute the approximate GCD. multifac1:=proc(f,ord,e) local i,f0,fx,gcd1,LL,pp,SS,cofactor,FF,RES; f0:=evalf(expand(f)/norm(expand(f),infinity)); fx:=diff(f0,ord[1]); Digits:=2*Digits; gcd1:=multigcd1(f,fx,ord,e); #print(`gcd=`,gcd1); gcd1:=cutoff1(gcd1,ord,e): #print(`gcd=`,gcd1); Digits:=Digits/2; LL:=multifac(gcd1,ord); pp:=1;SS:=[]; for i from 1 to nops(L) do pp:=pp*LL[i][1]^(LL[i][2]+1); SS:=[op(SS),[LL[i][1],LL[i][2]+1]]; end do; cofactor:=mcofactor(f0,pp,ord,e); FF:=appfac(cofactor,ord); FF:=makelist(FF); RES:=[op(FF),op(SS)]; RETURN(RES); end proc; # This routine compute the factorization of an approximate squarefree multivariate polynomials. appfac:=proc(f,ord ) # f is a polynomial which we want to factor # ord is bivariate variables local vars,VARS # the variables of the polynomials ,g,h # the unknown polynomials ,aa,bb # coefficients of the polynomials g and h ,f0 # the normalized polynomial ||f0||=1 ,i,j,j1,k # the parameter in for loop statement ,j2,kk,l,jj # the parameter in for loop statement ,temp,ti # the parameter in while loop statement ,E,EQ,eqns # the list of equations ,PP # the list of equations ,e0 # integral,tolerance level 10^(-e0) ,fx # derivate for x variable of the polynomials ,fy # derivate for y variable of the polynomial ,tf # the total degree of f0 ,pp # Ruppert's equation ,linsys # the coefficients of linear system. ,var # the variables in the Ruppert's equation except x and y ,lineqs # the linear systems. ,A # the matrix given by the linear systems. ,b # the unknown variables ,SS,VVt # SVD of A, A= U.SS.VVt ,rn # random vector ,nd # the number of the factors ,coefs # coefficients of the equation ,dimR # row dimension of the matrix ,dimC # column dimension of the matrix ,ddd # list converted by a vector ,sv # the singular values at the biggest gap ,addv # the zero vector ,rats # the list of gaps for singular values ,tq # quotient of two polynomials ,GG # list of some polynomials ,vv # some row vector in VVt ,gg # the polynomial generated by rn and the polynomials in GG ,vf # random number satisfied some condition ,fx1 # the univariate polynomial by fx ,f1 # the univariate polynomial by f0 ,tt # the remainder of two polynomials ,c # the vector of unknowns ,rp,tg,ffg,ffg2 # the polynomials ,L1 # list of the distance of eigenvalues ,cd,cd0 # the distance of eigenvalues ,rq0,rq # the eigenvalues of the matrix ,M1,M # the matrix ,ff # the polynomial which is factor , t # sets of the polynomials , v # sets of the least square solution ,S1 # the vector by degree of every variable. ,Q,P,V,time1 ,FF1, FF2,FF3,FF4 # factors of f ; time1:=time(); Digits:=10; vars:=mainvariable(f,ord); userinfo(1,appfac,print('vars'=vars)); f0:=evalf(f/norm(expand(f),2)); g:=0; tf:=degree(f0); if tf<=1 then RETURN([f]); fi; S1:=[]; for j from 1 to nops(ord) do S1:=[op(S1),degree(f0,vars[j])]; end do; #construct polynomial g[i] satisfies degree(g[i],ord[i])<=degree(f,ord[i])-1 and deg(g[i]) < deg(f); eqns:={}; for i from 1 to nops(vars) do Q[i]:=subtvec(S1,i); P[i]:=vecc(Q[i],tf); g[i]:=0; V:=Vector(nops(P[i])); for jj from 1 to LinearAlgebra:-Dimension(V) do V[jj]:=a[i,jj]; end do; pp:=convert(V,list); for kk from 1 to nops(pp) do g[i]:=g[i]+pp[kk]*vecmo(P[i][kk],vars); end do; eqns:=[op(eqns),g[i]]; end do; EQ:={};VARS:={}; for i from 2 to nops(vars) do PP[i]:=evalf(expand(g[i]*diff(f0,vars[1])-f0*diff(g[i],vars[1])+f0*diff(g[1],vars[i])-g[1]*diff(f0,vars[i]))); EQ:= EQ union {coeffs(PP[i],ord)}; VARS:= VARS union {op(indets(PP[i]))} minus {op(ord)}; end do; var:=[op(VARS)]; # construct the Ruppert matrix A lineqs:=['EQ[i]=0' $ 'i'=1..nops(EQ)]; A,b:=LinearAlgebra[GenerateMatrix](lineqs,var); userinfo(3,appfac,print('A'=A)); SS, VVt := LinearAlgebra[SingularValues](A, output=['S','Vt']): dimR:= LinearAlgebra[RowDimension](VVt); ddd:=convert(SS,list); nd:=1; dimC:=LinearAlgebra[ColumnDimension](A); # to ensure dimension(ddd)=dimC if nops(ddd)< dimC then addv:=convert(Vector(dimC-nops(ddd)),list); ddd:=[op(ddd),op(addv)]; fi; # compute the biggest gap in the singular values rats := seq(ddd[i]/ddd[i+1], i=(dimC-degree(f))..dimC-2); member(max(rats), [rats], 'nd'); userinfo(1,appfac,print(`max(rats)`, max(rats))); if max(rats) <2.0 then print(`the gap is too small, we think the polynomial is irreducible.`); return [f]; fi; nd := nops([rats])+2 - nd; userinfo(1,appfac,print('nd'=nd)); if nd=1 then print(`the gap is too small, we think the polynomial is irreducible`); return [f]; fi; sv:=evalf(ddd[dimC-nd+1]); userinfo(1,appfac,print('sv'=sv)); # obtain the tolerance using in multigcd program from the singular value sv or input . if nargs=3 then e0:=args[3]; else e0:=-ceil(log10(sv)); fi; userinfo(1,appfac,print('e0'=e0)); lprint(` the biggest gap, the last r-th singular values and the number of factors `, max(rats), evalf(sv/ddd[1]),nd); userinfo(1,appfac, print('maxrate'=max(rats))); userinfo(1,appfac,print('nd'=nd)); lprint(` The time for computing the number of factors*****`, time()-time1); #compute the basis GG of the null space of the matrix A. GG:=[]; for i from 1 to nd do # vv[i]:=LinearAlgebra[Row](VVt, dimR-nd+i); vv[i]:=LinearAlgebra[Column](LinearAlgebra:-HermitianTranspose(VVt), dimR-nd+i); gg[i]:=g[1]; for j from 1 to nops(var) do gg[i]:=subs(var[j]=vv[i][j], gg[i]); od; GG:=[op(GG),gg[i]]; od; tg:=0; cd:=0; rq0:=[]; for kk from 1 to 3 do # First,we choose 'rn' and 'vf' randomly for three times, then # we compute the smallest distance between eigenvalues. # We choose the 'rn' and 'vf' which result the largest distance of eigenvalues. rn:=zrandom3(nd); userinfo(2,appfac,print('rn'=rn)); gg:=eval(sum('GG[i]*rn[i]', 'i'=1..nd)); g:=0; g:=gg; fx:=diff(f0,vars[1]); temp:=false; while temp=false do vf:=zrandom1(nops(ord)-1); userinfo(2,appfac,print('vf'=vf)); # choose y=vf which keeps the degrees of fx and f1 in x # fx1:=subs(y=vf,fx); fx1:=expand(subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)},expand(fx))); if degree(fx,vars[1])=degree(fx1,vars[1]) then #f1:= subs(y=vf,f0); f1:=expand(subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)},expand(f0))); if degree(f0,vars[1])=degree(f1,vars[1]) then for i from 1 to nd do if degree(subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)},expand(GG[i]*g)))<>degree(expand(GG[i]*g),vars[1]) then temp:=false; break; elif degree(subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)},expand(GG[i]*fx1)))<>degree(expand(GG[i]*fx1),vars[1]) then temp:=false; break; fi; od; temp:=true; break; fi; fi; od; c:=[seq(c[i],i=1..nd)]; M1:=[]; #compute the characteristic polynomial Eg for i from 1 to nd do tt:= subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)},expand(GG[i]*g)); tq:=polydiv(tt,f1,vars[1]); tt:=expand(tt-tq*f1); for j from 1 to nd do t[j]:=subs({'vars[j+1]=vf[j]' $ 'j' = 1..nops(vf)}, expand(GG[j]*fx1)); tq:=polydiv(t[j],f1,vars[1]); t[j]:=expand(t[j]-tq*f1); od; rp:=expand( tt-expand(sum('c[k]*t[k]', 'k'=1..nd))); coefs:=[coeffs(rp,vars[1])]; E:={}; for l from 1 to nops(coefs) do # We get all the coefficients including var if degree(coefs[l])>0 then E:={op(E), coefs[l]=0} fi; od: v[i]:=subs(LinearAlgebra[LeastSquares](E,{op(c)}),c); M1:=[op(M1),v[i]]; od; M:=convert(M1,Matrix); rq:=LinearAlgebra[Eigenvalues](M,output=list); userinfo(2,appfac,print('rq'=rq)); # we compute the distance:=min{|rq[i]-rq[j]|,1<=i0 then if min(op(L1))>cd then cd:=min(op(L1)); tg:=g; rq0:=rq; fi; fi; od; userinfo(2,appfac,print('cd'=cd)); g:=tg; rq:=rq0; userinfo(2,appfac,print('rq'=rq)); FF1:=[]; FF2:=[]; FF3:=[]; fy:=diff(f0,y); # compute GCD to get the irreducible factors for i from 1 to nd do ff[i]:=multigcd(f0,g-expand(rq[i]*fx),vars,e0); if ff[i][1]=1 then FF1:=[op(FF1),ff[i][2]]; elif ff[i][1]=2 then FF3:=[op(FF3),ff[i][2]]; else FF2:=[op(FF2), ff[i][2]]; fi; od; lprint(`The time for the entire factorization******`, time()-time1); # check if sum of the degrees of all factors equal to the degree of f k:='k'; if nops(FF2)=0 then if sum('degree(FF1[i])','i'=1..nops(FF1))+sum('degree(FF3[i])','i'=1..nops(FF3))=degree(f0) then userinfo(1,appfac,`the factorization is correct`); RETURN([op(FF1), op(FF3)]); elif nops(FF1)=1 and nops(FF3)=1 then userinfo(1,appfac,`the factorization is wrong, special case as nops(FF1)=1,nops(FF3)=1`); RETURN([FF1[1], mcofactor(f0, FF1[1], ord, e0)]); elif nops(FF1)=1 then userinfo(1,appfac,`the factorization is not correct, check the factors in FF3`,nops(FF3)); if nops(FF3)=1 then ffg2:=mcofactor(f0, product(FF1[k], k=1..nops(FF1)),ord,e0); if degree(ffg2)>0 then RETURN([op(FF1),ffg2]); else RETURN(FF1); fi; else ffg2:=mcofactor(f0, product(FF1[k], k=1..nops(FF1)),ord,e0); if degree(ffg2)>0 then RETURN([op(FF1), op(appfac(mcofactor(f0, product(FF1[k], k=1..nops(FF1)), ord, e0),ord))]); else RETURN(FF1); fi; fi; elif nops(FF1)=nd then if nops(FF1)=2 then userinfo(1,appfac,`the factorization is wrong, try again for the rest factors except the first one`); RETURN([FF1[1], mcofactor(f0, product(FF1[k], k=2..nops(FF1)), ord, e0)]); else RETURN([FF1[1], op(appfac(mcofactor(f0, FF1[1], ord, e0),ord))]); fi; elif nops(FF1)<1 then if nops(FF3)>2 then userinfo(1,appfac,`the factorization is wrong, try again for the rest factors except the first one`); RETURN([FF3[1], op(appfac(mcofactor(f0, FF3[1], ord, e0),ord))]); else userinfo(1,appfac,`the factorization is wrong, try again for the rest factors except the first one`); RETURN([FF3[1], mcofactor(f0, product(FF3[k], k=2..nops(FF3)), ord, e0)]); fi; fi; elif nops(FF2)=1 then userinfo(1,appfac,`there is only one bad factor`); if nops(FF3)>0 and nops(FF1)>0 then if degree(product(FF1[k],k=1..nops(FF1))*product(FF3[k],k=1..nops(FF3)))<=degree(f0) then ffg:=mcofactor(f0,product(FF1[k],k=1..nops(FF1))*product(FF3[k],k=1..nops(FF3)),ord,e0); if sum('degree(FF1[i])','i'=1..nops(FF1))+sum('degree(FF3[i])','i'=1..nops(FF3))+deg(ffg)=degree(f0) then RETURN([op(FF1),op(FF3),ffg]); else RETURN([op(FF1),op(appfac(mcofactor(f0,product(FF1[k],k=1..nops(FF1)),ord,e0),ord,e0))]); fi; else RETURN([op(FF1),op(appfac(mcofactor(f0,product(FF1[k],k=1..nops(FF1)),ord,e0),ord,e0))]); fi; elif nops(FF3)=0 then RETURN([op(FF1), mcofactor(f0, product(FF1[k], k=1..nops(FF1)), ord, e0)]); elif nops(FF1)=0 then RETURN([op(FF3), mcofactor(f0, product(FF3[k], k=1..nops(FF3)), ord, e0)]); fi; elif nops(FF2)=nd then userinfo(1,appfac,`the factorization is totally wrong, try again for small e`); RETURN(appfac(f0,ord,e0-1)); else userinfo(1,appfac,`we get partial correct factorization,try the routine for the cofactor of FF1`,nops(FF2),nops(FF3)); if nops(FF3)>0 and nops(FF1)>0 then ffg:=cofactor(f0,product(FF1[k],k=1..nops(FF1))*product(FF3[k],k=1..nops(FF3)),ord,e0); FF4:=appfac(ffg,ord); if sum('degree(FF1[i])','i'=1..nops(FF1))+sum('degree(FF3[i])','i'=1..nops(FF3))+sum('degree(FF4[i])','i'=1..nops(FF4))=degree(f0) then RETURN([op(FF1),op(FF3),op(FF4)]); else RETURN([op(FF1),op(appfac(mcofactor(f0, product(FF1[k], k=1..nops(FF1)), ord, e0), ord))]); fi; elif nops(FF3)=0 then RETURN([op(FF1), op( appfac(mcofactor(f0, product(FF1[k], k=1..nops(FF1)), ord, e0), ord))]); elif nops(FF1)=0 then RETURN([op(FF3), op( appfac(mcofactor(f0, product(FF3[k], k=1..nops(FF3)), ord, e0), ord))]); fi; fi; end proc: # vecc code is to choose S[i] which satisfy degree(S[i],i=1..nops(S))< degree(f). vecc:=proc(S,m) local V,i,k,V1; V1:=vecc1(S); k:=nops(V1); V:=[]; for i from 1 to k do if numvecc(V1[i])< m then V:=[op(V),V1[i]]; fi; end do; RETURN(V); end proc; # vecc1 code is to compute all the vector which denotes the degree and every degree satisfy v[i] degree(g) and degree(da) > degree(f)) then userinfo(1,multigcd,`Something went wrong. Degree of the GCD is too big.`); RETURN([0,L1]); else if er < 10^(-e) or er<10^(-6) then userinfo(1,multigcd,`The estimated degree of GCD is correct, we get a good gcd.`); RETURN([1,da*L1]); elif er<10^(-e+1) then userinfo(1,multigcd,`The estimated degree of GCD may be correct`); RETURN([2,da*L1]); elif er<10^(-e+2) and er<10^(-2) then userinfo(1,multigcd,`The estimated degree of GCD may be wrong`); RETURN([3,da*L1]); else userinfo(1,multigcd,` The estimated degree of GCD maybe wrong, we get a bad gcd.`); RETURN([0, da*L1]); fi; fi; end proc: # Check if the tested total degree $k$ is the correct one. # Return the correct total degree of bivariate GCD and cofactors. testk:=proc(k,f,g,ord,e) local t,v0,i,epsilon,x0,y0,KK,VV,tt,QQ,RR,PP,QQ1,RR1,jj,v0R,v0R2,jmp,m0; t:=k; userinfo(1,testk,print('t'=t)); (PP,QQ,RR):=solve_x0(t,f,g,ord,e); epsilon:=PP[1]; userinfo(3,testk,print('sv'=epsilon)); x0:=PP[2]; jmp:=1; if epsilon > 10^(-e) then if k<2 then RETURN([k,x0]); fi; #The degree of the factor is at least 1. v0:=epsilon; if round(log10(v0))<=(-e+1) then # k is possibly small. v0R2:=solve_x02(k,f,g,QQ,RR,x0,ord,e); userinfo(3,testk,print('sv2'=v0R2)); if v0R2/v0 > 50 then RETURN([k,x0]); else jmp:=v0R2/v0; fi; m0:=x0; if k max(jmp,10) then RETURN([k+1,y0]); fi; fi; fi; fi; if k=1 then RETURN([k,x0]); fi; #The degree of the factor is at least 1. for jj from k-1 by -1 to 1 do # k is possible larger, do loop statement until we find rank deficiency 1 tt:=jj; (VV,QQ1,RR1):=solve_x0(tt,f,g,ord,e); v0R:=VV[1];y0:=VV[2]; userinfo(1,testk,`tt=`,tt); userinfo(3,testk,print('sv'=v0R)); if v0R < 10^(-e+1) and v0R > 10^(-e) then v0R2:=solve_x02(tt,f,g,QQ1,RR1,y0,ord,e); if v0R2/v0R > 50 then RETURN([tt,y0]); elif v0R2/v0R < jmp and jmp>10 then RETURN([tt+1,m0]); else m0:=y0; jmp:= v0R2/v0R; userinfo(3,testk,print('jmp'=jmp)); fi; elif v0R< 10^(-e) then RETURN([tt,y0]); fi; end do; RETURN([tt,y0]); # end if v0<10^(-e+1) then else # k is small possible. v0:=epsilon; KK:=solve_x01(k,f,g,x0,QQ,RR,v0,ord,e); RETURN(KK); fi; end proc; #Polynomial division, check the exact divisibility. polydiv:=proc(f::polynom, g::polynom,x) local n,m,A,i,j,sol,cf,q; n:=degree(f); m:=degree(g); A:=Matrix(n+1,n-m+1); if n=0 then A[i,j]:=coeff(g,x,m+j-i); fi; od; od; cf:=Vector(n+1,i->coeff(f,x,n-i+1)); sol:=LinearAlgebra[LeastSquares](A,cf); q:=add(sol[i]*x^(n-m-i+1),i=1..n-m+1); RETURN(q); end proc: # The 'ugcdtest' program is used to acquire the degree of GCD of two univariate polynomials. # First we compute SVD for Sylvester matrix, choose the biggest gap of the singular values between [10^(-e+1),10^(-e-4)]. # If the gap is large enough, we can get the degree of GCD and 1, here 1 shows we need not check. # Otherwise we correct the degree of GCD and see if the k-th Sylvester matrix is of rank deficiency 1. ugcdtest:=proc(f,g,x,e) local i,j,jj # loop variables ,m,n # the degree of polynomial ,f0,g0 #polynomials ,MM1,M #Sylvester matrix ,cc # singular values list ,r # norm of matrix ,addv # zero vector ,len1 # integral ,len2 # column dimension of matrix ,len3 # the number of sets ,rats # sets ,t # the degree of GCD got by singular value ,sv; # the singular value at the largest gap f0:=f; g0:=g; userinfo(3,ugcdtest,print('f0'=f0)); userinfo(3,ugcdtest,print('g0'=g0)); n:=degree(f0,x); m:=degree(g0,x); userinfo(2,ugcdtest,print('n'=n)); userinfo(2,ugcdtest,print('m'=m)); if m < 1 or n<1 then RETURN(0); fi; # constant MM1:=LinearAlgebra[SylvesterMatrix](f0,g0,x); M:=LinearAlgebra[Transpose](MM1); len2 := LinearAlgebra[ColumnDimension](M); cc:=LinearAlgebra[SingularValues](M,output='list'); #let dimension of cc equal to column dimension of M if nops(cc)< LinearAlgebra[ColumnDimension](M) then addv:=convert(Vector(LinearAlgebra[ColumnDimension](M)-nops(cc)),list); cc:=[op(cc),op(addv)]; fi; userinfo(2,ugcdtest,print('cc'=cc)); t:=0; if cc[-1]>10^(-e+1) then RETURN([0,1]); fi; len3:=nops(cc); for jj from nops(cc) by -1 to 1 do if cc[jj]>10^(-e+1) then # ensure the degree of GCD is not lower than len3 len3:=jj; break; fi; end do; if len3> len2-2 then len1:=len3; else for jj from len2-1 by -1 to len3 do if cc[jj]> 10^(-e-6) then len1:=jj; break; fi; end do; fi; # end if len3> len2-2 if len1-len3>0 then # see if there are singular values between 10^(-e+1) and 10^(-e-4) rats := seq(cc[i]/cc[i+1], i=len3..len1-1); member(max(rats), [rats], 't'); t:= nops([rats])+len2-len1 - t+1; elif len10 sv:=cc[len2-t+1]; userinfo(1,ugcdtest, print('maxrate'=max(rats))); userinfo(1, ugcdtest, print('sv'=sv)); if t=0 then RETURN([0,0]); fi; if max(rats)>10^3 then # we need not test the degree of GCD RETURN([t,1]); elif t>min(degree(f), degree(g)) then RETURN([min(degree(f), degree(g)),0]); else # we need to check if the degree of gcd is correct. t:=utest(t,f0,g0,e); for jj from len2-t[1]+1 to nops(cc)-1 do sv:=cc[jj]; #correct the degree if evalf(sv)<1.1*10^(-e) then #correct the degree userinfo(1, ugcdtest, print('KK'=len2-jj+1)); RETURN([len2-jj+1,t[2]]); fi; od; RETURN([1,t[2]]); fi; end proc; # Subroutine of ugcdtest which is used to correct the degree of univariate GCD and x0 utest:=proc(k,f,g,e) local t,tt # degree of polynomial ,i,jj # loop variable ,epsilon # the least singular value ,v0 # the least singular vale ,QQ,RR # QR decomposition ,QQ1,RR1 # QR decomposition ,PP # the least singular value and the null space ,VV # the least singular value and the null space ,KK # the least singular value and the null space ,v0R # the least singular value ,v0R2 # the second singular value ,jmp # the gap ,m0 # the vector about the least singular value ,x0,y0; # the vector about the least singular value t:=k; (PP,QQ,RR):=unsolve_x0(t,f,g,e); epsilon:=PP[1]; userinfo(1,utest,`t`,t); userinfo(3,utest,print('sv'=epsilon)); x0:=PP[2]; jmp:=1; if epsilon > 10^(-e) then if k=1 then RETURN([k,0]); fi; #The degree of the gcd is at least 1. v0:=epsilon; if v0<10^(-e+1) then #check if there is a big gap v0R2:=unsolve_x02(t,f,g,QQ,RR,x0,e); if v0R2/v0 > 10^5 then # k maybe correct userinfo(3,utest,print('KK'=[t,1])); RETURN([t,1]); elif v0R2/v0 > 50 then # k is possibly correct userinfo(3,utest,print('KK'=[t,0])); RETURN([t,0]); else jmp:=v0R2/v0; fi; # increase k to check if there is a bigger jump (VV,QQ1,RR1):=unsolve_x0(t+1,f,g,e); v0R:=VV[1];y0:=VV[2]; if v0R < 10^(-e+1) then v0R2:=unsolve_x02(t,f,g,QQ1,RR1,y0,e); if v0R2/v0R > 10^5 then # k+1 is the correct degree of gcd userinfo(3,utest,print('KK'=[t+1,1])); RETURN([t+1,1]); elif v0R2/v0R > max(jmp,50) then # k+1 is better userinfo(3,utest,print('KK'=[t+1,0])); RETURN([t+1,0]); fi; fi; fi; #decrease k to check if there is a bigger jump for jj from t-1 by -1 to 1 do # do loop statement until we find rank deficient 1 tt:=jj; if tt=1 then RETURN([tt,0]) fi; # the degree of the gcd is at least 1 (VV,QQ1,RR1):=unsolve_x0(tt,f,g,e); v0R:=VV[1];y0:=VV[2]; if v0R < 10^(-e+1) and v0R > 10^(-e) then # see if the rank deficient is 1 v0R2:=unsolve_x02(tt,f,g,QQ1,RR1,y0,e); m0:=y0; jmp:= v0R2/v0R; userinfo(3,utest,print('jmp'=jmp)); if v0R2/v0R > 10^5 then userinfo(1,utest,print('KK'=[tt,1])); RETURN([tt,1]); elif tt 10^5 then userinfo(1,utest,print('KK'=[tt+1,1])); RETURN([tt+1,1]); elif v0R2/v0R > jmp and v0R<10^(-e+1) then userinfo(1,utest,print('KK'=[tt+1,0])); RETURN([tt+1,0]); fi; else RETURN([tt,0]); fi; elif v0R< 10^(-e) then v0R2:=unsolve_x02(tt,f,g,QQ1,RR1,y0,e); if v0R2/v0R > 10^5 then userinfo(1,utest,print('KK'=[tt,1])); RETURN([tt,1]); elif v0R2/v0R > jmp then userinfo(1,utest,print('KK'=[tt,0])); RETURN([tt,0]); else userinfo(1,utest,print('KK'=[tt+1,0])); RETURN([tt+1,0]); fi; fi; end do; RETURN([1,0]); else v0:=epsilon; KK:=unsolve_x01(k,f,g,x0,QQ,RR,v0,e); userinfo(1,utest,print('KK'=KK)); RETURN(KK); fi; end proc; # Compute the singular value and singular vector of the $k$-th Sylvester matrix of the univariate polynomial f and g. # $k$ is the tested degree of the univariate GCD of f and g. # It adopted zeng's routine. unsolve_x0:=proc(k,f,g,e) local t,M,pp,QQ,RR; t:=k; M:=sych(f,g,x,t); (pp,QQ,RR):=leastsv(M,e); return pp,QQ,RR; end proc; # When the least singular values is small enough, we check if the $k$-th Sylvester matrix is of rank deficiency 1. # If the rank deficiency is bigger than 1 then we correct the degree by k=k+1 loop. unsolve_x01:=proc(k,f,g,x0,Q,R,v0,e) local t,M,len1,tmp,len2,tt,QQ,RR,kk,i,DD,dimen,bb,PP,QQ1,RR1,zz,epsilon,pp, r,bb1,bb2,v02,jj,y0,MM,MM1,z0,SS,QQ3,RR3,epsilon2,jmp; tmp:=degree(g)-k: len1:=((tmp+1)*(tmp+2))/2; tmp:=min(degree(f)-k,degree(g)-k): len2:=len1 + ((tmp+1)*(tmp+2))/2; if tmp=0 then RETURN([k,0]); fi; epsilon:=unsolve_x02(k,f,g,Q,R,x0,e); jmp:=epsilon/v0; if epsilon/v0>1000 then # k is right RETURN([k,1]); fi; t:=k; for jj from 1 to tmp do # increase k to right degree. t:=k+jj; M:=sych(f,g,x,t); (pp,QQ,RR):=leastsv(M,e); epsilon:=pp[1]; z0:=pp[2]; if epsilon1000 then RETURN([t,1]); fi; else RETURN([t-1,0]); fi; end do; RETURN([t,0]): end proc: #Compute the second least singular value by using zeng's routine. unsolve_x02:=proc(k,f,g,Q,R,x0,e) local pp,len1,kk,r,MM,QQ,RR,epsilon,QQ1,RR1,tmp,len2; tmp:=degree(g)-k: len1:=((tmp+1)*(tmp+2))/2; tmp:=degree(f)-k: len2:=len1 + ((tmp+1)*(tmp+2))/2; QQ1:=Q;RR1:=R; kk:=LinearAlgebra[RowDimension](RR1); r:=evalf(LinearAlgebra[MatrixNorm](RR1,infinity)); MM:=Matrix(LinearAlgebra[RowDimension](RR1)+1,LinearAlgebra[ColumnDimension](RR1),[[LinearAlgebra[Transpose](r*x0)],[RR1]]); (pp,QQ,RR):=leastsv(MM,e); epsilon:=pp[1]; RETURN(epsilon); end proc: # For the computed degree of the GCD of f and g, we form the bivariate Sylvester matrix. # Applying zeng's method to compute the least singular value and the singular vector about it. solve_x0:=proc(k,f,g,ord,e) local M,QQ,RR,pp; if nops(ord)=2 then M:=bvsylo(f,g,ord,k); else M:=mvsylo(f,g,ord,k); fi; (pp,QQ,RR):= leastsv(M,e); RETURN(pp,QQ,RR); end proc: # When the least singular values is small enough, we test if the $k$-th bivariate Sylvester matrix is rank deficiency one. # If the rank deficiency is bigger than 1 then we correct $k$ to $k+1$. solve_x01:=proc(k,f,g,x0,Q,R,v0,ord,e) local t,M,len1,tmp,len2,epsilon,jj,pp,QQ,RR,z0,y0,epsilon2,v02,jmp; tmp:=degree(g,ord)-k: len1:=((tmp+1)*(tmp+2))/2; tmp:=min(degree(f)-k,degree(g)-k): len2:=len1 + ((tmp+1)*(tmp+2))/2; if tmp=0 then RETURN([k,x0]); fi; epsilon:=solve_x02(k,f,g,Q,R,x0,ord,e); jmp:=epsilon/v0; userinfo(2, solve_x01, print('jmp'=jmp)); if epsilon/v0>10 then # k is right RETURN([k,x0]); fi; for jj from 1 to tmp do # increase k to get the right degree. t:=k+jj; userinfo(2,solve_x01,print('t'=t)); M:=mvsylo(f,g,ord,t); (pp,QQ,RR):=leastsv(M,e); epsilon:=pp[1]; z0:=pp[2]; if epsilon50 then RETURN([t,y0]); fi; fi; userinfo(3,solve_x01,print('KK'=t)); userinfo(3,solve_x01,print('sv'=epsilon)); if epsilon>10^(-e)then if epsilon>evalf(10^(-e+1)) then if t=k+1 then RETURN([k,x0]); fi; RETURN([t-1,y0]); else v02:=solve_x02(t,f,g,QQ,RR,z0,ord,e); userinfo(1,solve_x01,`jump=`, jmp, v02/epsilon); if v02/epsilon> max(5,jmp) then RETURN([t,z0]); else if t=k+1 then RETURN([k,x0]); else RETURN([t-1,y0]); fi; fi; fi; fi; end do; RETURN([k,x0]); end proc: # This routine computes approximately the degree of the GCD of two bivariate polynomials # Linear maps are used to convert bivariate polynomials to univariate polynomials. # We compute the degree of univariate polynomials GCD to estimate the total degree of the bivariate GCD. testtdeg:=proc(f,g,ord,e) local nx,t,a,mm,k,nx1,c,uc,b,A1,D1,pw,f1,g1,e0,tt,kk,pww,ww,nx2,nx3,nx4,jj; nx:=[]; t:=1; pw:=[]; mm:=0; k:=nops(ord); kk:=1; ww:=0; tt:=0; nx3:=[]; nx4:=[]; # Repeat several times of the linear maps to test the possible degree of GCD. while t<10 and mm<10 do b:=zrandom1(k); a:=zrandom1(k); # do the linear maps if degree(f)<10 or max(norm(expand(subs({'ord[j]=0' $ 'j' = 1..k},expand(f))),2), norm(expand(subs({'ord[j]=0' $ 'j' = 1..k},expand(g))),2))<10^(-e) then f1:=expand(subs({'ord[j]=((a[j]+b[j])/2*ord[1]-(-1)^(j+t)*b[j])' $ 'j' = 1..k},expand(f))); g1:=expand(subs({'ord[j]=( (a[j]+b[j])/2*ord[1]-(-1)^(j+t)*b[j])' $ 'j' =1..k},expand(g))); e0:=floor(log10(min(norm(f1,infinity), norm(g1,infinity)))); userinfo(1,testtdeg,`e,e0`,e,e0); else f1:=expand(subs(ord[k]=ord[1],subs({'ord[j]=((a[j]+b[j])/2*ord[k]-(-1)^(j+t)*b[j])' $ 'j' = 1..k-1},expand(f)))); g1:=expand(subs(ord[k]=ord[1], subs({'ord[j]=( (a[j]+b[j])/2*ord[k]-(-1)^(j+t)*b[j])' $ 'j' =1..k-1},expand(g)))); e0:=floor(log10(min(norm(f1,infinity), norm(g1,infinity)))); userinfo(1,testtdeg,`check the tolerance level of ugcdtest $$$$$$`, e, e0); fi; if e0>e-3 or e0<0 then if e0>=3 then if e>3 then nx1:=ugcdtest(f1,g1,ord[1],e-2); else nx1:=ugcdtest(f1,g1,ord[1],e); fi; #print(`nx1=`,nx1); else nx1:=ugcdtest(f1,g1,ord[1],e); fi; else nx1:=ugcdtest(f1,g1,ord[1],e-e0); fi; # stop the routine if the possible degree appears at least 3 times. nx4:=[nx1,op(nx4)]; tt:=nx1[1]; nx3:=[tt,op(nx3)]; # print(`nx3=`,nx3); if numboccur(nx3,tt)=3 then userinfo(1,testtdeg,`stop the ugcdtest early`, nx3); for jj from 1 to nops(nx4) do if nx4[jj][1]=tt then pww:=[nx4[jj],op(pww)]; fi; od; if min(pww[1][2],pww[2][2],pww[3][2])=1 then # there is big jump, the estimated degree is good RETURN([tt,1]); else RETURN([tt,0]); # there is no big jump, we need to check the estimated degree fi; fi; if nx1[1]=0 then mm:=mm+1; else nx:=[op(nx), (-1)^(nx1[2]+1)*nx1[1]]; t:=t+1; fi; od; userinfo(2,testtdeg,'nx'=nx ); if mm>10 then # the gcd is a constant RETURN([0,0]); fi; # we choose a least degree which appears in nx. #nx:=max(op(nx)); nx:=max(op(map(abs,nx))); RETURN([abs(nx),0]); #RETURN([abs(nx),(1+sign(nx))/2]); end proc: # We compute the second least singular value to check if the rank deficiency of the matrix is one. solve_x02:=proc(k,f,g,Q,R,x0,ord,e) local t,QQ1,RR1,kk,r,MM,pp,QQ,RR,epsilon; QQ1:=Q;RR1:=R; kk:=LinearAlgebra[RowDimension](RR1); r:=evalf(LinearAlgebra[MatrixNorm](RR1,infinity)); MM:=Matrix(LinearAlgebra[RowDimension](RR1)+1,LinearAlgebra[ColumnDimension](RR1),[[LinearAlgebra[Transpose](r*x0)],[RR1]]); (pp,QQ,RR):=leastsv(MM,e); epsilon:=pp[1]; RETURN(epsilon); end proc: # This routine computes the least singular value of the matrix and the associated singular vector leastsv:=proc(M,e) local QQ,RR,kk,x0,r,DD,bb1,dimen,QQ1,RR1,epsilon,i,bb2,bb,zz,rk,pp,sigma,dimC,M1,M2,b1,b2,z0; (QQ,RR,rk):=LinearAlgebra[QRDecomposition](M, output=['Q', 'R', 'rank'],fullspan); userinfo(1,leastsv,print('rk'=rk)); dimC:=LinearAlgebra[ColumnDimension](M); if rk degree(g) and degree(da) > degree(f)) then userinfo(1,multigcd,`Something went wrong. Degree of the GCD is too big.`); RETURN(1); else RETURN(da/norm(expand(da),2)); fi; end proc: # The routine computes a basis of the multivariate polynomial. numbve:=proc(n,m) local i,S,v,N; S:=[]; N:=(n+m)!/(n!*m!); for i from 0 to N-1 do v:=combinat[inttovec](i,m); S:=[op(S),v]; end do; RETURN(S); end proc: # The routine converts the coefficients of a multivariate polynomial to a vector. convec:=proc(f,ord) local i,m,n,vba,ve,N; m:=degree(f,{op(ord)}); n:=nops(ord); N:=(n+m)!/(n!*m!); vba:=numbve(m,n); ve:=Vector(N); for i from 1 to N do ve[N-i+1]:=convs(f,ord,vba[i]); end do; RETURN(ve); end proc; # The routine computes the coefficients w.r.t. given degrees of variables. convs:=proc(f,ord,V) local i,n,f0; n:=nops(ord); f0:=f; for i from 1 to n do; f0:=coeff(f0,ord[i],V[i]); end do; RETURN(f0); end proc; # This routine converts a vector to a multivariate polynomial. veccoe:=proc(v,ord) local i,n,N,f0; n:=nops(ord); N:=LinearAlgebra[Dimension](v); f0:=0; for i from 1 to N do f0:=f0+v[i]*vecmo(combinat[inttovec](N-i,n),ord); end do; end proc; # This routine converts a vector to a monomial. vecmo:=proc(v,ord) local i,n,s; n:=nops(ord); s:=1; for i from 1 to n do s:=s*ord[i]^(v[i]); end do; RETURN(s); end proc: # This routine computes the matrix which represents the multiplication of # f by a polynomial of total degree k. This is called the Cauchy matrix # by some. # This routine generalize John's definition for bivariate Cauchy matrix mvcauchyo := proc(f, ord,deg) local N,td,i, unknowns, g, eqns, pp,V ,n,A; # Protect the unknowns td := degree(f,{op(ord)}); n:=nops(ord); N:=(n+deg)!/(n!*deg!); g:=0; V:=Vector(N); for i from 1 to N do V[i]:=a[i]; end do; pp:=convert(V,list); g:=veccoe(V,ord); eqns := expand(f*g); eqns := convert(convec(eqns,ord),list); eqns := map(x->x=0, eqns); A:=LinearAlgebra[GenerateMatrix](eqns, pp)[1]; RETURN(A); end proc; # This routine call the Cauchy matrix routine to generate the multivariate version # of the Sylvester matrix of f and g. # The parameter m decreases the size of the matrix. m=1 gives the full # Sylvester-style matrix. # Generalize the algorithm bvsyo written by John P May mvsylo:= proc(f,g,ord,m) local M3 ; if nops(ord)=2 then M3:=Matrix([bvcauchyo(f,ord,degree(g,{op(ord)})-m),bvcauchyo(g,ord,degree(f,{op(ord)})-m)]); else M3:=Matrix([mvcauchyo(f,ord,degree(g,{op(ord)})-m),mvcauchyo(g,ord,degree(f,{op(ord)})-m)]); fi; RETURN(M3); end proc; # Obtain the epsilon for checking the approximate squarefree of the polynomial. epsilonsqrf:=proc(f,ord ) # f is a polynomial which we want to factor # ord is bivariate variables local S1,eqns,Q,P,V,jj,kk,EQ,VARS,PP,vars,f0,g,tf,i,j,h,pp,var,linsys,lineqs,A,b,SS,VVt,dimR,ddd,nd,dimC,addv,rats,sv,e0; vars:=ord; f0:=evalf(f/norm(expand(f),2)); g:=0; tf:=degree(f0); if tf<=0 then RETURN([]); fi; S1:=[]; for j from 1 to nops(ord) do S1:=[op(S1),degree(f0,vars[j])]; end do; #construct polynomial g[i] satisfies deg(g[i])<=(m-1,n) and deg(g[i]) < deg(f); eqns:={}; for i from 1 to nops(vars) do Q[i]:=subtvec(S1,i); P[i]:=vecc(Q[i],tf); g[i]:=0; V:=Vector(nops(P[i])); for jj from 1 to LinearAlgebra:-Dimension(V) do V[jj]:=a[i,jj]; end do; pp:=convert(V,list); for kk from 1 to nops(pp) do g[i]:=g[i]+pp[kk]*vecmo(P[i][kk],vars); end do; eqns:=[op(eqns),g[i]]; end do; EQ:={};VARS:={}; for i from 2 to nops(vars) do PP[i]:=evalf(expand(g[i]*diff(f0,ord[1])-f0*diff(g[i],ord[1])+f0*diff(g[1],ord[i])-g[1]*diff(f0,ord[i]))); EQ:= EQ union {coeffs(PP[i],ord)}; VARS:= VARS union {op(indets(PP[i]))} minus {op(ord)}; end do; var:=[op(VARS)]; # construct the Ruppert matrix A lineqs:=['EQ[i]=0' $ 'i'=1..nops(EQ)]; A,b:=LinearAlgebra[GenerateMatrix](lineqs,var); SS, VVt := LinearAlgebra[SingularValues](A, output=['S','Vt']): dimR:= LinearAlgebra[RowDimension](VVt); ddd:=convert(SS,list); nd:=1; dimC:=LinearAlgebra[ColumnDimension](A); # to ensure dimension(ddd)=dimC if nops(ddd)< dimC then addv:=convert(Vector(dimC-nops(ddd)),list); ddd:=[op(ddd),op(addv)]; fi; # compute the biggest gap in the singular values rats := seq(ddd[i]/ddd[i+1], i=(1..dimC-2)); member(max(rats), [rats], 'nd'); nd := nops([rats])+2 - nd; sv:=5*ddd[dimC-nd+1]; e0:=-ceil(log10(sv))-2; return max(e0-1,2),ddd; end proc; #This routine tests if the polynomial g is a factor of the polynomial f. testcofactor:=proc(f, g, ord) local k # degree of g ,Mf # bvcauchyo matrix ,da # the quotient ,tmp # the degree of the quotient ,fv # the vector converted by polynomial ,dv # the solution of Mg.dv=gv ,er; k:=degree(g); tmp:=degree(f)-k: if tmp<0 then RETURN(10); fi; if nops(ord)=2 then Mf:=bvcauchyo(g,ord,degree(f)-degree(g)); else Mf := mvcauchyo(g,ord,degree(f)-degree(g)); fi; fv := convec(f,ord); dv := LinearAlgebra[LeastSquares](Mf, fv); da := veccoe(dv,ord); userinfo(3,multigcd,print('da'=da)); userinfo(2,multigcd,print('deg_gcd'=degree(da), 'deg_g'=degree(g), 'deg_h'=degree(f)) ); if (degree(da) > degree(g) and degree(da) > degree(f)) then userinfo(1,multigcd,`Something went wrong. Degree of the GCD is too big.`); RETURN(10); else er:=LinearAlgebra[Norm](Mf.dv-fv,infinity); RETURN(er); fi; end proc: # The routine compute the multiplicity of the factors multiplicityoffactor:=proc(f,S,ord,e) local i,num,cofa,def,SS1,SS2,j,kk,jj,fy,er; cofa:=1; def:=degree(f); SS1:=[]; SS2:=[]; for i from 1 to nops(S) do if degree(f)>=degree(S[i]) then SS1:=[S[i],op(SS1)]; else SS2:=[[S[i],1],op(SS2)]; fi; end do; num:=nops(SS1); for j from 1 to num do if def-degree(cofa) >= degree(SS1[j]) then kk:=floor((def-degree(cofa))/degree(SS1[j])); userinfo(1,multiplicityoffactor,print('kk'=kk)); for jj from 1 to kk do fy:=expand(SS1[j]^jj*cofa); er:=testcofactor(f,fy,ord); userinfo(1,multiplicityoffactor,print('er'=er)); if er>10^(-e+2) and er>10^(-1) then SS2:=[[SS1[j],jj],op(SS2)]; cofa:=expand(SS1[j]^(jj-1)*cofa); break; else if jj=kk then SS2:=[[SS1[j],kk+1],op(SS2)]; cofa:=expand(SS1[j]^(kk)*cofa); fi; fi; end do; else SS2:=[[SS1[j],1],op(SS2)]; fi; end do; RETURN(SS2); end proc; # makelist assigns every polynomial whose multiple is 1. makelist:=proc(F) local n,i,SS; n:=nops(F); SS:=[]; for i from 1 to n do SS:=[[F[i],1],op(SS)]; end do; SS; end proc: #The routine arranges the factors according to the multiplicities orderlist:=proc(L) local n,i,j,S,S1,S2,kk,nm; if nops(L)=1 then RETURN(L);fi; S:=[L[1]]; n:=nops(L); for i from 2 to n do kk:=nops(S); if L[i][2]<=S[1][2] then S:=[L[i],op(S)]; elif L[i][2]>=S[kk][2] then S:=[op(S),L[i]]; else for j from 1 to kk-1 do if L[i][2]>=S[j][2] and L[i][2]<=S[j+1][2] then nm:=j; fi; end do; S1:=[seq(S[k],k=1..nm)]; S2:=[seq(S[k],k=nm+1..nops(S))]; S1:=[op(S1),L[i]]; S:=[op(S1),op(S2)]; fi; end do; end proc: # The routine chooses a main variable for multivariate polynomial. mainvariable:= proc(f::polynom,ord) local n,var,VAR,CC,x,i,b,f0; n:=nops(ord); if n=2 then RETURN(bivmainvariable(f,ord)); fi; b:=zrandom1(n-2); x:=ord[1]; for i from 2 to n do CC:=zrandom1(n-2); var:=[op({op(ord)} minus {x,ord[i]})]; f0:=expand(subs({'var[j]=CC[j]' $ 'j' = 1..n-2},expand(f))); VAR:=bivmainvariable(f0,[x,ord[i]]); x:=VAR[1]; end do; RETURN([x,op({op(ord)} minus {x})]); end proc; # This routine select the main variable of f w.r.t x, y # If there is a univariate factor in x only, then x should be chosen as the main variable # We still can not deal with the case where there exist univariate factors g(x) and h(y). bivmainvariable:= proc(f::polynom,ord) local L1, L2,i,u1,v1,sv1,M1,u2,v2,sv2,M2,j,pp1,pp2,QQ1,QQ2,RR1,RR2; L1:=[]; L2:=[]; for i from degree(f,ord[1]) by -1 to 0 do if evalf(norm(coeff(f,ord[1],i),2))>10^(-2) then if degree(coeff(f,ord[1],i))=0 then RETURN(ord); elif degree(coeff(f,ord[1],i))> 0 then L1:=[op(L1), coeff(f,ord[1],i)]; fi; fi; od; if nops(L1) <2 then RETURN([ord[2], ord[1]]); fi; for j from degree(f,ord[2]) by -1 to 0 do if evalf(norm(coeff(f,ord[2],j),2))>10^(-2)then if degree(coeff(f,ord[2],j))=0 then RETURN([ord[2], ord[1]]); elif degree(coeff(f,ord[2],j))> 0 then L2:=[op(L2), coeff(f,ord[2],j)]; fi; fi; od; if nops(L2) <2 then RETURN([ord[1], ord[2]]); fi; i:=2; while i evalf(sv2) then RETURN([ord[1],ord[2]]); else RETURN([ord[2],ord[1]]); fi; end proc: # cut off small terms may cause trouble. #example: input f=1+10^(-16)*I+(2+2*10^(-15)*I)*x*y; ord=[x,y],e=10 then #cutoff1(f,ord,e)=1+2*x*y. cutoff1:=proc(f,ord,e) local l,i,g,v,vv; #lprint(`f=`,f,ord,e); l:=[coeffs(expand(f),ord,'v')]; #lprint(`l=,v=`,l,v); g:=0; vv:=[v]; for i to nops(l) do if abs(Re(l[i]))> 10^(-e-2) then if abs(Im(l[i]))>10^(-e) then g:=g+l[i]*vv[i]; #lprint(`g=`,g); else g:=g+Re(l[i])*vv[i]; fi; fi; od; v:='v'; RETURN(g); end: # This routine converts an input bivariate polynomial into a vector # Written by John P May bv2veco := proc(f,ord) local i, j, tmpcoeff, td, fvec,x,y; x:=ord[1]; y:=ord[2]; tmpcoeff:=map(PolynomialTools[CoefficientList], PolynomialTools[CoefficientList](f,x) , y); td := degree(f, {x,y}); fvec := []; for i from td to 0 by -1 do for j from 0 to i do # Insert coeff of x^j*y^(i-j) try fvec := [op(fvec), tmpcoeff[j+1][i-j+1]]; catch "invalid subscript selector": fvec := [op(fvec), 0]; end try; 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 vec2bvo := proc(v,ord) local i, j, l, deg, f,x,y; x:=ord[1]; y:=ord[2]; l := LinearAlgebra[Dimension](v); f := 0; i:=-1; for deg from 0 to l do for j from deg to 0 by -1 do try f := f + v[i]*x^(j)*y^(deg-j); i := i - 1; catch: break; end try; end do; end do; f; end proc; # This routine computes the matrix which represents the multiplication of # f by a polynomial of total degree k. This is called the Cauchy matrix # by some. # Note: GenerateMatrix is probably not the most efficient way to do this. # Written by John P May bvcauchyo := proc(f, ord,deg) local i, j, td, unknowns, g, eqns, x,y, __v; # Protect the unknowns x:=ord[1]; y:=ord[2]; td := degree(f, {x,y}); unknowns:=[]: g:=0; for i from deg to 0 by -1 do for j from 0 to i do g := g+__v[j,i-j]*x^j*y^(i-j); unknowns := [op(unknowns), __v[j,i-j]]; end do; end do; eqns := expand(f*g); eqns := convert(bv2veco(eqns,ord),list); eqns := map(x->x=0, eqns); LinearAlgebra[GenerateMatrix](eqns, unknowns)[1]; end proc; # This routine call the Cauchy matrix routine to generate the bivariate version # of the Sylvester matrix of f and g. This also called the Macauley Matrix? # The parameter m decreases the size of the matrix. m=1 gives the full # Sylvester-style matrix. # Written by John P May bvsylo := proc(f,g,ord,m) local M, M1, M2, r, c,x,y; x:=ord[1]; y:=ord[2]; Matrix([bvcauchyo(f,ord,degree(g,{x,y})-m), bvcauchyo(g,ord,degree(f,{x,y})-m)]); end proc; # This routine computes the approximate GCD of univariate polynomials. ugcd:=proc(f,g,x,e) local n2,m2,f0,g0,n,m,MM1,M,cc,t,i,MM,PP,QQ,RR,x0,xp,p0,p1,q1,q2; n2:=degree(f,x); m2:=degree(g,x); if n2< m2 then f0:=g; g0:=f; else f0:=f; g0:=g; end if; n:=degree(f0,x); m:=degree(g0,x); MM1:=LinearAlgebra:-SylvesterMatrix(f0,g0,x); M:=LinearAlgebra:-Transpose(MM1); cc:=LinearAlgebra:-SingularValues(M,output='list'); t:=0; for i from 1 to nops(cc) do if evalf(cc[i]/cc[1])<10^(-e) then t:=t+1; fi: end do; if t=0 then return(1); fi; t:=m+1-t; MM:=LinearAlgebra:-DeleteRow(LinearAlgebra:-DeleteColumn(LinearAlgebra:-DeleteColumn(M,[n+t+1..n+m]),[t+1..m]),[n+t+1..n+m]); (PP,QQ,RR):=leastsv(MM,e); x0:=PP[2]; xp:=LinearAlgebra:-SubVector(x0,t+1..n-m+2*t); p0:=convert([seq(x^(n-m+t-i),i=1..n-m+t)],Vector); p1:=LinearAlgebra:-DotProduct(p0,xp,conjugate=false); q1:=polydiv(f0,p1,x); q2:=q1; if degree(q2,x)=max(n2,m2) then RETURN(g0); fi; #return(cutoff(q2,x,e)); return(q2); end proc; # This routine computes approximate GCD of multivariate polynomials. # It is based on Wen-shine Lee's interpolation algorithm multigcd1:=proc(f,g,ord,e) local j,n,S,b,var,ni,pi,si,num,VAR,ww,a,totn,TT,die,w,lcov,coe,f0,g0, fj,gj,cj,SS,PP,VV,res; f0:=evalf(expand(f/norm(expand(f),infinity))); g0:=evalf(expand(g/norm(expand(g),infinity))); n:=nops(ord); S:=[];PP:=[];w:=[];SS:=[]; die:=rand(1..1); for j from 1 to n do b:=zrandom1(n-1); var:=[op({op(ord)} minus {ord[j]}),ord[j]]; fj:=subs({'var[k]=b[k]' $ 'k' = 1..n-1},expand(f0)); gj:=subs({'var[k]=b[k]' $ 'k'=1..n-1},expand(g0)); cj:=ugcd(fj,gj,ord[-1],e); ni:=degree(cj,ord[-1])+3; S:=[op(S),ni]; pi:=nextprime(max(op(S),op(PP))); si:=die(); SS:=[op(SS),si]; w:=[op(w),evalf(exp(si*2*Pi*I/pi))]; PP:=[op(PP),pi]; end do; totn:=testtdeg(f,g,ord,e)[1]; #print(`totn=`,totn); TT:=2*(totn+3)!/(totn!*6)+3; #print(`TT=`,TT); lcov:=mainvar(f,g,ord,e); num:=lcov[1]; VAR:=[seq(ord[k],k=1..num-1),seq(ord[k],k=num+1..n),ord[num]]; PP:=[seq(PP[k],k=1..num-1),seq(PP[k],k=num+1..n),PP[num]]; SS:=[seq(SS[k],k=1..num-1),seq(SS[k],k=num+1..n),SS[num]]; ww:=[seq(w[k],k=1..num-1),seq(w[k],k=num+1..n),w[num]]; if lcov[2]=1 then for j from 0 to TT do fj:=subs({'VAR[k]=ww[k]^j' $ 'k'=1..n-1},expand(f0)); gj:=subs({'VAR[k]=ww[k]^j' $ 'k'=1..n-1},expand(g0)); cj:=ugcd(fj,gj,VAR[n],e); cj:=expand(evalf(cj/lcoeff(cj,VAR[n]))); a[j]:=subs(VAR[n]=ww[n]^j,cj); end do; else for j from 0 to TT do fj:=subs({'VAR[k]=ww[k]^j' $ 'k'=1..n-1},expand(f0)); gj:=subs({'VAR[k]=ww[k]^j' $ 'k'=1..n-1},expand(g0)); cj:=ugcd(fj,gj,VAR[n],e); cj:=expand(evalf(cj/lcoeff(cj,VAR[n]))); coe:=subs({'VAR[k]=ww[k]^j' $ 'k'=1..n-1},lcov[2]); cj:=expand(coe*cj); a[j]:=subs(VAR[n]=ww[n]^j,cj); end do; fi; VV:=[seq(a[j],j=0..TT)]; res:=GEPrimB(SS,PP,VV,VAR); RETURN(res); end proc: # Cut off the terms with coefficient less than 10^(-e) # Input: a polynomial f # Output: a polynomial g without small terms cutoff:=proc(f,ord,e) local l,i,g,v,vv; #lprint(`f=`,f,ord,e); l:=[coeffs(expand(f),ord,'v')]; #lprint(`l=,v=`,l,v); g:=0; vv:=[v]; for i to nops(l) do if abs(l[i])>5*10^(-e) then g:=g+l[i]*vv[i]; #lprint(`g=`,g); fi; od; v:='v'; RETURN(g); end: # The following code choose a main variable which satisfies whose GCD # of the coefficient of the largest degree is 1. mainvar:=proc(f,g,ord,e) local i,var,SS,pp,n; n:=nops(ord); for i from 1 to n do if degree(lcoeff(f,ord[i]))<=0 or degree(lcoeff(g,ord[i]))<=0 then RETURN([i,1]); fi; end do; for i from 1 to n do var:=[op({op(ord)} minus {ord[i]})]; SS:=[]; pp:=multigcd(lcoeff(f,ord[i]),lcoeff(g,ord[i]),var,e); if pp[1]<>1 or degree(pp[2])<=0 then RETURN([i,1]); fi; SS:=[op(SS),[i,pp[2]]]; end do; RETURN(SS[1]); end proc: # The following code is to test the polynomials whether is squarefree. testsquarefree:=proc(f,ord) local i, t,k,fx,a,b,f1,g1,M,er,PP,QQ,RR,e; k:=nops(ord); fx:=diff(f,x); for i from 1 to 3 do b:=zrandom1(k); a:=zrandom1(k); # do the linear maps f1:=expand(subs({'ord[j]=((a[j]+b[j])/2*ord[1]-(-1)^(j)*b[j])' $ 'j' = 1..k},expand(f))); g1:=expand(subs({'ord[j]=( (a[j]+b[j])/2*ord[1]-(-1)^(j)*b[j])' $ 'j' =1..k},expand(fx))); M:=LinearAlgebra:-SylvesterMatrix(f1,g1,ord[1]); if nargs=3 then e:=args[3]; (PP,QQ,RR):=leastsv(M,e); er:=PP[1]; if er>10^(-e) then RETURN(0); fi; fi; (PP,QQ,RR):=leastsv(M,3); er:=PP[1]; if er>10^(-2) then RETURN(0); fi; end do; RETURN(1); end proc; ###################################################################### # The following code is Wen-shin Lee's code #################################################################### ## Filename: GEPrimB.mm ## ## Procedure: ## ## Calling sequences: ## GEPrimB(start_list, prime_list, polyvalue_list, variable_list); ## ## Input: ## start_list - list of integers, determines the starting points of the PRU ## for each variable ## prime_list - list of primes, used for the PRU of each variable ## polyvalue_list - the evaluation list of the target polynomial ## variable_list - an ordered list of variable in the target polynomial ## ## Returns: ## p - a polynomial ## ## Description: ## Compute and round the integer exponents from a Vandermonde system with ## respect to the computed exponents to compute the coefficients. (Round the ## exponents before computing the coefficients.) ## ## See Giesbrecht, Labahn, and Lee. ## ## Author: Wen-shin Lee ## ## Last Update: 7/20/2003 ################################################################################ GEPrimB := proc( start_list, # a list of integers that defines the starting point of each # primitive root of unity prime_list, # a list of integers that are relatively prime polyvalue_list, # a list of target polynomial evaluated at the corresponding # powers of PRU variable_list # the list of the variables in the target polynomial ) local i # index variable ,j # index variable ,t # number of elements in termlist ,termlist # the term value list computed by GEShiftedHankel ,n # number of variables in target polynomial ,m # the product of elements in prime_list, used in the # reversed Chinese remainder theorem ,power_m_list # list of powers of term values in terms of the m-th # primitive root of unity ,term_exp_adjust # the exponent in each variable computed via the # reversed Chinese remainder theorem ,term_exp # a list represent the exponent of a multivariate term # in target polynomial. the exponent in each variable # is listed as the order in variable_list ,term_exp_list # list of exponents of terms in the target polynomial ,vlist # the list describe the Vandermonde system ,V # Vandermonde system ,VT # transposed Vandermonde system ,blist # the list of values for constructing B ,B # the evaluations of target polynomial used for # solving coefficients in Vx = B ,Coeff # the coefficients in the target polynomial ,p # store the returned target polynomial ,term # store a term in the target polynomial : #----------------------------------------------------------------------------- # make sure the input is a list #----------------------------------------------------------------------------- if not(type(start_list,list)) or not(type(prime_list,list)) or not(type(polyvalue_list,list)) or not(type(variable_list,list)) then error "each input argument must be a list"; end if; #----------------------------------------------------------------------------- # check whether start_list, prime_list, and variable_list have the same number # of elements #----------------------------------------------------------------------------- if not(nops(variable_list)=nops(start_list)) then error "first and fourth arguments must have the same number of elements"; else if not(nops(variable_list)=nops(prime_list)) then error "second and fourth arguments must have the same number of elements"; fi; fi; #----------------------------------------------------------------------------- # via generalized eigenvalue procedure, obtain the values of terms evaluated # at the corresponding PRU #----------------------------------------------------------------------------- termlist := GEShiftedHankel(polyvalue_list): #----------------------------------------------------------------------------- # determine the exponent of each term #----------------------------------------------------------------------------- n := nops(prime_list): m := 1: for i from 1 to n do m := m*prime_list[i]: od: t := nops(termlist): power_m_list:=[]: for i from 1 to t do power_m_list := [op(power_m_list),round(ExpLog(termlist[i],m))]: od; term_exp_list:=[]: for i from 1 to t do term_exp := []: for j from 1 to n do term_exp_adjust := RevChRem(power_m_list[i],prime_list[j],m): term_exp_adjust := term_exp_adjust/start_list[j] mod prime_list[j]: term_exp := [op(term_exp),term_exp_adjust]: od: term_exp_list := [op(term_exp_list), term_exp]: od: #----------------------------------------------------------------------------- # obtain the coefficients of the resulting polynomial by solving a transpose # Vandermonde system. This Vandermonde system is based on the evaluations of # the corresponding primitive roots of unity to the integer powers that has # been computed and rounded. vlist provides the list that describes such a # Vandermonde matrix. #----------------------------------------------------------------------------- vlist := []: for i from 1 to t do term := 1: for j from 1 to n do term := evalf(term*exp(2*Pi*I*start_list[j]/prime_list[j])^(term_exp_list[i,j])): od: vlist := [op(vlist),term]: od: V := LinearAlgebra[VandermondeMatrix](vlist): VT := LinearAlgebra[Transpose](V): blist := [op(1..t,polyvalue_list)]: B := convert(blist, Vector): Coeff := LinearAlgebra[LinearSolve](VT,B): #----------------------------------------------------------------------------- # reconstruct the resulting polynomial as p #----------------------------------------------------------------------------- p := 0: for i from 1 to t do term := Coeff[i]: for j from 1 to n do term := term*variable_list[j]^term_exp_list[i,j]: od: p := p + term: od: #----------------------------------------------------------------------------- # userinfo of this procedure #----------------------------------------------------------------------------- userinfo(3, GEPrimB, `entered with:`, start_list, prime_list, polyvalue_list, variable_list, `RETURN:`, p); return(p); end proc: # ################################################################################ ## ## Filename: GEShiftedHankel.mm ## ## Calling sequences: ## GEShiftedHankel(L) ## ## Input: ## L - a list, which provides the information of the generalized eigenvalue ## problem for shifted Hankel systems to be solved. ## ## Returns: ## terms - a list of the solutions of the generalized eigenvalues of the ## given shifted Hankel systems. ## ## Description: ## This procedure is used in the symbolic numeric sparse interpolations that ## use generalized eigenvalue procedure to determine the values of terms (the ## monomials without coefficients) in the target polynomial. If the input list ## has even number (say 2m) elements: the first 2m-1 elements are used to ## form a Hankel matrix H_0, and the last 2n-1 elements a Hankel matrix H_1. ## The procedure returns a list, whose elements are solutions to x in ## H_1 v = x H_0 v. ## The values in the output list, named terms, are actually the term values ## in the target polynomial evaluated at x_1=p_1, ..., x_n=p_m when the input ## list L=[f(1,...,1), f(p_1^1,...,p_n^1), ... f(p_1^{2m-1},...,p_n^{2m-1}]. ## In order to set up two Hankel system for the generalized eigen value ## problem, L needs to be a list of even number elements. If the input L is a ## list of odd number (2m+1) elements, this procedure will only use the first ## 2m elements, and warns that the last element is not used. ## ## Author: Wen-shin Lee ## ## Last Update: 7/7/2003 ################################################################################ GEShiftedHankel := proc( L # a list, provides the generalized eigenvalue problem for two shifted # Hankel systems that is to be solved. ) local t # store the number of elements in L ,hL # store the adjusted list from the input L ,h_0 # a list of values for Hankel matrix H_0 ,h_1 # a list of values for Hankel matrix H_1 ,H_0 # Hankel matrix from list h_0 ,H_1 # Hankel matrix from list h_1 ,terms # the list to be returned : #----------------------------------------------------------------------------- # make sure the input is a list #----------------------------------------------------------------------------- if not(type(L,list)) then error "input must be a list" end if; #----------------------------------------------------------------------------- # adjust the input list L and store the adjusted list to hL: # if there are even number elements in L, then L = hL # if there are odd number elements in L, then remove the last element in L and # store all the remainings in hL #----------------------------------------------------------------------------- t := nops(L): if not(type(t,even)) then print(L[t], `not used`): hL := subsop( t=NULL, L): t := t-1: else hL := L: end if; #----------------------------------------------------------------------------- # from hL, form two lists h_0 and h_1 # from h_0 and h_1, form the corresponding Halkel matrices H_0 and H_1 # then solve for the generalized eigenvalue problem for H_1 v = x H_0 v, the # solutions for x are returned as a list of all solutions. #----------------------------------------------------------------------------- h_0 := subsop( t=NULL, hL ): h_1 := subsop( 1=NULL, hL ): H_0:=LinearAlgebra[HankelMatrix](h_0): H_1:=LinearAlgebra[HankelMatrix](h_1): terms := LinearAlgebra[Eigenvalues](H_1, H_0, output=list); #----------------------------------------------------------------------------- # userinfo of this procedure #----------------------------------------------------------------------------- userinfo(3, GEShiftedHankel, `entered with:`, L, `RETURN:`, terms); return(terms): end proc: float_bb := proc( poly, # a maple polynomial variable_list # a list of the variables in poly ) local blackbox, # neutral form of procedure (black box of poly) numvars; # length of the list of indeterminants numvars := nops(variable_list); blackbox := `&proc`([pntsnf_i, maple_digits], [poly_b,i_b,numvars_b,variable_list_b],[], `&statseq`( `&:=`(`&local`[1], poly), `&:=`(`&local`[3], numvars), `&:=`(`&local`[4], variable_list), `&:=`(`&local`[1], 'evalf(subs({seq(`&local`[4] [`&local`[2]]=evalf(`&args`[1][`&local`[2]], `&args`[2]), `&local`[2]=1..`&local`[3])},`&local`[1]), `&args`[2])'), `&RETURN`(`&local`[1]) ) ); return(procmake(blackbox)); end: RevChRem := proc( f # integer, the input sum in mod M ,m_i # integer, ,M # integer, ) local M_i # M_i = m_1 * m_2 * ... * m_{i-1} * m_{i+1} *... * m_n ,r # store the result of igcdex(M_i, m_i, 's', 't') ,s # coefficient in 1 = s*M_i + t*m_i ,t # coefficient in 1 = s*M_i + t*m_i ,v # f mod m_i ,c_i # the integer coefficient of M_i to be returned : if not(type(m_i,posint)) or not(type(M,posint)) then error "second and third arguments must be positive integers"; fi; v := f mod m_i: M_i := M/m_i: if not(type(M_i, posint)) then error "m_i must divides M"; fi; r := igcdex(M_i, m_i, 's', 't'): if not(r=1) then error "m_i and M/m_i must be relatively prime"; fi; c_i := v*s mod m_i: #----------------------------------------------------------------------------- # userinfo of this procedure #----------------------------------------------------------------------------- userinfo(3, RevChRem, `entered with:`, f, m_i, M, `RETURN:`, c_i); return(c_i): end proc: ExpLog := proc( omega # the input complex number ,m # a number for the m-th primitive root of unity ) local explog # the quotient to be returned ,arg_omega # the argument of the input omega in floating point : #----------------------------------------------------------------------------- # compute the argument of the input omega in floating point #----------------------------------------------------------------------------- arg_omega := evalf(argument(omega)): #----------------------------------------------------------------------------- # determine the ratio of the argument of omega and 2*Pi/m #----------------------------------------------------------------------------- if arg_omega >= 0 then explog := evalf(arg_omega/(2*Pi/m)) else explog := evalf((2*Pi+arg_omega)/(2*Pi/m)): fi; #----------------------------------------------------------------------------- # userinfo of this procedure #----------------------------------------------------------------------------- userinfo(3, ExpLog, `entered with:`, omega, m, `RETURN:`, explog); return(explog); end proc: # The following code is used to compute the content of the polynomial. #e.g. f:=expand((x^3+1)*(y+1));We choose x is main # variable, e is tolerance. we can get output: result:=cont(f,x,[x,y],e)=y+1. cont_mul:=proc( f,ord, x,e) local i,j,k,LL,coeff_i,a,die,tt,fj,gj,var,y,cj,t1,f0,minus_var,min_deg,min_term; LL:=[]; f0:=cutoff(f,ord,e); f0:=evalf(expand(f0/norm(f0,infinity))); minus_var:=[op({op(ord)} minus {x})]; LL:=[coeffs(f0,minus_var,'v')]; min_deg:=degree(f0); min_term:=1; for i from 1 to nops(LL) do if degree(LL[i])<= 0 then if abs(coeff_i)>10^(-e) then RETURN(1); fi; fi; if degree(LL[i])<=min_deg then min_deg:=degree(LL[i]); min_term:=LL[i]; fi; end do; a:=zrandom1(nops(LL)); fj:=min_term; for j from 1 to 3 do gj:=sum('a[k]*LL[k]','k'=1..nops(LL)); #print(`gj=`,gj); cj:=ugcd(fj,gj,x,e); if cj=1 then RETURN(1); fi; fj:=cj; end do; RETURN(cj); end proc; cont:=proc( f,ord, x,e) local i,j,k,LL,coeff_i,a,die,tt,fj,gj,var,y,cj,t1,f0; LL:=[]; f0:=cutoff(f,ord,e); f0:=evalf(expand(f0/norm(f0,infinity))); t1:=degree(f,x); for i from 0 to t1 do coeff_i:=coeff(f0,x,i); if degree(coeff_i)<= 0 then if abs(coeff_i)>10^(-e) then RETURN(1); fi; fi; LL:=[op(LL),coeff_i]; end do; #print(`LL=`,LL); a:=zrandom1(nops(LL)); die:=rand(1..nops(LL)); tt:=die(); fj:=LL[nops(LL)]; for j from 1 to 3 do gj:=sum('a[k]*LL[k]','k'=1..nops(LL)); #print(`gj=`,gj); var:={op(ord)} minus {x}; y:=op(var); #print(`fj=,gj=`,fj); cj:=ugcd(fj,gj,y,e); if cj=1 then RETURN(1); fi; fj:=cj; end do; RETURN(cj); end proc; multifac_cont:=proc(f,ord,e) local f0,x,cont1,prim_f,FF1,FF,VARS; f0:=evalf(expand(f/norm(f,infinity))); x:=ord[2]; cont1:=cont_mul(f0,ord,x,e); prim_f:=mcofactor(f0,cont1,ord,e); prim_f:=cutoff(prim_f,ord,e); #print(`prim_f=`,prim_f); VARS:=[op(indets(prim_f))]; #print(`VARS=`,VARS); if nops({op(VARS)} minus {x}) < nops(VARS) then VARS:=[op({op(VARS)} minus {x}),x]; fi; #print(`VARS=`,VARS); FF1:=multifac(prim_f,VARS); FF:=[[cont1,1],op(FF1)]; RETURN(FF); end proc; # The following code is written by John May. # The code is to compute the error of two polynomails. # Cleaner computation of the square of the 2-norm of a real polynomial norm2 := proc(P::polynom, vars::list) local vec; vec := [coeffs(expand(P), vars)]; vec := map(_x->expand(_x^2), vec); convert(vec, `+`); end proc: # find the coefficient of P that has the largest abs value max_abs_coeff := proc(P::polynom(complex)) local vars, tmp, tmp2, bc, loc; vars := convert(indets(P), list); tmp := [coeffs(expand(P), vars)]; tmp2 := map(abs, tmp); bc := max(op(tmp2)); member(bc, tmp2, 'loc'); tmp[loc]; end proc: # Compute the constant that minimizes ||Porg - c*Papp||_2 minimum_scaling_real := proc(Porg::polynom, Papp::polynom, vars::list) local nm, _c, bc: #_c:='_c'; bc := max_abs_coeff(Papp); nm:=norm2(Porg - _c*1/bc*Papp, vars); 1/bc * solve(diff(nm, _c)=0, _c); end proc: # Compute the constant that minimizes ||Porg - c*Papp||_2 for complex P's minimum_scaling_complex := proc(Porg::polynom, Papp::polynom, vars::list) local nm, _c, _d, bc, tmp: bc := max_abs_coeff(Papp); assume(_c, real); assume(_d, real); nm:=expand(norm(expand(Porg - (_c + _d*I)*1./bc*Papp), 2, vars)^2); tmp:=solve({diff(nm, _c)=0, diff(nm, _d)=0}, {_c,_d}); 1/bc * (subs(tmp,_c) + subs(tmp,_d)*I); end proc: # Get rid of .0*I terms from the coefficients of P strip_ext_complex := proc(P::polynom(complex)) local F, newF, vars, trm, cof, t; vars:=convert(indets(P), list); F:=expand(1/max_abs_coeff(P) * P); F:=convert(F,list); newF:=0; for trm in F do cof:=lcoeff(trm, vars, 't'); if Im(cof) = 0 then newF := newF + Re(cof)*t; else newF := newF + trm; end if; end do; end proc: # Find a common monomial in the support of two polynomials common_monomial := proc(P1::polynom, P2::polynom) local vars, terms1, terms2; vars:=indets(P1); vars:=indets(P1) union indets(P2); coeffs(expand(P1), vars, 'terms1'); coeffs(expand(P2), vars, 'terms2'); {terms1} intersect {terms2}; end proc: # Find a constant so that c*Papp has a monomial in common with Porg # returns 0 if no common monomials commmon_mon_scaling := proc(P1::polynom, P2::polynom, vars::list) local mon, v, Fapp, Forg; mon := common_monomial(Forg, Fapp); if mon={} then return(0) end if; mon:=RandomTools[Generate](choose(mon)); print('common_monomial'=mon); Forg:=expand(Porg); Fapp:=expand(Papp); for v in vars do Fapp := coeff(Fapp, v, degree(mon, v)); Forg := coeff(Forg, v, degree(mon, v)); end do; Forg/Fapp; end proc: # Compute the (2-norm) relative backward error of the approximate factorization Papp backward_error := proc(Porg::polynom, Papp::polynom, vars::list) local Forg, Fapp, c, mon, v; Forg:=expand(Porg); Fapp:=expand(Papp); if not common_monomial(Forg, Fapp)={} then if type(Forg, polynom(numeric)) and (type(Fapp, polynom(numeric)) or type(strip_ext_complex(Fapp), polynom(numeric))) then c:=minimum_scaling_real(Porg, Papp, vars); else printf("Approximate factorization is not a real polynomial"); c:=minimum_scaling_real(Porg, Papp, vars); end if; else printf(" and "); c := 0; end if; if c = 0 then printf("there are no common monomials"); return( min( norm(expand(Forg-norm(Forg,2,vars)/norm(Fapp,2,vars)*Fapp),2,vars), norm(expand(Forg+norm(Forg,2,vars)/norm(Fapp,2,vars)*Fapp),2,vars) )/norm(Forg,2,vars) ); end if; print(('c') = c); #print(('abs_error') = norm(expand(Forg-c*Fapp),2,vars)); return norm(expand(Forg-c*Fapp), 2, vars) / norm(1.*Forg, 2, vars); end proc: homogenous:= proc( f, # the input polynomial ord, # the variables in f x # the variable ) local g # the homogenous polynomial to be returned , d # the total degree of f , i # the number of terms in f , L # the list of coefficients in f : #----------------------------------------------------------------------------- # Homogenize the polynomial f w.r.t. x #----------------------------------------------------------------------------- g:=0; d:=degree(f); L:=[coeffs(expand(f), ord, 't')]; for i from 1 to nops(L) do g:=g+L[i]*x^(d-degree(t[i]))*t[i]; od; return(g); end proc: