(original) (raw)

# 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.<a,b> = 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]=Tgen[i].coefficient(a[j//m]).coefficient(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 = cartesian_product([list(A),list(B)]).list() AB = [x[0]*x[1] for x in 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] = AB[i].coefficient(a[da]).coefficient(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([reduce(h, [g]) for h in 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,T0,V,k): global Nsols K = G[0,0].parent() # first put T0 in row-echelon form T = copy(T0) # to avoid modifying T0 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 = [reduce(g,T) for g in 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([postprocess(w,G) for w in Wsols], max_level=1) print ("Number of solutions:", len(sols)) # interpret solutions sols = [interpret(s,V) for s in 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); list(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); list(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); list(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.<a,b,c,u,v,w> = 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 = cartesian_product([A,U]).list() # AU = map(lambda x: x[0]*x[1], AU) AU = [x[0]*x[1] for x in 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] = AU[i].coefficient(A[da]).coefficient(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]=Tgen[i].coefficient(A[j//3]).coefficient(U[j%3]) return G, T, V </a,b,c,u,v,w></a,b>