Question: Improved the Speed of a Procedure in Modulo Two

Let A=[V₁, V₂ ,..., V_n ] be a list of binary vectors such that the length of V_i for i=1...n is the positive integer number m. For instance, in the following example we have n=6 and m=5.

A := [[1, 1, 1, 0, 0], [0, 1, 0, 1, 1], [1, 0, 1, 1, 1], [0, 1, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 1, 1, 1]]

Suppose that e_j with j=1...m are vectors of size m such that all entries of e_j 's are zero except the jth entry where is equal to 1. Now set S=[e₁,e₂,...,e_m]. For example by m=5 we get 

S:=[[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]]

First, we choose i,j in [1...m] such that i<>j and then we update S as follows S=[e₁,e₂,...,e_m,e_i+e_j mod 2]. For example, it follows from i=1 and j=2 that 

S:=[[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [1, 1, 0, 0, 0]]

Now consider the kth entry of A which is called A_k. Then we check that what is the minimum number of entries of S so that the summation of these entries mod 2 is equal to A_k. For example, for  S=[e₁,e₂,...,e_m,e₁+e₂ mod 2] we get [2, 3, 4, 3, 2, 4] which means the minimum number of entries of the S that are required to be added in order to obtain A₁ is 2 and so on. Therefore, for the given example we get 

E:=[[1, 2, [2, 3, 4, 3, 2, 4], 18, 7.615773106],
[1, 3, [2, 3, 3, 3, 3, 4], 18, 7.483314774],
[1, 4, [3, 3, 3, 3, 2, 4], 18, 7.483314774],
[1, 5, [3, 3, 3, 3, 3, 4], 19, 7.810249676],
[2, 3, [2, 3, 4, 2, 3, 3], 17, 7.141428429],
[2, 4, [3, 2, 4, 2, 2, 3], 16, 6.782329983],
[2, 5, [3, 2, 4, 3, 3, 3], 18, 7.483314774], 
[3, 4, [3, 3, 3, 2, 3, 3], 17, 7.], 
[3, 5, [3, 3, 3, 3, 3, 3], 18, 7.348469229],
[4, 5, [3, 2, 3, 3, 3, 3], 17, 7.]];

For instance, the interpretation of [3, 5, [3, 3, 3, 3, 3, 3], 18, 7.348469229] is that if S be [e₁,e₂,e₃,e_4,e₅,e₃+e₅ mod 2] then we get [3, 3, 3, 3, 3, 3] where 18 is the summation of [3, 3, 3, 3, 3, 3] and also the number 7.348469229 is obtained from the following command: 

MatrixNorm(convert([3, 3, 3, 3, 3, 3], Matrix), Frobenius);

Now we choose E_k for k in [1..nops(E)] such that E_k[4] be minimum over all E_k[4]'s. For example we choose [2, 4, [3, 2, 4, 2, 2, 3], 16, 6.782329983] since 16 is minimum between E_k[4]'s.

There are two points: First one is that if we have E_i and E_j such that E_i[4]=E_j[4] then we choose E_i if  E_i[5]>E_j[5] . The second point is that if E_i[4]=E_j[4] and also E_i[5]=E_j[5] then we choose one of them such as the first one. Finally we update the set S from the first two entries of E_k that we have obtained. For instance, the updated S in our example is:

S=[e[1],e[2],e[3],e[4],e[5],S[E[6][1]]+S[E[6][2]] mod 2]=[e[1],e[2],e[3],e[4],e[5],e[2]+e[4] mod 2]

Now we repeat this procedure for the updated S until that in one of the entries of E such as E_k we get E_k[3]=[1,1,..,1].

I have written a procedure in Maple for the mentioned question. But my procedure takes long time to compute when I run it over a list such as A with parameters n=m=64.

I want to kindly request you please modify the following code or suggest another fast procedure for this question.

 
 restart;

 with(LinearAlgebra):
 with(ListTools): 
 with(combinat):
 
 BP := proc (A::list)
    local n, m, r, S, U, tt, P, E, t, Q, Z, R, k, T, j, PP, i, QQ; 
    n := nops(A); 
    m := nops(A[1]);
    r := [seq(0, i = 1 .. m)]; 
    r[1] := 1; 
    S := [seq(Rotate(r, m-i+1), i = 1 .. m)]; 
    unassign('r'); 
    U := []; 
    tt := 1;
 while 0 < tt do
     P := choose(nops(S), 2);
     E := [];
     for t to nops(P) do
          Q := P[t];
          Z := []; 
          R := [S[], `mod`(S[Q[1]]+S[Q[2]], 2)];
          for k to n do
               T := A[k];
               for j to nops(R) do 
                    PP := choose(nops(R), j);
                    for i to nops(PP) do 
                         QQ := PP[i]; r := `mod`(add(R[QQ[i]], i = 1 .. nops(QQ)), 2);
                         if Occurrences(0, r-T) = m then
                            Z := [op(Z), j];
                            i := nops(PP)+1; 
                            j := nops(R)+1 
                         end if; 
                      end do; 
                      unassign('i, QQ, r, PP'): 
                  end do;
              end do;
              E := [op(E), [Q[], Z, add(Z[i], i = 1 .. nops(Z)), evalf(MatrixNorm(convert(Z, Matrix), Frobenius))]]; 
              unassign('k, Z, Q, R'): 
          end do;
          r := FindMinimalElement([seq(E[i][4], i = 1 .. nops(E))]);
          R := []; 
         for i to nops(E) do 
             if E[i][4] = r then R := [op(R), E[i]] end if 
         end do;
         T := [FindMaximalElement([seq(R[i][5], i = 1 .. nops(R))], position)];
         S := [S[], `mod`(S[R[T[2]][1]]+S[R[T[2]][2]], 2)]; 
         U := [op(U), [R[T[2]][1], R[T[2]][2]]];
         if Occurrences(1, R[T[2]][3]) = n then tt := 0 end if;
         unassign('r, i, R, T, E')
     end do;
     return U;
 end proc:
 
 A := [[1, 1, 1, 0, 0], [0, 1, 0, 1, 1], [1, 0, 1, 1, 1], [0, 1, 1, 1, 0], [1, 1, 0, 1, 0], [0, 1, 1, 1, 1]];

 BP(A);
      [[2, 4], [3, 6], [5, 7], [1, 2], [1, 6], [3, 8], [3, 9], [8, 9]]
 

For more information please see Section 2.1 of the following paper: https://eprint.iacr.org/2019/856.pdf

Thanks in advance for your consideration of this request.