#include <NTL/LLL.h>
#include <NTL/fileio.h>

#include <NTL/new.h>

NTL_START_IMPL




static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
// x = x - y*MU
{
   static ZZ T, MU;
   long k;

   long n = A.length();
   long i;

   MU = MU1;

   if (MU == 1) {
      for (i = 1; i <= n; i++)
         sub(A(i), A(i), B(i));

      return;
   }

   if (MU == -1) {
      for (i = 1; i <= n; i++)
         add(A(i), A(i), B(i));

      return;
   }

   if (MU == 0) return;

   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
      k = MakeOdd(MU);
   else
      k = 0;


   if (MU.WideSinglePrecision()) {
      long mu1;
      conv(mu1, MU);

      for (i = 1; i <= n; i++) {
         mul(T, B(i), mu1);
         if (k > 0) LeftShift(T, T, k);
         sub(A(i), A(i), T);
      }
   }
   else {
      for (i = 1; i <= n; i++) {
         mul(T, B(i), MU);
         if (k > 0) LeftShift(T, T, k);
         sub(A(i), A(i), T);
      }
   }
}

static void RowTransform2(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
// x = x + y*MU
{
   static ZZ T, MU;
   long k;

   long n = A.length();
   long i;

   MU = MU1;

   if (MU == 1) {
      for (i = 1; i <= n; i++)
         add(A(i), A(i), B(i));

      return;
   }

   if (MU == -1) {
      for (i = 1; i <= n; i++)
         sub(A(i), A(i), B(i));

      return;
   }

   if (MU == 0) return;

   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
      k = MakeOdd(MU);
   else
      k = 0;

   if (MU.WideSinglePrecision()) {
      long mu1;
      conv(mu1, MU);

      for (i = 1; i <= n; i++) {
         mul(T, B(i), mu1);
         if (k > 0) LeftShift(T, T, k);
         add(A(i), A(i), T);
      }
   }
   else {
      for (i = 1; i <= n; i++) {
         mul(T, B(i), MU);
         if (k > 0) LeftShift(T, T, k);
         add(A(i), A(i), T);
      }
   }
}

class GivensCache_RR {
public:
   GivensCache_RR(long m, long n);
   ~GivensCache_RR();

   void flush();
   void selective_flush(long l);
   void swap(long l);
   void swap();
   void touch();
   void incr();

   long sz;

   mat_RR buf;

   long *bl;
   long *bv;
   long bp;
};


GivensCache_RR::GivensCache_RR(long m, long n)
{
   sz = min(m, n)/10;
   if (sz < 2) 
      sz = 2;
   else if (sz > 20)
      sz = 20;

   typedef double *doubleptr;

   long i;

   buf.SetDims(sz, n);

   bl = NTL_NEW_OP long[sz];
   if (!bl) Error("out of memory");
   for (i = 0; i < sz; i++) bl[0] = 0;

   bv = NTL_NEW_OP long[sz];
   if (!bv) Error("out of memory");
   for (i = 0; i < sz; i++) bv[0] = 0;

   bp = 0;
}

GivensCache_RR::~GivensCache_RR()
{
   delete [] bl;
   delete [] bv;
}

void GivensCache_RR::flush()
{
   long i;
   for (i = 0; i < sz; i++) bl[i] = 0;
}

void GivensCache_RR::selective_flush(long l)
{
   long i;

   for (i = 0; i < sz; i++)
      if (bl[i] && bv[i] >= l)
         bl[i] = 0;
}

void GivensCache_RR::swap(long l)
{
   long k = bl[bp];
   long i;

   i = 0;
   while (i < sz && bl[i] != l)
      i++;

   if (i < sz) {
      bl[bp] = l;
      bl[i] = k;
   }
   else
      bl[bp] = l;

   selective_flush(l);
}

void GivensCache_RR::swap()
{
   swap(bl[bp] - 1);
}

void GivensCache_RR::touch()
{
   long k = bl[bp];
   bl[bp] = 0;
   selective_flush(k);
}

void GivensCache_RR::incr()
{
   long k = bl[bp];
   long k1 = k+1;
   long i;

   i = 0;
   while (i < sz && bl[i] != k1)
      i++;

   if (i < sz) {
      bp = i;
      return;
   }

   i = 0; 
   while (i < sz && bl[i] != 0)
      i++;

   if (i < sz) {
      bp = i;
      return;
   }

   long max_val = 0;
   long max_index = 0;
   for (i = 0; i < sz; i++) {
      long t = labs(bl[i]-k1);
      if (t > max_val) {
         max_val = t;
         max_index = i;
      }
   }

   bp = max_index;
   bl[max_index] = 0;
}


static
void GivensComputeGS(mat_RR& B1, mat_RR& mu, mat_RR& aux, long k, long n,
                     GivensCache_RR& cache)
{
   long i, j;

   RR c, s, a, b, t;
   RR T1, T2;

   vec_RR& p = mu(k);

   vec_RR& pp = cache.buf[cache.bp];

   if (!cache.bl[cache.bp]) {
      for (j = 1; j <= n; j++)
         pp(j) = B1(k,j);

      long backoff;
      backoff = k/4;
      if (backoff < 2)
         backoff = 2;
      else if (backoff > cache.sz + 2)
         backoff = cache.sz + 2; 

      long ub = k-(backoff-1);

      for (i = 1; i < ub; i++) {
         vec_RR& cptr = mu(i);
         vec_RR& sptr = aux(i);
   
         for (j = n; j > i; j--) {
            c = cptr(j);
            s = sptr(j);
   
            // a = c*pp(j-1) - s*pp(j);
            mul(T1, c, pp(j-1));
            mul(T2, s, pp(j));
            sub(a, T1, T2);

            // b = s*pp(j-1) + c*pp(j);
            mul(T1, s, pp(j-1));
            mul(T2, c, pp(j));
            add(b, T1, T2);
   
            pp(j-1) = a;
            pp(j) = b;
         }
   
         div(pp(i), pp(i), mu(i,i));
      }

      cache.bl[cache.bp] = k;
      cache.bv[cache.bp] = k-backoff;
   }

   for (j = 1; j <= n; j++)
      p(j) = pp(j);

   for (i = max(cache.bv[cache.bp]+1, 1); i < k; i++) {
      vec_RR& cptr = mu(i);
      vec_RR& sptr = aux(i);
  
      for (j = n; j > i; j--) {
         c = cptr(j);
         s = sptr(j);
  
         // a = c*p(j-1) - s*p(j);
         mul(T1, c, p(j-1));
         mul(T2, s, p(j));
         sub(a, T1, T2);

         // b = s*p(j-1) + c*p(j);
         mul(T1, s, p(j-1));
         mul(T2, c, p(j));
         add(b, T1, T2);
  
         p(j-1) = a;
         p(j) = b;
      }
  
      div(p(i), p(i), mu(i,i));
   }

   for (j = n; j > k; j--) {
      a = p(j-1);
      b = p(j);

      if (b == 0) {
         c = 1;
         s = 0;
      }
      else {
         abs(T1, b);
         abs(T2, a);

         if (T1 > T2) {
            // t = -a/b;
            div(T1, a, b);
            negate(t, T1);
   
            // s = 1/sqrt(1 + t*t);
            sqr(T1, t);
            add(T1, T1, 1);
            SqrRoot(T1, T1);
            inv(s, T1);
            
            // c = s*t;
            mul(c, s, t);
         }
         else {
            // t = -b/a;
            div(T1, b, a);
            negate(t, T1);
   
            // c = 1/sqrt(1 + t*t);
            sqr(T1, t);
            add(T1, T1, 1);
            SqrRoot(T1, T1);
            inv(c, T1);
   
            // s = c*t;
            mul(s, c, t);
         }
      }
   
      // p(j-1) = c*a - s*b;
      mul(T1, c, a);
      mul(T2, s, b);
      sub(p(j-1), T1, T2);

      p(j) = c;
      aux(k,j) = s;
   }

   if (k > n+1) Error("G_LLL_RR: internal error");
   if (k > n) p(k) = 0;

}

static RR red_fudge;
static long log_red = 0;

static void init_red_fudge()
{
   log_red = long(0.50*RR::precision());

   power2(red_fudge, -log_red);
}

static void inc_red_fudge()
{

   mul(red_fudge, red_fudge, 2);
   log_red--;

   cerr << "G_LLL_RR: warning--relaxing reduction (" << log_red << ")\n";

   if (log_red < 4)
      Error("G_LLL_RR: can not continue...sorry");
}




static long verbose = 0;

static unsigned long NumSwaps = 0;
static double StartTime = 0;
static double LastTime = 0;



static void G_LLLStatus(long max_k, double t, long m, const mat_ZZ& B)
{
   cerr << "---- G_LLL_RR status ----\n";
   cerr << "elapsed time: ";
   PrintTime(cerr, t-StartTime);
   cerr << ", stage: " << max_k;
   cerr << ", rank: " << m;
   cerr << ", swaps: " << NumSwaps << "\n";

   ZZ t1;
   long i;
   double prodlen = 0;

   for (i = 1; i <= m; i++) {
      InnerProduct(t1, B(i), B(i));
      if (!IsZero(t1))
         prodlen += log(t1);
   }

   cerr << "log of prod of lengths: " << prodlen/(2.0*log(2.0)) << "\n";

   if (LLLDumpFile) {
      cerr << "dumping to " << LLLDumpFile << "...";

      ofstream f;
      OpenWrite(f, LLLDumpFile);
      
      f << "[";
      for (i = 1; i <= m; i++) {
         f << B(i) << "\n";
      }
      f << "]\n";

      f.close();

      cerr << "\n";
   }

   LastTime = t;
   
}



static
long ll_G_LLL_RR(mat_ZZ& B, mat_ZZ* U, const RR& delta, long deep, 
           LLLCheckFct check, mat_RR& B1, mat_RR& mu, 
           mat_RR& aux, long m, long init_k, long &quit,
           GivensCache_RR& cache)
{
   long n = B.NumCols();

   long i, j, k, Fc1;
   ZZ MU;
   RR mu1, t1, t2, cc;
   ZZ T1;


   quit = 0;
   k = init_k;

   long counter;

   long trigger_index;
   long small_trigger;
   long cnt;

   RR half;
   conv(half,  0.5);
   RR half_plus_fudge;
   add(half_plus_fudge, half, red_fudge);

   long max_k = 0;
   double tt;

   cache.flush();

   while (k <= m) {

      if (k > max_k) {
         max_k = k;
      }

      if (verbose) {
         tt = GetTime();

         if (tt > LastTime + LLLStatusInterval)
            G_LLLStatus(max_k, tt, m, B);
      }

      GivensComputeGS(B1, mu, aux, k, n, cache);

      counter = 0;
      trigger_index = k;
      small_trigger = 0;
      cnt = 0;

      do {
         // size reduction

         counter++;
         if (counter > 10000) {
            cerr << "G_LLL_XD: warning--possible infinite loop\n";
            counter = 0;
         }


         Fc1 = 0;

         for (j = k-1; j >= 1; j--) {
            abs(t1, mu(k,j));
            if (t1 > half_plus_fudge) {

               if (!Fc1) {
                  if (j > trigger_index ||
                      (j == trigger_index && small_trigger)) {

                     cnt++;

                     if (cnt > 10) {
                        inc_red_fudge();
                        add(half_plus_fudge, half, red_fudge);
                        cnt = 0;
                     }
                  }

                  trigger_index = j;
                  small_trigger = (t1 < 4);
               }

               Fc1 = 1;
   
               mu1 = mu(k,j);
               if (sign(mu1) >= 0) {
                  sub(mu1, mu1, half);
                  ceil(mu1, mu1);
               }
               else {
                  add(mu1, mu1, half);
                  floor(mu1, mu1);
               }

               if (mu1 == 1) {
                  for (i = 1; i <= j-1; i++)
                     sub(mu(k,i), mu(k,i), mu(j,i));
               }
               else if (mu1 == -1) {
                  for (i = 1; i <= j-1; i++)
                     add(mu(k,i), mu(k,i), mu(j,i));
               }
               else {
                  for (i = 1; i <= j-1; i++) {
                     mul(t2, mu1, mu(j,i));
                     sub(mu(k,i), mu(k,i), t2);
                  }
               }

   
               conv(MU, mu1);

               sub(mu(k,j), mu(k,j), mu1);
   
               RowTransform(B(k), B(j), MU);
               if (U) RowTransform((*U)(k), (*U)(j), MU);
            }
         }

         if (Fc1) {
            for (i = 1; i <= n; i++)
               conv(B1(k, i), B(k, i));
            cache.touch();
            GivensComputeGS(B1, mu, aux, k, n, cache);
         }
      } while (Fc1);

      if (check && (*check)(B(k))) 
         quit = 1;

      if (IsZero(B(k))) {
         for (i = k; i < m; i++) {
            // swap i, i+1
            swap(B(i), B(i+1));
            swap(B1(i), B1(i+1));
            if (U) swap((*U)(i), (*U)(i+1));
         }

         cache.flush();

         m--;
         if (quit) break;
         continue;
      }

      if (quit) break;

      if (deep > 0) {
         // deep insertions
   
         Error("sorry...deep insertions not implemented");

      } // end deep insertions

      // test G_LLL reduction condition

      if (k <= 1) {
         cache.incr();
         k++;
      }
      else {
         sqr(t1, mu(k,k-1));
         sub(t1, delta, t1);
         sqr(t2, mu(k-1,k-1));
         mul(t1, t1, t2);
         sqr(t2, mu(k, k));
         if (t1 > t2) {
            // swap rows k, k-1
            swap(B(k), B(k-1));
            swap(B1(k), B1(k-1));
            if (U) swap((*U)(k), (*U)(k-1));

            cache.swap();
   
            k--;
            NumSwaps++;
         }
         else {
            cache.incr();
            k++;
         }
      }
   }

   if (verbose) {
      G_LLLStatus(m+1, GetTime(), m, B);
   }


   return m;
}

static
long G_LLL_RR(mat_ZZ& B, mat_ZZ* U, const RR& delta, long deep, 
           LLLCheckFct check)
{
   long m = B.NumRows();
   long n = B.NumCols();

   long i, j;
   long new_m, dep, quit;
   RR s;
   ZZ MU;
   RR mu1;

   RR t1;
   ZZ T1;

   init_red_fudge();