## ## Title: mroot - Multiple Roots and primal-dual structure ## Authors: Angelos Mantzaflaris ## Bernard Mourrain ## ## 2010-2011 ## See A.M. & B.M., `Deflation and Certified Isolation of Singular ## Zeros of Polynomial Systems` ## ---------------------------------------------------------------------- # # The module mroot # mroot := module() export #local # Deflation deflate, deflate_lin, deflate_incr, deflate_only, new_system,new_system1,new_system2, pert_poly, pert_system, #Combinatorics first_comb, next_comb, first_mon, next_mon, mac_dim, #Interval arithmetic operations mat_add,mat_sub,mat_prod, eval_on_ints,istart,iend, X_xs,i_box,mbox,dom_size, inclusion, # Operations on polynomials coeffof,lstcoefof, lstterms, lstmonof, lstmonsof,exp_of, dual_monomials, add_noise, diff_poly, diff_system, #Operations on matrices dmatrix2diffs, matrix2polys, matrix2diffs, rmaximal_minor, cmaximal_minor, nkernel, ncorank, # Operations on differentials has_term, add_df, sub_df, remove_zero,to_function, to_polynomial,to_monomial, refine_basis, symb_dual, int_df, diff_df, scalar_df, coefmatrix, # lex-ordering on differentials rlex_c,grlex_c,grlex_df,lex_c, lex_df, rlex_min,total_deg, leading_df, trailing_df, # Operations in the quotient normal_form, mult_table, # Dual Basis computation macaulay_block, cand_set, macaulay_step, mourrain_step, # Auxiliary check_list, var, rump, sign_of, app_funct, #Tests run_for, run_tests, cyclicn, #Experimental make_basis_pair, dual_constraints_1 ; export # rev-lex compare on differentials rlex_df, # Topological degee computation topological_degree, # Implicit curve branches computation curve_branches, # Certify root in domain certify_root, # Compute primal-dual orthogonal pair basis_pair, # Apply Rump's test rump_test, #Macaulay's method macaulay0, macaulay, # Integration method mourrain0, mourrain ; option package; use LinearAlgebra in ## ## Deflation by perturbation of the system ## deflate := proc( f::list, var::list, z::list, tol:=0, my_ind:=[] ) local n, AA, DD,N, J1,J2, ind, per ; n:=nops(var); per:=[cat(e, 1 .. nops(f))]; DD,AA:= basis_pair( expand(f),var,z,tol ) ; # compute dual basis #print( to_polynomial(DD), AA ) ; if nops(DD)=1 then lprint(`Warning: Simple root detected.`); fi; N:= new_system( expand(f),var, z, per, DD, AA ); #make perturbed system J2:=VectorCalculus:-Jacobian( N[1], var ) ; # compute Jacobian J2:= subs( seq(N[2][i]=z[i],i=1..n), seq(N[2][i]=0,i=n+1..nops(N[2]) ), J2); #print(N[1]); #print(Matrix(),N[2]); #print(`J=`, J2 ); #print(rmaximal_minor(J2), n ); if my_ind<>[] then #allow to choose which columbns to remove ind:= my_ind + [seq(n,i=1..n)] ; else ind:= rmaximal_minor(J2)+ [seq(n,i=1..n)]; # reduntant columns fi; #print(`Rows:`,ind-[seq(n,i=1..n)]); print(`Removing vars:`, N[2][ind] ); J1:= VectorCalculus:-Jacobian( N[1], N[2] );# compute augmented jacobian #print(J1); J1:= VectorCalculus:-Jacobian( subs(seq(N[2][ind[j]]=0,j=1..n), N[1]), subsop(seq( ind[j]=NULL,j=1..n),N[2]) );# compute jacobian print(`det`, subs( seq(N[2][i]=z[i],i=1..n), seq(N[2][i]=0,i=n+1..nops(N[2]) ), Determinant(J1) )); subs(seq(N[2][ind[j]]=0,j=1..n), N[1]), subsop(seq( ind[j]=NULL,j=1..n),N[2]); end: ## ## Deflation by incremental dual computation ## (under construction) deflate_incr := proc( f::list, var::list, z::list, tol:=0, my_ind:=[] ) local t, t1, i, n, BB, DD,N, J1,J2, ind, per ; n:=nops(var); per:=[cat(a, 1 .. nops(f))]; # Compute D_1 DD:=Matrix(n+1,1); DD[n+1,1]:=1; BB:=[1]; DD:= mourrain_step( expand(f),var,z,[DD], tol,BB ) ; # compute dual t:= coefmatrix(DD,var); t1:= cmaximal_minor( t[1], tol ); BB:= t[2][t1] ; # corresponding primal print(`Computed:`, to_polynomial(DD), BB ) ; if nops(DD)=1 then lprint(`Warning: Simple root detected.`); fi; for i to nops(DD) do print(to_polynomial( symb_dual(DD[i], cat(`a`,i) ) ) ); od; DD; end: # # Generates order 1 constraints as polynomial equations # dual_constraints_1:= proc( f::list, var::list, z::list ) local n, un, c, M, C, p; n:=nops(var); un:=[cat(`a`,1..n)]; # Condition (i) #c2:=combinat:-choose(n,2); c:=first_comb(n,2); M:=Matrix(nops(f),s); C:=NULL; while c <> NULL do p:= un[c[1]] - un[c[2]] ; C:= C, p; c:=next_comb(c,n); od; #for c in c2 do # p:= un[c[1]] - un[c[2]] ; # C:= C, p; #od; [C]; end: # # Deflate by adding linear perturbation on the equations # deflate_lin := proc( f::list, var::list, z::list, tol:=0, my_ind:=[] ) local n, AA, DD,N, J1,J2, ind, per ; n:=nops(var); per:=[cat(e, 1 .. nops(f))]; DD,AA:= basis_pair( expand(f),var,z,tol ) ; # compute dual basis print( to_polynomial(DD), AA ) ; if nops(DD)=1 then lprint(`Warning: Simple root detected.`); fi; N:= new_system1( expand(f),var, DD ); #make perturbed system J2:=VectorCalculus:-Jacobian( N[1], var ) ; # compute Jacobian J2:= subs( seq(N[2][i]=z[i],i=1..n), seq(N[2][i]=0,i=n+1..nops(N[2]) ), J2); if my_ind<>[] then #allow to choose which columbns to remove ind:= my_ind + [seq(n,i=1..n)] ; else ind:= rmaximal_minor(J2)+ [seq(n,i=1..n)]; # reduntant columns fi; #print(`Rows:`,ind-[seq(n,i=1..n)]); #print(`Removing vars:`, N[2][ind] ); J1:= VectorCalculus:-Jacobian( N[1], N[2] );# compute augmented jacobian print(`det`, subs( seq(N[2][i]=z[i],i=1..n), seq(N[2][i]=0,i=n+1..nops(N[2]) ), J1) ); subs(seq(N[2][ind[j]]=0,j=1..n), N[1]), subsop(seq( ind[j]=NULL,j=1..n),N[2]); end: # # Deflate without adding any new variables # deflate_only := proc( f::list, var::list, z::list, tol:=0, my_ind:=[] ) local n, AA, DD,N, J1,J2, ind, per ; n:=nops(var); per:=[cat(e, 1 .. nops(f))]; #DD,AA := basis_pair ( expand(f) ,var,z, tol ) ; DD,AA := mourrain ( expand(f) ,var,z, tol ) ; if nops(DD)=1 then lprint(`Warning: Simple root detected.`); fi; N:= new_system2(expand(f),var, DD ); #make augmented system J2:=VectorCalculus:-Jacobian( N, var ) ; # compute Jacobian J2:= subs(seq(var[i]=z[i],i=1..n) , J2); print( N, nops(N) ); print( J2, Dimension(J2) ); if my_ind<>[] then #allow to choose which columbns to remove ind:= my_ind; else ind:= rmaximal_minor(J2); # reduntant columns fi; #print("Poly,Diff=", seq([ind[k] mod nops(NN) +1, iquo(ind[k],nops(NN))], k=1..nops(ind)) ); #print(`Rows:`, ind ); [ seq( N[ind[j]],j=1..n) ]; end: # # Maximum length of a rectangular domain # dom_size:= proc(dom::list) local t,s,i; s:=dom[1][2]-dom[1][1]; for i from 2 to nops(dom) do t:=dom[i][2]-dom[i][1]; if t>s then s:=t; fi; od; s; end: # # Certify a multiple root of the approximate system f inside dom # certify_root := proc( f::list,var::list, dom::list ) local tol, DF, n, z; n:=nops(dom); tol:= dom_size(dom); z := [seq( (dom[k][1]+dom[k][2])/2, k=1..n) ]; print(`Certify tol:`,tol,`point:`,z ) ; DF:= deflate ( f , var , z , tol , []); #DF:= deflate_lin( f , var , z , tol , []); print(`New system size:`,nops(DF[2]) ) ; print( rump_test( DF[1], DF[2], i_box(dom, nops(DF[2]), tol) ) ); DF; end: # # Apply Rump's test for a simple root close to xs in a given range XX # rump_test :=proc ( f::list,var::list, xs::list, XX::list ) local i,E; E:= rump(f,var,xs,XX);# #print(`Inclusion:`, E ); #print(xs,XX); for i to nops(var) do if not inclusion(E[i,1],XX[i]) then print(`At position`,i,`:`, evalr(E[i,1]+xs[i]),evalr(XX[i]+xs[i]) ); RETURN(false); fi; E[i,1]:= evalr(E[i,1]+xs[i]); od: print(`Inclusion made.`); E; end: # # Compute Rump's inclusion formula # rump := proc ( f::list,var::list, xs::list, X::list) local Id, J, R, M, fev; J:= VectorCalculus:-Jacobian(f, var); Id:= Matrix(nops(var),shape=identity) ; #Jacobian evaluated on X+xs M:= eval_on_ints( var, X_xs(X,xs), J ) ; #print(`M=`,M); #print(f,var, J); #print( subs( seq( var[j]=xs[j],j=1..nops(var)),J ) ) ; print(`det`, #subs( seq( var[j]=xs[j],j=1..nops(var)),J ), ` = `, Determinant(subs( seq( var[j]=xs[j],j=1..nops(var)),J )) ) ; # Inverse of Jacobian evaluated on initial approximation #print(`J=`, J); #print( subs( seq( var[j]=xs[j],j=1..nops(var)),J ) ) ; #print( Determinant(subs( seq( var[j]=xs[j],j=1..nops(var)),J )) ) ; R:= MatrixInverse( subs( seq( var[j]=xs[j],j=1..nops(var)),J ) ) ; # Evaluation of f over xs fev:= Transpose(Matrix( subs(seq(var[j]=xs[j],j=1..nops(var)), f ) )) ; #print(`f_ev =`, fev); #print(`R.M =`, mat_prod(R,M)); # Rump's Formula: -R.fev + (Id-R.M).X mat_add( -R.fev , mat_prod( mat_sub(Id,mat_prod(R,M)), Transpose(Matrix(X))) ) ; end: ##################################### # Test Function ##################################### # f polynomial system # var variables # z exact root # ze root perturbation # bsize radius of the box with center z+ze # tol tolerance when computing numerical kernel run_for := proc(f::list, var::list, z, ze, bsize, tol:=0 ) local DS1, DS2; Digits:=16: #interface(displayprecision = 3): #interface(rtablesize=10): # Exact Deflated System #DS0:=deflate( f , vars , z ): # Rump's Test for DS0 #rump_test( DS0[1], DS0[2], i_box( mbox(z+ze,bsize), nops(DS0[2])) ) ; # Approximate deflated system DS1:=deflate( f , var , z + ze , tol): # Rump's Test for DS1 rump_test( DS1[1], DS1[2], i_box( mbox(z+ze,bsize), nops(DS1[2])) ) ; end: #mult := proc (f,var,z, tol:=0, comp:=rlex) #nops( get_dbasis(f,var,z,tol,rlex) ); #end: # # Sets coeff of primal monomial to 1 # (not used) # refine_basis:= proc(A::list, dual_basis::list, var::list ) local RB,i,j, n,m,c,k; n:= nops(var); RB:=dual_basis; for i to nops(A) do j:=1; while j<= Dimension(RB[i])[2] do if Column(RB[i],j)[1..n]=A[i][1..n] then RB[i]:= scalar_df( RB[i], 1/RB[i][n+1,j] ) ; #break; fi; for k to i-1 do if Column(RB[i],j)[1..n]=A[k][1..n] then RB[i]:= sub_df( RB[i],scalar_df(RB[k],RB[i][n+1,j])); fi; od; j:=j+1; od; od; RB; end: # # Given a list of monomials, computes their dual elements. # (not used) # dual_monomials := proc(var::list, A::list) local k,n,i,j,t,Dset; n:=nops(var); Dset:=NULL; for i from 1 to nops(A) do t:=Matrix(n+1,1); t[n+1]:=1; for j from 1 to nops(var) do t[j,1] := degree (A[i], var[j]); od; #print( A[i], t ); Dset:= Dset, t; od; [Dset]; end: # # Replaces coefficient of a dual element with variables ai # symb_dual := proc ( df::Matrix, aa:= `a` ) local n,m,un,j,sd ; n,m:= Dimension(df); un:=[cat(aa,1..m)]; sd:= copy(df); for j to m do sd[n,j]:= un[j]; od; sd; end: # # Based on new ideas # new_system:= proc ( f::list, var::list, z::list, per::list, dual_basis::list, pr_basis) local pp,all, n,g,i,j,ff,t, fu; all:=op(var); n:=nops(var); pp,all:= pert_system( f,z, var, per , pr_basis ) ; ff:=NULL; for j to n do for i to nops(dual_basis) do ff:= ff, diff_poly(pp[j],var,dual_basis[i]) ; od; od; [ff],all; end: # # apply linear perturbation e[j,i] on equations # new_system1:= proc ( f::list, var::list, dual_basis::list ) local pp,all, n,g,i,j,ff,t, fu; all:=op(var); n:=nops(var); t:=nops(f); ff:=NULL; for j to t do #ff:= ff, f[j]; for i from 2 to nops(dual_basis) do ff:= ff, diff_poly(f[j],var,dual_basis[i]) + e[j,i] ; all:=all , e[j,i]; od; od; ff:=[ff]; print(`New system is:`, Matrix() ); ff,[all]; end: # # Just add equations, not variables # new_system2:= proc ( f::list,var::list,dual_basis::list ) local i,j,ff ; ff:=NULL; for j to nops(f) do for i to nops(dual_basis) do print(to_polynomial(dual_basis[i]), f[j]); ff:= ff, diff_poly(f[j],var,dual_basis[i]) ; od; od; [ff]; end: # #Compute perturbed polynomial system wrt the elements of A # pert_system:= proc( f::list, z::list, var::list,per::list, A::list) local pp, i, j, all; all:=op(var); pp:=NULL; for i to nops(f) do pp:= pp, expand( pert_poly(f[i],z,var,A,per[i])) ; for j to nops(A) do all:= all, per[i][j]; od; od: [pp],[all]; end: # # Perturb polynomial wrt the elements of A # pert_poly := proc (p, z::list, var::list, A::list, k ) local n,i, g; #print("pert_poly input", p, z ); n:=nops(var); g:=p; # add new variables for i to nops(A) do g:= g+ subs(seq(var[j]=var[j]-z[j],j=1..nops(z)), A[i])* k[i]; od; #print("pert_poly output", g ); expand(g); end: # # Add noise to coefficients of a polynomial system # add_noise := proc ( f::list, ns:=1e-3 ) local Pt,m,i, g, j; g:=f; for i to nops(f) do #random numbers Pt:=Array(RandomTools:-Generate( list(float(range=-0..ns,'method=uniform'),nops(f[i]))), 'datatype=float[8]'); j:=1; for m in f[i] do # perturbe f_i g[i]:= g[i] + Pt[j]* m ; j:=j+1; od; od; g; end: # # Differentiate polynomial (normalized) # diff_poly := proc (p, var::list, df::Matrix ) local s,j,i,n, g, res; n:=nops(var); s:= Dimension(df)[2]; res:=0; for j to s do g:=p; for i to n do # normalized differentiation g:= diff( g, [ x[n-i+1]$df[n-i+1,j] ]) / df[n-i+1,j]! ; od; res:=res + df[n+1,j]*g; #coeff times g #res:= res + g; od; res; end: # # Apply (normalized) functionals on system f # diff_system := proc ( f::list,var::list,dual_basis::list ) local i,j,ff ; ff:=NULL; for j to nops(f) do for i to nops(dual_basis) do ff:= ff, diff_poly(f[j],var,dual_basis[i]) ; od; od; [ff]; end: # # Return the indices of n independent rows of an mxn Matrix # rmaximal_minor:= proc( M::Matrix, tol:=0.05) cmaximal_minor(Transpose(M), tol ); end: # # Returns the indices of n independent columns of an nxm Matrix # cmaximal_minor:= proc( M::Matrix, tol:=0.05) local RR, i,j,n, cc; n:= Dimension(M)[1]; RR:= QRDecomposition( M , output = ['R'] ); #Find principal elements in RR-matrix cc:=NULL; i:=1; j:=1; while i<=n do if evalf(abs(RR[i,j]))>tol then # evalf(.) added for Maple2017 cc:= cc, j; i:=i+1; fi; j:=j+1; od; return [cc]; end: ############################## # Rump's Test help functions ############################## # Matrix Multiplication with evalr mat_prod := proc ( A::Matrix, B::Matrix ) local P,i,j,k,n1,n2,m1,m2; n1,m1:= Dimension(A); n2,m2:= Dimension(B); if m1<>n2 then print("Wrong Dim");RETURN(-1); fi; P:= Matrix(n1,m2); for i to n1 do for j to m2 do for k to n2 do P[i,j]:= evalr(P[i,j] + evalr(A[i,k]*B[k,j])); od; od; od; P; end: # Matrix Addition with evalr mat_add := proc ( A::Matrix, B::Matrix ) local P,i,j,k,n1,n2,m1,m2; n1,m1:= Dimension(A); n2,m2:= Dimension(B); if n1<>n2 or m1<>m2 then print("Wrong Dim");RETURN(-1); fi; P:= Matrix(n1,m1); for i to n1 do for j to m1 do P[i,j]:= evalr(A[i,j] + B[i,j]); od; od; P; end: # Matrix Substruction with evalr mat_sub := proc ( A::Matrix, B::Matrix ) local P,i,j,k,n1,n2,m1,m2; n1,m1:= Dimension(A); n2,m2:= Dimension(B); if n1<>n2 or m1<>m2 then print("Wrong Dim");RETURN(-1); fi; P:= Matrix(n1,m1); for i to n1 do for j to m1 do P[i,j]:= evalr(A[i,j] - B[i,j]); od; od; P; end: # Evaluate Polynomial Matrix on intervals eval_on_ints:= proc( var::list ,int::list , M::Matrix ) local E, n,m,i,j; #print(`Eval_on_ints IN:`, var, int, M ); n,m:= Dimension(M); E:= Matrix(n,m); for i to n do for j to m do E[i,j]:= evalr( subs( seq( var[j]=int[j], j=1..nops(var)) , M[i,j] ) ) ; od: od: E; end: # Compute X+xs with X interval and xs value X_xs:= proc( X::list, xs::list) local i, E; E:=X; for i to nops(X) do E[i]:= evalr( X[i]+ xs[i]); od; E; end: #Interval box for use in rump's test i_box:= proc( bx::list, nvars::integer, tol:=0.01 ) local i,zs,n, X,cc ; n:= nops(bx); zs:= [ seq( (bx[i][1]+bx[i][2])/2, i=1..nops(bx) ), seq( 0, i=1..nvars-nops(bx) ) ]; X:=NULL; for i to nops(bx) do cc:=bx[i]- [zs[i], zs[i]]; X:= X, INTERVAL(cc[1]..cc[2]); od: X:= X, seq( INTERVAL(-tol..tol),j= nops(bx)+1..nops(zs) ) ; #print(`i_box: `, zs ); zs ,[X] ; end: # Make a box around z with certain tolerance mbox:= proc (z::list, tol:=0.02) [seq( [z[j]-tol,z[j]+tol], j=1..nops(z))]; end: # Return starting value of interval istart:= proc ( a) if whattype(a)=function then RETURN(op( op(a) )[1]); else RETURN(a); fi: end: # Return ending value of interval iend:= proc ( a ) if whattype(a)=function then RETURN(op( op(a) )[2]); else RETURN(a); fi: end: #Checks if a is included in b inclusion:= proc ( a, b) #a:function or number #print(`inclusion IN:`,a,b); istart(a)>istart(b) and iend(a)0 then if df[i] > 0 then return true; else return false; fi; fi; od; fi; false; end: lex_df := proc ( ddf1::Matrix, ddf2::Matrix ) lex_c(leading_df(ddf1, lex_c), leading_df(ddf2, lex_c) ); end: # # Graded lexicographic order from the end # rlex_c := proc ( df1::Vector, df2::Vector ) # df1, df2 are nx1 Matrix # df1 0 then if df[i] > 0 then return true; else return false; fi; fi; od; fi; false; end: rlex_df := proc ( ddf1::Matrix, ddf2::Matrix ) rlex_c(leading_df(ddf1, rlex_c), leading_df(ddf2, rlex_c) ); end: rlex_min := proc ( p , var::list) Groebner:-TrailingTerm(p, grlex(op(var) ))[2] ; end: # # Returns the total degree (order) of a differential # total_deg := proc(df1::Vector) local s,i; s:=0: for i from 1 to Dimension(df1)-1 do s:=s+ df1[i]; od; s; end: # # Graded order, not used anymore here. # grlex_c := proc ( df1::Vector, df2::Vector ) # df1, df2 are nx1 Matrix local s1,s2,n,i, df; n:= Dimension(df1); s1:=0: s2:=0: for i from 1 to n-1 do s1:=s1+ df1[i]; s2:=s2+ df2[i]; od; if s1< s2 then return true; fi; # go to rlex if s1=s2 then return rlex_c(df1,df2); fi; false; end: grlex_df := proc ( ddf1::Matrix, ddf2::Matrix ) grlex_c(leading_df(ddf1, grlex_c), leading_df(ddf2, grlex_c) ); end: # # leading term of differential wrt an ordering # leading_df:= proc( df::Matrix, compare:=rlex_c ) local t,lt,i,n,m; n,m:= Dimension(df); lt:= Column(df,1); for i from 2 to m do t:= Column(df,i); if not compare(t,lt) then lt:=t; fi; od; lt; end: # # Trailing term of differential wrt an ordering # trailing_df:= proc( df::Matrix, compare:=rlex_c ) local t,lt,i,n,m; n,m:= Dimension(df); lt:= Column(df,1); for i from 2 to m do t:= Column(df,i); if compare(t,lt) then lt:=t; fi; od; lt; end: # #Return the coeff. of p in variables var of the monomial m #(from multires) coeffof := proc(m,p,vvar) local i,zs,lm,lc,k; lc := [coeffs(p,vvar,'lm')]; lm := [lm]; if member(m,lm,'k') then lc[k] else 0 fi; end: # # Return the exponents of a monomial list # exp_of := proc(mm::list,var::list) local i,j, E , L; L:=NULL; for j to nops(mm) do E:=NULL; for i to nops(var) do E:= E, degree(mm[j],var[i]); od; L:= L,[E]; od; [L]; end: # # List coefficients of polynomial # lstcoefof := proc(p,var) local lm,r,c; lm := NULL; r := sort(p,var,tdeg); while r <> 0 do c := tcoeff(r,var,'m'); lm := lm, c; r := expand(r-c*m); od; lm; end: # #List the monomials of a polynomial p in the variables var: #(from multires) lstmonof := proc(p,var:=indets(PP)) local lm,r,c; lm := NULL; r := sort(p,var,tdeg); while r <> 0 do c := tcoeff(r,var,'m'); lm := lm, m; r := expand(r-c*m); od; #ListTools:-MakeUnique([lm]); #print(lm); op(sort([lm], proc (a, b) options operator, arrow;Groebner:-TestOrder(a, b, grlex(op(var))) end proc)); #lm; end: # #List the monomials of a polynomial system in the variables var: # lstmonsof := proc(PP::list,var:=indets(PP) ) local lm,p; lm := NULL; for p in PP do lm:= lm, lstmonof(p,var); od; ListTools:-MakeUnique([lm]); sort(%, proc (a, b) options operator, arrow;Groebner:-TestOrder(a, b, grlex(op(var))) end proc); end: # # List monomials of a differential # lstterms := proc(df::Matrix ) local t,n1,i,j,A; n1:= Dimension(df)[1]; A:=NULL; #for df in DD do for j to Dimension(df)[2] do t:=Column(df,j); t[n1]:=1; A:= A, t; od; #od; #sort(ListTools:-MakeUnique([A]), cat(comp, _c) ); #sort([A], cat(comp, _c) ); sort([A], rlex_c ); end: # # Check for dublicates in a sorted list # (not used) check_list := proc ( l::list ) local i; for i to nops(l)-1 do if l[i]=l[i+1] then lprint(`Warning: Dublicates in the list`); return false; fi; od; true; end: # # Return a differential as an operator (to be applied on a polynomial) # to_function := overload([ proc( db::Matrix,var::list) option overload; local n,m; #n is number of variables-1, m number of terms n,m:= Dimension(db); unassign(`h`); unapply( add( db[n,k]* mul(1/factorial(db[i,k]), i=1..n-1) * diff(h(op(var)), [seq(var[i]$(db[i,k]), i=1..n-1)]) #diff , k=1..m) #add , h); #unapply end, proc( db::list,var::list) option overload; local t, FF; FF:= []; for t in db do FF:= [op(FF), to_function(t,var) ]; od; FF; end ]); # # Apply Df on function t # app_funct := proc( Df, t ) Df( unapply(t) ); end: # # Return a differential as a polynomial # to_polynomial := overload([ proc( db::Matrix, var :=[seq(d[i],i=1..(Dimension(db)[1]-1))]) option overload; local t,j,n,m; #print("to_polynomial_input:",db,var , Dimension(db[1])); n:=nops(var); #m= number of terms m:= Dimension(db)[2]; t:=0; for j to m do t:= t+ db[n+1,j]*mul( var[s]^db[s,j],s=1..n); od; #print("to_polynomial_Output:",t); t; end, proc( db::list, var::list :=[seq(d[i],i=1..Dimension(db[1])[1]-1)] ) option overload; local t, FF; FF:= NULL; for t in db do FF:= FF, to_polynomial(t,var); od; [FF]; end ]); # # Return a differential with one terms as a monomial # to_monomial := overload([ proc( db::Vector, var :=[seq(d[i],i=1..(Dimension(db) -1))]) option overload; local n; #print("to_monomial_input:",db,var , Dimension(db)); n:=nops(var); #db[n+1]* mul( var[s]^db[s],s=1..n); end, proc( db::list, var::list :=[seq(d[i],i=1..Dimension(db[1])[1]-1)] ) option overload; local t, FF; #print( var); FF:= NULL; for t in db do FF:= FF, to_monomial(t,var); od; [FF]; end ]); # # Input: coefficient matrix A, indexing diffs NN # Output: Differentials induced by A's rows # dmatrix2diffs := proc(R::Matrix, NN::list, tol:=0) local n,p,M, j, i, t ,s ; #print("to_diffs_Input:", R, NN ); n:=nops(NN[1]); s:=nops(NN); p:=NULL; for j to Dimension(R)[2] do #M:=NULL; M:=Matrix(n+1,1); for i to s do if abs(R[i,j])>tol then #print(i, NN[i], R[i,j]); t:= Matrix(); M:=add_df(M,t,tol); fi; od; p:= p, Matrix([M]); od; #print("to_diffs_Output:", p ); [p]; end: # # Create variable list [ x[1],...,x[n] ] # var := proc (n) [seq(x[i],i=1..n)]; end: # # Run tests ! # run_tests := proc(test::integer ) local Db,A, dd,i,j,z0,z1,z2,G,v,p ; Digits:= 32: interface(displayprecision = 5): interface(rtablesize=30): ############################# # Systems ############################# G:=[seq(0,k=1..12) ]; p:=[seq(0,k=1..12) ]; v:=[seq(0,k=1..12) ]; #1. cbms1: G[1]:=[ x[1]^3 - x[2]*x[3], x[2]^3 - x[1]*x[3], x[3]^3 - x[1]*x[2] ]; p[1]:= [ [0, 0, 0] ]; v[1]:= var(nops(%[1])); #2. cbms2: G[2]:=[ x[1]^3 - 3*x[1]^2*x[2] + 3*x[1]*x[2]^2 - x[2]^3 - x[3]^2, x[3]^3 - 3*x[3]^2*x[1] + 3*x[3]*x[1]^2 - x[1]^3 - x[2]^2, x[2]^3 - 3*x[2]^2*x[3] + 3*x[2]*x[3]^2 - x[3]^3 - x[1]^2 ]; p[2]:= [ [0, 0, 0] ]; v[2]:= var(nops(%[1])); #3. mth191: G[3]:=[ x[1]^3 + x[2]^2 + x[3]^2 - 1, x[1]^2 + x[2]^3 + x[3]^2 - 1, x[1]^2 + x[2]^2 + x[3]^3 - 1 ]; p[3]:=[ [0, 1, 0] ]; v[3]:= var(nops(%[1])); #4. decker2: G[4]:=[ x[1] + x[2]^3, x[1]^2*x[2] - x[2]^4 ]; p[4]:= [ [0, 0] ]; v[4]:= var(nops(%[1])); #5. Ojika2: G[5]:=[ x[1]^2 + x[2] + x[3] - 1, x[1] + x[2]^2 + x[3] - 1, x[1] + x[2] + x[3]^2 - 1 ]; p[5]:=[ #[0, 0, 1], [1, 0, 0] ]; v[5]:= var(nops(%[1])); #6. Ojika3: G[6]:=[ x[1] + x[2] + x[3] - 1, 2*x[1]^3 + 5*x[2]^2 - 10*x[3] + 5*x[3]^3 + 5, 2*x[1] + 2*x[2] + x[3]^2 - 1 ]; p[6]:=[ [0, 0, 1] #,[-5/2, 5/2, 1] ]; v[6]:= var(3); #7. KSS: G[7]:=[ x[1]^2-2*x[1] + x[1]+x[2]+x[3]+x[4]+x[5]-4, x[2]^2-2*x[2] + x[1]+x[2]+x[3]+x[4]+x[5]-4, x[3]^2-2*x[3] + x[1]+x[2]+x[3]+x[4]+x[5]-4, x[4]^2-2*x[4] + x[1]+x[2]+x[3]+x[4]+x[5]-4, x[5]^2-2*x[5] + x[1]+x[2]+x[3]+x[4]+x[5]-4 ]; p[7]:= [ [1,1,1,1,1] ]; v[7]:= var(nops(%[1])); #8. Caprasse: G[8]:=[ -x[1]^3*x[3] +4*x[1]*x[2]^2*x[3] +4*x[1]^2*x[2]*x[4] +2*x[2]^3*x[4] +4*x[1]^2 -10*x[2]^2 +4*x[1]*x[3] -10*x[2]*x[4] +2, -x[1]*x[3]^3+4*x[2]*x[3]^2*x[4]+4*x[1]*x[3]*x[4]^2+2*x[2]*x[4]^3+4*x[1]*x[3]+4*x[3]^2- 10*x[2]*x[4]- 10*x[4]^2+ 2, x[2]^2*x[3] +2*x[1]*x[2]*x[4] -2*x[1] -x[3], x[1]*x[4]^2 +2*x[2]*x[3]*x[4] -x[1] -2*x[3] ]; p[8]:= [ [ 2, -I*sqrt(3), 2, I*sqrt(3) ] ]; v[8]:= var(nops(%[1])); #9. Cyclic nine: G[9]:= [ x[1]+x[2]+x[3]+x[4]+x[5]+x[6]+x[7]+x[8]+x[9], x[1]*x[2]+x[2]*x[3]+x[3]*x[4]+x[4]*x[5]+x[5]*x[6]+x[6]*x[7]+x[7]*x[8]+x[8]*x[9]+x[9]*x[1], x[1]*x[2]*x[3]+x[2]*x[3]*x[4]+x[3]*x[4]*x[5]+x[4]*x[5]*x[6]+x[5]*x[6]*x[7]+x[6]*x[7]*x[8]+x[7]*x[8]*x[9]+x[8]*x[9]*x[1]+x[9]*x[1]*x[2], x[1]*x[2]*x[3]*x[4]+x[2]*x[3]*x[4]*x[5]+x[3]*x[4]*x[5]*x[6]+x[4]*x[5]*x[6]*x[7]+x[5]*x[6]*x[7]*x[8]+x[6]*x[7]*x[8]*x[9]+x[7]*x[8]*x[9]*x[1]+x[8]*x[9]*x[1]*x[2]+x[9]*x[1]*x[2]*x[3], x[1]*x[2]*x[3]*x[4]*x[5]+x[2]*x[3]*x[4]*x[5]*x[6]+x[3]*x[4]*x[5]*x[6]*x[7]+x[4]*x[5]*x[6]*x[7]*x[8]+x[5]*x[6]*x[7]*x[8]*x[9]+x[6]*x[7]*x[8]*x[9]*x[1]+x[7]*x[8]*x[9]*x[1]*x[2]+x[8]*x[9]*x[1]*x[2]*x[3]+x[9]*x[1]*x[2]*x[3]*x[4], x[1]*x[2]*x[3]*x[4]*x[5]*x[6]+x[2]*x[3]*x[4]*x[5]*x[6]*x[7]+x[3]*x[4]*x[5]*x[6]*x[7]*x[8]+x[4]*x[5]*x[6]*x[7]*x[8]*x[9]+x[5]*x[6]*x[7]*x[8]*x[9]*x[1]+x[6]*x[7]*x[8]*x[9]*x[1]*x[2]+x[7]*x[8]*x[9]*x[1]*x[2]*x[3]+x[8]*x[9]*x[1]*x[2]*x[3]*x[4]+x[9]*x[1]*x[2]*x[3]*x[4]*x[5], x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]+x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]+x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]+x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[1]+x[5]*x[6]*x[7]*x[8]*x[9]*x[1]*x[2]+x[6]*x[7]*x[8]*x[9]*x[1]*x[2]*x[3]+x[7]*x[8]*x[9]*x[1]*x[2]*x[3]*x[4]+x[8]*x[9]*x[1]*x[2]*x[3]*x[4]*x[5]+x[9]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6], x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]+x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]+x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[1]+x[4]*x[5]*x[6]*x[7]*x[8]*x[9]*x[1]*x[2]+x[5]*x[6]*x[7]*x[8]*x[9]*x[1]*x[2]*x[3]+x[6]*x[7]*x[8]*x[9]*x[1]*x[2]*x[3]*x[4]+x[7]*x[8]*x[9]*x[1]*x[2]*x[3]*x[4]*x[5]+x[8]*x[9]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6]+x[9]*x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7], 1-x[1]*x[2]*x[3]*x[4]*x[5]*x[6]*x[7]*x[8]*x[9] ]; z0:=-0.9396926-0.3520201*I: z1:=-2.4601472-0.8954204*I: z2:=-0.3589306-0.1306401*I: p[9]:=[ [z0,z1,z2,z0,-z2,-z1,z0,-z2,-z1] ]; v[9]:= var(nops(%[1])); #10. DZ1: G[10]:= [ x[1]^4 -x[2]*x[3]*x[4], x[2]^4 -x[1]*x[3]*x[4], x[3]^4 -x[1]*x[2]*x[4], x[4]^4 -x[1]*x[2]*x[3] ]; p[10]:= [ [0, 0, 0, 0] ]; v[10]:= var(nops(%[1])); #11. DZ2: G[11]:= [ x[1]^4, x[1]^2*x[2] + x[2]^4, x[3] + x[3]^2 - 7*x[1]^3 - 8*x[1]^2 ]; p[11]:= [ [0, 0, -1] ]; v[11]:= var(nops(%[1])); #12. DZ3: G[12]:= [ 14*x[1] + 33*x[2] - 3*sqrt(5)* (x[1]^2 + 4*x[1]*x[2] + 4*x[2]^2 + 2) + sqrt(7) + x[1]^3 + 6*x[1]^2*x[2] + 12*x[1]*x[2]^2 + 8*x[2]^3, 41*x[1] - 18*x[2] - sqrt(5) + 8*x[1]^3 - 12*x[1]^2*x[2] + 6*x[1]*x[2]^2 - x[2]^3 + 3*sqrt(7)*(4*x[1]*x[2] - 4*x[1]^2 - x[2]^2 - 2) ]; # with coeffcients rounded to 5 digits, at approximate zero dd:=Digits; Digits:=5; G[12]:= evalf(G[12]); Digits:=dd; p[12]:= [ [1.5055, 0.36528] ]; v[12]:= var(nops(%[1])); ############################# # Compute Dual Basis test ############################# if test=1 then ############################# for i in [seq(4..4)] do for j to nops(p[i]) do print("--",i,"--", p[i][j], "--"); #Db, A := mourrain(G[i],v[i],p[i][j] ); #Db := mourrain0(G[i],v[i],p[i][j] ); Db, A := macaulay(G[i],v[i],p[i][j] ); print( nops( Db ), Db ); print("Mac_Dim:", mac_dim(nops(G[i]),nops(v[i]), degree(to_polynomial(Db[-1], v[i]))+1 ) ); od; od; fi; ############################# # Certify root test ############################# if test=2 then ############################# for i in [9] do for j to nops(p[i]) do print("--",i,"--", p[i][j], "--"); Db, A := mourrain(G[i],v[i],p[i][j]); print( nops( Db ) ); print("Mac_Dim:", mac_dim(nops(G[i]),nops(v[i]), degree(to_polynomial(Db[-1], v[i]))+1 ) ); od; od; fi; ## ## More instances ## # Ojika87 #f1 := x[1]^2+x[2] - 3 ; #f2 := x[1] + 1/8*x[2]^2 -3/2 ; #z:=[1,2]: #ze := [0.002,0.001]: #print(`Runing for`,[f1,f2], z+ze ); #DS:=deflate( [f1,f2] , vars , z + ze , 0.03 ); #rump_test( DS[1], DS[2], i_box( mbox(z+ze,0.02), nops(DS[2])) ) ; #RumpGr-success #f1 := x[1]^2 - x[2]^2 ; #f2 := x[1] - x[2]^2 ; #z:=[0,0]: #ze := [0.002,0.001]: #print(`Runing for`,[f1,f2], z+ze ); #DS:=deflate( [f1,f2] , vars , z + ze , 0.03 ); #rump_test( DS[1], DS[2], i_box( mbox(z+ze,0.02), nops(DS[2])) ) ; #RumpGr-fail #f1 := x[1]^2*x[2] - x[1]*x[2]^2 ; #f2 := x[1] - x[2]^2 ; #z:=[0,0]: #ze := [0.002,0.001]: #print(`Runing for`,[f1,f2], z+ze ); #DS:=deflate( [f1,f2] , vars , z + ze , 0.03 ); #rump_test( DS[1], DS[2], i_box( mbox(z+ze,0.02), nops(DS[2])) ) ; end: # # Makes Macaulay's matrix block of depth dp # (not used) # macaulay_block := proc(f::list, var::list, z::list, dp::integer:=1 ) local n, r, c, M, ci, mr, ri, i, t; n:= nops(var); if dp<0 then return [0]; fi; if dp=0 then return Matrix([0]); fi; r:= nops(f)*binomial(n+dp-2,dp-1); c:= binomial(n+dp-1,dp); M:= Matrix(r,c); r:=0: c:=0: ci:= first_mon( dp , n); while ci<>NULL do c:=c+1; mr:=0; ri:= first_mon(dp-1, n); while ri<>NULL do r:=r+1; for i to nops(f) do# f_i t:=diff_poly( mul(var[s]^ri[s],s=1..n)*f[i], var, Matrix() ); M[mr + i,c]:= eval(t,[seq(var[s]=z[s],s=1..n)]); od; mr := mr+ i-1; ri:= next_mon( ri , dp-1); od; ci:= next_mon( ci , dp); od; M; end: # # `Integrate` a differential wrt i-th variable # (normalized) # int_df := proc( df::Matrix, i::integer) local n,t, m, ci, j, k; #print("int_df_input:",df,"var",i); t:= copy(df); n:= Dimension(t)[1]-1; if i>n then lprint("Integration out of range"); return df; fi; m:=Dimension(t)[2]; for j to m do # Integrate t[i,j]:= t[i,j]+1; od; #print("int_df_output:",t); t; end: # # Remove terms with near-to-zero coefficient # from a differential # remove_zero :=proc (df::Matrix, tol:=0 ) local n1,m,i,c; n1,m:= Dimension(df); c:= NULL; for i to m do if evalf(abs(df[n1,i]))<=tol then c:=c,i; fi; od; if nops([c])=0 then return df; fi; if nops([c])=m then return Matrix(n1,1); fi; DeleteColumn(df,[c]); end: # # Add two differentials, df1 + df2 # add_df := proc (df1::Matrix, df2::Matrix, tol:=0) local ex, flag, t,n,i,j; #print(" add_df_input",df1,df2); n:= Dimension(df1)[1]-1; if Equal(df1,Matrix(n+1,1)) then return copy(df2); fi; if Equal(df2,Matrix(n+1,1)) then return copy(df1); fi; t:=copy(df1); ex:= 1..n; for i to Dimension(df2)[2] do flag:=false; j:=1; while j<= Dimension(t)[2] do #print( t[ex,j],df2[ex,i] ); if Equal(t[ex,j],df2[ex,i]) then t[n+1,j]:= t[n+1,j] + df2[n+1,i]; if abs(evalf(t[n+1,j]))<=tol*.1 then t:=DeleteColumn(t,j);j:=j-1; fi; flag:=true; fi; j:=j+1; od; if flag=false then t:=Matrix([t,Column(df2,i)]); fi; od; #print(" add_df_Output",t); t; end: # # Add two differentials, df1 - df2 # sub_df := proc (df1::Matrix, df2::Matrix, tol:=0) local ex, flag, t,n,i,j; #print(" sub_df_input",df1,df2); n:= Dimension(df1)[1]-1; if Equal(df2,Matrix(n+1,1)) then return df1; fi; t:=copy(df1); ex:= 1..n; for i to Dimension(df2)[2] do flag:=false; j:=1; while j<= Dimension(t)[2] do #print( i,j ); if Equal(t[ex,j],df2[ex,i]) then t[n+1,j]:= t[n+1,j] - df2[n+1,i]; if abs(t[n+1,j])<=tol*.1 then t:=DeleteColumn(t,j);j:=j-1; fi; flag:=true; fi; j:=j+1; od; if flag=false then t:=Matrix([Column(df2,i),t]); t[n+1,1]:= - t[n+1,1]; fi; od; #print(" sub_df_Output",t); t; end: # # Multiply a differential by a scalar value # scalar_df := proc(df::Matrix, s ) local n1,k, t; n1:= Dimension(df)[1]; t:=df; for k to Dimension(df)[2] do t[n1,k]:= s * t[n1,k]; od; t; end: # # `Differentiate` a differential wrt variable i # diff_df := proc( df::Matrix, i::integer) local n,t, m, ci, j, k; #print("diff_df_input:",df,i); t:= Matrix(df); n:= Dimension(t)[1]-1; if i>n then lprint("Differentiation out of range"); return df; fi; m:=Dimension(t)[2]; # nullify terms of deg=0 wrt i ci:=NULL; for j to m do if t[i,j]=0 then ci:= ci, j; fi; od;#nullif t:=DeleteColumn(t,[ci]); m:=Dimension(t)[2]; if m=0 then return Matrix(n+1,1); fi; for j to m do # Differentiate t[i,j]:= t[i,j]-1; od; # print("diff_df_output:",df,i); t; end: # # Check of the leading term of df appears as a term in LT. # (not used) has_term:= proc( df, LT::list) local d, i, m; #print("has_term_Input:", df, LT ); d:= Dimension(df); for m to d[2] do for i in LT do if convert([seq(i[s]=df[s,m],s=1..(d[1]-1))], `and` ) then return true; fi; od; od; #print("has_term_Output:", false); false; end: # # Produce candidate set for the Integration Method # cand_set:= proc( DD::list ) local n, df, t, m, ci,i, j, k, L; n:= Dimension(DD[1])[1]-1; #print(" cand_set In:", Matrix()); L:=NULL; for i to n do ## CAUTION: It does not work the other way!!! for df in DD do t:= copy(df); m:=Dimension(df)[2]; #print(`next integration: `, d[i]); ci:=NULL; for k from i+1 to n do # set i+1..n diffs to zero for j to m do if t[k,j]>0 then ci:= ci, j; fi; od; od;#set ci:=ListTools:-MakeUnique([ci]); if nops(ci)=m then t:=Matrix(n+1,1); else t:=DeleteColumn(t,ci); fi; L:= L, int_df(t,i); od: od; # print(" cand_set In:", Matrix(), # "Out:", Matrix()); [L]; end: # # Compute corank by SVD # ncorank := proc( M::Matrix, tol:=1e-3) nops(nkernel(M,tol) ); end: ## ## Compute Numerical Kernel using SVD ## nkernel := proc(M::Matrix, tol:=1e-3) local m,n,U,Sv,V,K,i,j; #print(`nkernel_Input:`, M ); if tol=0 then return NullSpace(M);fi; U,Sv,V:= MTM:-svd(M); n,m:= Dimension(M); if n=1 then if Sv<=tol then return [V]; else return {}; fi; fi; #print( Sv, V, n, m ) ; K:=NULL; for i to min(n,m) do if abs(Sv[i,i])<=tol then K:= K, V[..,i]; fi; od; if min(n,m) < m then for i from min(n,m)+1 to m do K:= K, V[..,i]; od; fi; #print("nkernel_Output:",K, "Maple:",NullSpace(evalf(M)) ); if nops([K])=0 then return {}; fi; [K]; end: ##################################################################### # # Macaulay's matrix dimension for s eq in n vars in depth dp # mac_dim:= proc(s::integer,n::integer,dp::integer ) s*binomial(dp-1+n,dp-1), binomial(dp+n,dp); end: # # Step in depth dp of Macaulay's method for the dual basis # macaulay_step := proc(f::list, var::list, z::list, dp::integer, tol:=0, BB:=[]) local E, n,NN,R,p,row, r, c, M, mc, j, ci, mr, k, ri, i, t; n:= nops(var); if dp<0 then return [0]; fi; if dp=0 then return [Matrix()]; fi; r:= nops(f)*binomial(dp-1+n,dp-1); c:= binomial(dp+n,dp); NN:=[seq(0,s=1..n)]; mc:=1; if BB<>[] then #MM2011 c:= c-nops(BB); mc:=0; NN:=NULL; fi; M:= Matrix(r,c); E:= exp_of(BB,var); for j to dp do#depth ci:= first_mon( j , n); c:=0: while ci<>NULL do #order of diff(column) if BB<>[] and member(ci,E) ## MM2011 then ci:= next_mon( ci , j); next; fi; c:=c+1; NN:=NN,ci;# keep column indices mr:=0; for k to j do#multiplier degree+1 ri:= first_mon(k-1, n); r:=0: while ri<>NULL do#multiplier monomials r:=r+1; for i to nops(f) do# f_i t:=diff_poly( mul((var[s]-z[s])^ri[s],s=1..n)*f[i], var, Matrix() ); #print(`at`,mr + i,mc + c); M[mr + i,mc + c]:= eval(t,[seq(var[s]=z[s],s=1..n)]); od; mr := mr+ i-1; ri:= next_mon( ri , k-1); od; od; ci:= next_mon( ci , j); od; mc:= mc + c; #loop ends in c-value Ok od; c:= Dimension(M)[2]; ## Condition (iii) MM2011 c:=NULL; if BB <> [] then R:=NULL; c:=NULL; fi; print("Dimension:",Dimension(M)); R:= nkernel( M , tol); p:=NULL; for row in R do #M:=NULL; M:=Matrix(n+1,1); for i to nops([NN]) do if abs(row[i])>tol then t:= Matrix(); #M:=M, t; M:=add_df(M,t, tol); fi; od; p:= p, Matrix([M]); od; #[p]; sort([p], rlex_df ); end: # # Step of the Integration Method given partial dual basis DD # and (optionally) partial standard basis BB of the quotient # mourrain_step := proc( f::list, var::list, z::list,DD::list, tol:=0, BB:=[]) local sd,un,c,t1,t2,tt1,tt2,p, k, row,R,rows,rr,n, M, NN, DIF, s, ro, i,j; n:= nops(var); sd:=nops(DD); #print(" mourrain_step_Input",DD); # Vector Lambda NN:= cand_set( DD ); # print( Matrix() ); #print("NN=", NN); s:=nops(NN); ###### Condition (i) ro:=0;# number of rows added #c2:=combinat:-choose(n,2); M:=Matrix(nops(f),s); un:=[cat(a,1..s)]; #cnt:=0; c:=first_comb(n,2); #for c in c2 do #all combinations while c <> NULL do #cnt:=cnt+1; #print( "diffing", seq(diff_df(DD[j], c[1]),j=1..sd)); #print( "diffing", seq(diff_df(DD[j], c[2]),j=1..sd)); t1:=[ seq(diff_df(DD[j], c[2]),j=1..sd)]; t2:=[ seq(diff_df(DD[j], c[1]),j=1..sd)]; #print( t1,t2); tt1:= to_polynomial(t1,var); tt2:= to_polynomial(t2,var); #print( Matrix(),Matrix() ); #print("OK", s, d); p:=0; for k to sd do #print("--", k, p); p:= p + un[(c[1]-1)*sd + k] * tt1[k] - un[(c[2]-1)*sd + k] * tt2[k]; # if not has_term(t1[k] , LT) then # p:= p + un[(c[1]-1)*d + k] * tt1[k]; # fi; # if not has_term(t2[k] , LT) then # p:= p - un[(c[2]-1)*d + k] * tt2[k]; # fi; od; #print("p=", p); #print(expand(p,var)); # if p<>0 then rows:= Array([lstcoefof(p,var)]); # fi; #print(rows); R:=NULL; for row in rows do #rr:= [seq(0,u=1..s)]; rr:= Array([seq(0,u=1..s)]); #lists don't work for big dimension for i to s do rr[i]:= coeffof(cat(a,i),row,un); ### #rr:=subsop(rr[i]=coeffof(cat(a,i),row,un),rr); od; #R:=R,rr; R:=R, convert(rr,`list`); od; if R<>NULL then M:= Matrix( ); ro:= ro + nops([R]); #print("M=",M); fi; c:=next_comb(c,n); od;#all combinations ##### Condition (ii), nops(f) constraints DIF:= to_function(NN,var); for i to nops(f) do # f[i] for j to s do M[ro+i,j]:= eval( app_funct( DIF[j], f[i]) , [seq(var[k]=z[k],k=1..n)]); od; od; ##### Condition (iii), MM2011 # ADD CONDITION: DUALS OF BB MUST VANISH, nops(DD)-1 constaints c:=NULL; if BB <> [] then R:=NULL; c:=NULL; t1:= expand(add( un[k]*to_polynomial(NN[k],var), k=1..s)); for i from 2 to nops(BB) do tt1:= coeffof(BB[i],t1,var); if nops(tt1)=1 then for j to s do if un[j]=tt1 then c:=c, j; break; fi; od; else rr:= Array([seq(0,u=1..s)]); for j to s do rr[j]:= coeffof(cat(a,j),tt1,un); od; R:=R,convert(rr,`list`); fi; od; if R<>NULL then M:= Matrix( ); #print(`Condition III added`, nops([R]),`rows` ); fi; if c<>NULL then c:= [c]; M:= DeleteColumn(M,c); NN:=subsop(seq(c[k]=NULL,k=1..nops(c)),NN); #print(`Condition III removed`,nops(c),`cols` ); fi; fi; #print("Dimension:",Dimension(M)); R:= nkernel( M, tol ); #print(`Found`, nops(R). `vectors`); #### Computing new elements p:=NULL; for row in R do M:=Matrix(n+1,1); for i to nops(NN) do if abs(evalf(row[i]))>=tol then t1:= copy(NN[i]); for j to Dimension(t1)[2] do t1[n+1,j]:= t1[n+1,j]*row[i]; od; M:=add_df(M,t1,tol); fi; od; p:= p, remove_zero(M); #p:= p, M; od; [p]; end: # # Macaulay' s 1916 algorithm to compute dual basis # macaulay0:= proc( f::list, var::list, z::list, tol:=0 ) local mu, DD,t,ind, n,j,i, flag; n:= nops(var); DD:= [ ]; i:=0; mu:=0; while true do t:= macaulay_step(f,var,z,i, tol); #print(mu,nops(t), "step=",t); if mu=nops(t) then #print("Yes!", mu,nops(DD)- 1 ); break; else DD:= [ op(t) ]; mu:=nops(t); i:=i+1; fi; od; sort( DD, rlex_df ); end: # # Macaulay's 1916 Algorithm with MM11 improvement, # Output is the primal-dual structure # macaulay:= proc( f::list, var::list, z::list, tol:=0 ) local mu, BB, DD,t, t1, n,i ; n:= nops(var); DD:= [ ]; BB:= []; i:=0; mu:=0; while true do t:= macaulay_step(f,var,z,i, tol, BB); #print(mu,nops(t), "step=",t); if 0=nops(t) then break; else DD:= [ op(DD), op(t) ]; ##### update partial standard basis t:= coefmatrix(t,var); t1:= cmaximal_minor( t[1], tol ); BB:= [ op(BB), op(t[2][t1]) ]; i:=i+1; #print(`B=`,BB); print("Added:", nops(t1) ); fi; od; DD,BB; end: # # Primal-Dual structure by Mourrain's 97 algorithm, with MM11 improvement # mourrain := proc( f::list, var::list, z::list, tol:=0 ) local BB, DD,t,t1, n; n:= nops(var); DD:=Matrix(n+1,1); DD[n+1,1]:=1; DD:= [ DD ]; BB:= [1]; while true do t:= mourrain_step(f,var,z,DD, tol, BB); #print(mu,nops(t)); #print("step=",t); if 0=nops(t) then break; else #print(`new`, nops(DD),nops(t)); DD:= [ op(DD), op(t) ]; ##### update Partial Standard Basis #print( Matrix() ); t:= coefmatrix(t,var); t1:= cmaximal_minor( t[1], tol ); BB:= [ op(BB), op(t[2][t1]) ]; #print("Added:", nops(t1) ); fi; od; DD, BB; end: # # Dual basis by classic Mourrain's '97 `Integration Method` # mourrain0 := proc( f::list, var::list, z::list, tol:=0 ) local mu, DD,t, n ; n:= nops(var); DD:=Matrix(n+1,1); DD[n+1,1]:=1; DD:= [ DD ]; mu:=0; while true do t:= mourrain_step(f,var,z,DD, tol); #print(mu,nops(t)); #print("step=",t); if mu=nops(t) then #print("Yes!", mu,nops(DD)- 1 ); break; else DD:= [DD[1], op(t) ]; mu:=nops(t); print("mult. is now:", mu ); fi; #cnt:=cnt+1;if cnt>3 then break; fi; od; #sort(DD, cat(comp,_df) ); sort(DD, rlex_df ); end: ######################################################## ######################################################## #### Applications ######################################################## ######################################################## # # Computes Normal form of p, given the primal-dual pair (BB,AA) # normal_form := proc ( p, var::list, z::list, BB::list, AA::list) local Sb, nf, i, m, dfs; m:=nops(BB); nf:=0; Sb:= seq(var[k]=z[k], k=1..nops(z)); for i to m do nf:= nf + eval( app_funct( to_function( BB[i],var ) , p ), [Sb]) * AA[i] ; od; nf; end: # # Computes the matrix of multiplication in of R/Q_z # mult_table := proc ( f::list, var::list, z::list) local BB,AA, m, MM,j,i ; BB, AA := basis_pair(f,var,z); m:= nops(BB); MM:= Matrix(m); for i to m do for j to m do MM[i,j]:= normal_form( AA[i]*AA[j], var, z, BB,AA ); od; od; MM; end: # # Return sign of number # sign_of := proc (k) if k=0 then return 0; fi; sign(k); end: # # Return the number of half-branches of real curve f=0 at point z # curve_branches := proc ( f::list, var::list, z::list ) local n, Sb,m,i,MM,j,JJ, ff, BB, AA, ev, Phi; n:= nops(var); Sb:= [seq(var[k]=z[k],k=1..n)]; ff:= [ op(f), Determinant(VectorCalculus:-Jacobian([op(f),add((var[k]-z[k])^2,k=1..n)] , var ))]; #print(`Map is ff= `, ff); BB, AA := basis_pair(ff,var,z ); m:= nops(BB); #print(to_polynomial(BB) ); JJ:= Determinant ( VectorCalculus:-Jacobian( ff , var ) ); #print(`J= `, JJ); # find positive quad form for i to m do ev:= eval( app_funct( to_function( BB[i],var ), JJ ), Sb) ; if ev<>0 then if ev>0 then Phi:= BB[i]; else print(`negative`); fi; break; fi; od; print(`Phi=`, to_polynomial(Phi) ) ; # Compute Matrix representation of quadratic form Phi Phi:= to_function(Phi,var); MM:= Matrix(m); for i to m do for j to m do MM[i,j]:= eval( app_funct( Phi , normal_form( AA[i]*AA[j], var, z, BB,AA ) ), Sb ); od; od; #print(MM); #Compute signature MM:=Eigenvalues( evalf(MM) ); #print(`Eigenvalues=`, MM); 2*add( sign_of(MM[i]), i=1..m ); end: # # Compute topological degree of mapping ff at point z # topological_degree := proc ( ff::list, var::list, z::list, tol:=0) local n, Sb,m,i,MM,j,JJ, BB, AA, ev, Phi; n:= nops(var); Sb:= [seq(var[k]=z[k],k=1..n)]; BB, AA := basis_pair(ff,var,z, tol ); m:= nops(BB); JJ:= Determinant ( VectorCalculus:-Jacobian( ff , var ) ); # find positive quad form for i to m do ev:= eval( app_funct( to_function( BB[i],var ), JJ ), Sb) ; if abs(ev)> tol then Phi:= BB[i]; break; fi; od; print(`got Phi=`, to_polynomial(Phi) ) ; # Compute Matrix representation of quadratic form Phi Phi:= to_function(Phi,var); MM:= Matrix(m); for i to m do for j to m do MM[i,j]:= eval( app_funct( Phi , normal_form( AA[i]*AA[j], var, z, BB,AA ) ), Sb ); od; od; #print(MM); #Compute signature and take absolute value Eigenvalues( evalf(MM) ); abs( add( sign_of(MM[i]), i=1..m ) ); end: # # Create coefficient matrix of differentials BB # coefmatrix := proc (BB,v:=[] ) local SBB,i,j,A,c,n,var; n:=Dimension(BB[1])[1]-1; if v=[] then var:=[seq(d[s],s=1..n)]; else var:=v; fi; c:=to_polynomial( BB,var); #print("coefmat",c); SBB:= lstmonsof(c): A:= Matrix(nops(BB), nops(SBB) ) : for i to nops(c) do for j to nops(SBB) do A[i,j]:= coeffof(SBB[j],c[i], var ): od: od: [ A, SBB ] ; end: # # Input: coefficient matrix A, indexing polys SBB # Output: Polynomials induced by A's rows # matrix2polys:= proc ( A, SBB ) local u,v,i,j,t,L; u,v := Dimension(A); L:=NULL; for i to u do t:=0; for j to v do t:= t + A[i,j]* SBB[j]; od; L:=L, t; od; [L]; end: # # Input: coefficient matrix A, indexing polys SBB # Output: Differentials induced by A's rows # matrix2diffs:= proc ( A::Matrix, SBB::list, var::list, tol:=0 ) local n,k,df,u,v,i,j,t,L; #print("matrix2diffs_IN:", A, SBB, var); u,v := Dimension(A); n:=nops(var); L:=NULL; t:=Matrix(n+1,1); for k to u do df:= Matrix(n+1,1); for i to v do if (not abs(A[k,i]) <= tol ) then t[n+1,1]:=A[k,i]; for j to n do t[j,1] := degree (SBB[i], var[j]); od; df:= add_df(df,t,tol); fi; od; L:=L, copy(df); od; #print("matrix2diffs_OUT:", L ); [L]; end: # # Compute a primal-dual orthogonal pair # basis_pair:= proc( f::list, var::list, z::list, tol:=0 ) local Db,A,i,C,ind; Db,A:= mourrain(f,var,z, tol); C:= coefmatrix( Db, var ); ind:=NULL; for i to nops(A) do ind:=ind, ListTools:-Search(A[i], C[2] ) ; od; ind:=[ind]; C[1]:= ReducedRowEchelonForm( Matrix([ C[1][..,ind], DeleteColumn(C[1],ind) ])); A:=C[2][ind]; C[2]:=[op(A),op(subsop(seq( ind[j]=NULL,j=1..nops(ind)),C[2]))]; matrix2diffs( C[1], C[2],var, tol), A; end: # # Given dual basis BB, find AA. (under constuction) # make_basis_pair:= proc( BB::list, var::list) local Db,A,B,i,C,ind; C:= coefmatrix( BB, var ); ind:=NULL; for i to nops(A) do ind:=ind, ListTools:-Search(A[i], C[2] ) ; od; ind:=[ind]; #print( Matrix([ C[1][..,ind], DeleteColumn(C[1],ind) ]) ); C[1]:= ReducedRowEchelonForm( Matrix([ C[1][..,ind], DeleteColumn(C[1],ind) ])); A:=C[2][ind]; C[2]:=[op(A),op(subsop(seq( ind[j]=NULL,j=1..nops(ind)),C[2]))]; #print(C[2]); #print(C[1]); #matrix2polys( ReducedRowEchelonForm(C[1]), subs(seq(var[k]=d[k],k=1..nops(var)), C[2])), matrix2diffs( C[1], C[2],var, tol), A; end: #################################################### #################################################### ### Auxiliary functions #################################################### #################################################### # # Computes the `Cyclic-n` polynmial system in variables x[1]..x[n] # cyclicn := proc( n ) local m,i,j, t, CN; CN:=NULL; m:=n-1; for i from 0 to m-1 do # t=f_i t:=0; for j from 0 to m do # sum t:= t+ mul( x[(k mod n)+1], k= j..(i+j)); od; CN:= CN, t; od; CN:= CN, 1- mul( x[k], k=1..n); [CN]; end: # # Next r-combination of [1..n] after c in lex-order # nops(c)=r # next_comb := proc( c::list, n::integer) local nc,r,i,j; r:= nops(c); if c=[seq(n-r+1..n)] then return NULL; fi; #finished nc:=c; i:=r; while c[i]=n-r+i do i:=i-1; od; nc[i]:=c[i]+1; for j from i+1 to r do nc[j]:= nc[i]+j-i; od; nc; end: # # Lex-first r-combination of [1..n] # first_comb:= proc(n::integer,r::integer) if r<=n then return [seq(1..r)]; else print(`Warning: r>n combination requested`,r,n); return NULL; fi; end: # # RevLex-first degree r monomial in n variables # first_mon:= proc(r::integer,n::integer) return [r,seq(0,k=1..n-1)]; end: # # RevLex-Next degree r monomial in n variables after c # next_mon := proc( c::list, r::integer) local nc,n,i,j; n:= nops(c); if c[n]=r then return NULL; fi; #finished nc:=c; i:=1; while c[i]=0 do i:=i+1; od; nc[i+1]:=nc[i+1]+1; #complete first non-zero element from the right with nc[i]-1 if i=1 then nc[i]:=nc[i]-1; else nc[1]:=nc[i]-1; nc[i]:=0; fi; nc; end: ##################################### end use; lprint(`Multiple Roots package, A. Mantzaflaris & B. Mourrain, 2011`); end: ## ## procedure analyzer ## mmint := proc( f, {level::posint:=2}, {toggle::list(posint):=[]}, {header::truefalse:=false}, {extras::string:=""} ) local fname,k,oldiface,executable; fname := cat(kernelopts('homedir'),kernelopts('dirsep'),"minttemp.txt"); #fname := "/tmp/minttemp.txt"; # or some other temp location oldiface := interface('verboseproc'); interface('verboseproc'=3); writeto(fname); printf("%a := %a:\n", f, eval(f)); writeto(terminal): fclose(fname); interface('verboseproc'=oldiface); executable:=`if`(kernelopts('platform')="windows", cat("\"",kernelopts('bindir'),"\\mint","\""), cat(kernelopts('mapledir'),kernelopts('dirsep'),"bin/mint")); k:=ssystem(cat(executable, cat(" -i ",min(level,4)," "), seq(cat(" -t ",x," "), x in remove(t->t>32,toggle)), `if`(not(header)," -q ",NULL), " ",extras," ","\"",fname,"\"")); if k[1]>=0 then printf("%s",k[2]); end if; NULL; end proc: