#include <NTL/image.h>
#include <image.c>

NTL_CLIENT

int main()
{
   mat_ZZ B, L;
   long s, s1;
   long r;

   long i=0;
   long j=0;
   long n=0;
   long m=0;
#if 0
   cin >> B;
   m = B.NumRows();
   n = B.NumCols();
#else
   cout << "\n n is the number of rows.\n "; 
   cout << "m is the number of columns.\n "; 
   cout << "bit_in is the maximum number of bits in entry of input.\n ";
   cout << "t_im is the time taken by image to compute the kernel.\n ";
   cout << "t_HNF is the time taken by HNF_image to compute the kernel.\n ";
   cout << "t_HNF/im is their ratio.\n ";
   cout << "bit_im is the maximum number of bits in kernel from image.\n ";
   cout << "bit_HNF is the maximum number of bits in kernel from HNF_image.\n ";
   cout << "bit_HNF/im  is their ratio.\n ";
   cout << "\n";
   
   cerr << "\nn: ";
   cin >> n;

while ( n ){
   cerr << "m: ";
   cin >> m;

   long k;
   cerr << "k: ";
   cin >> k;

   B.SetDims(n, m);
   for (i = 1; i <= n; i++)
      for (j = 1; j <= m; j++) {
         RandomLen(B(i,j), k);
         if (RandomBnd(2)) negate(B(i,j), B(i,j));
      }

   L.SetDims(m, m);
   for (i = 1; i <= m; i++)
      for (j = 1; j <= m; j++) {
         RandomLen(L(i,j), k);
         if (RandomBnd(2)) negate(L(i,j), L(i,j));
      }

#endif

   mat_ZZ U, U1, B0, B1, B2;

   B0 = B;

   double t, t0, t1, t2, ra;
   ZZ d, d1, max1, max2, max3, co;

//   cerr << "B0 = " << B0 << "\n";   

B = B0;



for (r = 1; r <= 2 ; r++){
   B = B0;
   t = GetTime();
   s = image(d,B,U);
   U.SetDims(n - s, n);
   t1 = (GetTime()-t);
//   cerr <<  "time = " << t1 << "\n";

   max1=0;
      for (i = 1; i <= n - s; i++)
      for (j = 1; j <= n; j++) {
         abs(co,U(i,j));
         if (co>max1) max1=co;
      }


   B1=B;
   d1=d;
   U1=U;
   
   B = B0;
//   cerr << "HNF_image...\n";
   t = GetTime();
   s1 = HNF_image(B, U);
   U.SetDims(n - s1, n);
   t2 = GetTime()-t;
   if (s1 != s) Error("Other answer for rank");
//   cout << "B = " << B << "\n";
   mul(B1, U, B0);
   B.SetDims(0,0);
   B.SetDims(n - s1,m);
   if (B1 != B) Error("bad kernel");

   max2=0;
      for (i = 1; i <= n - s1; i++)
      for (j = 1; j <= n; j++) {
         abs(co,U(i,j));
         if (co>max2) max2=co;
      }
   max3=0;
      for (i = 1; i <= n; i++)
      for (j = 1; j <= m; j++) {
         abs(co,B0(i,j));
         if (co>max3) max3=co;
      }
   ra=NumBits(max2);
   ra=ra/NumBits(max1);
   cout << "m = " << m << ", n = " << n << ", bit_in = " << NumBits(max3);
   cout <<", t_im = "<<t1<< ", t_HNF = "<< t2 ;
   cout << "\n";
   cout << "     t_HNF/im = "<< t2/t1; 
   cout <<"  bit_im = " << NumBits(max1) ;
   cout <<"  bit_HNF = " << NumBits(max2) ;
   cout <<"  bit_HNF/im = " << ra ;
   cout << "\n";
   
   mul(B0,B0,L);
   
   } ;  
   cerr << "\nn: ";
   cin >> n;
   }
   cout <<"bye!\n";
}