#include <NTL/LLL.h>#include <NTL/fileio.h>#include <NTL/vec_double.h>#include <NTL/new.h>NTL_START_IMPLstatic 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 double max_abs(double *v, long n){   long i;   double res, t;   res = 0;   for (i = 1; i <= n; i++) {      t = fabs(v[i]);      if (t > res) res = t;   }   return res;}static void RowTransformStart(double *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] < TR_BND && a[i] > -TR_BND);      inf = inf & in_a[i];   }   in_float = inf;}static void RowTransformFinish(vec_ZZ& A, double *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]);      }      else {         conv(a[i], A(i));      }   }}static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1,                          double mu, double *a, double *b, long *in_a,                         double& max_a, double max_b, long& in_float)// x = x - y*MU{   static ZZ T, MU;   long k;   long n = A.length();   long i;   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] -= b[i];         return;      }      if (mu == -1) {         for (i = 1; i <= n; i++)            a[i] += b[i];         return;      }      if (mu == 0) return;      for (i = 1; i <= n; i++)         a[i] -= mu*b[i];      return;   }   MU = MU1;   if (MU == 1) {      for (i = 1; i <= n; i++) {         if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND &&             b[i] < TR_BND && b[i] > -TR_BND) {            a[i] -= b[i];         }         else {            if (in_a[i]) {               conv(A(i), a[i]);               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] < TR_BND && a[i] > -TR_BND &&             b[i] < TR_BND && b[i] > -TR_BND) {            a[i] += b[i];         }         else {            if (in_a[i]) {               conv(A(i), a[i]);               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]);               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] < TR_BND && a[i] > -TR_BND &&                b[i] < b_bnd && b[i] > -b_bnd) {                  a[i] -= b[i]*mu;            }            else {               if (in_a[i]) {                  conv(A(i), a[i]);                  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]);            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);      }   }}class GivensCache_FP {public:   GivensCache_FP(long m, long n);   ~GivensCache_FP();   void flush();   void selective_flush(long l);   void swap(long l);   void swap();   void touch();   void incr();   long sz;   double **buf;   long *bl;   long *bv;   long bp;};GivensCache_FP::GivensCache_FP(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 = NTL_NEW_OP doubleptr[sz];    if (!buf) Error("out of memory");   for (i = 0; i < sz; i++)      if (!(buf[i] = NTL_NEW_OP double[n+1])) Error("out of memory");   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_FP::~GivensCache_FP(){   long i;   for (i = 0; i < sz; i++) delete [] buf[i];   delete [] buf;   delete [] bl;   delete [] bv;}void GivensCache_FP::flush(){   long i;   for (i = 0; i < sz; i++) bl[i] = 0;}void GivensCache_FP::selective_flush(long l){   long i;   for (i = 0; i < sz; i++)      if (bl[i] && bv[i] >= l)         bl[i] = 0;}void GivensCache_FP::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_FP::swap(){   swap(bl[bp] - 1);}void GivensCache_FP::touch(){   long k = bl[bp];   bl[bp] = 0;   selective_flush(k);}void GivensCache_FP::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;}staticvoid GivensComputeGS(double **B1, double **mu, double **aux, long k, long n,                     GivensCache_FP& cache){   long i, j;   double c, s, a, b, t;   double *p = mu[k];   double *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++) {         double *cptr = mu[i];         double *sptr = aux[i];            for (j = n; j > i; j--) {            c = cptr[j];            s = sptr[j];               a = c*pp[j-1] - s*pp[j];            b = s*pp[j-1] + c*pp[j];               pp[j-1] = a;            pp[j] = b;         }            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++) {      double *cptr = mu[i];      double *sptr = aux[i];        for (j = n; j > i; j--) {         c = cptr[j];         s = sptr[j];           a = c*p[j-1] - s*p[j];         b = s*p[j-1] + c*p[j];           p[j-1] = a;         p[j] = b;      }        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 if (fabs(b) > fabs(a)) {         t = -a/b;         s = 1/sqrt(1 + t*t);         c = s*t;      }      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_FP: internal error");   if (k > n) p[k] = 0;