long G_BKZ_QP(mat_ZZ& BB, double delta, 
         long beta, long prune, LLLCheckFct check, long verb)
{
   verbose = verb;
   NumSwaps = 0;
   if (verbose) {
      StartTime = GetTime();
      LastTime = StartTime;
   }



   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
   if (beta < 2) Error("G_BKZ_QP: bad block size");

   return G_BKZ_QP(BB, 0, to_quad_float(delta), beta, prune, check);
}



static
long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ* UU, quad_float delta, 
         long beta, long prune, LLLCheckFct check)
{

   long m = BB.NumRows();
   long n = BB.NumCols();
   long m_orig = m;
   
   long i, j;
   ZZ MU;

   ZZ T1;
   quad_float *tp;

   init_red_fudge();

   mat_ZZ B;
   B = BB;

   B.SetDims(m+1, n);


   quad_float **B1;  // approximates B

   typedef quad_float *quad_floatptr;

   B1 = NTL_NEW_OP quad_floatptr[m+2];
   if (!B1) Error("G_BKZ_QP: out of memory");

   for (i = 1; i <= m+1; i++) {
      B1[i] = NTL_NEW_OP quad_float[n+1];
      if (!B1[i]) Error("G_BKZ_QP: out of memory");
   }

   quad_float **mu;
   mu = NTL_NEW_OP quad_floatptr[m+2];
   if (!mu) Error("G_BKZ_QP: out of memory");

   for (i = 1; i <= m+1; i++) {
      mu[i] = NTL_NEW_OP quad_float[n+2];
      if (!mu[i]) Error("G_BKZ_QP: out of memory");
   }

   quad_float **aux;
   aux = NTL_NEW_OP quad_floatptr[m+2];
   if (!aux) Error("G_BKZ_QP: out of memory");

   for (i = 1; i <= m+1; i++) {
      aux[i] = NTL_NEW_OP quad_float[n+1];
      if (!aux[i]) Error("G_BKZ_QP: out of memory");
   }

   quad_float *c; // squared lengths of Gramm-Schmidt basis vectors

   c = NTL_NEW_OP quad_float[m+2];
   if (!c) Error("G_BKZ_QP: out of memory");

   double cbar;

   double *ctilda;
   ctilda = NTL_NEW_OP double[m+2];
   if (!ctilda) Error("G_BKZ_QP: out of memory");

   double *vvec;
   vvec = NTL_NEW_OP double[m+2];
   if (!vvec) Error("G_BKZ_QP: out of memory");

   double *yvec;
   yvec = NTL_NEW_OP double[m+2];
   if (!yvec) Error("G_BKZ_QP: out of memory");

   double *uvec;
   uvec = NTL_NEW_OP double[m+2];
   if (!uvec) Error("G_BKZ_QP: out of memory");

   double *utildavec;
   utildavec = NTL_NEW_OP double[m+2];
   if (!utildavec) Error("G_BKZ_QP: out of memory");


   long *Deltavec;
   Deltavec = NTL_NEW_OP long[m+2];
   if (!Deltavec) Error("G_BKZ_QP: out of memory");

   long *deltavec;
   deltavec = NTL_NEW_OP long[m+2];
   if (!deltavec) Error("G_BKZ_QP: out of memory");

   mat_ZZ Ulocal;
   mat_ZZ *U;

   if (UU) {
      Ulocal.SetDims(m+1, m);
      for (i = 1; i <= m; i++)
         conv(Ulocal(i, i), 1);
      U = &Ulocal;
   }
   else
      U = 0;

   long quit;
   long new_m;
   long z, jj, kk;
   long s, t;
   long h;

   double eta;

   for (i = 1; i <=m; i++)
      for (j = 1; j <= n; j++) {
         conv(B1[i][j], B(i, j));
         CheckFinite(&B1[i][j]);
      }


   GivensCache_QP cache(m, n);

   m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);



   double tt;

   double enum_time = 0;
   unsigned long NumIterations = 0;
   unsigned long NumTrivial = 0;
   unsigned long NumNonTrivial = 0;
   unsigned long NumNoOps = 0;

   long verb = verbose;

   verbose = 0;

   long clean = 1;

   if (m < m_orig) {
      for (i = m_orig+1; i >= m+2; i--) {
         // swap i, i-1

         swap(B(i), B(i-1));
         if (U) swap((*U)(i), (*U)(i-1));
      }
   }

   if (!quit && m > 1) {
      // cerr << "continuing\n";
      if (beta > m) beta = m;

      if (prune > 0)
         ComputeG_BKZConstant(beta, prune);

      z = 0;
      jj = 0;
   
      while (z < m-1) {
         jj++;
         kk = min(jj+beta-1, m);
   
         if (jj == m) {
            jj = 1;
            kk = beta;
            clean = 1;
         }

         if (verb) {
            tt = GetTime();
            if (tt > LastTime + LLLStatusInterval)
               G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
                         NumNonTrivial, NumNoOps, m, B);
         }



   
         // ENUM

         double tt1;

         if (verb) {
            tt1 = GetTime();
         }

         for (i = jj; i <= kk; i++) {
            c[i] = mu[i][i]*mu[i][i];
            CheckFinite(&c[i]);
         }

         if (prune > 0)
            ComputeG_BKZThresh(&c[jj], kk-jj+1);

   
         cbar = to_double(c[jj]);
         utildavec[jj] = uvec[jj] = 1;
   
         yvec[jj] = vvec[jj] = 0;
         Deltavec[jj] = 0;
   
   
         s = t = jj;
         deltavec[jj] = 1;
   
         for (i = jj+1; i <= kk+1; i++) {
            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
            Deltavec[i] = 0;
            vvec[i] = 0;
            deltavec[i] = 1;
         }

         long enum_cnt = 0;
   
         while (t <= kk) {
            if (verb) {
               enum_cnt++;
               if (enum_cnt > 100000) {
                  enum_cnt = 0;
                  tt = GetTime();
                  if (tt > LastTime + LLLStatusInterval) {
                     enum_time += tt - tt1;
                     tt1 = tt;
                     G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
                               NumNonTrivial, NumNoOps, m, B);
                  }
               }
            }

            ctilda[t] = ctilda[t+1] + 
               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*to_double(c[t]);

            ForceToMem(&ctilda[t]); // prevents an infinite loop
   
            if (prune > 0 && t > jj) {
               eta = to_double(G_BKZThresh(t-jj));
            }
            else
               eta = 0;

            if (ctilda[t] < cbar - eta) {
               if (t > jj) {
                  double t1;

                  t--;
                  t1 = 0;
                  for (i = t+1; i <= s; i++) {
                     t1 += utildavec[i]*to_double(mu[i][t]);
                  }


                  yvec[t] = t1;
                  t1 = -t1;
                  if (t1 >= 0)
                     t1 = ceil(t1-0.5);
                  else
                     t1 = floor(t1+0.5);

                  utildavec[t] = vvec[t] = t1;
                  Deltavec[t] = 0;
                  if (utildavec[t] > -yvec[t]) 
                     deltavec[t] = -1;
                  else
                     deltavec[t] = 1;
               }
               else {
                  cbar = ctilda[jj];
                  for (i = jj; i <= kk; i++) {
                     uvec[i] = utildavec[i];
                  }
               }
            }
            else {
               t++;
               s = max(s, t);
               if (t < s) Deltavec[t] = -Deltavec[t];
               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
               utildavec[t] = vvec[t] + Deltavec[t];
            }
         }

         if (verb) {
            tt1 = GetTime() - tt1;
            enum_time += tt1;
         }

         NumIterations++;

         h = min(kk+1, m);

         quad_float t1;
   
         if ((delta-8*red_fudge)*c[jj] > cbar*(1+64/NTL_FDOUBLE_PRECISION)) {

            clean = 0;

            // we treat the case that the new vector is b_s (jj < s <= kk)
            // as a special case that appears to occur most of the time.
   
            s = 0;
            for (i = jj+1; i <= kk; i++) {
               if (uvec[i] != 0) {
                  if (s == 0)
                     s = i;
                  else
                     s = -1;
               }
            }
   
            if (s == 0) Error("G_BKZ_QP: internal error");
   
            if (s > 0) {
               // special case

               NumTrivial++;
   
               for (i = s; i > jj; i--) {
                  // swap i, i-1
                  swap(B(i-1), B(i));
                  if (U) swap((*U)(i-1), (*U)(i));
                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
               }
   
               // cerr << "special case\n";
               new_m = ll_G_LLL_QP(B, U, delta, 0, check,
                                B1, mu, aux, h, jj, quit, cache);
               if (new_m != h) Error("G_BKZ_QP: internal error");
               if (quit) break;
            }
            else {
               // the general case

               NumNonTrivial++;

   
               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);

               if (U) {
                  for (i = 1; i <= m_orig; i++)
                     conv((*U)(m+1, i), 0);
               }

               for (i = jj; i <= kk; i++) {
                  if (uvec[i] == 0) continue;
                  conv(MU, uvec[i]);
                  RowTransform2(B(m+1), B(i), MU);
                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
               }
      
               for (i = m+1; i >= jj+1; i--) {
                  // swap i, i-1
                  swap(B(i-1), B(i));
                  if (U) swap((*U)(i-1), (*U)(i));
                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
               }
      
               for (i = 1; i <= n; i++) {
                  conv(B1[jj][i], B(jj, i));
                  CheckFinite(&B1[jj][i]);
               }
      
               if (IsZero(B(jj))) Error("G_BKZ_QP: internal error"); 
      
               // remove linear dependencies
   
               // cerr << "general case\n";
               new_m = ll_G_LLL_QP(B, U, delta, 0, 0, B1, mu, aux,
                                  kk+1, jj, quit, cache);
              
               if (new_m != kk) Error("G_BKZ_QP: internal error"); 

               // remove zero vector
      
               for (i = kk+2; i <= m+1; i++) {
                  // swap i, i-1
                  swap(B(i-1), B(i));
                  if (U) swap((*U)(i-1), (*U)(i));
                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
               }
      
               quit = 0;
               if (check) {
                  for (i = 1; i <= kk; i++)
                     if ((*check)(B(i))) {
                        quit = 1;
                        break;
                     }
               }

               if (quit) break;
   
               if (h > kk) {
                  // extend reduced basis
   
                  new_m = ll_G_LLL_QP(B, U, delta, 0, check,
                                   B1, mu, aux, h, h, quit, cache);
   
                  if (new_m != h) Error("G_BKZ_QP: internal error");
                  if (quit) break;
               }
            }
   
            z = 0;
         }
         else {
            // G_LLL_QP
            // cerr << "progress\n";

            NumNoOps++;


            if (!clean) {
               new_m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux,
                                   h, h, quit, cache);

               if (new_m != h) Error("G_BKZ_QP: internal error");
               if (quit) break;
            }
   
            z++;
         }
      }
   }


   if (verb) {
      G_BKZStatus(GetTime(), enum_time, NumIterations, NumTrivial, NumNonTrivial,
                NumNoOps, m, B);
   }


   // clean up


   if (m_orig > m) {
      // for consistency, we move zero vectors to the front

      for (i = m+1; i <= m_orig; i++) {
         swap(B(i), B(i+1));
         if (U) swap((*U)(i), (*U)(i+1));
      }

      for (i = 0; i < m; i++) {
         swap(B(m_orig-i), B(m-i));
         if (U) swap((*U)(m_orig-i), (*U)(m-i));
      }
   }

   B.SetDims(m_orig, n);
   BB = B;

   if (U) {
      U->SetDims(m_orig, m_orig);
      *UU = *U;
   }

   for (i = 1; i <= m_orig+1; i++) {
      delete [] B1[i];
   }

   delete [] B1;

   for (i = 1; i <= m_orig+1; i++) {
      delete [] mu[i];
   }

   delete [] mu;

   for (i = 1; i <= m_orig+1; i++) {
      delete [] aux[i];
   }

   delete [] aux;


   delete [] c;
   delete [] ctilda;
   delete [] vvec;
   delete [] yvec;
   delete [] uvec;
   delete [] utildavec;
   delete [] Deltavec;
   delete [] deltavec;

   return m;
}

long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ& UU, double delta, 
         long beta, long prune, LLLCheckFct check, long verb)
{
   verbose = verb;
   NumSwaps = 0;
   if (verbose) {
      StartTime = GetTime();
      LastTime = StartTime;
   }


   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
   if (beta < 2) Error("G_BKZ_QP: bad block size");

   return G_BKZ_QP1(BB, &UU, to_quad_float(delta), beta, prune, check);
}

long G_BKZ_QP1(mat_ZZ& BB, double delta, 
         long beta, long prune, LLLCheckFct check, long verb)
{
   verbose = verb;
   NumSwaps = 0;
   if (verbose) {
      StartTime = GetTime();
      LastTime = StartTime;
   }



   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
   if (beta < 2) Error("G_BKZ_QP: bad block size");

   return G_BKZ_QP1(BB, 0, to_quad_float(delta), beta, prune, check);
}

NTL_END_IMPL