T0 = []
T1 = []
T2 = []
X = []

n = ZZ(10000)

for i in range(1,100,1):
  X.append(i)
  T0.append(0.0)
  T1.append(0.0)
  T2.append(0.0)
  magma._start()
  ctr = 3
  for j in range(ctr):
    A = random_matrix(GF(2),n,n,density=i/n)

    t0 = cputime()
    E0 = A.copy().echelonize(algorithm='pluq')
    t0 = cputime(t0)
    T0[-1] += t0

    t1 = cputime()
    E1 = A.copy().echelonize(algorithm='m4ri')
    t1 = cputime(t1)
    T1[-1] += t1

    if E0 != E1:
      print "FUCK"
      break

#     AM = magma(A)
#     t2 = magma.cputime()
#     _ = AM.EchelonForm()
#     t2 = magma.cputime(t2)
#     T2[-1] += t2

  T0[-1] = T0[-1]/float(ctr)
  T1[-1] = T1[-1]/float(ctr)
#   T2[-1] = T2[-1]/float(ctr)
  
  print i,T0[-1], T1[-1], T2[-1]

X0 = [i/2.0 for i in X]

import pylab
pylab.clf() # clear the figure first
pylab.figure(1)

# plot some data and add a legend
pylab.plot(X0,T0,label="PLUQ") 
pylab.plot(X0,T1,label="M4RI") 
#pylab.plot(X0,T2,label="Magma") 

pylab.legend() # print the legend

pylab.title("Sensitivity to Density (%s x %s)"%(str(n),str(n)))
pylab.ylabel("execution time $t$") # label the axes
pylab.xlabel("non-zero entries per row on average")

pylab.savefig('pluq-m4ri-magma-%s.png'%str(n),dpi=int(150)) # fire!