title: Optimital algorithm for solving XOR equation systems
body: <p>I am trying to solve the XOR equation system. For example:</p>
  
  <pre><code>A = [[1, 1, 1, 0, 0], [0, 1, 1, 1, 0], [0, 0, 1, 1, 1], [0, 1, 1, 0, 1], [0, 1, 0, 1, 1]]
  s = [3, 14, 13, 5, 2]
  m = 5 # len(s)
  Ax = s =&gt; x = [12, 9, 6, 1, 10]
  </code></pre>
  
  <p>I tried 2 ways: </p>
  
  <ul>
  <li>The first way is Gaussian elimination (~2.5 second) which was showed <a href="https://stackoverflow.com/questions/19350769/solving-systems-of-xor-equations">here</a></li>
  <li>The second way to invert modular matrix A (with modulo 2) and then, XOR multiply with A_invert and s. (~7.5 second)</li>
  </ul>
  
  <p>Could you please show me is there a way or a python library to speed up. Even I tried to use gmpy2 library, but it cannot reduce much. Below I described python code so that you can easily follow.</p>
  
  <p>Using Gaussian elimination:</p>
  
  <pre><code>def SolveLinearSystem (A, B, N):
      for K in range (0, N):
          if (A[K][K] == 0):
              for i in range (K+1, N):
                  if (A[i][K]!=0):
                      for L in range (0, N):
                          s = A[K][L]
                          A[K][L] = A[i][L]
                          A[i][L] = s
                      s = B[i]
                      B[i] = B[K]
                      B[K] = s
                      break
          for I in range (0, N):
              if (I!=K):
                  if (A[I][K]):
                      #M = 0
                      for M in range (K, N):
                          A[I][M] = A[I][M] ^ A[K][M]
                      B[I] = B[I] ^ B[K]
  
  SolveLinearSystem (A, s, 5)
  </code></pre>
  
  <p>Using Inversion</p>
  
  <pre><code>def identitymatrix(n):
      return [[long(x == y) for x in range(0, n)] for y in range(0, n)]
  
  def multiply_vector_scalar (vector, scalar, q):
      kq = []
      for i in range (0, len(vector)):
          kq.append (vector[i] * scalar %q)
      return kq
  
  def minus_vector_scalar(vector1, scalar, vector2, q):
      kq = []
      for i in range (0, len(vector1)):
          kq.append ((vector1[i] - scalar * vector2[i]) %q)
      return kq
  
  def inversematrix(matrix, q):
      n = len(matrix)
      A =[]
      for j in range (0, n):
          temp = []
          for i in range (0, n):
              temp.append (matrix[j][i])
          A.append(temp)
  
      Ainv = identitymatrix(n)
  
      for i in range(0, n):
          factor = gmpy2.invert(A[i][i], q) #invert mod q
          A[i] = multiply_vector_scalar(A[i],factor,q)
          Ainv[i] = multiply_vector_scalar(Ainv[i],factor,q)
          for j in range(0, n):
              if (i != j):
                  factor = A[j][i]
                  A[j] = minus_vector_scalar(A[j], factor, A[i], q)
                  Ainv[j] = minus_vector_scalar(Ainv[j], factor, Ainv[i], q)
      return Ainv
  
  def solve_equation (A, y):
      result = []
      for i in range (0, m):
          temp = 0
          for j in range (0, m):
              temp = (temp ^ A[i][j]* y[j])
          result.append(temp)
      return result
  
  A_invert = inversematrix(A, 2)
  print solve_equation (A_invert, s)
  </code></pre>