lll_qp.cpp

数值算法库for Windows

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

#include <NTL/new.h>

NTL_START_IMPL


static quad_float InnerProduct(quad_float *a, quad_float *b, long n)
{
   quad_float s;
   long i;

   s = 0;
   for (i = 1; i <= n; i++) 
      s += a[i]*b[i];

   return s;
}


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);
      }
   }
}

#define TR_BND (NTL_FDOUBLE_PRECISION/2.0)
// Just to be safe!!

static quad_float max_abs(quad_float *v, long n)
{
   long i;
   quad_float res, t;

   res = 0;

   for (i = 1; i <= n; i++) {
      t = fabs(v[i]);
      if (t > res) res = t;
   }

   return res;
}


static void RowTransformStart(quad_float *a, long *in_a, long& in_float, long n)
{
   long i;
   long inf = 1;

   for (i = 1; i <= n; i++) {
      in_a[i] = (a[i].hi < TR_BND && a[i].hi > -TR_BND);
      inf = inf & in_a[i];
   }

   in_float = inf;
}


static void RowTransformFinish(vec_ZZ& A, quad_float *a, long *in_a)
{
   long n = A.length();
   long i;

   for (i = 1; i <= n; i++) {
      if (in_a[i])  {
         conv(A(i), a[i].hi);
      }
      else {
         conv(a[i], A(i));
      }
   }
}



   



static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1, 
                         quad_float *a, quad_float *b, long *in_a, 
                         quad_float& max_a, quad_float max_b, long& in_float)
// x = x - y*MU
{
   static ZZ T, MU;
   long k;
   double mu;


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

   conv(mu, MU1);

   if (in_float) {
      max_a += fabs(mu)*max_b; 
      if (max_a >= TR_BND) {
         in_float = 0;
      }
   }

   if (in_float) {
      if (mu == 1) {
         for (i = 1; i <= n; i++)
            a[i].hi -= b[i].hi;

         return;
      }

      if (mu == -1) {
         for (i = 1; i <= n; i++)
            a[i].hi += b[i].hi;

         return;
      }

      if (mu == 0) return;

      for (i = 1; i <= n; i++)
         a[i].hi -= mu*b[i].hi;


      return;
   }

   MU = MU1;

   if (MU == 1) {
      for (i = 1; i <= n; i++) {
         if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
             b[i].hi < TR_BND && b[i].hi > -TR_BND) {

            a[i].hi -= b[i].hi;
         }
         else {
            if (in_a[i]) {
               conv(A(i), a[i].hi);
               in_a[i] = 0;
            }

            sub(A(i), A(i), B(i));
         }
      }

      return;
   }

   if (MU == -1) {
      for (i = 1; i <= n; i++) {
         if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
             b[i].hi < TR_BND && b[i].hi > -TR_BND) {

            a[i].hi += b[i].hi;
         }
         else {
            if (in_a[i]) {
               conv(A(i), a[i].hi);
               in_a[i] = 0;
            }

            add(A(i), A(i), B(i));
         }
      }

      return;
   }

   if (MU == 0) return;

   double b_bnd = fabs(TR_BND/mu) - 1;
   if (b_bnd < 0) b_bnd = 0;


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


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

      if (k > 0) {
         for (i = 1; i <= n; i++) {
            if (in_a[i]) {
               conv(A(i), a[i].hi);
               in_a[i] = 0;
            }

            mul(T, B(i), mu1);
            LeftShift(T, T, k);
            sub(A(i), A(i), T);
         }
      }
      else {
         for (i = 1; i <= n; i++) {
            if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
                b[i].hi < b_bnd && b[i].hi > -b_bnd) {
  
               a[i].hi -= b[i].hi*mu;
            }
            else {
               if (in_a[i]) {
                  conv(A(i), a[i].hi);
                  in_a[i] = 0;
               }
               mul(T, B(i), mu1);
               sub(A(i), A(i), T);
           }
         }
      }
   }
   else {
      for (i = 1; i <= n; i++) {
         if (in_a[i]) {
            conv(A(i), a[i].hi);
            in_a[i] = 0;
         }
         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);
      }
   }
}

static
void ComputeGS(mat_ZZ& B, quad_float **B1, quad_float **mu, quad_float *b, 
               quad_float *c, long k, double bound, long st, quad_float *buf)
{
   long n = B.NumCols();
   long i, j;
   quad_float s, t1, y, t;
   ZZ T1;
   long test;

   quad_float *mu_k = mu[k];

   if (st < k) {
      for (i = 1; i < st; i++)
         buf[i] = mu_k[i]*c[i];
   }

   for (j = st; j <= k-1; j++) {
      if (b[k].hi/NTL_FDOUBLE_PRECISION < NTL_FDOUBLE_PRECISION/b[j].hi) {
         // we can compute inner product exactly in double precision

         double z = 0;
         quad_float *B1_k = B1[k];
         quad_float *B1_j = B1[j];

         for (i = 1; i <= n; i++) 
            z += B1_k[i].hi * B1_j[i].hi;

         s = z;
      }
      else {
         s = InnerProduct(B1[k], B1[j], n);
   
         y = fabs(s);
         if (y.hi == 0)
            test = (b[k].hi != 0);
         else {
            double t = y.hi/b[j].hi;
            double t1 = b[k].hi/y.hi;
            if (t <= 1)
               test = (t*bound <= t1);
            else if (t1 >= 1)
               test = (t <= t1/bound);
            else
               test = 0;
         }

         if (test) {
            InnerProduct(T1, B(k), B(j));
            conv(s, T1);
         }
      }


      quad_float *mu_j = mu[j];

      t1 = 0;
      for (i = 1; i <= j-1; i++)
         t1 += mu_j[i]*buf[i];

      mu_k[j] = (buf[j] = (s - t1))/c[j];
   }

   s = 0;
   for (j = 1; j <= k-1; j++)
      s += mu_k[j]*buf[j];

   c[k] = b[k] - s;
}

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 LLLStatus(long max_k, double t, long m, const mat_ZZ& B)
{
   cerr << "---- 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 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 << "LLL_QP: warning--relaxing reduction (" << log_red << ")\n";

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


static
long ll_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
           LLLCheckFct check, quad_float **B1, quad_float **mu, 
           quad_float *b, quad_float *c,
           long m, long init_k, long &quit)
{
   long n = B.NumCols();

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

   quad_float t1;
   ZZ T1;
   quad_float *tp;


   static double bound = 0;


   if (bound == 0) {
      // we tolerate a 15% loss of precision in computing
      // inner products in ComputeGS.

      bound = 1;
      for (i = 2*long(0.15*2*NTL_DOUBLE_PRECISION); i > 0; i--) {
         bound = bound * 2;
      }
   }


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