#include <NTL/mat_ZZ_p.h>
#include <flint/fmpz.h>
#include <flint/fmpz_mat.h>
#include <cstdio>

NTL_CLIENT


struct FlintZZ {
   fmpz_t value;

   explicit
   FlintZZ(const ZZ& a)
   {
      long n = NumBits(a);
   
      fmpz_init(value);
   
      for (long i = 0; i < n; i++)
         if (bit(a, i)) fmpz_setbit(value, i);
   
      if (a < 0)
         fmpz_neg(value, value);

   }

   ~FlintZZ() { fmpz_clear(value); }

};



struct Flint_mat_ZZ {

   fmpz_mat_t value;

   explicit
   Flint_mat_ZZ(const mat_ZZ_p& a)
   {
      long n = a.NumRows();
      long m = a.NumCols();
      fmpz_mat_init(value, n, m);

      for (long i = 0; i < n; i++) {
         for (long j = 0; j < m; j++) {
            FlintZZ aij(rep(a[i][j]));
            fmpz_set(fmpz_mat_entry (value, i, j), aij.value);
         }
      }
   }

   explicit
   Flint_mat_ZZ(const ZZ& p, long n)
   {
      fmpz_mat_init(value, n, n);
      FlintZZ pp(p*p*n);

      for (long i = 0; i < n; i++) {
         for (long j = 0; j < n; j++) {
            fmpz_set(fmpz_mat_entry (value, i, j), pp.value);
         }
      }
   }

   ~Flint_mat_ZZ() { fmpz_mat_clear(value); }
};


void reduce(Flint_mat_ZZ& X, const Flint_mat_ZZ& A, const FlintZZ& p, long n)
{
   for (long i = 0; i < n; i++)
      for (long j = 0; j < n; j++)
         fmpz_mod(fmpz_mat_entry(X.value, i, j), 
                  fmpz_mat_entry(A.value, i, j), p.value);
}



void FillRandom(Mat<ZZ_p>& A)
{
   long n = A.NumRows();
   long m = A.NumCols();
   for (long i = 0; i < n; i++)
      for (long j = 0; j < m; j++)
         random(A[i][j]);
}


int main()
{
   setbuf(stdout, NULL);

   // similar format as flint-factor, one less row, two less columns
   for (long l = 256; l <= 2048; l *= 2) {
      for (long idx = 0; idx < 7; idx ++) {
         long n  = 256*(1L << idx/2);
         if (idx & 1) n += n/2;

	 SetSeed((ZZ(l) << 64) + ZZ(n));

	 ZZ p;

	 RandomLen(p, l);
	 if (p % 2 == 0) p++;
	 ZZ_p::init(p);

	 mat_ZZ_p A, B, C;

	 A.SetDims(n, n);
	 FillRandom(A);

	 B.SetDims(n, n);
	 FillRandom(B);

	 C.SetDims(n, n);



	 double t;

	 long iter = 1;
	 do {
	    t = GetTime();
	    for (long i = 0; i < iter; i++) mul(C, A, B);
	    t = GetTime() - t;
	    iter *= 2;
	 } while (t < 3);
	 iter /= 2;

	 double NTLTime = t;

	 Flint_mat_ZZ fA(A), fB(B), fC(C), fCC(p, n);
         FlintZZ fp(p);


	 t = GetTime();
	 for (long i = 0; i < iter; i++) {
            fmpz_mat_mul(fCC.value, fA.value, fB.value);
            reduce(fC, fCC, fp, n);
         }
	 t = GetTime()-t;
	 double FlintTime = t;

	 printf("%8.2f", FlintTime/NTLTime);
      }

      printf("\n");
   }
}


