#------------------------------------------------------------------------------------------------------------------------- # author : Zhengfeng Yang # date : Jan. 15, 2006. #-------------------------------------------------------------------------------------------------------------------------- ########################################################################################################################### ## Description: ## Given the polynomial list L:=[f1,f2,....fr]\in C[x1,x2...xn] and given the integer k, we want to compute the $tilde{L}$=[\tilde{f1},\tilde{f2},....\tilde{fn}]$ ## such that degree($\tilde{f1},\tilde{f2},....\tilde{fn})$=k,and $\tilde{f1},\tilde{f2},....\tilde{fn} satisfy the linear constraints. ## This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ ## generated by $\tilde{f1},\tilde{f2},....\tilde{fn}]$, and rank deficiency of $tilde(M)$ is 1. ## Here, $\tilde{f1},\tilde{f2},....\tilde{fn}]$ in C field. ## In C_con_multi, M=[B,A], we want to compute the matrix res=[f+B,A+E] such that $f+B \in Range(A+E)$. ## M is a $k$th- multiple polynomials Sylvester matrix by f1,f2,....fn ## B is a column of M ## A is a remaining columns of M ## Input: ## L=[f1,f2,....fn] ---- the list of the polynomials ## k ---- degree of gcd of the polynomials; ## e ---- the tolerance ## Con_L ----- the matrix of linear constraints ## DD ---- the weight matrix ## ord ---- the variables ## Output ## pp ---- gcd of res_L ## prim_list ---- the list of the cofactors ############################################################################################### C_con_multi:=proc(L,ord,k,e,Con_L,DD) local M # the Sylvester matrix ,zz # the vectors ,A # the submatrices of the matrix M ,b # the column of M ,r # the number of the polynomials ,deg_L # the list of the degree of the polynomials ,LL # the list of the polynomials ,La # the list of the perturbed polynomials ,var_dim # the number of the free variables ,res_L # the result of the polynomials list ,prim_list # the list of the cofactors ,num_org # the number of the vectors ,m # the degree of the polynomials ,N # the integer ,V # the vector ,int_vec # the vector ,gcd_res # the gcd of res_L ,X # the least squares solution of A.X=b ,ww # the weight value ,P0_Matrix # the matrix such that $ P0_Matrix.zz$ is the column of M ,XY # the matrix ,delta_x,f_A # the vectors ,Aa,A_R,A_I # the submatrices of the matrix M ,ba # the column of M ,fa,ga # the polynomials ,Ma # the k-th Sylvester matrix ,zz_R,zz_I # the real part and the image part of zz ,Con_M # the matrix of the linear constraints ,dd # the vector of the linear constraints ,zeta # the vectors ,zeta_R,zeta_I # the real part and the image part of zeta ,Con_M_R,Con_M_I # the real part and the image part of Con_M ,D_R,D_I # the real part and the image part of the matrix DD ,dd_R,dd_I # the real part and the image part of ,f_var # the dimension of the free vector ,zeta_A_R,zeta_A_I # the real part and the image part of the vector zeta_A ,b_R,b_I # the real part and the image part of b ,EE_R,EE_I #the real part and the image part of EE ,A_A_R,A_A_I #the real part and the image part of A_A ,R_R,R_I # the real part and the image part of R ,A_A_sub1,A_A_sub2 # the submatrix of the matrix A_A ,A_A_sub # the submatrix of the matrix A_A ,XY_sub1,XY_sub2 # the submatrix of the matrix XY ,XY_sub # the submatrix of the matrix XY ,MM_sub1,MM_sub2 # the submatrix of the matrix MM ,MM_sub # the submatrix of the matrix MM ,zz_A_R,zz_A_I # the real part and the image part of zz_A ,X_A_R,X_A_I #the real part and the image part of X_A ,X_R,X_I # the least squares solution of A.X=b ,XY_R,XY_I # the matrix ,f_A_R,f_A_I # the vector ,num # the number of the iterations ,EE,E # the matrices ,i,kk,uu,u,v,vv,j,jj # the parameter in for loop statement ,XYY # the matrix ,M_M # the whole matrix ,M_sub,M_sub2 # the sub matrix of M_M ,row_dim # the row dimension of the matrix ,col_dim # the column dimension of the matrix ,ff,gg # the result we get which satisfy GCD(ff,gg)=k ,var # the list of the variables ,m1,m2,n1,n2 # the degrees of the polynomials. ,n # the number of the variables ,error_orgin # the error of $b-A.X$. ,pp # the null vector of M ,vvv,v_sub # the subvector of pp ,vec_org # the vector by the coefficients of f and g ,vec_res # the vector by the coefficients of ff and gg ,UU,_SS,_VV # the singular values decomposition ,num_S # rank deficiency of M_M ,M_new_inverse # the inverse matrix of M_M ,mm,nn # binomial number ,I_matrix # the identity matrix ,zero_matrix # the zero matrix ,X_int # the matrix ,M_tem # the matrix ,elta_s # the vector ,R # the residual i.e. b-A.X ,A_A # the matrix which denotes the change of the matrix A ,D_elta # the vector ,X_A_Y # the solution of LS problem ,X_A # the vector which denotes the change of the vector X ,zz_A # the vector which denotes the change of the vector zz ,res_M # the matrix of the k-th sylvester matrix of ff and gg ,_U,_S,_V; # the singular values decomposition of the matrix n:=nops(ord); r:=nops(L); deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(L[i],ord)]; end do; member(min(op(deg_L)),deg_L, 'nd'); if nd=r then LL:=L; elif nd=1 then LL:=[seq(L[i],i=2..r),L[1]]; else LL:=[seq(L[i],i=1..nd-1),seq(L[i],i=nd+1..r),L[nd]]; fi; LL:=map(x->evalf(x),LL); deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(LL[i])]; end do; vec_org:=con_vec_fg(LL,ord); # we compute the matrix of linear constraints if args[5]='default' then Con_M:=LinearAlgebra:-IdentityMatrix(nops(vec_org)); Con_M_R:=convert(map(x->Re(x),convert(Con_M,matrix)),Matrix); Con_M_I:=convert(map(x->Im(x),convert(Con_M,matrix)),Matrix); dd:=Vector(nops(vec_org)); else Con_M:=LinearAlgebra:-SubMatrix(Con_L,1..-1,1..-2); Con_M_R:=map(x->Re(x),Con_M); Con_M_I:=map(x->Im(x),Con_M); dd:=LinearAlgebra:-Column(Con_L,-1); fi; ## construct the $k$-th Sylvster matrix, and choose the column $b$ such that $A.X appxo = b$. M:=multi_generate_syl(LL,ord,k); pp:= LinearAlgebra[SingularValues](M, output=['Vt']); vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1); v_sub:=LinearAlgebra:-SubVector(vv,1..(deg_L[1]+n-k)!/((deg_L[1]-k)!*n!)): v_sub:=convert(v_sub,list): member(max(op(map(x ->abs(x),v_sub))), map(x->abs(x),v_sub),'choo_k'): A:=LinearAlgebra:-DeleteColumn(M,choo_k); b:=LinearAlgebra:-Column(M,choo_k); error_orgin:=LinearAlgebra:-Norm(M.vv,infinity); # We add complex randomized preconditioning of the inputs if LinearAlgebra:-Norm(map(x->Im(x),convert(vec_org,Vector)+dd),infinity)<=10^(-4) or error_orgin> 10^(-3) then int_vec:=0.1*error_orgin*I*Vector([seq(LinearAlgebra[RandomVector][row](LinearAlgebra:-ColumnDimension(Con_M),generator=-0.5..0.5)[j],j=1..LinearAlgebra:-ColumnDimension(Con_M))]); La:=con_pol_fg(LL,ord,convert(dd,list)+vec_org+convert(Con_M.int_vec,list)); Ma:=multi_generate_syl(La,ord,k); else La:=con_pol_fg(LL,ord,dd+convert(vec_org,Vector)); Ma:=multi_generate_syl(La,ord,k); fi; pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']); vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1); # we compute the initial vector $X$ via the Lagrangian multiplier method, then do the iteration X_int:=Int_point_matrix(vv,LL,k,ord).Con_M; M_tem:=X_int.LinearAlgebra:-HermitianTranspose(X_int); UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_tem,output=['U','S', 'Vt']); num_S:=0; for i from 1 to LinearAlgebra:-ColumnDimension(M_tem) do if _SS[LinearAlgebra:-ColumnDimension(M_tem)-i+1]1/x,LinearAlgebra:-SubVector(_SS,1..LinearAlgebra:-ColumnDimension(M_tem)-num_S)),Vector(num_S)]); _SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_tem),LinearAlgebra:-RowDimension(M_tem)); M_new_inverse:=LinearAlgebra:-HermitianTranspose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-HermitianTranspose(UU),1..LinearAlgebra:-ColumnDimension(M_tem),1..-1); elta_s:=-LinearAlgebra:-HermitianTranspose(X_int).M_new_inverse.Ma.vv; X:=-vv/vv[choo_k]; X:=convert(LinearAlgebra:-DeleteRow(convert(X,Matrix),choo_k),Vector); error_orgin:=LinearAlgebra:-Norm(b-A.X,infinity); userinfo(1,C_con_multi,print('error_orgin'=error_orgin)); row_dim:=LinearAlgebra:-RowDimension(A); col_dim:=LinearAlgebra:-ColumnDimension(A); var_dim:=nops(vec_org); # we get the weight matrix. if args[6]='default' then D_R:=LinearAlgebra:-IdentityMatrix(var_dim); D_I:=Matrix(var_dim,var_dim); else D_R:=map(x->Re(x),DD); D_I:=convert(map(x->Im(x),convert(DD,matrix)),Matrix); fi; f_var:=LinearAlgebra:-ColumnDimension(Con_M); zeta:=elta_s; zeta_R:=map(x->Re(x),zeta); zeta_I:=map(x->Im(x),zeta); X_R:=map(x->Re(x),X); X_I:=map(x->Im(x),X); A_R:=map(x->Re(x),A); A_I:=map(x->Im(x),A); b_R:=map(x->Re(x),b); b_I:=map(x->Im(x),b); ww:=evalf(10^(min(4-ceil(log10(error_orgin)),6))); P0_Matrix:=multi_p_matrix(LL,ord,k,choo_k); delta_x:=1; I_matrix:=Matrix([D_R.Con_M_R-D_I.Con_M_I,-D_R.Con_M_I-D_I.Con_M_R]); I_matrix:=Matrix([[I_matrix],[Matrix([D_R.Con_M_I+D_I.Con_M_R,D_R.Con_M_R-D_I.Con_M_I])]]); zero_matrix:=Matrix(2*(var_dim),2*col_dim); zz:=Con_M.(zeta_R+I*zeta_I)+dd; zz_R:=map(x->Re(x),zz); zz_I:=map(x->Im(x),zz); dd_R:=map(x->Re(x),dd); dd_I:=map(x->Im(x),dd); f_A:=P0_Matrix.zz; f_A_R:=map(x->Re(x),f_A); f_A_I:=map(x->Im(x),f_A); num:=0; # the following, we compute the $\delta f$ and $\delta E$ such that # $b+\delta f \in Range(A+\delta E)$ while delta_x>10^(-e) do num:=num+1; if num>300 then break; fi; ## we construct the matrix $E$ and $Y_{k}$ in our paper, according to X and zz. EE_R:=mul_construct_E(zz_R,LL,k,choo_k,ord); EE_I:=mul_construct_E(zz_I,LL,k,choo_k,ord); XY_R:=mul_construct_XK(X_R,LL,k,choo_k,ord); XY_I:=mul_construct_XK(X_I,LL,k,choo_k,ord); A_A_R:=A_R+EE_R; A_A_I:=A_I+EE_I; R_R:=b_R+f_A_R-A_A_R.X_R+A_A_I.X_I; R_I:=b_I+f_A_I-A_A_R.X_I-A_A_I.X_R; ## we construct the large matrix which is used to do STLN userinfo(3,C_con_multi,print('num'=num)); userinfo(3,C_con_multi,print(`norm of R=`,LinearAlgebra:-Norm(R_R+I*R_I,2))); if LinearAlgebra:-Norm(R_R+I*R_I,2)<10^(-10) then break; fi; A_A_sub1:=Matrix([A_A_R,-A_A_I]); A_A_sub2:=Matrix([A_A_I,A_A_R]); A_A_sub:=Matrix([[A_A_sub1],[A_A_sub2]]); XY_sub1:=Matrix([XY_R.Con_M_R-XY_I.Con_M_I-P0_Matrix.Con_M_R,-(XY_R.Con_M_I+XY_I.Con_M_R-P0_Matrix.Con_M_I)]); XY_sub2:=Matrix([XY_I.Con_M_R+XY_R.Con_M_I-P0_Matrix.Con_M_I,XY_R.Con_M_R-XY_I.Con_M_I-P0_Matrix.Con_M_R]); XY_sub:=Matrix([[XY_sub1],[XY_sub2]]); M_sub:=ww*Matrix([A_A_sub,XY_sub]); M_sub2:=Matrix([zero_matrix,I_matrix]); M_M:=Matrix([[M_sub],[M_sub2]]); D_elta:=Vector([ww*R_R,ww*R_I,-(D_R.zz_R-D_I.zz_I),-(D_R.zz_I+D_I.zz_R)]); X_A_Y:=LinearAlgebra:-LeastSquares(M_M,D_elta,optimize); delta_x:=LinearAlgebra:-Norm(X_A_Y,infinity); userinfo(3,C_con_multi,print('delta_x'=delta_x)); zeta_A_R:=LinearAlgebra:-SubVector(X_A_Y,2*col_dim+1..2*col_dim+f_var); zeta_A_I:=LinearAlgebra:-SubVector(X_A_Y,2*col_dim+f_var+1..2*col_dim+2*f_var); X_A_R:=LinearAlgebra:-SubVector(X_A_Y,1..col_dim); X_A_I:=LinearAlgebra:-SubVector(X_A_Y,col_dim+1..2*col_dim); X_R:=X_R+X_A_R; X_I:=X_I+X_A_I; zeta_R:=zeta_R+zeta_A_R; zeta_I:=zeta_I+zeta_A_I; zz_R:=Con_M_R.zeta_R-Con_M_I.zeta_I+dd_R; zz_I:=Con_M_R.zeta_I+Con_M_I.zeta_R+dd_I; ## compute the first column of E. f_A_R:=P0_Matrix.zz_R; f_A_I:=P0_Matrix.zz_I; end do; # end the iterative method userinfo(1,C_con_multi,print(`number of the iterations` =num)); # compute the matrix E,res_M, ff,gg. E:=mul_construct_E(zz_R+I*zz_I,LL,k,choo_k,ord); E:=E+A; X:=X_R+I*X_I; vec_res:=vec_org+convert(zz_R+I*zz_I,list); res_L:=con_pol_fg(LL,ord,vec_res); # we compute the cofactor vector. if choo_k=1 then X:=Vector([-1,X]); else X:=Vector([LinearAlgebra:-SubVector(X,1..choo_k-1),-1,LinearAlgebra:-SubVector(X,choo_k..-1)]); fi; prim_list:=[]; num_org:=0; for j from 1 to r do ff:=0; m:=degree(LL[j],ord)-k; N:=(m+n)!/(m!*n!); for i from num_org+1 to num_org+N do V:=combinat[inttovec](N-i+num_org,n); ff:=ff+X[i]*mul(ord[kkk]^V[kkk],kkk=1..n); end do; if j=r then ff:=-ff; fi; prim_list:=[op(prim_list),ff]; num_org:=num_org+N; end do; gcd_res:=multi_polydiv(res_L,prim_list,k,ord); if nd=r then res_L:=res_L; prim_list:=prim_list; elif nd=1 then res_L:=[res_L[r],seq(res_L[i],i=1..r-1)]; prim_list:=[prim_list[r],seq(prim_list[i],i=1..r-1)]; else res_L:=[seq(res_L[i],i=1..nd-1),res_L[r],seq(res_L[i],i=nd..r-1)]; prim_list:=[seq(prim_list[i],i=1..nd-1),prim_list[r],seq(prim_list[i],i=nd..r-1)]; fi; if nargs= 7 then if args[7]='cofa' then return gcd_res,prim_list; else return gcd_res; fi; else return gcd_res; fi; end proc; ############################################################################################################################# ## Description: ## Given the polynomial list L:=[f1,f2,....fr]\in R[x1,x2...xn] and given the integer k, we want to compute the $tilde{L}$=[\tilde{f1},\tilde{f2},....\tilde{fn}]$ ## such that degree($\tilde{f1},\tilde{f2},....\tilde{fn})$=k,and $\tilde{f1},\tilde{f2},....\tilde{fn} satisfy the linear constraints. ## This procedure is to compute the $k$-th multiple polynomials Sylvester matrix $tilde(M)$ ## generated by $\tilde{f1},\tilde{f2},....\tilde{fn}]$, and rank deficiency of $tilde(M)$ is 1. ## Here, $\tilde{f1},\tilde{f2},....\tilde{fn}]$ in R field. ## In R_con_multi, M=[B,A], we want to compute the matrix res=[f+B,A+E] such that $f+B \in Range(A+E)$. ## M is a $k$th- multiple polynomials Sylvester matrix by f1,f2,....fn ## B is a column of M ## A is a remaining columns of M ## Input: ## L=[f1,f2,....fn] ---- the list of the polynomials ## k ---- degree of gcd of the polynomials; ## e ---- the tolerance ## Con_L ----- the matrix of linear constraints ## DD ---- the weight matrix ## ord ---- the variables ## Output ## pp ---- gcd of res_L ## prim_list ---- the list of the cofactors ############################################################################################### R_con_multi:=proc(L,ord,k,e,Con_L,DD) local M # the Sylvester matrix ,zz # the vectors ,A # the submatrices of the matrix M ,b # the column of M ,r # the number of the polynomials ,deg_L # the list of the degree of the polynomials ,LL # the list of the polynomials ,var_dim # the number of the free variables ,res_L # the result of the polynomials list ,X_int # the matrix ,M_tem # the matrix ,elta_s # the vector ,prim_list # the list of the cofactors ,num_org # the number of the vectors ,m # the degree of the polynomials ,N # the integer ,V # the vector ,gcd_res # the gcd of res_L ,X # the least squares solution of A.X=b ,ww # the weight value ,DD_test # the weight matrix ,Con_M # the matrix of the linear constraints ,dd # the vector of the linear constraints ,zeta # the vector ,P0_Matrix # the matrix such that $ P0_Matrix.zz$ is the column of M ,XY # the matrix ,delta_x,f_A # the vectors ,A_A_sub1,A_A_sub2 # the submatrix of the matrix A_A ,A_A_sub # the submatrix of the matrix A_A ,XY_sub1,XY_sub2 # the submatrix of the matrix XY ,XY_sub # the submatrix of the matrix XY ,MM_sub1,MM_sub2 # the submatrix of the matrix MM ,MM_sub # the submatrix of the matrix MM ,num # the number of the iterations ,EE,E # the matrices ,i,kk,uu,u,v,vv,j,jj # the parameter in for loop statement ,XYY # the matrix ,M_M # the whole matrix ,M_sub,M_sub2 # the sub matrix of M_M ,row_dim # the row dimension of the matix ,Aa,Ma # the matrices ,ba # the column of the matrix Ma ,La # the list of the polynomials ,f_var # the free variables ,zeta_A # the vector ,col_dim # the column dimension of the matrix ,ff,gg # the result we get which satisfy GCD(ff,gg)=k ,var # the list of the variables ,m1,m2,n1,n2 # the degrees of the polynomials. ,n # the number of the variables ,error_orgin # the error of $b-A.X$. ,pp # the null vector of M ,vvv,v_sub # the subvector of pp ,vec_org # the vector by the coefficients of f and g ,vec_res # the vector by the coefficients of ff and gg ,UU,_SS,_VV # the singular values decomposition ,num_S # rank deficiency of M_M ,M_new_inverse # the inverse matrix of M_M ,mm,nn # binomial number ,I_matrix # the identity matrix ,zero_matrix # the zero matrix ,R # the residual i.e. b-A.X ,A_A # the matrix which denotes the change of the matrix A ,D_elta # the vector ,X_A_Y # the solution of LS problem ,X_A # the vector which denotes the change of the vector X ,zz_A # the vector which denotes the change of the vector zz ,res_M # the matrix of the k-th sylvester matrix of ff and gg ,_U,_S,_V; # the singular values decomposition of the matrix n:=nops(ord); r:=nops(L); deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(L[i],ord)]; end do; member(min(op(deg_L)),deg_L, 'nd'); if nd=r then LL:=L; elif nd=1 then LL:=[seq(L[i],i=2..r),L[1]]; else LL:=[seq(L[i],i=1..nd-1),seq(L[i],i=nd+1..r),L[nd]]; fi; LL:=map(x->evalf(x),LL); deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(LL[i])]; end do; vec_org:=con_vec_fg(LL,ord); # we get the matrix of linear constraints. if args[5]='default' then Con_M:=LinearAlgebra:-IdentityMatrix(nops(vec_org)); dd:=Vector(nops(vec_org)); else Con_M:=LinearAlgebra:-SubMatrix(Con_L,1..-1,1..-2); dd:=LinearAlgebra:-Column(Con_L,-1); fi; ## construct the $k$-th Sylvster matrix, and choose the column $b$ such that $A.X appxo = b$. M:=multi_generate_syl(LL,ord,k); La:=con_pol_fg(LL,ord,dd+convert(vec_org,Vector)); Ma:=multi_generate_syl(La,ord,k); pp:= LinearAlgebra[SingularValues](Ma, output=['Vt']); vv:=LinearAlgebra:-Column(LinearAlgebra:-HermitianTranspose(pp),-1); v_sub:=LinearAlgebra:-SubVector(vv,1..(deg_L[1]+n-k)!/((deg_L[1]-k)!*n!)): v_sub:=convert(v_sub,list): member(max(op(map(x ->abs(x),v_sub))), map(x->abs(x),v_sub),'choo_k'): error_orgin:=LinearAlgebra:-Norm(Ma.vv,infinity); # we compute the initial vector $X$ via the Lagrangian multiplier method, then do the iteration X_int:=Int_point_matrix(vv,LL,k,ord).Con_M; M_tem:=X_int.LinearAlgebra:-HermitianTranspose(X_int); UU,_SS,_VV:=LinearAlgebra:-SingularValues(M_tem,output=['U','S', 'Vt']); num_S:=0; for i from 1 to LinearAlgebra:-ColumnDimension(M_tem) do if _SS[LinearAlgebra:-ColumnDimension(M_tem)-i+1]1/x,LinearAlgebra:-SubVector(_SS,1..LinearAlgebra:-ColumnDimension(M_tem)-num_S)),Vector(num_S)]); _SS:=LinearAlgebra:-DiagonalMatrix(_SS,LinearAlgebra:-RowDimension(M_tem),LinearAlgebra:-RowDimension(M_tem)); M_new_inverse:=LinearAlgebra:-HermitianTranspose(_VV).LinearAlgebra:-SubMatrix(_SS.LinearAlgebra:-HermitianTranspose(UU),1..LinearAlgebra:-ColumnDimension(M_tem),1..-1); elta_s:=-LinearAlgebra:-HermitianTranspose(X_int).M_new_inverse.Ma.vv; X:=-vv/vv[choo_k]; X:=convert(LinearAlgebra:-DeleteRow(convert(X,Matrix),choo_k),Vector); A:=LinearAlgebra:-DeleteColumn(M,choo_k); b:=LinearAlgebra:-Column(M,choo_k); error_orgin:=LinearAlgebra:-Norm(b-A.X,infinity); userinfo(1,R_con_multi,print('error_orgin'=error_orgin)); row_dim:=LinearAlgebra:-RowDimension(A); col_dim:=LinearAlgebra:-ColumnDimension(A); var_dim:=nops(vec_org); f_var:=LinearAlgebra:-ColumnDimension(Con_M); zeta:=elta_s; P0_Matrix:=multi_p_matrix(LL,ord,k,choo_k); ## get the weight matrix. if args[6]='default' then DD_test:=LinearAlgebra:-IdentityMatrix(var_dim); else DD_test:=DD; fi; I_matrix:=DD_test.Con_M; zero_matrix:=Matrix(var_dim,col_dim); ww:=evalf(10^(min(4-ceil(log10(error_orgin)),6))); delta_x:=1; zz:=Con_M.zeta+dd; f_A:=P0_Matrix.zz; num:=0; # the following, we compute the $\delta f$ and $\delta E$ such that # $b+\delta f \in Range(A+\delta E)$ while delta_x>10^(-e) do num:=num+1; if num>100 then break; fi; ## we construct the matrix $E$ and $Y_{k}$ in our paper, according to X and zz. EE:=mul_construct_E(zz,LL,k,choo_k,ord); XY:=mul_construct_XK(X,LL,k,choo_k,ord); XY:=XY.Con_M; XY:=XY-P0_Matrix.Con_M; ## we construct the large matrix which is used to do STLN A_A:=A+EE; R:=b+f_A-A_A.X; userinfo(3,R_con_multi,print('num'=num)); userinfo(2,R_con_multi,print('R'=LinearAlgebra:-Norm(R,infinity))); M_sub:=ww*Matrix([A_A,XY]); M_sub2:=Matrix([zero_matrix,I_matrix]); M_M:=Matrix([[M_sub],[M_sub2]]); D_elta:=Vector([ww*R,-DD_test.zz]); X_A_Y:=LinearAlgebra:-LeastSquares(M_M,D_elta); delta_x:=LinearAlgebra:-Norm(X_A_Y,infinity); userinfo(2,R_con_mulpoly,print('delta_x'=delta_x)); zeta_A:=LinearAlgebra:-SubVector(X_A_Y,col_dim+1..-1); X_A:=LinearAlgebra:-SubVector(X_A_Y,1..col_dim); X:=X+X_A; zeta:=zeta+zeta_A; zz:=Con_M.zeta+dd; # compute the num_b-th column of E. f_A:=P0_Matrix.zz; end do; # end the iterative method userinfo(1,R_con_multi,print(`number of the iterations` =num)); E:=mul_construct_E(zz,LL,k,choo_k,ord); E:=E+A; vec_res:=vec_org+convert(zz,list); res_L:=con_pol_fg(LL,ord,vec_res); # compute the cofactor vector X. if choo_k=1 then X:=Vector([-1,X]); else X:=Vector([LinearAlgebra:-SubVector(X,1..choo_k-1),-1,LinearAlgebra:-SubVector(X,choo_k..-1)]); fi; prim_list:=[]; num_org:=0; for j from 1 to r do ff:=0; m:=degree(LL[j],ord)-k; N:=(m+n)!/(m!*n!); for i from num_org+1 to num_org+N do V:=combinat[inttovec](N-i+num_org,n); ff:=ff+X[i]*mul(ord[kkk]^V[kkk],kkk=1..n); end do; if j=r then ff:=-ff; fi; prim_list:=[op(prim_list),ff]; num_org:=num_org+N; end do; gcd_res:=multi_polydiv(res_L,prim_list,k,ord); if nd=r then res_L:=res_L; prim_list:=prim_list; elif nd=1 then res_L:=[res_L[r],seq(res_L[i],i=1..r-1)]; prim_list:=[prim_list[r],seq(prim_list[i],i=1..r-1)]; else res_L:=[seq(res_L[i],i=1..nd-1),res_L[r],seq(res_L[i],i=nd..r-1)]; prim_list:=[seq(prim_list[i],i=1..nd-1),prim_list[r],seq(prim_list[i],i=nd..r-1)]; fi; if nargs= 7 then if args[7]='cofa' then return gcd_res,prim_list; else return gcd_res; fi; else return gcd_res; fi; end proc; ##################################################################################################################################### ## The following code is to compute the approximate GCD of two multivariate polynomials ## Input-- ## F,G -- the polynomials ## ord -- the list of variables ## e -- the tolerance ## X -- the vector of the cofactors ## Output-- ## a list= [v,da*L1,fbara,gbara] ## where v is an integer describes if the backward error is big. ## da*L1 is the GCD ## fbara is the prime part of the polynomial F ## gbara is the prime part of the polynomial G ################################################################################################################################### multigcd_M:= proc(F::polynom, G::polynom,ord,k,X ) local n # the number of variables ,QQ,RR # QR decomposition ,SS,PP,KK # the least singular value and the vector ,f,g,f0,g0 # polynomials expanded ,i,kk,j # loop variables ,Mg # bvcauchyo matrix ,Mf # bvcauchyo matrix ,Mfg # the matrix constructed by Mf and Mg ,x0 # the null space corresponding to the least singular values ,gbara # subvector of x0 ,fbara # subvector of x0 ,tmp # degree of GCD ,da # GCD ,len1 # the number terms in a polynomial ,gv # the vector converted by polynomial g ,fv # the vector converted by polynomial f ,u_vec # the vector converted by the prime part of f ,v_vec # the vector converted by the prime part of g ,fv_leading # the leading coefficient of the polynomial f ,gv_leading # the leading coefficient of the polynomial g ,u_leading # the leading coefficient of the polynomial u ,v_leading # the leading coefficient of the polynomial v ,dv # the solution of Mg.dv=gv , er # the backward error ,L1 # special factors x, y ,fgv; # the vector by fv and gv f:=expand(F); g:=expand(G); n:=nops(ord); x0:=X; tmp:=degree(g)-k: fv:=convec(f,ord); gv := convec(g,ord); len1:=(n+tmp)!/(n!*tmp!); v_vec:=x0[1..len1]; u_vec:=x0[len1+1..-1]; tmp:=degree(f)-k: fbara:=-veccoe(u_vec,ord); gbara:=veccoe(v_vec,ord); if n=2 then Mg:=bvcauchyo(gbara,ord,degree(g)-degree(gbara)); else Mg := mvcauchyo(gbara,ord,degree(g)-degree(gbara)); fi; if n=2 then Mf:=bvcauchyo(fbara,ord,degree(f)-degree(fbara)); else Mf:=mvcauchyo(fbara,ord,degree(f)-degree(fbara)); fi; Mfg:=Matrix([[Mf],[Mg]]); fgv:=Vector([fv,gv]); dv:=LinearAlgebra:-LeastSquares(Mfg,fgv); da := veccoe(dv,ord); er:=LinearAlgebra[Norm](Mfg.dv-fgv,2); userinfo(1,multigcd_M,print(`check the error`, er)); return da,fbara,gbara; end proc: ############################################################################################################################# ## The following code is to compute the matrix $P_{k}$ as in the paper, when the column $b$ is not the first ## column of the k-th Sylvester matrix. In this code the number of variables >=2 ## Input-- ## f,g --the polynomials ## num_b-- the number of the column in the matrix. ## k-- the degree of the GCD of polynomials ## Output-- ## res_M-- the matrix $P_{k}$ such that $P_{k}.\zz=b$ #################################################################################################################################### pk_matrix:=proc(f,g,ord,k,num_b) local n,m1,m2,row_dim,col_dim,S,SS,i,SK,M,N1,col_num,deg_i; n:=nops(ord); if n=2 then return pk_matrix_2(f,g,ord,k,num_b); fi; m1:=degree(f,ord); m2:=degree(g,ord); row_dim:=(m1+m2-k+n)!/((m1+m2-k)!*n!); col_dim:=(m1+n)!/(m1!*n!)+(m2+n)!/(m2!*n!); M:=Matrix(row_dim,col_dim); N1:=(m2-k+n)!/((m2-k)!*n!); S:=combinat[inttovec](N1-num_b,n); deg_i:=add( i, i=S); for i from row_dim-(m1+deg_i+n)!/((m1+deg_i)!*n!)+1 to row_dim do SS:=combinat[inttovec](row_dim-i,n); if min(op(SS-S))>=0 then SK:=SS-S; col_num:=(m1+n)!/(m1!*n!)-combinat[vectoint](SK); M[i,col_num]:=1; fi; end do; return M; end proc; ############################################################################################################################# ## The following code is to compute the matrix $P_{k}$ as in the paper, when the column $b$ is not the first ## column of the k-th Sylvester matrix, and the number of variables is 2 ## Input-- ## f,g --the polynomials ## num_b-- the number of the column in the matrix. ## k-- the degree of the GCD of polynomials ## Output-- ## res_M-- the matrix $P_{k}$ such that $P_{k}.\zz=b$ #################################################################################################################################### pk_matrix_2:=proc(f,g,ord,k,num_b) local M,M1,M2,res_M,m1,m2,row_dim,col_dim,N1,S,i,SS,SK,col_num,MM,deg_i; M:=p0_matrix(f,g,ord,k); if num_b=1 then return M; fi; m1:=degree(f,ord); m2:=degree(g,ord); if num_b<=(m2+2-k)!/((m2-k)!*2)-(m2+2-k-1)!/((m2-k-1)!*2) then M1:=Matrix(num_b-1,LinearAlgebra:-ColumnDimension(M)); M2:=LinearAlgebra:-SubMatrix(M,1..LinearAlgebra:-RowDimension(M)-num_b+1,1..-1); res_M:=Matrix([[M1],[M2]]); return res_M; fi; row_dim:=(m1+m2-k+2)!/((m1+m2-k)!*2); col_dim:=(m1+2)!/(m1!*2)+(m2+2)!/(m2!*2); MM:=Matrix(row_dim,col_dim); N1:=(m2-k+2)!/((m2-k)!*2); S:=combinat[inttovec](N1-num_b,2); deg_i:=add( i, i=S); for i from row_dim-(m1+deg_i+2)!/((m1+deg_i)!*2)+1 to row_dim do SS:=combinat[inttovec](row_dim-i,2); if min(op(SS-S))>=0 then SK:=SS-S; col_num:=(m1+2)!/(m1!*2)-combinat[vectoint](SK); MM[i,col_num]:=1; fi; end do; return MM; end proc; ########################################################################################################################################## ## The following code is to compute the matrix M such that M.zz is the first column of k-th sylvester matrix ## Input-- ## f,g -- the polynomials ## ord -- the list of the variables ## k -- the given integer ## Output-- ## M -- the matrix ################################################################################################################################### p0_matrix:=proc(f,g,ord,k) local i,j,m1,m2,f_N,g_N,row_dim,col_dim,M,n,N; m1:=degree(f,ord); m2:=degree(g,ord); n:=nops(ord); f_N:=(m1+n)!/(m1!*n!); g_N:=(m2+n)!/(m2!*n!): row_dim:=(m1+m2-k+n)!/((m1+m2-k)!*n!); col_dim:=f_N+g_N; M:=Matrix(row_dim,col_dim); N:=0; for i from 1 to m1+1 do for j from 1 to m1+2-i do M[N+j+(m2-k)*(i-1),N+j]:=1; end do; N:=N+m1+2-i; end do; return M; end proc; ############################################################################################################################# ## The following code is to compute the matrix $Y_{k}$ satisfies the conditions: $Y_{k}\zz$=$E.X$. ## Input-- ## L- denotes the list $X$ defined in the paper. ## LL- a list of the polynomials. ## k- degree of GCD of the polynomials ## num_b- b is the $num_b$-th column of the matrix Sylvester matrix. ## ord- variables. ## Output-- ## M-the matrix $Y_{k}$ in the paper. ######################################################################################################################################### mul_construct_XK:=proc(L,LL,k,num_b,ord) local n,r,res_L,num_org,j,ff,m,N,i,V,M_list1,M_list2,MM_sub1 ,v,M_sub2,M_sub1,MM,vec_add; r:=nops(LL); n:=nops(ord); #num_dim:=nops(con_vec_fg(LL,ord)); res_L:=[]; num_org:=0; for j from 1 to r do ff:=0; m:=degree(LL[j],ord)-k; N:=(m+n)!/(m!*n!); for i from num_org+1 to num_org+N do V:=combinat[inttovec](N-i+num_org,n); ff:=ff+aa[i]*mul(ord[kkk]^V[kkk],kkk=1..n); end do; res_L:=[op(res_L),ff]; num_org:=num_org+N; end do; M_list1:=[]:M_list2:=[]: for j from 1 to r-1 do ff:=res_L[j]; if n=2 then MM_sub1:=bvcauchyo(ff,ord,degree(LL[r],{op(ord)})); else MM_sub1:=mvcauchyo(ff,ord,degree(LL[r],{op(ord)})); fi; M_list1:=[op(M_list1),MM_sub1]; end do; M_list2:=[]; for v from 1 to r-1 do ff:=res_L[r]; if n=2 then MM_sub1:=bvcauchyo(ff,ord,degree(LL[v],{op(ord)})); else MM_sub1:=mvcauchyo(ff,ord,degree(LL[v],{op(ord)})); fi; M_list2:=[op(M_list2),MM_sub1]; end do; M_sub2:=LinearAlgebra:-DiagonalMatrix(M_list2); if r=2 then M_sub1:=M_list1[1]; else M_sub1:=M_list1[1]; for j from 2 to r-1 do M_sub1:=Matrix([[M_sub1],[M_list1[j]]]); end do; fi; MM:=Matrix([M_sub2,M_sub1]); if num_b=1 then vec_add:=Vector([0,L]); else vec_add:=Vector([LinearAlgebra:-SubVector(L,1..num_b-1),0,LinearAlgebra:-SubVector(L,num_b..-1)]); fi; MM:=subs({'aa[j]=vec_add[j]'$'j'=1..LinearAlgebra:-Dimension(L)+1},MM); return MM; end proc; #################################################################################################################################### ## The following code is to compute the matrix $E$ defined in the paper. ## Input- ## L-- a list represents the matrix $E$. ## LL-- a list including the degrees of f, g and GCD ## num_b-- the number of the column $b$. ## output-- ## M- the matrix $E$. ############################################################################################################################# mul_construct_E:=proc(L,LL,kk,num_b,ord) local m,n,k,i,j,N,num_org,num_dim,r,ff,V,M,aa,res_L; r:=nops(LL); n:=nops(ord); num_dim:=nops(con_vec_fg(LL,ord)); res_L:=[]; num_org:=0; for j from 1 to r do ff:=0; m:=degree(LL[j],ord); N:=(m+n)!/(m!*n!); for i from num_org+1 to num_org+N do V:=combinat[inttovec](N-i+num_org,n); ff:=ff+aa[i]*mul(ord[k]^V[k],k=1..n); end do; res_L:=[op(res_L),ff]; num_org:=num_org+N; end do; M:=LinearAlgebra:-DeleteColumn(multi_generate_syl(res_L,ord,kk),num_b); M:=subs({'aa[j]=L[j]'$'j'=1..num_dim},M); return M; end proc; ################################################################################################################################# ## The following code is to change the polynomials f,g into a list includes the coefficients of f and g ## Input-- ## L -- the polynomials ## ord -- the list of the variables ## Output-- ## S --the list of the coefficients of f and g ############################################################################################################################################# con_vec_fg:=proc(L,ord) local m,n,i,j,N1,N2,a,LL,S,r,N; r:=nops(L); S:=[]; for j from 1 to r do m:=degree(L[j],ord); n:=nops(ord); N:=(m+n)!/(m!*n!); for i from 1 to N do LL:=combinat[inttovec](N-i,n); a:=coeff_f(L[j],LL,ord); S:=[op(S),a]; end do; end do; return S; end proc; ################################################################################################################################### ## The following code is to change a list into two polynomial,here the coefficients of the polynomials ## is the element of the vector ## Input-- ## LL -- the polynomials ## ord -- the list of the variables ## L -- the list which is used to obtain the polynomials ## Output-- ## f0,g0-- the polynomials ##################################################################################################################################### con_pol_fg:=proc(LL,ord,V) local m,n,N1,f,i,r,NN,res_L,N,j,LLL,a,N2,f0; n:=nops(ord); r:=nops(LL); N1:=(m1+n)!/(m1!*n!); N2:=(m2+n)!/(m2!*n!); f:=0; NN:=0; res_L:=[]; for i from 1 to r do f:=LL[i]; m:=degree(f,ord); N:=(m+n)!/(m!*n!); f0:=0; for j from 1 to N do LLL:=combinat[inttovec](N-j,n); a:=V[j+NN]*mul(ord[k]^LLL[k],k=1..nops(LLL)); f0:=f0+a; end do; NN:=NN+N; res_L:=[op(res_L),f0]; end do; return res_L; end proc; # the following code is to compute the coefficients of monomials. coeff_f:=proc(f,L,ord) local m,n,a; if nops(ord)=1 then a:=coeff(f,ord[1],L[1]); else a:=coeff(coeff_f(f,[seq(L[i],i=2..nops(L))],[seq(ord[i],i=2..nops(ord))]),ord[1],L[1]); fi; return a; end proc; #The algorithm generate a random polynomial with the same support as the input polynomial. mul_randp:=proc(f,ord) local d,i,norm_f,n,V,N,r,m,res_f; n:=nops(ord); m:=degree(f,ord); N:=(m+n)!/(m!*n!); norm_f:=norm(f,2); V:=LinearAlgebra[RandomVector][row](N,generator=-0.02*norm_f..0.02*norm_f); res_f:=op(con_pol_fg([f],ord,V)); return res_f; end proc: ##################################################################################################################################3 ## The following procedure compute the k-th Sylvester matrix. ## Input-- ## L -- the list of the polynomials ## x -- the variable ## k -- the integer ## Output-- ## MM -- the k-th Sylvester matrix of f and g ################################################################################################################################### multi_generate_syl:=proc(L,ord,k) local r,n,M_list1,M_list2,i,j,MM1,M,m,t,MM,MM_sub1,MM_sub2,M_sub1,M_sub2; r:=nops(L); n:=nops(ord); m:=degree(L[r],ord); M_list1:=[]; M_list2:=[]; for i from 1 to r-1 do MM1:=LinearAlgebra:-SylvesterMatrix(L[r],L[i],x); if n=2 then MM_sub1:=bvcauchyo(L[r],ord,degree(L[i],{op(ord)})-k); MM_sub2:=bvcauchyo(L[i],ord,degree(L[r],{op(ord)})-k); else MM_sub1:=mvcauchyo(L[r],ord,degree(L[i],{op(ord)})-k); MM_sub2:=mvcauchyo(L[i],ord,degree(L[r],{op(ord)})-k); fi; M_list1:=[op(M_list1),MM_sub1]; M_list2:=[op(M_list2),MM_sub2]; end do; M_sub1:=LinearAlgebra:-DiagonalMatrix(M_list1); if r=2 then M_sub2:=M_list2[1]; else M_sub2:=M_list2[1]; for j from 2 to r-1 do M_sub2:=Matrix([[M_sub2],[M_list2[j]]]); end do; fi; MM:=Matrix([M_sub1,M_sub2]); return(MM); end proc; ############################################################################################################################# ## The following code is to compute the matrix $P_{k}$ as in the paper, when the column $b$ is not the first ## column of the k-th Sylvester matrix. In this code the number of variables >=2 ## Input-- ## L -- the polynomials ## num_b -- the number of the column in the matrix. ## k -- the degree of the GCD of polynomials ## Output-- ## res_M-- the matrix $P_{k}$ such that $P_{k}.\zz=b$ #################################################################################################################################### multi_p_matrix:=proc(L,ord,k,num_b) local r,n,M1,M2,num1_aux,num2_aux,num3_aux,m,M3,M,row_dim,i; r:=nops(L); n:=nops(ord); m:=degree(L[r]); row_dim:=add((degree(L[r])+degree(L[i])+n-k)!/((degree(L[r])+degree(L[i])-k)!*n!),i=1..r-1); M1:=pk_matrix(L[r],L[1],ord,k,num_b); num1_aux:=nops(con_vec_fg([L[r],L[1]],ord)); num2_aux:=nops(con_vec_fg(L,ord)); num3_aux:=nops(con_vec_fg([seq(L[i],i=1..r-1)],ord)); M2:=Matrix(num1_aux,num2_aux); for i from 1 to (m+n)!/(m!*n!) do M2[i,i+num3_aux]:=1; end do; M3:=Matrix([[LinearAlgebra:-IdentityMatrix(LinearAlgebra:-RowDimension(M1))],[Matrix(row_dim-LinearAlgebra:-RowDimension(M1),LinearAlgebra:-RowDimension(M1))]]); M:=M3.M1.M2; return M; end proc; ################################################################################################################################## ## According to the primary part LL, and the list of the polynomials L we compute the polynomial $p$. ## i.e. we compute the polynomial $p$ such that $f=pu$ and $f=pv$. ## Input-- ## L -- the list of the polynomials, where the polynomials in C[x1,x2....xn]. ## LL -- the list of the prime part of L ## Output-- ## gcd_res-- gcd of L ############################################################################################################################# multi_polydiv:=proc(L,LL,k,ord) local f_vec,n,M_sub,M,N,r,i,X,gcd_res,V; A_1,A_2,i,j,M,V,X,xp,gcd_res; r:=nops(L); n:=nops(ord); f_vec:=convert(con_vec_fg(L,ord),Vector); if n=2 then M:=bvcauchyo(LL[1],ord,k); else M:=mvcauchyo(LL[1],ord,k); fi; for i from 1 to r-1 do if n=2 then M_sub:=bvcauchyo(LL[i+1],ord,k); else M_sub:=mvcauchyo(LL[i+1],ord,k); fi; M:=Matrix([[M],[M_sub]]); end do; X:=LinearAlgebra:-LeastSquares(M,f_vec); #print(`error??=`,LinearAlgebra:-Norm(f_vec-M.X,infinity)); N:=(k+n)!/(k!*n!); gcd_res:=0; for i from 1 to N do V:=combinat[inttovec](N-i,n); gcd_res:=gcd_res+X[i]*mul(ord[jj]^V[jj],jj=1..n); end do; return gcd_res; end proc; ########################################################################################### ## The subroutine transforms linear constraints into a matrix about many polynomials(r>=2) mulpoly_multi_trans_mat:=proc(LL,L,var,ord) local i,j,num_dim,b,v,r,deg_L,res_L,res_var,va,var_dim,n,N,num_org,L_sub,vec_org,A,fg_org,L_sol,bb,XX,res_M,res_b,vv; r:=nops(LL); n:=nops(ord); deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(LL[i],ord)]; end do; member(min(op(deg_L)),deg_L, 'nd'); if nd=r then res_L:=LL; res_var:=var; elif nd=1 then res_L:=[seq(LL[i],i=2..r),LL[1]]; res_var:=[seq(var[i],i=2..r),var[1]]; else res_L:=[seq(LL[i],i=1..nd-1),seq(LL[i],i=nd+1..r),LL[nd]]; res_var:=[seq(var[i],i=1..nd-1),seq(var[i],i=nd+1..r),var[nd]]; fi; deg_L:=[]; for i from 1 to r do deg_L:=[op(deg_L),degree(res_L[i],ord)]; end do; for i from 1 to r do va[i]:=res_var[i]; end do; vec_org:=convert(con_vec_fg(res_L,ord),Vector); var_dim:=LinearAlgebra:-Dimension(vec_org); v:=[seq(vv[i],i=1..var_dim)]; num_org:=0; L_sub:=L; for i from 1 to r do N:=(deg_L[i]+n)!/((deg_L[i])!*n!); L_sub:=subs({'va[i][op(combinat[inttovec](N-j,n))]=vv[j+num_org]'$'j'=1..N},L_sub); num_org:=num_org+N; end do; A,b:=LinearAlgebra:-GenerateMatrix(L_sub,v); vec_org:=convert(con_vec_fg(res_L,ord),Vector); b:=A.vec_org-b; L_sol:=LinearAlgebra:-LinearSolve(A,b,free='s'); bb:=[seq(L_sol[i]=0,i=1..var_dim)]; res_M,res_b:=LinearAlgebra:-GenerateMatrix(bb,[op(indets(bb))]); XX:=LinearAlgebra:-LeastSquares(res_M,res_b); res_b:=res_b-res_M.XX; res_M:=Matrix([res_M,res_b]); end proc; ############################################################################################################################# ## The following code is to compute the matrix $H$ satisfies the conditions: $H\zz$=$S.X$. ## Input-- ## L-denotes the list $X$ defined in the paper. ## LL- a list of the polynomials. ## k- degree of GCD of the polynomials ## ord- variables. ## Output-- ## M-the matrix $H$ in the paper. ######################################################################################################################################### Int_point_matrix:=proc(L,LL,k,ord) local n,r,res_L,num_org,j,ff,m,N,i,V,M_list1,M_list2,MM_sub1 ,v,M_sub2,M_sub1,MM,vec_add; r:=nops(LL); n:=nops(ord); res_L:=[]; num_org:=0; for j from 1 to r do ff:=0; m:=degree(LL[j],ord)-k; N:=(m+n)!/(m!*n!); for i from num_org+1 to num_org+N do V:=combinat[inttovec](N-i+num_org,n); ff:=ff+aa[i]*mul(ord[kkk]^V[kkk],kkk=1..n); end do; res_L:=[op(res_L),ff]; num_org:=num_org+N; end do; M_list1:=[]:M_list2:=[]: for j from 1 to r-1 do ff:=res_L[j]; if n=2 then MM_sub1:=bvcauchyo(ff,ord,degree(LL[r],{op(ord)})); else MM_sub1:=mvcauchyo(ff,ord,degree(LL[r],{op(ord)})); fi; M_list1:=[op(M_list1),MM_sub1]; end do; M_list2:=[]; for v from 1 to r-1 do ff:=res_L[r]; if n=2 then MM_sub1:=bvcauchyo(ff,ord,degree(LL[v],{op(ord)})); else MM_sub1:=mvcauchyo(ff,ord,degree(LL[v],{op(ord)})); fi; M_list2:=[op(M_list2),MM_sub1]; end do; M_sub2:=LinearAlgebra:-DiagonalMatrix(M_list2); if r=2 then M_sub1:=M_list1[1]; else M_sub1:=M_list1[1]; for j from 2 to r-1 do M_sub1:=Matrix([[M_sub1],[M_list1[j]]]); end do; fi; MM:=Matrix([M_sub2,M_sub1]); vec_add:=L; MM:=subs({'aa[j]=vec_add[j]'$'j'=1..LinearAlgebra:-Dimension(L)},MM); return MM; end proc;