source: git/ntl/src/G_LLL_QP.c @ 33a041

spielwiese
Last change on this file since 33a041 was 311902, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: 5.3.2-C++-fix git-svn-id: file:///usr/local/Singular/svn/trunk@7474 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 43.5 KB
Line 
1
2#include <NTL/LLL.h>
3#include <NTL/vec_quad_float.h>
4#include <NTL/fileio.h>
5
6
7#include <NTL/new.h>
8
9NTL_START_IMPL
10
11
12static inline
13void CheckFinite(double *p)
14{
15   if (!IsFinite(p)) Error("G_LLL_QP: numbers too big...use G_LLL_XD");
16}
17
18
19static inline
20void CheckFinite(quad_float *p)
21{
22   if (!IsFinite(p)) Error("G_LLL_QP: numbers too big...use G_LLL_XD");
23}
24
25
26
27static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
28// x = x - y*MU
29{
30   static ZZ T, MU;
31   long k;
32
33   long n = A.length();
34   long i;
35
36   MU = MU1;
37
38   if (MU == 1) {
39      for (i = 1; i <= n; i++)
40         sub(A(i), A(i), B(i));
41
42      return;
43   }
44
45   if (MU == -1) {
46      for (i = 1; i <= n; i++)
47         add(A(i), A(i), B(i));
48
49      return;
50   }
51
52   if (MU == 0) return;
53
54   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
55      k = MakeOdd(MU);
56   else
57      k = 0;
58
59
60   if (MU.WideSinglePrecision()) {
61      long mu1;
62      conv(mu1, MU);
63
64      for (i = 1; i <= n; i++) {
65         mul(T, B(i), mu1);
66         if (k > 0) LeftShift(T, T, k);
67         sub(A(i), A(i), T);
68      }
69   }
70   else {
71      for (i = 1; i <= n; i++) {
72         mul(T, B(i), MU);
73         if (k > 0) LeftShift(T, T, k);
74         sub(A(i), A(i), T);
75      }
76   }
77}
78
79#define TR_BND (NTL_FDOUBLE_PRECISION/2.0)
80// Just to be safe!!
81
82static double max_abs(quad_float *v, long n)
83{
84   long i;
85   double res, t;
86
87   res = 0;
88
89   for (i = 1; i <= n; i++) {
90      t = fabs(v[i].hi);
91      if (t > res) res = t;
92   }
93
94   return res;
95}
96
97
98static void RowTransformStart(quad_float *a, long *in_a, long& in_float, long n)
99{
100   long i;
101   long inf = 1;
102
103   for (i = 1; i <= n; i++) {
104      in_a[i] = (a[i].hi < TR_BND && a[i].hi > -TR_BND);
105      inf = inf & in_a[i];
106   }
107
108   in_float = inf;
109}
110
111
112static void RowTransformFinish(vec_ZZ& A, quad_float *a, long *in_a)
113{
114   long n = A.length();
115   long i;
116
117   for (i = 1; i <= n; i++) {
118      if (in_a[i])  {
119         conv(A(i), a[i].hi);
120      }
121      else {
122         conv(a[i], A(i));
123         CheckFinite(&a[i]);
124      }
125   }
126}
127
128
129
130   
131
132
133
134static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1, 
135                         quad_float *a, quad_float *b, long *in_a, 
136                         double& max_a, double max_b, long& in_float)
137// x = x - y*MU
138{
139   static ZZ T, MU;
140   long k;
141   double mu;
142
143
144   long n = A.length();
145   long i;
146
147   conv(mu, MU1);
148   CheckFinite(&mu);
149
150   if (in_float) {
151      double mu_abs = fabs(mu);
152      if (mu_abs > 0 && max_b > 0 && (mu_abs >= TR_BND || max_b >= TR_BND)) {
153         in_float = 0;
154      }
155      else {
156         max_a += mu_abs*max_b;
157         if (max_a >= TR_BND)
158            in_float = 0;
159      }
160   }
161
162   if (in_float) {
163      if (mu == 1) {
164         for (i = 1; i <= n; i++)
165            a[i].hi -= b[i].hi;
166
167         return;
168      }
169
170      if (mu == -1) {
171         for (i = 1; i <= n; i++)
172            a[i].hi += b[i].hi;
173
174         return;
175      }
176
177      if (mu == 0) return;
178
179      for (i = 1; i <= n; i++)
180         a[i].hi -= mu*b[i].hi;
181
182
183      return;
184   }
185
186   MU = MU1;
187
188   if (MU == 1) {
189      for (i = 1; i <= n; i++) {
190         if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
191             b[i].hi < TR_BND && b[i].hi > -TR_BND) {
192
193            a[i].hi -= b[i].hi;
194         }
195         else {
196            if (in_a[i]) {
197               conv(A(i), a[i].hi);
198               in_a[i] = 0;
199            }
200
201            sub(A(i), A(i), B(i));
202         }
203      }
204
205      return;
206   }
207
208   if (MU == -1) {
209      for (i = 1; i <= n; i++) {
210         if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
211             b[i].hi < TR_BND && b[i].hi > -TR_BND) {
212
213            a[i].hi += b[i].hi;
214         }
215         else {
216            if (in_a[i]) {
217               conv(A(i), a[i].hi);
218               in_a[i] = 0;
219            }
220
221            add(A(i), A(i), B(i));
222         }
223      }
224
225      return;
226   }
227
228   if (MU == 0) return;
229
230   double b_bnd = fabs(TR_BND/mu) - 1;
231   if (b_bnd < 0) b_bnd = 0;
232
233
234   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
235      k = MakeOdd(MU);
236   else
237      k = 0;
238
239
240   if (MU.WideSinglePrecision()) {
241      long mu1;
242      conv(mu1, MU);
243
244      if (k > 0) {
245         for (i = 1; i <= n; i++) {
246            if (in_a[i]) {
247               conv(A(i), a[i].hi);
248               in_a[i] = 0;
249            }
250
251            mul(T, B(i), mu1);
252            LeftShift(T, T, k);
253            sub(A(i), A(i), T);
254         }
255      }
256      else {
257         for (i = 1; i <= n; i++) {
258            if (in_a[i] && a[i].hi < TR_BND && a[i].hi > -TR_BND &&
259                b[i].hi < b_bnd && b[i].hi > -b_bnd) {
260 
261               a[i].hi -= b[i].hi*mu;
262            }
263            else {
264               if (in_a[i]) {
265                  conv(A(i), a[i].hi);
266                  in_a[i] = 0;
267               }
268               mul(T, B(i), mu1);
269               sub(A(i), A(i), T);
270           }
271         }
272      }
273   }
274   else {
275      for (i = 1; i <= n; i++) {
276         if (in_a[i]) {
277            conv(A(i), a[i].hi);
278            in_a[i] = 0;
279         }
280         mul(T, B(i), MU);
281         if (k > 0) LeftShift(T, T, k);
282         sub(A(i), A(i), T);
283      }
284   }
285}
286
287static void RowTransform2(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
288// x = x + y*MU
289{
290   static ZZ T, MU;
291   long k;
292
293   long n = A.length();
294   long i;
295
296   MU = MU1;
297
298   if (MU == 1) {
299      for (i = 1; i <= n; i++)
300         add(A(i), A(i), B(i));
301
302      return;
303   }
304
305   if (MU == -1) {
306      for (i = 1; i <= n; i++)
307         sub(A(i), A(i), B(i));
308
309      return;
310   }
311
312   if (MU == 0) return;
313
314   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
315      k = MakeOdd(MU);
316   else
317      k = 0;
318
319   if (MU.WideSinglePrecision()) {
320      long mu1;
321      conv(mu1, MU);
322
323      for (i = 1; i <= n; i++) {
324         mul(T, B(i), mu1);
325         if (k > 0) LeftShift(T, T, k);
326         add(A(i), A(i), T);
327      }
328   }
329   else {
330      for (i = 1; i <= n; i++) {
331         mul(T, B(i), MU);
332         if (k > 0) LeftShift(T, T, k);
333         add(A(i), A(i), T);
334      }
335   }
336}
337class GivensCache_QP {
338public:
339   GivensCache_QP(long m, long n);
340   ~GivensCache_QP();
341
342   void flush();
343   void selective_flush(long l);
344   void swap(long l);
345   void swap();
346   void touch();
347   void incr();
348
349   long sz;
350
351   quad_float **buf;
352   long *bl;
353   long *bv;
354   long bp;
355};
356
357
358GivensCache_QP::GivensCache_QP(long m, long n)
359{
360   sz = min(m, n)/10;
361   if (sz < 2) 
362      sz = 2;
363   else if (sz > 20)
364      sz = 20;
365
366   typedef quad_float *quad_floatptr;
367
368   long i;
369   buf = NTL_NEW_OP quad_floatptr[sz]; 
370   if (!buf) Error("out of memory");
371   for (i = 0; i < sz; i++)
372      if (!(buf[i] = NTL_NEW_OP quad_float[n+1])) Error("out of memory");
373
374   bl = NTL_NEW_OP long[sz];
375   if (!bl) Error("out of memory");
376   for (i = 0; i < sz; i++) bl[0] = 0;
377
378   bv = NTL_NEW_OP long[sz];
379   if (!bv) Error("out of memory");
380   for (i = 0; i < sz; i++) bv[0] = 0;
381
382   bp = 0;
383}
384
385GivensCache_QP::~GivensCache_QP()
386{
387   long i;
388
389   for (i = 0; i < sz; i++) delete [] buf[i];
390   delete [] buf;
391   delete [] bl;
392   delete [] bv;
393}
394
395void GivensCache_QP::flush()
396{
397   long i;
398   for (i = 0; i < sz; i++) bl[i] = 0;
399}
400
401void GivensCache_QP::selective_flush(long l)
402{
403   long i;
404
405   for (i = 0; i < sz; i++)
406      if (bl[i] && bv[i] >= l)
407         bl[i] = 0;
408}
409
410void GivensCache_QP::swap(long l)
411{
412   long k = bl[bp];
413   long i;
414
415   i = 0;
416   while (i < sz && bl[i] != l)
417      i++;
418
419   if (i < sz) {
420      bl[bp] = l;
421      bl[i] = k;
422   }
423   else
424      bl[bp] = l;
425
426   selective_flush(l);
427}
428
429void GivensCache_QP::swap()
430{
431   swap(bl[bp] - 1);
432}
433
434void GivensCache_QP::touch()
435{
436   long k = bl[bp];
437   bl[bp] = 0;
438   selective_flush(k);
439}
440
441void GivensCache_QP::incr()
442{
443   long k = bl[bp];
444   long k1 = k+1;
445   long i;
446
447   i = 0;
448   while (i < sz && bl[i] != k1)
449      i++;
450
451   if (i < sz) {
452      bp = i;
453      return;
454   }
455
456   i = 0; 
457   while (i < sz && bl[i] != 0)
458      i++;
459
460   if (i < sz) {
461      bp = i;
462      return;
463   }
464
465   long max_val = 0;
466   long max_index = 0;
467   for (i = 0; i < sz; i++) {
468      long t = labs(bl[i]-k1);
469      if (t > max_val) {
470         max_val = t;
471         max_index = i;
472      }
473   }
474
475   bp = max_index;
476   bl[max_index] = 0;
477}
478
479
480static
481void GivensComputeGS(quad_float **B1, quad_float **mu, quad_float **aux, long k, long n,
482                     GivensCache_QP& cache)
483{
484   long i, j;
485
486   quad_float c, s, a, b, t;
487
488   quad_float *p = mu[k];
489
490   quad_float *pp = cache.buf[cache.bp];
491
492   if (!cache.bl[cache.bp]) {
493      for (j = 1; j <= n; j++)
494         pp[j] = B1[k][j];
495
496      long backoff;
497      backoff = k/4;
498      if (backoff < 2)
499         backoff = 2;
500      else if (backoff > cache.sz + 2)
501         backoff = cache.sz + 2; 
502
503      long ub = k-(backoff-1);
504
505      for (i = 1; i < ub; i++) {
506         quad_float *cptr = mu[i];
507         quad_float *sptr = aux[i];
508   
509         for (j = n; j > i; j--) {
510            c = cptr[j];
511            s = sptr[j];
512   
513            a = c*pp[j-1] - s*pp[j];
514            b = s*pp[j-1] + c*pp[j];
515   
516            pp[j-1] = a;
517            pp[j] = b;
518         }
519   
520         pp[i] = pp[i]/mu[i][i]; 
521      }
522
523      cache.bl[cache.bp] = k;
524      cache.bv[cache.bp] = k-backoff;
525   }
526
527   for (j = 1; j <= n; j++)
528      p[j] = pp[j];
529
530   for (i = max(cache.bv[cache.bp]+1, 1); i < k; i++) {
531      quad_float *cptr = mu[i];
532      quad_float *sptr = aux[i];
533 
534      for (j = n; j > i; j--) {
535         c = cptr[j];
536         s = sptr[j];
537 
538         a = c*p[j-1] - s*p[j];
539         b = s*p[j-1] + c*p[j];
540 
541         p[j-1] = a;
542         p[j] = b;
543      }
544 
545      p[i] = p[i]/mu[i][i];
546   }
547
548   for (j = n; j > k; j--) {
549      a = p[j-1];
550      b = p[j];
551
552      if (b == 0) {
553         c = 1;
554         s = 0;
555      }
556      else if (fabs(b) > fabs(a)) {
557         t = -a/b;
558         s = 1/sqrt(1 + t*t);
559         c = s*t;
560      }
561      else {
562         t = -b/a;
563         c = 1/sqrt(1 + t*t);
564         s = c*t;
565      }
566   
567      p[j-1] = c*a - s*b;
568      p[j] = c;
569      aux[k][j] = s;
570   }
571
572   if (k > n+1) Error("G_LLL_QP: internal error");
573   if (k > n) p[k] = 0;
574
575   for (i = 1; i <= k; i++)
576      CheckFinite(&p[i]);
577}
578
579static quad_float red_fudge = to_quad_float(0);
580static long log_red = 0;
581
582static long verbose = 0;
583
584static unsigned long NumSwaps = 0;
585static double StartTime = 0;
586static double LastTime = 0;
587
588
589
590static void G_LLLStatus(long max_k, double t, long m, const mat_ZZ& B)
591{
592
593   ZZ t1;
594   long i;
595   double prodlen = 0;
596
597   for (i = 1; i <= m; i++) {
598      InnerProduct(t1, B(i), B(i));
599      if (!IsZero(t1))
600         prodlen += log(t1);
601   }
602
603   LastTime = t;
604   
605}
606
607
608static void init_red_fudge()
609{
610   long i;
611
612   // initial log_red should be <= NTL_DOUBLE_PRECISION-2,
613   // to help ensure stability in G_BKZ_QP1
614
615   log_red = NTL_DOUBLE_PRECISION-2; 
616
617   red_fudge = 1;
618
619   for (i = log_red; i > 0; i--)
620      red_fudge = red_fudge*0.5;
621}
622
623static void inc_red_fudge()
624{
625
626   red_fudge = red_fudge * 2;
627   log_red--;
628
629   if (log_red < 4)
630      Error("G_LLL_QP: too much loss of precision...stop!");
631}
632
633
634static
635long ll_G_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
636           LLLCheckFct check, quad_float **B1, quad_float **mu, 
637           quad_float **aux,
638           long m, long init_k, long &quit, GivensCache_QP& cache)
639{
640   long n = B.NumCols();
641
642   long i, j, k, Fc1;
643   ZZ MU;
644   quad_float mu1;
645
646   quad_float t1;
647   double dt1;
648   ZZ T1;
649   quad_float *tp;
650
651   quad_float half = to_quad_float(0.5);
652   quad_float half_plus_fudge = 0.5 + red_fudge;
653
654   quit = 0;
655   k = init_k;
656
657   vec_long in_vec_mem;
658   in_vec_mem.SetLength(n+1);
659   long *in_vec = in_vec_mem.elts();
660
661   double *max_b;
662   max_b = NTL_NEW_OP double [m+1];
663   if (!max_b) Error("out of memory in lll_G_LLL_QP");
664
665   for (i = 1; i <= m; i++)
666      max_b[i] = max_abs(B1[i], n);
667
668   long in_float;
669
670
671   long counter;
672
673   long trigger_index;
674   long small_trigger;
675   long cnt;
676
677   long max_k = 0;
678
679   double tt;
680
681   cache.flush();
682
683   while (k <= m) {
684
685      if (k > max_k) {
686         max_k = k;
687      }
688
689      if (verbose) {
690         tt = GetTime();
691
692         if (tt > LastTime + LLLStatusInterval)
693            G_LLLStatus(max_k, tt, m, B);
694      }
695
696
697      GivensComputeGS(B1, mu, aux, k, n, cache);
698
699      counter = 0;
700      trigger_index = k;
701      small_trigger = 0;
702      cnt = 0;
703
704      do {
705         // size reduction
706
707         counter++;
708         if (counter > 10000) {
709            //cerr << "G_LLL_QP: warning--possible infinite loop\n";
710            counter = 0;
711         }
712
713
714         Fc1 = 0;
715   
716         for (j = k-1; j >= 1; j--) {
717            t1 = fabs(mu[k][j]);
718            if (t1 > half_plus_fudge) {
719
720               if (!Fc1) {
721                  if (j > trigger_index ||
722                      (j == trigger_index && small_trigger)) {
723
724                     cnt++;
725
726                     if (cnt > 10) {
727                        inc_red_fudge();
728                        half_plus_fudge = 0.5 + red_fudge;
729                        cnt = 0;
730                     }
731                  }
732
733                  trigger_index = j;
734                  small_trigger = (t1 < 4);
735
736                  Fc1 = 1;
737                  RowTransformStart(B1[k], in_vec, in_float, n);
738               }
739
740
741   
742               mu1 = mu[k][j];
743               if (mu1 >= 0)
744                  mu1 = ceil(mu1-half);
745               else
746                  mu1 = floor(mu1+half);
747   
748   
749               quad_float *mu_k = mu[k];
750               quad_float *mu_j = mu[j];
751 
752               if (mu1 == 1) {
753                  for (i = 1; i <= j-1; i++)
754                     mu_k[i] -= mu_j[i];
755               }
756               else if (mu1 == -1) {
757                  for (i = 1; i <= j-1; i++)
758                     mu_k[i] += mu_j[i];
759               }
760               else {
761                  for (i = 1; i <= j-1; i++)
762                     mu_k[i] -= mu1*mu_j[i];
763               }
764
765               // cout << j << " " << mu[k][j] << " " << mu1 << "\n";
766 
767               mu_k[j] -= mu1;
768
769               conv(MU, mu1);
770
771   
772               RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec,
773                            max_b[k], max_b[j], in_float);
774
775               if (U) RowTransform((*U)(k), (*U)(j), MU);
776            }
777         }
778
779         if (Fc1) {
780            RowTransformFinish(B(k), B1[k], in_vec);
781            max_b[k] = max_abs(B1[k], n);
782            cache.touch();
783            GivensComputeGS(B1, mu, aux, k, n, cache);
784         }
785      } while (Fc1);
786
787      if (check && (*check)(B(k))) 
788         quit = 1;
789
790      if (IsZero(B(k))) {
791         for (i = k; i < m; i++) {
792            // swap i, i+1
793            swap(B(i), B(i+1));
794            tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp;
795            dt1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = dt1;
796            if (U) swap((*U)(i), (*U)(i+1));
797         }
798
799         cache.flush();
800
801         m--;
802         if (quit) break;
803         continue;
804      }
805
806      if (quit) break;
807
808      if (deep > 0) {
809         // deep insertions
810   
811         Error("sorry...deep insertions not implemented");
812      } // end deep insertions
813
814      // test LLL reduction condition
815
816      if (k > 1 &&
817         sqrt(delta - mu[k][k-1]*mu[k][k-1])*fabs(mu[k-1][k-1]) >
818         fabs(mu[k][k])) {
819
820         // swap rows k, k-1
821         swap(B(k), B(k-1));
822         tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp;
823         dt1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = dt1;
824         if (U) swap((*U)(k), (*U)(k-1));
825
826         cache.swap();
827
828         k--;
829         NumSwaps++;
830         // cout << "- " << k << "\n";
831      }
832      else {
833         cache.incr();
834         k++;
835         // cout << "+ " << k << "\n";
836      }
837   }
838
839   if (verbose) {
840      G_LLLStatus(m+1, GetTime(), m, B);
841   }
842
843
844   delete [] max_b;
845
846   return m;
847}
848
849static
850long G_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
851           LLLCheckFct check)
852{
853   long m = B.NumRows();
854   long n = B.NumCols();
855
856   long i, j;
857   long new_m, dep, quit;
858   quad_float s;
859   ZZ MU;
860   quad_float mu1;
861
862   quad_float t1;
863   ZZ T1;
864
865   init_red_fudge();
866
867   if (U) ident(*U, m);
868
869   quad_float **B1;  // approximates B
870
871   typedef quad_float *quad_floatptr;
872
873   B1 = NTL_NEW_OP quad_floatptr[m+1];
874   if (!B1) Error("G_LLL_QP: out of memory");
875
876   for (i = 1; i <= m; i++) {
877      B1[i] = NTL_NEW_OP quad_float[n+1];
878      if (!B1[i]) Error("G_LLL_QP: out of memory");
879   }
880
881   quad_float **mu;
882   mu = NTL_NEW_OP quad_floatptr[m+1];
883   if (!mu) Error("G_LLL_QP: out of memory");
884
885   for (i = 1; i <= m; i++) {
886      mu[i] = NTL_NEW_OP quad_float[n+2];
887      if (!mu[i]) Error("G_LLL_QP: out of memory");
888   }
889
890   quad_float **aux;
891   aux = NTL_NEW_OP quad_floatptr[m+1];
892   if (!aux) Error("G_LLL_QP: out of memory");
893
894   for (i = 1; i <= m; i++) {
895      aux[i] = NTL_NEW_OP quad_float[n+1];
896      if (!aux[i]) Error("G_LLL_QP: out of memory");
897   }
898
899   for (i = 1; i <=m; i++)
900      for (j = 1; j <= n; j++) {
901         conv(B1[i][j], B(i, j));
902         CheckFinite(&B1[i][j]);
903      }
904
905
906   GivensCache_QP cache(m, n);
907
908   new_m = 
909      ll_G_LLL_QP(B, U, delta, deep, check, B1, mu, aux, m, 1, quit, cache);
910
911
912
913   dep = m - new_m;
914   m = new_m;
915
916   if (dep > 0) {
917      // for consistency, we move all of the zero rows to the front
918
919      for (i = 0; i < m; i++) {
920         swap(B(m+dep-i), B(m-i));
921         if (U) swap((*U)(m+dep-i), (*U)(m-i));
922      }
923   }
924
925
926   // clean-up
927
928   for (i = 1; i <= m; i++) {
929      delete [] B1[i];
930   }
931
932   delete [] B1;
933
934   for (i = 1; i <= m; i++) {
935      delete [] mu[i];
936   }
937
938   delete [] mu;
939
940   for (i = 1; i <= m; i++) {
941      delete [] aux[i];
942   }
943
944   delete [] aux;
945
946   return m;
947}
948
949         
950
951long G_LLL_QP(mat_ZZ& B, double delta, long deep, LLLCheckFct check,
952            long verb)
953{
954   verbose = verb;
955   NumSwaps = 0;
956   if (verbose) {
957      StartTime = GetTime();
958      LastTime = StartTime;
959   }
960
961   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
962   if (deep < 0) Error("G_LLL_QP: bad deep");
963   return G_LLL_QP(B, 0, to_quad_float(delta), deep, check);
964}
965
966long G_LLL_QP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, 
967           LLLCheckFct check, long verb)
968{
969   verbose = verb;
970   NumSwaps = 0;
971   if (verbose) {
972      StartTime = GetTime();
973      LastTime = StartTime;
974   }
975
976
977   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
978   if (deep < 0) Error("G_LLL_QP: bad deep");
979   return G_LLL_QP(B, &U, to_quad_float(delta), deep, check);
980}
981
982
983
984static vec_quad_float G_BKZConstant;
985
986static
987void ComputeG_BKZConstant(long beta, long p)
988{
989   const quad_float c_PI = 
990      to_quad_float("3.141592653589793238462643383279502884197");
991   const quad_float LogPI = 
992      to_quad_float("1.144729885849400174143427351353058711647");
993
994   G_BKZConstant.SetLength(beta-1);
995
996   vec_quad_float Log;
997   Log.SetLength(beta);
998
999
1000   long i, j, k;
1001   quad_float x, y;
1002
1003   for (j = 1; j <= beta; j++)
1004      Log(j) = log(to_quad_float(j));
1005
1006   for (i = 1; i <= beta-1; i++) {
1007      // First, we compute x = gamma(i/2)^{2/i}
1008
1009      k = i/2;
1010
1011      if ((i & 1) == 0) { // i even
1012         x = 0;
1013         for (j = 1; j <= k; j++)
1014            x = x + Log(j);
1015         
1016         x = x * (1/to_quad_float(k));
1017
1018         x = exp(x);
1019      }
1020      else { // i odd
1021         x = 0;
1022         for (j = k + 2; j <= 2*k + 2; j++)
1023            x = x + Log(j);
1024
1025         x = 0.5*LogPI + x - 2*(k+1)*Log(2);
1026
1027         x = x * (2.0/to_quad_float(i));
1028
1029         x = exp(x);
1030      }
1031
1032      // Second, we compute y = 2^{2*p/i}
1033
1034      y = -(2*p/to_quad_float(i))*Log(2);
1035      y = exp(y);
1036
1037      G_BKZConstant(i) = x*y/c_PI;
1038   }
1039}
1040
1041static vec_quad_float G_BKZThresh;
1042
1043static 
1044void ComputeG_BKZThresh(quad_float *c, long beta)
1045{
1046   G_BKZThresh.SetLength(beta-1);
1047
1048   long i;
1049   quad_float x;
1050
1051   x = 0;
1052
1053   for (i = 1; i <= beta-1; i++) {
1054      x += log(c[i-1]);
1055      G_BKZThresh(i) = exp(x/to_quad_float(i))*G_BKZConstant(i);
1056      if (!IsFinite(&G_BKZThresh(i))) G_BKZThresh(i) = 0;
1057   }
1058}
1059
1060
1061static 
1062void G_BKZStatus(double tt, double enum_time, unsigned long NumIterations, 
1063               unsigned long NumTrivial, unsigned long NumNonTrivial, 
1064               unsigned long NumNoOps, long m, 
1065               const mat_ZZ& B)
1066{
1067   ZZ t1;
1068   long i;
1069   double prodlen = 0;
1070
1071   for (i = 1; i <= m; i++) {
1072      InnerProduct(t1, B(i), B(i));
1073      if (!IsZero(t1))
1074         prodlen += log(t1);
1075   }
1076
1077   LastTime = tt;
1078   
1079}
1080
1081
1082static
1083long G_BKZ_QP(mat_ZZ& BB, mat_ZZ* UU, quad_float delta, 
1084         long beta, long prune, LLLCheckFct check)
1085{
1086
1087   long m = BB.NumRows();
1088   long n = BB.NumCols();
1089   long m_orig = m;
1090   
1091   long i, j;
1092   ZZ MU;
1093
1094   quad_float t1;
1095   ZZ T1;
1096   quad_float *tp;
1097
1098   init_red_fudge();
1099
1100   mat_ZZ B;
1101   B = BB;
1102
1103   B.SetDims(m+1, n);
1104
1105
1106   quad_float **B1;  // approximates B
1107
1108   typedef quad_float *quad_floatptr;
1109
1110   B1 = NTL_NEW_OP quad_floatptr[m+2];
1111   if (!B1) Error("G_BKZ_QP: out of memory");
1112
1113   for (i = 1; i <= m+1; i++) {
1114      B1[i] = NTL_NEW_OP quad_float[n+1];
1115      if (!B1[i]) Error("G_BKZ_QP: out of memory");
1116   }
1117
1118   quad_float **mu;
1119   mu = NTL_NEW_OP quad_floatptr[m+2];
1120   if (!mu) Error("G_BKZ_QP: out of memory");
1121
1122   for (i = 1; i <= m+1; i++) {
1123      mu[i] = NTL_NEW_OP quad_float[n+2];
1124      if (!mu[i]) Error("G_BKZ_QP: out of memory");
1125   }
1126
1127   quad_float **aux;
1128   aux = NTL_NEW_OP quad_floatptr[m+2];
1129   if (!aux) Error("G_BKZ_QP: out of memory");
1130
1131   for (i = 1; i <= m+1; i++) {
1132      aux[i] = NTL_NEW_OP quad_float[n+1];
1133      if (!aux[i]) Error("G_BKZ_QP: out of memory");
1134   }
1135
1136   quad_float *c; // squared lengths of Gramm-Schmidt basis vectors
1137
1138   c = NTL_NEW_OP quad_float[m+2];
1139   if (!c) Error("G_BKZ_QP: out of memory");
1140
1141   quad_float cbar;
1142
1143   quad_float *ctilda;
1144   ctilda = NTL_NEW_OP quad_float[m+2];
1145   if (!ctilda) Error("G_BKZ_QP: out of memory");
1146
1147   quad_float *vvec;
1148   vvec = NTL_NEW_OP quad_float[m+2];
1149   if (!vvec) Error("G_BKZ_QP: out of memory");
1150
1151   quad_float *yvec;
1152   yvec = NTL_NEW_OP quad_float[m+2];
1153   if (!yvec) Error("G_BKZ_QP: out of memory");
1154
1155   quad_float *uvec;
1156   uvec = NTL_NEW_OP quad_float[m+2];
1157   if (!uvec) Error("G_BKZ_QP: out of memory");
1158
1159   quad_float *utildavec;
1160   utildavec = NTL_NEW_OP quad_float[m+2];
1161   if (!utildavec) Error("G_BKZ_QP: out of memory");
1162
1163
1164   long *Deltavec;
1165   Deltavec = NTL_NEW_OP long[m+2];
1166   if (!Deltavec) Error("G_BKZ_QP: out of memory");
1167
1168   long *deltavec;
1169   deltavec = NTL_NEW_OP long[m+2];
1170   if (!deltavec) Error("G_BKZ_QP: out of memory");
1171
1172   mat_ZZ Ulocal;
1173   mat_ZZ *U;
1174
1175   if (UU) {
1176      Ulocal.SetDims(m+1, m);
1177      for (i = 1; i <= m; i++)
1178         conv(Ulocal(i, i), 1);
1179      U = &Ulocal;
1180   }
1181   else
1182      U = 0;
1183
1184   long quit;
1185   long new_m;
1186   long z, jj, kk;
1187   long s, t;
1188   long h;
1189   quad_float eta;
1190
1191
1192   for (i = 1; i <=m; i++)
1193      for (j = 1; j <= n; j++) {
1194         conv(B1[i][j], B(i, j));
1195         CheckFinite(&B1[i][j]);
1196      }
1197         
1198
1199   GivensCache_QP cache(m, n);
1200
1201   m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1202
1203   double tt;
1204
1205   double enum_time = 0;
1206   unsigned long NumIterations = 0;
1207   unsigned long NumTrivial = 0;
1208   unsigned long NumNonTrivial = 0;
1209   unsigned long NumNoOps = 0;
1210
1211   long verb = verbose;
1212
1213   verbose = 0;
1214
1215   long clean = 1;
1216
1217   if (m < m_orig) {
1218      for (i = m_orig+1; i >= m+2; i--) {
1219         // swap i, i-1
1220
1221         swap(B(i), B(i-1));
1222         if (U) swap((*U)(i), (*U)(i-1));
1223      }
1224   }
1225
1226   if (!quit && m > 1) {
1227      // cerr << "continuing\n";
1228      if (beta > m) beta = m;
1229
1230      if (prune > 0)
1231         ComputeG_BKZConstant(beta, prune);
1232
1233      z = 0;
1234      jj = 0;
1235   
1236      while (z < m-1) {
1237         jj++;
1238         kk = min(jj+beta-1, m);
1239   
1240         if (jj == m) {
1241            jj = 1;
1242            kk = beta;
1243            clean = 1;
1244         }
1245
1246         if (verb) {
1247            tt = GetTime();
1248            if (tt > LastTime + LLLStatusInterval)
1249               G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1250                         NumNonTrivial, NumNoOps, m, B);
1251         }
1252
1253
1254
1255   
1256         // ENUM
1257
1258         double tt1;
1259
1260         if (verb) {
1261            tt1 = GetTime();
1262         }
1263
1264         for (i = jj; i <= kk; i++) {
1265            c[i] = mu[i][i]*mu[i][i];
1266            CheckFinite(&c[i]);
1267         }
1268
1269         if (prune > 0)
1270            ComputeG_BKZThresh(&c[jj], kk-jj+1);
1271
1272   
1273         cbar = c[jj];
1274         utildavec[jj] = uvec[jj] = 1;
1275   
1276         yvec[jj] = vvec[jj] = 0;
1277         Deltavec[jj] = 0;
1278   
1279   
1280         s = t = jj;
1281         deltavec[jj] = 1;
1282   
1283         for (i = jj+1; i <= kk+1; i++) {
1284            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
1285            Deltavec[i] = 0;
1286            vvec[i] = 0;
1287            deltavec[i] = 1;
1288         }
1289
1290         long enum_cnt = 0;
1291   
1292         while (t <= kk) {
1293            if (verb) {
1294               enum_cnt++;
1295               if (enum_cnt > 100000) {
1296                  enum_cnt = 0;
1297                  tt = GetTime();
1298                  if (tt > LastTime + LLLStatusInterval) {
1299                     enum_time += tt - tt1;
1300                     tt1 = tt;
1301                     G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1302                               NumNonTrivial, NumNoOps, m, B);
1303                  }
1304               }
1305            }
1306
1307            ctilda[t] = ctilda[t+1] + 
1308               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*c[t];
1309   
1310            if (prune > 0 && t > jj) {
1311               eta = G_BKZThresh(t-jj);
1312            }
1313            else
1314               eta = 0;
1315
1316            if (ctilda[t] < cbar - eta) {
1317               if (t > jj) {
1318                  t--;
1319                  t1 = 0;
1320                  for (i = t+1; i <= s; i++) {
1321                     t1 += utildavec[i]*mu[i][t];
1322                  }
1323
1324
1325                  yvec[t] = t1;
1326                  t1 = -t1;
1327                  if (t1 >= 0)
1328                     t1 = ceil(t1-0.5);
1329                  else
1330                     t1 = floor(t1+0.5);
1331
1332                  utildavec[t] = vvec[t] = t1;
1333                  Deltavec[t] = 0;
1334                  if (utildavec[t] > -yvec[t]) 
1335                     deltavec[t] = -1;
1336                  else
1337                     deltavec[t] = 1;
1338               }
1339               else {
1340                  cbar = ctilda[jj];
1341                  for (i = jj; i <= kk; i++) {
1342                     uvec[i] = utildavec[i];
1343                  }
1344               }
1345            }
1346            else {
1347               t++;
1348               s = max(s, t);
1349               if (t < s) Deltavec[t] = -Deltavec[t];
1350               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1351               utildavec[t] = vvec[t] + Deltavec[t];
1352            }
1353         }
1354
1355         if (verb) {
1356            tt1 = GetTime() - tt1;
1357            enum_time += tt1;
1358         }
1359
1360         NumIterations++;
1361
1362         h = min(kk+1, m);
1363   
1364         if ((delta-8*red_fudge)*c[jj] > cbar) {
1365
1366            clean = 0;
1367
1368            // we treat the case that the new vector is b_s (jj < s <= kk)
1369            // as a special case that appears to occur most of the time.
1370   
1371            s = 0;
1372            for (i = jj+1; i <= kk; i++) {
1373               if (uvec[i] != 0) {
1374                  if (s == 0)
1375                     s = i;
1376                  else
1377                     s = -1;
1378               }
1379            }
1380   
1381            if (s == 0) Error("G_BKZ_QP: internal error");
1382   
1383            if (s > 0) {
1384               // special case
1385
1386               NumTrivial++;
1387   
1388               for (i = s; i > jj; i--) {
1389                  // swap i, i-1
1390                  swap(B(i-1), B(i));
1391                  if (U) swap((*U)(i-1), (*U)(i));
1392                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1393               }
1394   
1395               // cerr << "special case\n";
1396               new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1397                                B1, mu, aux, h, jj, quit, cache);
1398               if (new_m != h) Error("G_BKZ_QP: internal error");
1399               if (quit) break;
1400            }
1401            else {
1402               // the general case
1403
1404               NumNonTrivial++;
1405
1406   
1407               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);
1408
1409               if (U) {
1410                  for (i = 1; i <= m_orig; i++)
1411                     conv((*U)(m+1, i), 0);
1412               }
1413
1414               for (i = jj; i <= kk; i++) {
1415                  if (uvec[i] == 0) continue;
1416                  conv(MU, uvec[i]);
1417                  RowTransform2(B(m+1), B(i), MU);
1418                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
1419               }
1420     
1421               for (i = m+1; i >= jj+1; i--) {
1422                  // swap i, i-1
1423                  swap(B(i-1), B(i));
1424                  if (U) swap((*U)(i-1), (*U)(i));
1425                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1426               }
1427     
1428               for (i = 1; i <= n; i++) {
1429                  conv(B1[jj][i], B(jj, i));
1430                  CheckFinite(&B1[jj][i]);
1431               }
1432     
1433               if (IsZero(B(jj))) Error("G_BKZ_QP: internal error"); 
1434     
1435               // remove linear dependencies
1436   
1437               // cerr << "general case\n";
1438               new_m = ll_G_LLL_QP(B, U, delta, 0, 0, B1, mu, aux,
1439                                  kk+1, jj, quit, cache);
1440             
1441               if (new_m != kk) Error("G_BKZ_QP: internal error"); 
1442
1443               // remove zero vector
1444     
1445               for (i = kk+2; i <= m+1; i++) {
1446                  // swap i, i-1
1447                  swap(B(i-1), B(i));
1448                  if (U) swap((*U)(i-1), (*U)(i));
1449                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1450               }
1451     
1452               quit = 0;
1453               if (check) {
1454                  for (i = 1; i <= kk; i++)
1455                     if ((*check)(B(i))) {
1456                        quit = 1;
1457                        break;
1458                     }
1459               }
1460
1461               if (quit) break;
1462   
1463               if (h > kk) {
1464                  // extend reduced basis
1465   
1466                  new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1467                                   B1, mu, aux, h, h, quit, cache);
1468   
1469                  if (new_m != h) Error("G_BKZ_QP: internal error");
1470                  if (quit) break;
1471               }
1472            }
1473   
1474            z = 0;
1475         }
1476         else {
1477            // G_LLL_QP
1478            // cerr << "progress\n";
1479
1480            NumNoOps++;
1481
1482
1483            if (!clean) {
1484               new_m =
1485                  ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, 
1486                              h, h, quit, cache);
1487               if (new_m != h) Error("G_BKZ_QP: internal error");
1488               if (quit) break;
1489            }
1490   
1491            z++;
1492         }
1493      }
1494   }
1495
1496
1497   if (verb) {
1498      G_BKZStatus(GetTime(), enum_time, NumIterations, NumTrivial, NumNonTrivial,
1499                NumNoOps, m, B);
1500   }
1501
1502
1503   // clean up
1504
1505
1506   if (m_orig > m) {
1507      // for consistency, we move zero vectors to the front
1508
1509      for (i = m+1; i <= m_orig; i++) {
1510         swap(B(i), B(i+1));
1511         if (U) swap((*U)(i), (*U)(i+1));
1512      }
1513
1514      for (i = 0; i < m; i++) {
1515         swap(B(m_orig-i), B(m-i));
1516         if (U) swap((*U)(m_orig-i), (*U)(m-i));
1517      }
1518   }
1519
1520   B.SetDims(m_orig, n);
1521   BB = B;
1522
1523   if (U) {
1524      U->SetDims(m_orig, m_orig);
1525      *UU = *U;
1526   }
1527
1528   for (i = 1; i <= m+1; i++) {
1529      delete [] B1[i];
1530   }
1531
1532   delete [] B1;
1533
1534   for (i = 1; i <= m+1; i++) {
1535      delete [] mu[i];
1536   }
1537
1538   delete [] mu;
1539
1540   for (i = 1; i <= m+1; i++) {
1541      delete [] aux[i];
1542   }
1543
1544   delete [] aux;
1545
1546
1547   delete [] c;
1548   delete [] ctilda;
1549   delete [] vvec;
1550   delete [] yvec;
1551   delete [] uvec;
1552   delete [] utildavec;
1553   delete [] Deltavec;
1554   delete [] deltavec;
1555
1556   return m;
1557}
1558
1559long G_BKZ_QP(mat_ZZ& BB, mat_ZZ& UU, double delta, 
1560         long beta, long prune, LLLCheckFct check, long verb)
1561{
1562   verbose = verb;
1563   NumSwaps = 0;
1564   if (verbose) {
1565      StartTime = GetTime();
1566      LastTime = StartTime;
1567   }
1568
1569
1570   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1571   if (beta < 2) Error("G_BKZ_QP: bad block size");
1572
1573   return G_BKZ_QP(BB, &UU, to_quad_float(delta), beta, prune, check);
1574}
1575
1576long G_BKZ_QP(mat_ZZ& BB, double delta, 
1577         long beta, long prune, LLLCheckFct check, long verb)
1578{
1579   verbose = verb;
1580   NumSwaps = 0;
1581   if (verbose) {
1582      StartTime = GetTime();
1583      LastTime = StartTime;
1584   }
1585
1586
1587
1588   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1589   if (beta < 2) Error("G_BKZ_QP: bad block size");
1590
1591   return G_BKZ_QP(BB, 0, to_quad_float(delta), beta, prune, check);
1592}
1593
1594
1595
1596static
1597long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ* UU, quad_float delta, 
1598         long beta, long prune, LLLCheckFct check)
1599{
1600
1601   long m = BB.NumRows();
1602   long n = BB.NumCols();
1603   long m_orig = m;
1604   
1605   long i, j;
1606   ZZ MU;
1607
1608   ZZ T1;
1609   quad_float *tp;
1610
1611   init_red_fudge();
1612
1613   mat_ZZ B;
1614   B = BB;
1615
1616   B.SetDims(m+1, n);
1617
1618
1619   quad_float **B1;  // approximates B
1620
1621   typedef quad_float *quad_floatptr;
1622
1623   B1 = NTL_NEW_OP quad_floatptr[m+2];
1624   if (!B1) Error("G_BKZ_QP: out of memory");
1625
1626   for (i = 1; i <= m+1; i++) {
1627      B1[i] = NTL_NEW_OP quad_float[n+1];
1628      if (!B1[i]) Error("G_BKZ_QP: out of memory");
1629   }
1630
1631   quad_float **mu;
1632   mu = NTL_NEW_OP quad_floatptr[m+2];
1633   if (!mu) Error("G_BKZ_QP: out of memory");
1634
1635   for (i = 1; i <= m+1; i++) {
1636      mu[i] = NTL_NEW_OP quad_float[n+2];
1637      if (!mu[i]) Error("G_BKZ_QP: out of memory");
1638   }
1639
1640   quad_float **aux;
1641   aux = NTL_NEW_OP quad_floatptr[m+2];
1642   if (!aux) Error("G_BKZ_QP: out of memory");
1643
1644   for (i = 1; i <= m+1; i++) {
1645      aux[i] = NTL_NEW_OP quad_float[n+1];
1646      if (!aux[i]) Error("G_BKZ_QP: out of memory");
1647   }
1648
1649   quad_float *c; // squared lengths of Gramm-Schmidt basis vectors
1650
1651   c = NTL_NEW_OP quad_float[m+2];
1652   if (!c) Error("G_BKZ_QP: out of memory");
1653
1654   double cbar;
1655
1656   double *ctilda;
1657   ctilda = NTL_NEW_OP double[m+2];
1658   if (!ctilda) Error("G_BKZ_QP: out of memory");
1659
1660   double *vvec;
1661   vvec = NTL_NEW_OP double[m+2];
1662   if (!vvec) Error("G_BKZ_QP: out of memory");
1663
1664   double *yvec;
1665   yvec = NTL_NEW_OP double[m+2];
1666   if (!yvec) Error("G_BKZ_QP: out of memory");
1667
1668   double *uvec;
1669   uvec = NTL_NEW_OP double[m+2];
1670   if (!uvec) Error("G_BKZ_QP: out of memory");
1671
1672   double *utildavec;
1673   utildavec = NTL_NEW_OP double[m+2];
1674   if (!utildavec) Error("G_BKZ_QP: out of memory");
1675
1676
1677   long *Deltavec;
1678   Deltavec = NTL_NEW_OP long[m+2];
1679   if (!Deltavec) Error("G_BKZ_QP: out of memory");
1680
1681   long *deltavec;
1682   deltavec = NTL_NEW_OP long[m+2];
1683   if (!deltavec) Error("G_BKZ_QP: out of memory");
1684
1685   mat_ZZ Ulocal;
1686   mat_ZZ *U;
1687
1688   if (UU) {
1689      Ulocal.SetDims(m+1, m);
1690      for (i = 1; i <= m; i++)
1691         conv(Ulocal(i, i), 1);
1692      U = &Ulocal;
1693   }
1694   else
1695      U = 0;
1696
1697   long quit;
1698   long new_m;
1699   long z, jj, kk;
1700   long s, t;
1701   long h;
1702
1703   double eta;
1704
1705   for (i = 1; i <=m; i++)
1706      for (j = 1; j <= n; j++) {
1707         conv(B1[i][j], B(i, j));
1708         CheckFinite(&B1[i][j]);
1709      }
1710
1711
1712   GivensCache_QP cache(m, n);
1713
1714   m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1715
1716
1717
1718   double tt;
1719
1720   double enum_time = 0;
1721   unsigned long NumIterations = 0;
1722   unsigned long NumTrivial = 0;
1723   unsigned long NumNonTrivial = 0;
1724   unsigned long NumNoOps = 0;
1725
1726   long verb = verbose;
1727
1728   verbose = 0;
1729
1730   long clean = 1;
1731
1732   if (m < m_orig) {
1733      for (i = m_orig+1; i >= m+2; i--) {
1734         // swap i, i-1
1735
1736         swap(B(i), B(i-1));
1737         if (U) swap((*U)(i), (*U)(i-1));
1738      }
1739   }
1740
1741   if (!quit && m > 1) {
1742      // cerr << "continuing\n";
1743      if (beta > m) beta = m;
1744
1745      if (prune > 0)
1746         ComputeG_BKZConstant(beta, prune);
1747
1748      z = 0;
1749      jj = 0;
1750   
1751      while (z < m-1) {
1752         jj++;
1753         kk = min(jj+beta-1, m);
1754   
1755         if (jj == m) {
1756            jj = 1;
1757            kk = beta;
1758            clean = 1;
1759         }
1760
1761         if (verb) {
1762            tt = GetTime();
1763            if (tt > LastTime + LLLStatusInterval)
1764               G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1765                         NumNonTrivial, NumNoOps, m, B);
1766         }
1767
1768
1769
1770   
1771         // ENUM
1772
1773         double tt1;
1774
1775         if (verb) {
1776            tt1 = GetTime();
1777         }
1778
1779         for (i = jj; i <= kk; i++) {
1780            c[i] = mu[i][i]*mu[i][i];
1781            CheckFinite(&c[i]);
1782         }
1783
1784         if (prune > 0)
1785            ComputeG_BKZThresh(&c[jj], kk-jj+1);
1786
1787   
1788         cbar = to_double(c[jj]);
1789         utildavec[jj] = uvec[jj] = 1;
1790   
1791         yvec[jj] = vvec[jj] = 0;
1792         Deltavec[jj] = 0;
1793   
1794   
1795         s = t = jj;
1796         deltavec[jj] = 1;
1797   
1798         for (i = jj+1; i <= kk+1; i++) {
1799            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
1800            Deltavec[i] = 0;
1801            vvec[i] = 0;
1802            deltavec[i] = 1;
1803         }
1804
1805         long enum_cnt = 0;
1806   
1807         while (t <= kk) {
1808            if (verb) {
1809               enum_cnt++;
1810               if (enum_cnt > 100000) {
1811                  enum_cnt = 0;
1812                  tt = GetTime();
1813                  if (tt > LastTime + LLLStatusInterval) {
1814                     enum_time += tt - tt1;
1815                     tt1 = tt;
1816                     G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1817                               NumNonTrivial, NumNoOps, m, B);
1818                  }
1819               }
1820            }
1821
1822            ctilda[t] = ctilda[t+1] + 
1823               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*to_double(c[t]);
1824   
1825            if (prune > 0 && t > jj) {
1826               eta = to_double(G_BKZThresh(t-jj));
1827            }
1828            else
1829               eta = 0;
1830
1831            if (ctilda[t] < cbar - eta) {
1832               if (t > jj) {
1833                  double t1;
1834
1835                  t--;
1836                  t1 = 0;
1837                  for (i = t+1; i <= s; i++) {
1838                     t1 += utildavec[i]*to_double(mu[i][t]);
1839                  }
1840
1841
1842                  yvec[t] = t1;
1843                  t1 = -t1;
1844                  if (t1 >= 0)
1845                     t1 = ceil(t1-0.5);
1846                  else
1847                     t1 = floor(t1+0.5);
1848
1849                  utildavec[t] = vvec[t] = t1;
1850                  Deltavec[t] = 0;
1851                  if (utildavec[t] > -yvec[t]) 
1852                     deltavec[t] = -1;
1853                  else
1854                     deltavec[t] = 1;
1855               }
1856               else {
1857                  cbar = ctilda[jj];
1858                  for (i = jj; i <= kk; i++) {
1859                     uvec[i] = utildavec[i];
1860                  }
1861               }
1862            }
1863            else {
1864               t++;
1865               s = max(s, t);
1866               if (t < s) Deltavec[t] = -Deltavec[t];
1867               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1868               utildavec[t] = vvec[t] + Deltavec[t];
1869            }
1870         }
1871
1872         if (verb) {
1873            tt1 = GetTime() - tt1;
1874            enum_time += tt1;
1875         }
1876
1877         NumIterations++;
1878
1879         h = min(kk+1, m);
1880
1881         quad_float t1;
1882   
1883         if ((delta-8*red_fudge)*c[jj] > cbar*(1+64/NTL_FDOUBLE_PRECISION)) {
1884
1885            clean = 0;
1886
1887            // we treat the case that the new vector is b_s (jj < s <= kk)
1888            // as a special case that appears to occur most of the time.
1889   
1890            s = 0;
1891            for (i = jj+1; i <= kk; i++) {
1892               if (uvec[i] != 0) {
1893                  if (s == 0)
1894                     s = i;
1895                  else
1896                     s = -1;
1897               }
1898            }
1899   
1900            if (s == 0) Error("G_BKZ_QP: internal error");
1901   
1902            if (s > 0) {
1903               // special case
1904
1905               NumTrivial++;
1906   
1907               for (i = s; i > jj; i--) {
1908                  // swap i, i-1
1909                  swap(B(i-1), B(i));
1910                  if (U) swap((*U)(i-1), (*U)(i));
1911                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1912               }
1913   
1914               // cerr << "special case\n";
1915               new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1916                                B1, mu, aux, h, jj, quit, cache);
1917               if (new_m != h) Error("G_BKZ_QP: internal error");
1918               if (quit) break;
1919            }
1920            else {
1921               // the general case
1922
1923               NumNonTrivial++;
1924
1925   
1926               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);
1927
1928               if (U) {
1929                  for (i = 1; i <= m_orig; i++)
1930                     conv((*U)(m+1, i), 0);
1931               }
1932
1933               for (i = jj; i <= kk; i++) {
1934                  if (uvec[i] == 0) continue;
1935                  conv(MU, uvec[i]);
1936                  RowTransform2(B(m+1), B(i), MU);
1937                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
1938               }
1939     
1940               for (i = m+1; i >= jj+1; i--) {
1941                  // swap i, i-1
1942                  swap(B(i-1), B(i));
1943                  if (U) swap((*U)(i-1), (*U)(i));
1944                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1945               }
1946     
1947               for (i = 1; i <= n; i++) {
1948                  conv(B1[jj][i], B(jj, i));
1949                  CheckFinite(&B1[jj][i]);
1950               }
1951     
1952               if (IsZero(B(jj))) Error("G_BKZ_QP: internal error"); 
1953     
1954               // remove linear dependencies
1955   
1956               // cerr << "general case\n";
1957               new_m = ll_G_LLL_QP(B, U, delta, 0, 0, B1, mu, aux,
1958                                  kk+1, jj, quit, cache);
1959             
1960               if (new_m != kk) Error("G_BKZ_QP: internal error"); 
1961
1962               // remove zero vector
1963     
1964               for (i = kk+2; i <= m+1; i++) {
1965                  // swap i, i-1
1966                  swap(B(i-1), B(i));
1967                  if (U) swap((*U)(i-1), (*U)(i));
1968                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1969               }
1970     
1971               quit = 0;
1972               if (check) {
1973                  for (i = 1; i <= kk; i++)
1974                     if ((*check)(B(i))) {
1975                        quit = 1;
1976                        break;
1977                     }
1978               }
1979
1980               if (quit) break;
1981   
1982               if (h > kk) {
1983                  // extend reduced basis
1984   
1985                  new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1986                                   B1, mu, aux, h, h, quit, cache);
1987   
1988                  if (new_m != h) Error("G_BKZ_QP: internal error");
1989                  if (quit) break;
1990               }
1991            }
1992   
1993            z = 0;
1994         }
1995         else {
1996            // G_LLL_QP
1997            // cerr << "progress\n";
1998
1999            NumNoOps++;
2000
2001
2002            if (!clean) {
2003               new_m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux,
2004                                   h, h, quit, cache);
2005
2006               if (new_m != h) Error("G_BKZ_QP: internal error");
2007               if (quit) break;
2008            }
2009   
2010            z++;
2011         }
2012      }
2013   }
2014
2015
2016   if (verb) {
2017      G_BKZStatus(GetTime(), enum_time, NumIterations, NumTrivial, NumNonTrivial,
2018                NumNoOps, m, B);
2019   }
2020
2021
2022   // clean up
2023
2024
2025   if (m_orig > m) {
2026      // for consistency, we move zero vectors to the front
2027
2028      for (i = m+1; i <= m_orig; i++) {
2029         swap(B(i), B(i+1));
2030         if (U) swap((*U)(i), (*U)(i+1));
2031      }
2032
2033      for (i = 0; i < m; i++) {
2034         swap(B(m_orig-i), B(m-i));
2035         if (U) swap((*U)(m_orig-i), (*U)(m-i));
2036      }
2037   }
2038
2039   B.SetDims(m_orig, n);
2040   BB = B;
2041
2042   if (U) {
2043      U->SetDims(m_orig, m_orig);
2044      *UU = *U;
2045   }
2046
2047   for (i = 1; i <= m+1; i++) {
2048      delete [] B1[i];
2049   }
2050
2051   delete [] B1;
2052
2053   for (i = 1; i <= m+1; i++) {
2054      delete [] mu[i];
2055   }
2056
2057   delete [] mu;
2058
2059   for (i = 1; i <= m+1; i++) {
2060      delete [] aux[i];
2061   }
2062
2063   delete [] aux;
2064
2065
2066   delete [] c;
2067   delete [] ctilda;
2068   delete [] vvec;
2069   delete [] yvec;
2070   delete [] uvec;
2071   delete [] utildavec;
2072   delete [] Deltavec;
2073   delete [] deltavec;
2074
2075   return m;
2076}
2077
2078long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ& UU, double delta, 
2079         long beta, long prune, LLLCheckFct check, long verb)
2080{
2081   verbose = verb;
2082   NumSwaps = 0;
2083   if (verbose) {
2084      StartTime = GetTime();
2085      LastTime = StartTime;
2086   }
2087
2088
2089   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
2090   if (beta < 2) Error("G_BKZ_QP: bad block size");
2091
2092   return G_BKZ_QP1(BB, &UU, to_quad_float(delta), beta, prune, check);
2093}
2094
2095long G_BKZ_QP1(mat_ZZ& BB, double delta, 
2096         long beta, long prune, LLLCheckFct check, long verb)
2097{
2098   verbose = verb;
2099   NumSwaps = 0;
2100   if (verbose) {
2101      StartTime = GetTime();
2102      LastTime = StartTime;
2103   }
2104
2105
2106
2107   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
2108   if (beta < 2) Error("G_BKZ_QP: bad block size");
2109
2110   return G_BKZ_QP1(BB, 0, to_quad_float(delta), beta, prune, check);
2111}
2112
2113NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.