Problem:
 Score : 1500 points    Problem Statement  Snuke received two matrices A and B as birthday presents.
Each of the matrices is an N by N matrix that consists of only 0 and 1 .  Then he computed the product of the two matrices, C = AB .
Since he performed all computations in modulo two, C was also an N by N matrix that consists of only 0 and 1 .
For each 1 ≤ i, j ≤ N , you are given c_{i, j} , the (i, j) -element of the matrix C .  However, Snuke accidentally ate the two matrices A and B , and now he only knows C .
Compute the number of possible (ordered) pairs of the two matrices A and B , modulo 10^9+7 .      Constraints  1 ≤ N ≤ 300  c_{i, j} is either 0 or 1 .         Input The input is given from Standard Input in the following format:  N  c_{1, 1}  ...  c_{1, N} : c_{N, 1}  ...  c_{N, N}       Output Print the number of possible (ordered) pairs of two matrices A and B (modulo 10^9+7 ).        Sample Input 1 2
0 1
1 0      Sample Output 1 6       Sample Input 2 10
1 0 0 1 1 1 0 0 1 0
0 0 0 1 1 0 0 0 1 0
0 0 1 1 1 1 1 1 1 1
0 1 0 1 0 0 0 1 1 0
0 0 1 0 1 1 1 1 1 1
1 0 0 0 0 1 0 0 0 0
1 1 1 0 1 0 0 0 0 1
0 0 0 1 0 0 1 0 1 0
0 0 0 1 1 1 0 0 0 0
1 0 1 0 0 1 1 1 1 1      Sample Output 2 741992411

Solution (Python):
import sys
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 7)

import numpy as np

MOD = 10**9 + 7

N = int(readline())
C = np.array([line.split() for line in readlines()],np.int8)

def rank(A):
    if (A==0).all():
        return 0
    i = np.nonzero(A[:,0])[0]
    if len(i) == 0:
        return rank(A[:,1:])
    i = i[0]
    temp = A[i].copy()
    A[i] = A[0]
    A[0] = temp
    A[1:] ^= A[1:,0][:,None] * A[0][None,:]
    return 1 + rank(A[1:,1:])

def cumprod(arr,MOD):
    L = len(arr); Lsq = int(L**.5+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    return arr.ravel()[:L]

def power_mod(a,n,MOD):
    if n == 0:
        return 1
    x = power_mod(a,n//2,MOD)
    x *= x
    x %= MOD
    return (a*x)%MOD if n&1 else x

r = rank(C)

x = np.full(N*N+100,2,np.int64); x[0] = 1
pow2 = cumprod(x,MOD)

G = np.zeros((N+1,N+1),np.int64)
for n in range(N+1):
    G[n,1:] = pow2[n] - pow2[:N]
G[:,0] = 1
for n in range(1,N+1):
    G[:,n] *= G[:,n-1]; G[:,n] %= MOD

Gd = np.diagonal(G)

D = G * power_mod(Gd,MOD-2,MOD)[None,:] % MOD
F = np.zeros((N+1,N+1),np.int64)
for n in range(N+1):
    F[n,:n+1] = Gd[:n+1] * D[n,n::-1] % MOD

B = D[N] * F[N] % MOD

C = D[N,r] * F[:,r] % MOD * pow2[N*N::-N] % MOD

A = (B[r:N+1] * C[r:N+1] % MOD).sum() % MOD

answer = A * pow(int(B[r]),MOD-2,MOD) % MOD
print(answer)

Solution (Python):
from numpy import *
P=10**9+7
N=int(input())
C=array([input().split() for _ in range(N)],int8)
def r(A):
  if(A==0).all():return 0
  i=A[:,0].nonzero()[0]
  if len(i) == 0:return r(A[:,1:])
  t=A[i[0]].copy();A[i[0]]=A[0];A[0]=t
  A[1:]^= A[1:,0][:,None]*A[0][None,:]
  return 1+r(A[1:,1:])
r=r(C)
p=ones(N+1,int64)
for n in range(1,N+1):p[n]=p[n-1]*2%P
d=zeros((N+1,N+1,N+1),int64);d[:,0,0]=1
for M in range(1,N+1):
  d[:,M,:M]+=d[:,M-1,:M]*p[:M]%P
  d[:,M,1:M+1]+=d[:,M-1,0:M]*(p[:,None]-p[None,0:M])%P
  d[:,M,:]%=P
print(sum(d[N,N,n]*d[N,n,r]%P*pow(2,N*(N-n),P)%P for n in range(r,N+1))%P*pow(int(d[N,N,r]),P-2,P)%P)

Solution (Python):
import sys
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 7)

import numpy as np

MOD = 10**9 + 7

N = int(readline())
C = np.array([line.split() for line in readlines()],np.int8)

def rank(A):
    if (A==0).all():
        return 0
    i = np.nonzero(A[:,0])[0]
    if len(i) == 0:
        return rank(A[:,1:])
    i = i[0]
    temp = A[i].copy()
    A[i] = A[0]
    A[0] = temp
    A[1:] ^= A[1:,0][:,None] * A[0][None,:]
    return 1 + rank(A[1:,1:])

r = rank(C)

pow2 = np.ones(301,dtype=np.int64)
for n in range(1,301):
    pow2[n] = pow2[n-1]*2%MOD


dp = np.zeros((301,301,301),dtype=np.int64)
dp[:,0,0] = 1
for M in range(1,301):
    dp[:,M,:M] += dp[:,M-1,:M] * pow2[:M] % MOD
    dp[:,M,1:M+1] += dp[:,M-1,0:M] * (pow2[:,None]-pow2[None,0:M]) % MOD
    dp[:,M,:] %= MOD


x = 0
for n in range(r,N+1):
    x += dp[N,N,n] * dp[N,n,r] % MOD * pow(2,N*(N-n),MOD) % MOD
x %= MOD

answer = x * pow(int(dp[N,N,r]),MOD-2,MOD) % MOD
print(answer)