reset()
import itertools, multiprocessing
#----------------------------------------------------------------------------------------------------------------
#function computing the structure constants and then the commutator matrix
#----------------------------------------------------------------------------------------------------------------
def CommMatrix(matrix_basis,K):
n=len(matrix_basis)
R=PolynomialRing(K,'x',n)
var=vector(R.gens())
matrix_basis = [Matrix(K, b) for b in matrix_basis]
A = Matrix(K, [b.list() for b in matrix_basis])
if A.rank() != n:
raise TypeError('matrices are linearly dependent')
mult = lambda x,y: x*y - y*x
return Matrix([[var*A.solve_left(vector(mult(b,c).list())) for c in matrix_basis] for b in matrix_basis])
#----------------------------------------------------------------------------------------------------------------
#Fixing basis. Here we choose the basis e_{ij} i \neq j with e_{11} - e{22}, e_{11} - e{33}, e_{11} - e{44}.
#----------------------------------------------------------------------------------------------------------------
E = [[matrix(ZZ,[[1 if (i,j) == (k,l) else 0 for i in range(4)] for j in range(4)]) for k in range(4)] for l in range(4)]
basis=[E[i][j] for i in range(4) for j in range(4) if i != j] + [E[0][0] - E[1][1], E[0][0] - E[2][2], E[0][0] - E[3][3]]
#----------------------------------------------------------------------------------------------------------------
#create the commutator matrix
#----------------------------------------------------------------------------------------------------------------
M=CommMatrix(basis,ZZ)
m=len(basis)
S=MatrixSpace(PolynomialRing(ZZ,'x',m),m,m)
R=MatrixSpace(ZZ,m,m)
M=S(M)
#----------------------------------------------------------------------------------------------------------------
# Examples (de-comment as needed)
#----------------------------------------------------------------------------------------------------------------
#fix the prime
p = 3
#----------------------------------------------------------------------------------------------------------------
#pre-print example
#----------------------------------------------------------------------------------------------------------------
u = [10, 21, 0, 0, 9, 21, 0, 9, 10, 0, 0, 0, 9, 0, 0]
# starting level
lev = 3
#----------------------------------------------------------------------------------------------------------------
# defining function X |-> M(u) + p^lev * M(X).
# Note: This is needed because using matrices over Polynomial ring and evaluating them is much slower
#----------------------------------------------------------------------------------------------------------------
c_matrix_original = M(u)
c_matrix_rows_x = M.rows()
# defining the string p^lev * M(X)
var('y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14')
c_matrix_list_y = [[p^(lev) * c_matrix_rows_x[i][j].subs(x0 = y0,x1 = y1, x2 = y2,x3 = y3,x4 = y4,x5 = y5,x6 = y6, x7 = y7, x8 = y8, x9 = y9, x10 = y10, x11 = y11, x12 = y12, x13 = y13, x14 = y14) for j in xrange(m)] for i in xrange(m)]
c_matrix_string = str(c_matrix_list_y)
del y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14
exec('''
def M1(list):
y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14 = list
return c_matrix_original + matrix(ZZ, {formula})''').format( formula = c_matrix_string)
#----------------------------------------------------------------------------------------------------------------
#print data
#----------------------------------------------------------------------------------------------------------------
print 'Prime: ', p
print 'Level: ', lev
print('Running through all {size} lifts of').format(size = p^m)
b = sum(basis[i] * u[i] for i in range(m))
print(b)
print ('Elementary divisors of the original commutator matrix over the integers:\n {el}').format( el = M(u).elementary_divisors())
el_div_orig = [e if e.mod(p^(lev))!=0 else 0 for e in M(u).elementary_divisors() ]
print ('Elementary divisors of the original commutator matrix over the integers modulo {prime}^{level}:\n {el}').format(
el = el_div_orig, prime = p, level = lev)
ncpus = multiprocessing.cpu_count()
print 'Number of cores in use: ', ncpus
#----------------------------------------------------------------------------------------------------------------
# defining parallel counting function
#----------------------------------------------------------------------------------------------------------------
@parallel(ncpus)
def f(sample):
h,i,j,k,l=0,0,0,0,0
for y in sample:
x = ZZ(y).digits(p, padto = m)
r=len([e for e in M1(x).elementary_divisors() if e.mod(p^(lev + 1))!=0])
if r==0:
h+=1
if r==6:
i+=1
if r==8:
j+=1
if r==10:
k+=1
if r==12:
l+=1
return vector([h,i,j,k,l])
#----------------------------------------------------------------------------------------------------------------
#main loop
#----------------------------------------------------------------------------------------------------------------
space=xrange(p^m)
inp=[[space[j] for j in xrange(i, len(space), ncpus)] for i in xrange(ncpus)]
res = vector([0,0,0,0,0])
for arg,out in f(inp):
if out == 'NO DATA':
out = f.func(arg[0])
res += out
#----------------------------------------------------------------------------------------------------------------
# Printing results
#----------------------------------------------------------------------------------------------------------------
e = 15 - el_div_orig.count(0)
if e==0:
pos = 0
if e==6:
pos = 1
if e==8:
pos = 2
if e==10:
pos = 3
if e==12:
pos = 4
print('\nRESULT:')
print('There are {number} lifts preserving the number of non-zero elementary divisors of the comm. matrix.').format(number = res[pos])
print('\nIn detail:')
print(' # non-zero el. div. # of lifts')
print(' 6 {lifts}').format(lifts = res[1])
print(' 8 {lifts}').format(lifts = res[2])
print(' 10 {lifts}').format(lifts = res[3])
print(' 12 {lifts}').format(lifts = res[4])
Like this:
Like Loading...
Michele Zordan's Webpage 