# This is a implementation in Sage of Algorithm 1 from "Finding Optimal
# Formulae for Bilinear Maps" by Razvan Barbulescu, Jérémie Detrey,
# Nicolas Estibals and Paul Zimmermann, http://hal.inria.fr/hal-00640165
# This function returns the generator set G, the target set T, and the
# basis of the vector space V, for n x m multiplication over GF(2)
# example: G, T, V = polymulgf2(2,3)
def polymulgf2(n,m):
K = GF(2)
S. = InfinitePolynomialRing(K)
V = [a[i]*b[j] for i in range(n) for j in range(m)]
Tgen = [sum([a[i]*b[d-i] for i in range(max(0,d-m+1),min(n,d+1))])
for d in range(n+m-1)]
T = matrix(K, n+m-1, n*m)
for i in range(n+m-1):
for j in range(n*m):
T[i,j]=SR(Tgen[i]).coefficient(SR(a[j//m])).coefficient(SR(b[j%m]))
A = [a[i] for i in range(n)]
A = Combinations(A).list()
if A[0] <> []:
raise ValueError, "A[0] <> []"
del A[0]
A = map(sum, A)
B = [b[j] for j in range(m)]
B = Combinations(B).list()
if B[0] <> []:
raise ValueError, "B[0] <> []"
del B[0]
B = map(sum, B)
AB = CartesianProduct(A,B).list()
AB = map(lambda x: x[0]*x[1], AB)
G = matrix(K,len(AB),n*m)
for i in range(len(AB)):
for j in range(n*m):
da = j // m
db = j % m
if V[j] <> a[da]*b[db]:
print V[j], a[da]*b[db]
raise ValueError, "V[j] <> a[da]*b[db]"
G[i,j] = SR(AB[i]).coefficient(SR(a[da])).coefficient(SR(b[db]))
return G, T, V
# return p, a list of swaps
def rowechelon(T):
m = T.ncols()
p = []
for i in range(T.nrows()):
# look for pivot in row i
for j in range(m):
if T[i,j] <> 0:
# found pivot
p.append([i,j])
T.swap_columns(i,j)
# reduce further rows
for k in range(i+1,T.nrows()):
T.add_multiple_of_row(k,i,-T[k,i]/T[i,i])
break
return T, p
def reduce(g,T):
nm = len(T[0])
for t in T:
# look for pivot
for j in range(nm):
if t[j] <> 0:
if g[j] <> 0:
g = g - g[j]/t[j] * t
break
return g
def pivot(g):
for j in range(len(g)):
if g[j] <> 0:
return j
# remove zero vectors from l and duplicate vectors
def uniq(l):
keep = []
for i in range(len(l)):
if not l[i].is_zero():
ok = True
for j in keep:
if l[i] == l[j]:
ok = False
break
if ok:
keep.append(i)
return [l[i] for i in keep]
Nsols = 0
# W is a list
# H is a matrix
# G is a list
def enlarge_to_dim_k(W,H,k,G):
global Nsols
sols = []
if len(W) == k:
WG = []
for g in G:
if reduce(g,W).is_zero():
WG.append(g)
r = matrix(WG).rank()
if r==k:
sols = [W]
Nsols = Nsols + 1
print Nsols, W
sys.stdout.flush()
else:
while H <> []:
# we take g with the leftmost pivot
imin = 0
jmin = pivot(H[0])
for i in range(1,len(H)):
j = pivot(H[i])
if j < jmin:
imin = i
jmin = j
g = H[imin]
H[imin] = H[0]
del H[0]
newH = uniq(map(lambda h:reduce(h, [g]), H))
sols = sols + enlarge_to_dim_k(W + [g], newH, k, G)
return sols
def postprocess(W,G):
sols = []
k = len(W)
# first reduce G by W
G = [g for g in G if reduce(g,W).is_zero()]
# then look for all k-subsets of G
for s in subsets(G):
if len(s) == k:
if matrix(s).rank() == k:
sols = sols + [s]
return sols
def interpret_aux(s,V):
return factor(sum([s[i]*V[i] for i in range(len(s))]),proof=False)
def interpret(sol,V):
return map(lambda s: interpret_aux(s,V), sol)
# This function implements Algorithm 1 from the paper. It takes as input the
# generator set G, the target set T (both should be matrices), the basis V
# of the ambient vector space (which gives the meaning for the columns of G
# and T), and the target bilinear rank k.
# It returns a set sols of solutions, where each "solution" is a list of k
# generators from G that span the solution space T.
# Example: G,T,V=polymulgf2(3,3); sols=Algorithm1(G,T,V,6); sols[0]
# [a_0*b_0, a_2*b_2, (b_0 + b_2)*(a_0 + a_1), (b_1 + b_2)*(a_0 + a_2), (b_0 + b_1)*(a_1 + a_2), (b_0 + b_1 + b_2)*(a_0 + a_1 + a_2)]
def Algorithm1(G,T,V,k):
global Nsols
K = G[0,0].parent()
# first put T in row-echelon form
T, p = rowechelon(T)
# swap columns of G according to p
for s in p:
G.swap_columns(s[0],s[1])
tmp=V[s[0]]; V[s[0]]=V[s[1]]; V[s[1]]=tmp
print "Permuted columns:", V
print "Initial number of generators:", G.nrows()
sys.stdout.flush()
G0 = map(lambda g: reduce(g,T), G)
G0 = uniq(G0)
print "After reduction by T, remains:", len(G0)
sys.stdout.flush()
Nsols = 0
Wsols = enlarge_to_dim_k(list(T), G0, k, G)
print "Number of W solutions:", len(Wsols)
sys.stdout.flush()
Wsols2 = []
for w in Wsols:
ok = True
for u in Wsols2:
if matrix(K,w).stack(matrix(K,u)).rank() == k:
ok = False
break
if ok:
Wsols2.append(w)
Wsols = Wsols2
print "Number of unique W solutions:", len(Wsols)
sols = flatten(map(lambda w: postprocess(w,G), Wsols), max_level=1)
print "Number of solutions:", len(sols)
# interpret solutions
sols = map(lambda s: interpret(s,V),sols)
return sols
# generates G, T, V for Exercise 1.18 from "Modern Computer Arithmetic"
# Richard Brent and Paul Zimmermann, Cambridge University Press, 2010.
# Produces Algorithm 1 input for Exercise1_18 over GF(p) for
# a*u+c*w,a*v+b*w,b*u+c*v if sign=1, and a*u-c*w,a*v-b*w,b*u-c*v if sign=-1
#
# For p=2, sign=1 or -1, k=5, there are 7 solution spaces W, and 42 solutions:
# G,T,V = Exercise1_18(2,1); sols=Algorithm1(G,T,V,5); sols[0]
# [(v + w) * a, (u + v) * c, u * (b + c), (u + v + w) * (a + c), w * (a + b)]
#
# For p=3, sign=1, k=4, there is 1 solution space W, and 1 solution:
# G,T,V = Exercise1_18(3,1); sols=Algorithm1(G,T,V,4); sols[0]
# [(u + v + w) * (a + b + c), (-u + v + w) * (-a - b + c), (-1) * (-u + v - w) * (a - b + c), (-1) * (u + v - w) * (-a + b + c)]
#
# For p=3, sign=-1, k=5, there are 234 solution spaces W, and 1404 solutions:
# G,T,V = Exercise1_18(3,-1); sols=Algorithm1(G,T,V,5); sols[0]
# [(v + w) * a, (u - v) * (b + c), (-1) * (u + v) * (-b + c), w * (a + b), (u + v - w) * (a - b + c)]
def Exercise1_18(p,sign):
R. = GF(p)[]
V = [a*u,a*v,a*w,b*u,b*v,b*w,c*u,c*v,c*w]
A = [a,b,c]
U = [u,v,w]
# we impose the first non-zero of alpha,beta,gamma to be 1
for alpha in range(p):
if alpha <> 0 and alpha <> 1:
continue
for beta in range(p):
if alpha == 0 and beta <> 0 and beta <> 1:
continue
for gamma in range(p):
if alpha == 0 and beta == 0 and gamma <> 0 and gamma <> 1:
continue
if alpha+beta+gamma > 1:
A.append(alpha*a+beta*b+gamma*c)
U.append(alpha*u+beta*v+gamma*w)
AU = CartesianProduct(A,U).list()
AU = map(lambda x: x[0]*x[1], AU)
G = matrix(GF(p),len(AU),len(V))
for i in range(len(AU)):
for j in range(len(V)):
da = j // 3
du = j % 3
if V[j] <> A[da]*U[du]:
print V[j], A[da]*U[du]
raise ValueError, "V[j] <> A[da]*U[du]"
G[i,j] = SR(AU[i]).coefficient(SR(A[da])).coefficient(SR(U[du]))
if sign==1:
Tgen = [a*u+c*w,a*v+b*w,b*u+c*v]
else:
Tgen = [a*u-c*w,a*v-b*w,b*u-c*v]
T = matrix(GF(p), len(Tgen), len(V))
for i in range(len(Tgen)):
for j in range(len(V)):
T[i,j]=SR(Tgen[i]).coefficient(SR(A[j//3])).coefficient(SR(U[j%3]))
return G, T, V