# Algorithms of the paper: # Symmetric Determinantal Representations in Characteristic 2 # Bruno Grenet, Thierry Monteil, Stéphan Thomassé # Organization: # - Routines # - I. Multilinear poly => List of couples (factor, ideal) # - II. List of couples => SDR # - III. SDR => List of factors modulo an ideal # - Algorithms of the paper ######################## # AUXILIARY ROUTINES # ######################## #====================# # Algs: Dico & Ideal # #====================# # To switch from two representations of ideals # Inputs: An ideal < x1^2+l1^2, x2^2+l2^2, ... > represented as (either) # - dictionary {x1:l1^2, x2:l2^2,...} # - list [ x1^2+l1^2, x2^2+l2^2, ... ] # - single element l^2 => represents < x1^2+l^2, x2^2+l^2, ... > # Outputs: # - Dico: A dictionary { x1:l1, x2:l2, ... } # - Ideal: A list [ x1^2+l1^2, x2^2+l2^2, ... ] def ListToDico(l): d={} for c in l: d[c[0]]=c[1] return d def Dico(arg): if isinstance(arg,dict): dico=ListToDico([(t,arg[t].sqrt()) for t in arg]) elif isinstance(arg,list): dico=ListToDico([(t.variable(),t.constant_coefficient().sqrt()) for t in arg]) else: dico=ListToDico([(t,arg.sqrt()) for t in R.gens()]) return dico def Ideal(arg): if isinstance(arg,dict): ideal=[t^2+arg[t] for t in arg] elif isinstance(arg,list): ideal=arg else: ideal=[t^2+arg for t in R.gens()] return ideal def SmallestNonMonomial(p): pv=p.variables() nv=len(pv) pm=p.monomials() #print nv for i in range(1,nv+1): #print i if len([m for m in pm if m.degree()==i])nv-1: ex[t]=False ex[t+1]=True else: trues=0 while ex[t]: t-=1 trues+=1 while not(ex[t]): t-=1 ex[t]=False for t1 in range(trues+1): ex[t+trues+1]=True ############################################################################ # I. From a multilinear polynomial to a list of couples (factor, ideal) # # ( Main algorithm ) # ############################################################################ #==================# # Alg: Preparation # #==================# # Input: multilinear polynomial p # list L of couples (l:linear poly,b={0,1}) # Ouput: multilinear polynomial p' # list L' of couples = L+[(l1,b1),...,(lk,bk)] # Prop: p' has valuation exactly 1 # p=(...(p'*lk mod ) * ...) * l1 mod # Ref: Algorithm 3 + Lemma 4.7 def Prep(p,l): R=p.parent() pv=p.variables() pm=p.monomials() # Linear poly: nothing to do if p.degree()<2: return p,l # "Full polynomial" (Lemma 4.4) elif len(pm)==2^(len(pv)): x=pv[0] px=p.monomial_coefficient(x) p0=p.constant_coefficient() p*=(px*x+p0) p=p.reduce([x^2]) # only x can be squared: p=p.reduce([t^2 for t in pv]) return Prep(p,l+[(R((px*x+p0)/p0^2),0)]) # Valuation = 0 (Lemma 4.5) elif p.constant_coefficient()<>0: #allMon=prod([1+t for t in pv]).monomials() #for i in range(len(allMon)): # if not(allMon[-i-1] in pm): # m=allMon[-i-1] # break m=SmallestNonMonomial(p) return ((m*p).reduce([t^2+1 for t in pv]),l+[(v,1) for v in m.variables()]) # Valuation > 1 (Lemma 4.5) elif p.monomials()[-1].degree()>1: m=R(pm[-1]//(pm[-1].variables()[0])) p*=m return (p.reduce([t^2+1 for t in pv]),l+[(v,1) for v in m.variables()]) # Valuation = 1: already "prepared" else: return (p,l) #========================================# # Main algorithm # # Common part to IsFactorizable & SymDet # #========================================# # Input: multilinear polynomial p # list L of couples (l:linear poly,b={0,1}) # Output: Boolean B # L+[(l1,b1),...,(lk,bk),(l,-1)] # Prop: If B==False, p is not factorizable # Else, p=(...(l*lk mod ) * ... ) * l1 mod # Refs: Algorithm 4 + Thm 4.8 + Thm 4.10 def Main(p,l): #print p,l (p,l)=Prep(p,l) #print p,l pv=p.variables() pm=p.monomials() # p is linear: OK if p.degree()<2: return (True,l+[(p,-1)]) # p is not linear (Lemmas 4.2 and 4.3) else: linPart=sum([t*p.monomial_coefficient(t) for t in pm if t.degree()<2]) alpha=linPart.lc() x=linPart.variables()[0] p0=p.coefficient(x) pTilde=(linPart/alpha*p0).reduce([t^2 for t in pv]) # p is factorizable iff p==pTilde and p0 is factorizable if p==pTilde: return Main(p0,l+[(linPart/alpha,0)]) else: return False,l+[(p,-1)] ########################################################### # II. From a list of couples (factor, ideal) to an SDR: # # (Second part of SymDet) # ########################################################### #=====================# # Alg: LinPoly2Matrix # #=====================# # Input: Linear polynomial p # Output: Matrix M with each entry: variable or constant # Prop: det(M)=p # Ref: Proposition 2.4 def LinPoly2Matrix(p): R=p.parent() pm=p.monomials() terms=[t*p.monomial_coefficient(t) for t in pm] if len(terms)==1: return matrix(R,[[p]]) elif len(terms)==2: return matrix(R,[[terms[0],1,0],[1,0,1],[0,1,terms[1]]]) size=2*len(pm)+1 M=matrix(R,size) for i in range(len(pm)): M[0,2*i+1]=M[2*i+1,0]=1; M[2*i+1,2*i+2]=M[2*i+2,2*i+1]=1; M[2*i+2,2*i+2]=terms[i] return M #================# # Alg: CleanDiag # #================# # Input: Matrix M, b={0,1} (diag entries of M have degree at most 1) # Output: Matrix M* # Prop: det(M)=det(M*) mod < x^2+b, ... > # Each variable appears at most once on the diagonal # Ref: Lemma 3.14 def CleanDiag(M,b): R=M.base_ring() n=M.ncols() for i in range(n-1): # If the entry has degree 1: if M[i,i].degree()>0: # Looking for another diagonal entry with the same variable for j in range(i+1,n): if M[j,j].degree()>0: vi=M[i,i].variable(0) vj=M[j,j].variable(0) if vi==vj: ci=M[i,i].lc() cj=M[j,j].lc() M=Add(M,i,j,(cj/ci).sqrt(),Dico(b)) return M #============# # Alg: Merge # #============# # Input: Matrices M & N, b={0,1} # Output: Matrix R # Prop: det(R)=det(M)*det(N) mod < x^2+b, ... > # Ref: Lemma 4.9 def Merge(M,N,b): R=block_diagonal_matrix(M,N) return CleanDiag(R,b) #==================# # Alg: List2Matrix # #==================# # Input: List L of couples (l:linear poly, b={0,1}) # L=[(l1,b1), ..., (lk,bk), (l,b)] # Output: Matrix M with entries of degree at most 1 # Prop: det(M) = ( ... ( l*lk mod ) * ... ) * l1 mod # Ref: Theorem 4.10 def List2Matrix(l): l.reverse() M=LinPoly2Matrix(l[0][0]) for i in range(1,len(l)): N=LinPoly2Matrix(l[i][0]) M=Merge(M,N,l[i][1]) return M ########################################### # III. From an SDR to a list of factors # ########################################### #=====================# # Algs: Add & Isolate # #=====================# # Input: Matrix M, index i, (index j, coeff alpha), dico {x:l, ... } # Output: Matrix M* # Prop: M* obtained after adding alpha*(Col i) to Col j and alpha*(Row i) to Row j (for Add) # isolating coefficient M[i,i] (for Isolate) # Operations performed modulo < x^2+l^2, ... > # Refs: Section 3.2: Necessary condition def Add(M,i,j,alpha,dico): R=M.base_ring() M.add_multiple_of_column(j,i,alpha) M.add_multiple_of_row(j,i,alpha) v=M[i,j].subs(dico) M[i,j]=M[j,i]=v return M def Isolate(M,i,dico): for j in range(M.ncols()): if j<>i and M[i,j]<>0: Mii=M[i,i].subs(dico) M=Add(M,i,j,M[i,j]/Mii,dico) return M #===========# # Alg: Diag # #===========# # Input: Matrix M, dico {x:l, ...} # Output: Diagonal matrix M* # Prop: det(M)=det(M*) mod < x^2+l^2, ... > # M* is diagonal, and diagonal entries have degree at most 1 # Ref: Theorem 3.8 def Diag(M,dico): R=M.base_ring() s=M.ncols() nonzeroes=[] # Looking for nonzero diag entry for i in range(s): if M[i,i]<>0: isolated=M.nonzero_positions_in_row(i)==[i] nonzeroes+=[(i,isolated)] # if M[i,i] is nonzero, with a nonzero square, isolate it! if not(isolated) and M[i,i].subs(dico)<>0: return Diag(Isolate(M,i,dico),dico) # All nonzero diag entries have a null square # If all diag entries are zero, det(M) is a constant if len(nonzeroes)==0: return matrix(R,[M.det()]) # Else, diag entries with nonzero squares are created (Lemma 3.13) for (i,b) in nonzeroes: if not(b): N=matrix(R,s+1) N[:s,:s]=M N[i,i]+=1 N[s,s]=N[i,s]=N[s,i]=1 N=Isolate(N,i,dico) return Diag(N,dico) return M #=========================# # Alg: DiagMatrix2Factors # #=========================# # Input: Diagonal matrix M, list [ x^2+l^2, ... ] # Output: List of factors of det(M) modulo < x^2+l^2, ... > # Prop: If the list of diagonal elements is returned, there can be some # duplicate factors. For instance, if x+1 appears twice, and the # ideal is , both occurences can be removed since (x+1)^2=1. def DiagMatrix2Factors(M,ideal): # First factor is the (only) constant factor factors=[1] # List all diag entries for f in M.diagonal(): # Test if this new factor annihilates another one for i in range(1,len(factors)): r=(f*factors[i]).reduce(ideal) if r.degree()<2: # det(M)=0 (should not happen!) if r==0: return [r] # If the product of 2 factors is a constant: # - update the constant factor # - remove the annihilated factor elif r.degree()<1: factors=[r*factors[0]]+factors[1:i]+factors[i+1:] break # Last case: prod of 2 factors has degree 1 else: factors=factors[:i]+[r]+factors[i+1:] break # If this is a really new factor, add it to the list else: factors+=[f] return factors ########################### # ALGORITHMS OF THE PAPER # ########################### #=============================================# # Algs: IsFactorizable, SymDet, Factorization # #=============================================# # Algorithms described in the paper: # # - IsFactorizable: Input: A multilinear polynomial p # Output: True/False # # - SymDet: Input: A multilinear polynomial p # Output: A matrix M s.t. p=det(M) # # - Factorization: # Input: A multilinear polynomial p # An ideal < x^2+l^2 ... > represented as # - { x:l^2, ... } # - [ x^2+l^2, ... ] # - l: < x1^2+l^2, x2^2+l^2, ... > # Output: List [L1, L2, ... ] s.t. p=L1*L2*... mod < x^2+l^2, ... > def IsFactorizable(p): b,l=Main(p,[]) return b def SymDet(p): b,l=Main(p,[]) if b: M=List2Matrix(l) return M else: return False def Factorization(p,arg): dico=Dico(arg) ideal=Ideal(arg) M=SymDet(p) if M==False: return False M=Diag(M,dico) return DiagMatrix2Factors(M,ideal)