g_lll_qp.cpp

数值算法库for Windows

      else {
         t = -b/a;
         c = 1/sqrt(1 + t*t);
         s = c*t;
      }
   
      p[j-1] = c*a - s*b;
      p[j] = c;
      aux[k][j] = s;
   }

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

   for (i = 1; i <= k; i++)
      if (!IsFinite(&p[i])) Error("G_LLL_QP...numbers too big");
}

static quad_float red_fudge = to_quad_float(0);
static long log_red = 0;

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_QP 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 void init_red_fudge()
{
   long i;

   // initial log_red should be <= NTL_DOUBLE_PRECISION-2,
   // to help ensure stability in G_BKZ_QP1

   log_red = NTL_DOUBLE_PRECISION-2; 

   red_fudge = 1;

   for (i = log_red; i > 0; i--)
      red_fudge = red_fudge*0.5;
}

static void inc_red_fudge()
{

   red_fudge = red_fudge * 2;
   log_red--;

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

   if (log_red < 4)
      Error("G_LLL_QP: too much loss of precision...stop!");
}


static
long ll_G_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
           LLLCheckFct check, quad_float **B1, quad_float **mu, 
           quad_float **aux,
           long m, long init_k, long &quit, GivensCache_QP& cache)
{
   long n = B.NumCols();

   long i, j, k, Fc1;
   ZZ MU;
   quad_float mu1;

   quad_float t1;
   ZZ T1;
   quad_float *tp;

   quad_float half = to_quad_float(0.5);
   quad_float half_plus_fudge = 0.5 + red_fudge;

   quit = 0;
   k = init_k;

   vec_long in_vec_mem;
   in_vec_mem.SetLength(n+1);
   long *in_vec = in_vec_mem.elts();

   quad_float *max_b;
   max_b = NTL_NEW_OP quad_float [m+1];
   if (!max_b) Error("out of memory in lll_G_LLL_QP");

   for (i = 1; i <= m; i++)
      max_b[i] = max_abs(B1[i], n);

   long in_float;


   long counter;

   long trigger_index;
   long small_trigger;
   long cnt;

   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_QP: warning--possible infinite loop\n";
            counter = 0;
         }


         Fc1 = 0;
   
         for (j = k-1; j >= 1; j--) {
            t1 = fabs(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();
                        half_plus_fudge = 0.5 + red_fudge;
                        cnt = 0;
                     }
                  }

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

                  Fc1 = 1;
                  RowTransformStart(B1[k], in_vec, in_float, n);
               }


   
               mu1 = mu[k][j];
               if (mu1 >= 0)
                  mu1 = ceil(mu1-half);
               else
                  mu1 = floor(mu1+half);
   
   
               quad_float *mu_k = mu[k];
               quad_float *mu_j = mu[j];
  
               if (mu1 == 1) {
                  for (i = 1; i <= j-1; i++)
                     mu_k[i] -= mu_j[i];
               }
               else if (mu1 == -1) {
                  for (i = 1; i <= j-1; i++)
                     mu_k[i] += mu_j[i];
               }
               else {
                  for (i = 1; i <= j-1; i++)
                     mu_k[i] -= mu1*mu_j[i];
               }

               // cout << j << " " << mu[k][j] << " " << mu1 << "\n";
  
               mu_k[j] -= mu1;

               conv(MU, mu1);

   
               RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec,
                            max_b[k], max_b[j], in_float);

               if (U) RowTransform((*U)(k), (*U)(j), MU);
            }
         }

         if (Fc1) {
            RowTransformFinish(B(k), B1[k], in_vec);
            max_b[k] = max_abs(B1[k], n);
            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));
            tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp;
            t1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = t1;
            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 LLL reduction condition

      if (k > 1 &&
         sqrt(delta - mu[k][k-1]*mu[k][k-1])*fabs(mu[k-1][k-1]) >
         fabs(mu[k][k])) {

         // swap rows k, k-1
         swap(B(k), B(k-1));
         tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp;
         t1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = t1;
         if (U) swap((*U)(k), (*U)(k-1));

         cache.swap();

         k--;
         NumSwaps++;
         // cout << "- " << k << "\n";
      }
      else {
         cache.incr();
         k++;
         // cout << "+ " << k << "\n";
      }
   }

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


   delete [] max_b;

   return m;
}

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

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

   quad_float t1;
   ZZ T1;

   init_red_fudge();

   if (U) ident(*U, m);

   quad_float **B1;  // approximates B

   typedef quad_float *quad_floatptr;

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

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

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

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

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

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

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

   for (i = 1; i <= m; i++)
      for (j = 1; j <= n; j++)
         if (!IsFinite(&B1[i][j]))
             Error("G_LLL_QP: numbers too big...use G_LLL_XD");

   GivensCache_QP cache(m, n);

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



   dep = m - new_m;
   m = new_m;

   if (dep > 0) {
      // for consistency, we move all of the zero rows to the front

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


   // clean-up

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

   delete [] B1;

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

   delete [] mu;

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

   delete [] aux;

   return m;
}

         

long G_LLL_QP(mat_ZZ& B, double delta, long deep, LLLCheckFct check,
            long verb)
{
   verbose = verb;
   NumSwaps = 0;
   if (verbose) {
      StartTime = GetTime();
      LastTime = StartTime;
   }

   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
   if (deep < 0) Error("G_LLL_QP: bad deep");
   return G_LLL_QP(B, 0, to_quad_float(delta), deep, check);
}

long G_LLL_QP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, 
           LLLCheckFct check, long verb)
{
   verbose = verb;
   NumSwaps = 0;
   if (verbose) {
      StartTime = GetTime();
      LastTime = StartTime;
   }


   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
   if (deep < 0) Error("G_LLL_QP: bad deep");
   return G_LLL_QP(B, &U, to_quad_float(delta), deep, check);
}



static vec_quad_float G_BKZConstant;

static
void ComputeG_BKZConstant(long beta, long p)
{
   const quad_float c_PI = 
      to_quad_float("3.141592653589793238462643383279502884197");
   const quad_float LogPI = 
      to_quad_float("1.144729885849400174143427351353058711647");

   G_BKZConstant.SetLength(beta-1);

   vec_quad_float Log;
   Log.SetLength(beta);


   long i, j, k;
   quad_float x, y;

   for (j = 1; j <= beta; j++)
      Log(j) = log(to_quad_float(j));

   for (i = 1; i <= beta-1; i++) {
      // First, we compute x = gamma(i/2)^{2/i}

      k = i/2;

      if ((i & 1) == 0) { // i even
         x = 0;
         for (j = 1; j <= k; j++)
            x = x + Log(j);
          
         x = x * (1/to_quad_float(k));

         x = exp(x);
      }
      else { // i odd
         x = 0;
         for (j = k + 2; j <= 2*k + 2; j++)
            x = x + Log(j);

         x = 0.5*LogPI + x - 2*(k+1)*Log(2);

         x = x * (2.0/to_quad_float(i));

         x = exp(x);
      }

      // Second, we compute y = 2^{2*p/i}

      y = -(2*p/to_quad_float(i))*Log(2);
      y = exp(y);

      G_BKZConstant(i) = x*y/c_PI;
   }
}

static vec_quad_float G_BKZThresh;

static 
void ComputeG_BKZThresh(quad_float *c, long beta)
{
   G_BKZThresh.SetLength(beta-1);

   long i;
   quad_float x;

   x = 0;

   for (i = 1; i <= beta-1; i++) {
      x += log(c[i-1]);
      G_BKZThresh(i) = exp(x/to_quad_float(i))*G_BKZConstant(i);
      if (!IsFinite(&G_BKZThresh(i))) G_BKZThresh(i) = 0;
   }
}


static 
void G_BKZStatus(double tt, double enum_time, long NumIterations, 
               long NumTrivial, long NumNonTrivial, long NumNoOps, long m, 
               const mat_ZZ& B)
{
   cerr << "---- G_BKZ_QP status ----\n";
   cerr << "elapsed time: ";
   PrintTime(cerr, tt-StartTime);
   cerr << ", enum time: ";