source: git/ntl/src/G_LLL_QP.c @ 45fefa

fieker-DuValspielwiese
Last change on this file since 45fefa was 447abc, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: NTL 5.4.1 git-svn-id: file:///usr/local/Singular/svn/trunk@10336 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 40.3 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 init_red_fudge()
591{
592   long i;
593
594   // initial log_red should be <= NTL_DOUBLE_PRECISION-2,
595   // to help ensure stability in G_BKZ_QP1
596
597   log_red = NTL_DOUBLE_PRECISION-2; 
598
599   red_fudge = 1;
600
601   for (i = log_red; i > 0; i--)
602      red_fudge = red_fudge*0.5;
603}
604
605static void inc_red_fudge()
606{
607
608   red_fudge = red_fudge * 2;
609   log_red--;
610
611   //cerr << "G_LLL_QP: warning--relaxing reduction (" << log_red << ")\n";
612
613   if (log_red < 4)
614      Error("G_LLL_QP: too much loss of precision...stop!");
615}
616
617
618static
619long ll_G_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
620           LLLCheckFct check, quad_float **B1, quad_float **mu, 
621           quad_float **aux,
622           long m, long init_k, long &quit, GivensCache_QP& cache)
623{
624   long n = B.NumCols();
625
626   long i, j, k, Fc1;
627   ZZ MU;
628   quad_float mu1;
629
630   quad_float t1;
631   double dt1;
632   ZZ T1;
633   quad_float *tp;
634
635   quad_float half = to_quad_float(0.5);
636   quad_float half_plus_fudge = 0.5 + red_fudge;
637
638   quit = 0;
639   k = init_k;
640
641   vec_long in_vec_mem;
642   in_vec_mem.SetLength(n+1);
643   long *in_vec = in_vec_mem.elts();
644
645   double *max_b;
646   max_b = NTL_NEW_OP double [m+1];
647   if (!max_b) Error("out of memory in lll_G_LLL_QP");
648
649   for (i = 1; i <= m; i++)
650      max_b[i] = max_abs(B1[i], n);
651
652   long in_float;
653
654
655   long counter;
656
657   long trigger_index;
658   long small_trigger;
659   long cnt;
660
661   long max_k = 0;
662
663   double tt;
664
665   cache.flush();
666
667   while (k <= m) {
668
669      if (k > max_k) {
670         max_k = k;
671      }
672
673      GivensComputeGS(B1, mu, aux, k, n, cache);
674
675      counter = 0;
676      trigger_index = k;
677      small_trigger = 0;
678      cnt = 0;
679
680      do {
681         // size reduction
682
683         counter++;
684         if (counter > 10000) {
685            Error("G_LLL_QP: warning--possible infinite loop\n");
686            counter = 0;
687         }
688
689
690         Fc1 = 0;
691   
692         for (j = k-1; j >= 1; j--) {
693            t1 = fabs(mu[k][j]);
694            if (t1 > half_plus_fudge) {
695
696               if (!Fc1) {
697                  if (j > trigger_index ||
698                      (j == trigger_index && small_trigger)) {
699
700                     cnt++;
701
702                     if (cnt > 10) {
703                        inc_red_fudge();
704                        half_plus_fudge = 0.5 + red_fudge;
705                        cnt = 0;
706                     }
707                  }
708
709                  trigger_index = j;
710                  small_trigger = (t1 < 4);
711
712                  Fc1 = 1;
713                  RowTransformStart(B1[k], in_vec, in_float, n);
714               }
715
716
717   
718               mu1 = mu[k][j];
719               if (mu1 >= 0)
720                  mu1 = ceil(mu1-half);
721               else
722                  mu1 = floor(mu1+half);
723   
724   
725               quad_float *mu_k = mu[k];
726               quad_float *mu_j = mu[j];
727 
728               if (mu1 == 1) {
729                  for (i = 1; i <= j-1; i++)
730                     mu_k[i] -= mu_j[i];
731               }
732               else if (mu1 == -1) {
733                  for (i = 1; i <= j-1; i++)
734                     mu_k[i] += mu_j[i];
735               }
736               else {
737                  for (i = 1; i <= j-1; i++)
738                     mu_k[i] -= mu1*mu_j[i];
739               }
740
741               // cout << j << " " << mu[k][j] << " " << mu1 << "\n";
742 
743               mu_k[j] -= mu1;
744
745               conv(MU, mu1);
746
747   
748               RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec,
749                            max_b[k], max_b[j], in_float);
750
751               if (U) RowTransform((*U)(k), (*U)(j), MU);
752            }
753         }
754
755         if (Fc1) {
756            RowTransformFinish(B(k), B1[k], in_vec);
757            max_b[k] = max_abs(B1[k], n);
758            cache.touch();
759            GivensComputeGS(B1, mu, aux, k, n, cache);
760         }
761      } while (Fc1);
762
763      if (check && (*check)(B(k))) 
764         quit = 1;
765
766      if (IsZero(B(k))) {
767         for (i = k; i < m; i++) {
768            // swap i, i+1
769            swap(B(i), B(i+1));
770            tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp;
771            dt1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = dt1;
772            if (U) swap((*U)(i), (*U)(i+1));
773         }
774
775         cache.flush();
776
777         m--;
778         if (quit) break;
779         continue;
780      }
781
782      if (quit) break;
783
784      if (deep > 0) {
785         // deep insertions
786   
787         Error("sorry...deep insertions not implemented");
788      } // end deep insertions
789
790      // test LLL reduction condition
791
792      if (k > 1 &&
793         sqrt(delta - mu[k][k-1]*mu[k][k-1])*fabs(mu[k-1][k-1]) >
794         fabs(mu[k][k])) {
795
796         // swap rows k, k-1
797         swap(B(k), B(k-1));
798         tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp;
799         dt1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = dt1;
800         if (U) swap((*U)(k), (*U)(k-1));
801
802         cache.swap();
803
804         k--;
805         NumSwaps++;
806         // cout << "- " << k << "\n";
807      }
808      else {
809         cache.incr();
810         k++;
811         // cout << "+ " << k << "\n";
812      }
813   }
814
815   delete [] max_b;
816
817   return m;
818}
819
820static
821long G_LLL_QP(mat_ZZ& B, mat_ZZ* U, quad_float delta, long deep, 
822           LLLCheckFct check)
823{
824   long m = B.NumRows();
825   long n = B.NumCols();
826
827   long i, j;
828   long new_m, dep, quit;
829   quad_float s;
830   ZZ MU;
831   quad_float mu1;
832
833   quad_float t1;
834   ZZ T1;
835
836   init_red_fudge();
837
838   if (U) ident(*U, m);
839
840   quad_float **B1;  // approximates B
841
842   typedef quad_float *quad_floatptr;
843
844   B1 = NTL_NEW_OP quad_floatptr[m+1];
845   if (!B1) Error("G_LLL_QP: out of memory");
846
847   for (i = 1; i <= m; i++) {
848      B1[i] = NTL_NEW_OP quad_float[n+1];
849      if (!B1[i]) Error("G_LLL_QP: out of memory");
850   }
851
852   quad_float **mu;
853   mu = NTL_NEW_OP quad_floatptr[m+1];
854   if (!mu) Error("G_LLL_QP: out of memory");
855
856   for (i = 1; i <= m; i++) {
857      mu[i] = NTL_NEW_OP quad_float[n+2];
858      if (!mu[i]) Error("G_LLL_QP: out of memory");
859   }
860
861   quad_float **aux;
862   aux = NTL_NEW_OP quad_floatptr[m+1];
863   if (!aux) Error("G_LLL_QP: out of memory");
864
865   for (i = 1; i <= m; i++) {
866      aux[i] = NTL_NEW_OP quad_float[n+1];
867      if (!aux[i]) Error("G_LLL_QP: out of memory");
868   }
869
870   for (i = 1; i <=m; i++)
871      for (j = 1; j <= n; j++) {
872         conv(B1[i][j], B(i, j));
873         CheckFinite(&B1[i][j]);
874      }
875
876
877   GivensCache_QP cache(m, n);
878
879   new_m = 
880      ll_G_LLL_QP(B, U, delta, deep, check, B1, mu, aux, m, 1, quit, cache);
881
882
883
884   dep = m - new_m;
885   m = new_m;
886
887   if (dep > 0) {
888      // for consistency, we move all of the zero rows to the front
889
890      for (i = 0; i < m; i++) {
891         swap(B(m+dep-i), B(m-i));
892         if (U) swap((*U)(m+dep-i), (*U)(m-i));
893      }
894   }
895
896
897   // clean-up
898
899   for (i = 1; i <= m+dep; i++) {
900      delete [] B1[i];
901   }
902
903   delete [] B1;
904
905   for (i = 1; i <= m+dep; i++) {
906      delete [] mu[i];
907   }
908
909   delete [] mu;
910
911   for (i = 1; i <= m+dep; i++) {
912      delete [] aux[i];
913   }
914
915   delete [] aux;
916
917   return m;
918}
919
920         
921
922long G_LLL_QP(mat_ZZ& B, double delta, long deep, LLLCheckFct check,
923            long verb)
924{
925   verbose = verb;
926   NumSwaps = 0;
927
928   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
929   if (deep < 0) Error("G_LLL_QP: bad deep");
930   return G_LLL_QP(B, 0, to_quad_float(delta), deep, check);
931}
932
933long G_LLL_QP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, 
934           LLLCheckFct check, long verb)
935{
936   verbose = verb;
937   NumSwaps = 0;
938
939   if (delta < 0.50 || delta >= 1) Error("G_LLL_QP: bad delta");
940   if (deep < 0) Error("G_LLL_QP: bad deep");
941   return G_LLL_QP(B, &U, to_quad_float(delta), deep, check);
942}
943
944
945
946static vec_quad_float G_BKZConstant;
947
948static
949void ComputeG_BKZConstant(long beta, long p)
950{
951   const quad_float c_PI = 
952      to_quad_float("3.141592653589793238462643383279502884197");
953   const quad_float LogPI = 
954      to_quad_float("1.144729885849400174143427351353058711647");
955
956   G_BKZConstant.SetLength(beta-1);
957
958   vec_quad_float Log;
959   Log.SetLength(beta);
960
961
962   long i, j, k;
963   quad_float x, y;
964
965   for (j = 1; j <= beta; j++)
966      Log(j) = log(to_quad_float(j));
967
968   for (i = 1; i <= beta-1; i++) {
969      // First, we compute x = gamma(i/2)^{2/i}
970
971      k = i/2;
972
973      if ((i & 1) == 0) { // i even
974         x = 0;
975         for (j = 1; j <= k; j++)
976            x = x + Log(j);
977         
978         x = x * (1/to_quad_float(k));
979
980         x = exp(x);
981      }
982      else { // i odd
983         x = 0;
984         for (j = k + 2; j <= 2*k + 2; j++)
985            x = x + Log(j);
986
987         x = 0.5*LogPI + x - 2*(k+1)*Log(2);
988
989         x = x * (2.0/to_quad_float(i));
990
991         x = exp(x);
992      }
993
994      // Second, we compute y = 2^{2*p/i}
995
996      y = -(2*p/to_quad_float(i))*Log(2);
997      y = exp(y);
998
999      G_BKZConstant(i) = x*y/c_PI;
1000   }
1001}
1002
1003static vec_quad_float G_BKZThresh;
1004
1005static 
1006void ComputeG_BKZThresh(quad_float *c, long beta)
1007{
1008   G_BKZThresh.SetLength(beta-1);
1009
1010   long i;
1011   quad_float x;
1012
1013   x = 0;
1014
1015   for (i = 1; i <= beta-1; i++) {
1016      x += log(c[i-1]);
1017      G_BKZThresh(i) = exp(x/to_quad_float(i))*G_BKZConstant(i);
1018      if (!IsFinite(&G_BKZThresh(i))) G_BKZThresh(i) = 0;
1019   }
1020}
1021
1022
1023static
1024long G_BKZ_QP(mat_ZZ& BB, mat_ZZ* UU, quad_float delta, 
1025         long beta, long prune, LLLCheckFct check)
1026{
1027
1028   long m = BB.NumRows();
1029   long n = BB.NumCols();
1030   long m_orig = m;
1031   
1032   long i, j;
1033   ZZ MU;
1034
1035   quad_float t1;
1036   ZZ T1;
1037   quad_float *tp;
1038
1039   init_red_fudge();
1040
1041   mat_ZZ B;
1042   B = BB;
1043
1044   B.SetDims(m+1, n);
1045
1046
1047   quad_float **B1;  // approximates B
1048
1049   typedef quad_float *quad_floatptr;
1050
1051   B1 = NTL_NEW_OP quad_floatptr[m+2];
1052   if (!B1) Error("G_BKZ_QP: out of memory");
1053
1054   for (i = 1; i <= m+1; i++) {
1055      B1[i] = NTL_NEW_OP quad_float[n+1];
1056      if (!B1[i]) Error("G_BKZ_QP: out of memory");
1057   }
1058
1059   quad_float **mu;
1060   mu = NTL_NEW_OP quad_floatptr[m+2];
1061   if (!mu) Error("G_BKZ_QP: out of memory");
1062
1063   for (i = 1; i <= m+1; i++) {
1064      mu[i] = NTL_NEW_OP quad_float[n+2];
1065      if (!mu[i]) Error("G_BKZ_QP: out of memory");
1066   }
1067
1068   quad_float **aux;
1069   aux = NTL_NEW_OP quad_floatptr[m+2];
1070   if (!aux) Error("G_BKZ_QP: out of memory");
1071
1072   for (i = 1; i <= m+1; i++) {
1073      aux[i] = NTL_NEW_OP quad_float[n+1];
1074      if (!aux[i]) Error("G_BKZ_QP: out of memory");
1075   }
1076
1077   quad_float *c; // squared lengths of Gramm-Schmidt basis vectors
1078
1079   c = NTL_NEW_OP quad_float[m+2];
1080   if (!c) Error("G_BKZ_QP: out of memory");
1081
1082   quad_float cbar;
1083
1084   quad_float *ctilda;
1085   ctilda = NTL_NEW_OP quad_float[m+2];
1086   if (!ctilda) Error("G_BKZ_QP: out of memory");
1087
1088   quad_float *vvec;
1089   vvec = NTL_NEW_OP quad_float[m+2];
1090   if (!vvec) Error("G_BKZ_QP: out of memory");
1091
1092   quad_float *yvec;
1093   yvec = NTL_NEW_OP quad_float[m+2];
1094   if (!yvec) Error("G_BKZ_QP: out of memory");
1095
1096   quad_float *uvec;
1097   uvec = NTL_NEW_OP quad_float[m+2];
1098   if (!uvec) Error("G_BKZ_QP: out of memory");
1099
1100   quad_float *utildavec;
1101   utildavec = NTL_NEW_OP quad_float[m+2];
1102   if (!utildavec) Error("G_BKZ_QP: out of memory");
1103
1104
1105   long *Deltavec;
1106   Deltavec = NTL_NEW_OP long[m+2];
1107   if (!Deltavec) Error("G_BKZ_QP: out of memory");
1108
1109   long *deltavec;
1110   deltavec = NTL_NEW_OP long[m+2];
1111   if (!deltavec) Error("G_BKZ_QP: out of memory");
1112
1113   mat_ZZ Ulocal;
1114   mat_ZZ *U;
1115
1116   if (UU) {
1117      Ulocal.SetDims(m+1, m);
1118      for (i = 1; i <= m; i++)
1119         conv(Ulocal(i, i), 1);
1120      U = &Ulocal;
1121   }
1122   else
1123      U = 0;
1124
1125   long quit;
1126   long new_m;
1127   long z, jj, kk;
1128   long s, t;
1129   long h;
1130   quad_float eta;
1131
1132
1133   for (i = 1; i <=m; i++)
1134      for (j = 1; j <= n; j++) {
1135         conv(B1[i][j], B(i, j));
1136         CheckFinite(&B1[i][j]);
1137      }
1138         
1139
1140   GivensCache_QP cache(m, n);
1141
1142   m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1143
1144   double tt;
1145
1146   double enum_time = 0;
1147   unsigned long NumIterations = 0;
1148   unsigned long NumTrivial = 0;
1149   unsigned long NumNonTrivial = 0;
1150   unsigned long NumNoOps = 0;
1151
1152   long verb = verbose;
1153
1154   verbose = 0;
1155
1156   long clean = 1;
1157
1158   if (m < m_orig) {
1159      for (i = m_orig+1; i >= m+2; i--) {
1160         // swap i, i-1
1161
1162         swap(B(i), B(i-1));
1163         if (U) swap((*U)(i), (*U)(i-1));
1164      }
1165   }
1166
1167   if (!quit && m > 1) {
1168      // cerr << "continuing\n";
1169      if (beta > m) beta = m;
1170
1171      if (prune > 0)
1172         ComputeG_BKZConstant(beta, prune);
1173
1174      z = 0;
1175      jj = 0;
1176   
1177      while (z < m-1) {
1178         jj++;
1179         kk = min(jj+beta-1, m);
1180   
1181         if (jj == m) {
1182            jj = 1;
1183            kk = beta;
1184            clean = 1;
1185         }
1186
1187         // ENUM
1188
1189         double tt1;
1190
1191         for (i = jj; i <= kk; i++) {
1192            c[i] = mu[i][i]*mu[i][i];
1193            CheckFinite(&c[i]);
1194         }
1195
1196         if (prune > 0)
1197            ComputeG_BKZThresh(&c[jj], kk-jj+1);
1198
1199   
1200         cbar = c[jj];
1201         utildavec[jj] = uvec[jj] = 1;
1202   
1203         yvec[jj] = vvec[jj] = 0;
1204         Deltavec[jj] = 0;
1205   
1206   
1207         s = t = jj;
1208         deltavec[jj] = 1;
1209   
1210         for (i = jj+1; i <= kk+1; i++) {
1211            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
1212            Deltavec[i] = 0;
1213            vvec[i] = 0;
1214            deltavec[i] = 1;
1215         }
1216
1217         long enum_cnt = 0;
1218   
1219         while (t <= kk) {
1220
1221            ctilda[t] = ctilda[t+1] + 
1222               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*c[t];
1223   
1224            if (prune > 0 && t > jj) {
1225               eta = G_BKZThresh(t-jj);
1226            }
1227            else
1228               eta = 0;
1229
1230            if (ctilda[t] < cbar - eta) {
1231               if (t > jj) {
1232                  t--;
1233                  t1 = 0;
1234                  for (i = t+1; i <= s; i++) {
1235                     t1 += utildavec[i]*mu[i][t];
1236                  }
1237
1238
1239                  yvec[t] = t1;
1240                  t1 = -t1;
1241                  if (t1 >= 0)
1242                     t1 = ceil(t1-0.5);
1243                  else
1244                     t1 = floor(t1+0.5);
1245
1246                  utildavec[t] = vvec[t] = t1;
1247                  Deltavec[t] = 0;
1248                  if (utildavec[t] > -yvec[t]) 
1249                     deltavec[t] = -1;
1250                  else
1251                     deltavec[t] = 1;
1252               }
1253               else {
1254                  cbar = ctilda[jj];
1255                  for (i = jj; i <= kk; i++) {
1256                     uvec[i] = utildavec[i];
1257                  }
1258               }
1259            }
1260            else {
1261               t++;
1262               s = max(s, t);
1263               if (t < s) Deltavec[t] = -Deltavec[t];
1264               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1265               utildavec[t] = vvec[t] + Deltavec[t];
1266            }
1267         }
1268
1269         NumIterations++;
1270
1271         h = min(kk+1, m);
1272   
1273         if ((delta-8*red_fudge)*c[jj] > cbar) {
1274
1275            clean = 0;
1276
1277            // we treat the case that the new vector is b_s (jj < s <= kk)
1278            // as a special case that appears to occur most of the time.
1279   
1280            s = 0;
1281            for (i = jj+1; i <= kk; i++) {
1282               if (uvec[i] != 0) {
1283                  if (s == 0)
1284                     s = i;
1285                  else
1286                     s = -1;
1287               }
1288            }
1289   
1290            if (s == 0) Error("G_BKZ_QP: internal error");
1291   
1292            if (s > 0) {
1293               // special case
1294
1295               NumTrivial++;
1296   
1297               for (i = s; i > jj; i--) {
1298                  // swap i, i-1
1299                  swap(B(i-1), B(i));
1300                  if (U) swap((*U)(i-1), (*U)(i));
1301                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1302               }
1303   
1304               // cerr << "special case\n";
1305               new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1306                                B1, mu, aux, h, jj, quit, cache);
1307               if (new_m != h) Error("G_BKZ_QP: internal error");
1308               if (quit) break;
1309            }
1310            else {
1311               // the general case
1312
1313               NumNonTrivial++;
1314
1315   
1316               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);
1317
1318               if (U) {
1319                  for (i = 1; i <= m_orig; i++)
1320                     conv((*U)(m+1, i), 0);
1321               }
1322
1323               for (i = jj; i <= kk; i++) {
1324                  if (uvec[i] == 0) continue;
1325                  conv(MU, uvec[i]);
1326                  RowTransform2(B(m+1), B(i), MU);
1327                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
1328               }
1329     
1330               for (i = m+1; i >= jj+1; i--) {
1331                  // swap i, i-1
1332                  swap(B(i-1), B(i));
1333                  if (U) swap((*U)(i-1), (*U)(i));
1334                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1335               }
1336     
1337               for (i = 1; i <= n; i++) {
1338                  conv(B1[jj][i], B(jj, i));
1339                  CheckFinite(&B1[jj][i]);
1340               }
1341     
1342               if (IsZero(B(jj))) Error("G_BKZ_QP: internal error"); 
1343     
1344               // remove linear dependencies
1345   
1346               // cerr << "general case\n";
1347               new_m = ll_G_LLL_QP(B, U, delta, 0, 0, B1, mu, aux,
1348                                  kk+1, jj, quit, cache);
1349             
1350               if (new_m != kk) Error("G_BKZ_QP: internal error"); 
1351
1352               // remove zero vector
1353     
1354               for (i = kk+2; i <= m+1; i++) {
1355                  // swap i, i-1
1356                  swap(B(i-1), B(i));
1357                  if (U) swap((*U)(i-1), (*U)(i));
1358                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1359               }
1360     
1361               quit = 0;
1362               if (check) {
1363                  for (i = 1; i <= kk; i++)
1364                     if ((*check)(B(i))) {
1365                        quit = 1;
1366                        break;
1367                     }
1368               }
1369
1370               if (quit) break;
1371   
1372               if (h > kk) {
1373                  // extend reduced basis
1374   
1375                  new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1376                                   B1, mu, aux, h, h, quit, cache);
1377   
1378                  if (new_m != h) Error("G_BKZ_QP: internal error");
1379                  if (quit) break;
1380               }
1381            }
1382   
1383            z = 0;
1384         }
1385         else {
1386            // G_LLL_QP
1387            // cerr << "progress\n";
1388
1389            NumNoOps++;
1390
1391
1392            if (!clean) {
1393               new_m =
1394                  ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, 
1395                              h, h, quit, cache);
1396               if (new_m != h) Error("G_BKZ_QP: internal error");
1397               if (quit) break;
1398            }
1399   
1400            z++;
1401         }
1402      }
1403   }
1404
1405
1406   // clean up
1407
1408
1409   if (m_orig > m) {
1410      // for consistency, we move zero vectors to the front
1411
1412      for (i = m+1; i <= m_orig; i++) {
1413         swap(B(i), B(i+1));
1414         if (U) swap((*U)(i), (*U)(i+1));
1415      }
1416
1417      for (i = 0; i < m; i++) {
1418         swap(B(m_orig-i), B(m-i));
1419         if (U) swap((*U)(m_orig-i), (*U)(m-i));
1420      }
1421   }
1422
1423   B.SetDims(m_orig, n);
1424   BB = B;
1425
1426   if (U) {
1427      U->SetDims(m_orig, m_orig);
1428      *UU = *U;
1429   }
1430
1431   for (i = 1; i <= m_orig+1; i++) {
1432      delete [] B1[i];
1433   }
1434
1435   delete [] B1;
1436
1437   for (i = 1; i <= m_orig+1; i++) {
1438      delete [] mu[i];
1439   }
1440
1441   delete [] mu;
1442
1443   for (i = 1; i <= m_orig+1; i++) {
1444      delete [] aux[i];
1445   }
1446
1447   delete [] aux;
1448
1449
1450   delete [] c;
1451   delete [] ctilda;
1452   delete [] vvec;
1453   delete [] yvec;
1454   delete [] uvec;
1455   delete [] utildavec;
1456   delete [] Deltavec;
1457   delete [] deltavec;
1458
1459   return m;
1460}
1461
1462long G_BKZ_QP(mat_ZZ& BB, mat_ZZ& UU, double delta, 
1463         long beta, long prune, LLLCheckFct check, long verb)
1464{
1465   verbose = verb;
1466   NumSwaps = 0;
1467
1468
1469   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1470   if (beta < 2) Error("G_BKZ_QP: bad block size");
1471
1472   return G_BKZ_QP(BB, &UU, to_quad_float(delta), beta, prune, check);
1473}
1474
1475long G_BKZ_QP(mat_ZZ& BB, double delta, 
1476         long beta, long prune, LLLCheckFct check, long verb)
1477{
1478   verbose = verb;
1479   NumSwaps = 0;
1480
1481   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1482   if (beta < 2) Error("G_BKZ_QP: bad block size");
1483
1484   return G_BKZ_QP(BB, 0, to_quad_float(delta), beta, prune, check);
1485}
1486
1487
1488
1489static
1490long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ* UU, quad_float delta, 
1491         long beta, long prune, LLLCheckFct check)
1492{
1493
1494   long m = BB.NumRows();
1495   long n = BB.NumCols();
1496   long m_orig = m;
1497   
1498   long i, j;
1499   ZZ MU;
1500
1501   ZZ T1;
1502   quad_float *tp;
1503
1504   init_red_fudge();
1505
1506   mat_ZZ B;
1507   B = BB;
1508
1509   B.SetDims(m+1, n);
1510
1511
1512   quad_float **B1;  // approximates B
1513
1514   typedef quad_float *quad_floatptr;
1515
1516   B1 = NTL_NEW_OP quad_floatptr[m+2];
1517   if (!B1) Error("G_BKZ_QP: out of memory");
1518
1519   for (i = 1; i <= m+1; i++) {
1520      B1[i] = NTL_NEW_OP quad_float[n+1];
1521      if (!B1[i]) Error("G_BKZ_QP: out of memory");
1522   }
1523
1524   quad_float **mu;
1525   mu = NTL_NEW_OP quad_floatptr[m+2];
1526   if (!mu) Error("G_BKZ_QP: out of memory");
1527
1528   for (i = 1; i <= m+1; i++) {
1529      mu[i] = NTL_NEW_OP quad_float[n+2];
1530      if (!mu[i]) Error("G_BKZ_QP: out of memory");
1531   }
1532
1533   quad_float **aux;
1534   aux = NTL_NEW_OP quad_floatptr[m+2];
1535   if (!aux) Error("G_BKZ_QP: out of memory");
1536
1537   for (i = 1; i <= m+1; i++) {
1538      aux[i] = NTL_NEW_OP quad_float[n+1];
1539      if (!aux[i]) Error("G_BKZ_QP: out of memory");
1540   }
1541
1542   quad_float *c; // squared lengths of Gramm-Schmidt basis vectors
1543
1544   c = NTL_NEW_OP quad_float[m+2];
1545   if (!c) Error("G_BKZ_QP: out of memory");
1546
1547   double cbar;
1548
1549   double *ctilda;
1550   ctilda = NTL_NEW_OP double[m+2];
1551   if (!ctilda) Error("G_BKZ_QP: out of memory");
1552
1553   double *vvec;
1554   vvec = NTL_NEW_OP double[m+2];
1555   if (!vvec) Error("G_BKZ_QP: out of memory");
1556
1557   double *yvec;
1558   yvec = NTL_NEW_OP double[m+2];
1559   if (!yvec) Error("G_BKZ_QP: out of memory");
1560
1561   double *uvec;
1562   uvec = NTL_NEW_OP double[m+2];
1563   if (!uvec) Error("G_BKZ_QP: out of memory");
1564
1565   double *utildavec;
1566   utildavec = NTL_NEW_OP double[m+2];
1567   if (!utildavec) Error("G_BKZ_QP: out of memory");
1568
1569
1570   long *Deltavec;
1571   Deltavec = NTL_NEW_OP long[m+2];
1572   if (!Deltavec) Error("G_BKZ_QP: out of memory");
1573
1574   long *deltavec;
1575   deltavec = NTL_NEW_OP long[m+2];
1576   if (!deltavec) Error("G_BKZ_QP: out of memory");
1577
1578   mat_ZZ Ulocal;
1579   mat_ZZ *U;
1580
1581   if (UU) {
1582      Ulocal.SetDims(m+1, m);
1583      for (i = 1; i <= m; i++)
1584         conv(Ulocal(i, i), 1);
1585      U = &Ulocal;
1586   }
1587   else
1588      U = 0;
1589
1590   long quit;
1591   long new_m;
1592   long z, jj, kk;
1593   long s, t;
1594   long h;
1595
1596   double eta;
1597
1598   for (i = 1; i <=m; i++)
1599      for (j = 1; j <= n; j++) {
1600         conv(B1[i][j], B(i, j));
1601         CheckFinite(&B1[i][j]);
1602      }
1603
1604
1605   GivensCache_QP cache(m, n);
1606
1607   m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1608
1609
1610
1611   double tt;
1612
1613   double enum_time = 0;
1614   unsigned long NumIterations = 0;
1615   unsigned long NumTrivial = 0;
1616   unsigned long NumNonTrivial = 0;
1617   unsigned long NumNoOps = 0;
1618
1619   long verb = verbose;
1620
1621   verbose = 0;
1622
1623   long clean = 1;
1624
1625   if (m < m_orig) {
1626      for (i = m_orig+1; i >= m+2; i--) {
1627         // swap i, i-1
1628
1629         swap(B(i), B(i-1));
1630         if (U) swap((*U)(i), (*U)(i-1));
1631      }
1632   }
1633
1634   if (!quit && m > 1) {
1635      // cerr << "continuing\n";
1636      if (beta > m) beta = m;
1637
1638      if (prune > 0)
1639         ComputeG_BKZConstant(beta, prune);
1640
1641      z = 0;
1642      jj = 0;
1643   
1644      while (z < m-1) {
1645         jj++;
1646         kk = min(jj+beta-1, m);
1647   
1648         if (jj == m) {
1649            jj = 1;
1650            kk = beta;
1651            clean = 1;
1652         }
1653
1654         // ENUM
1655
1656         double tt1;
1657
1658         for (i = jj; i <= kk; i++) {
1659            c[i] = mu[i][i]*mu[i][i];
1660            CheckFinite(&c[i]);
1661         }
1662
1663         if (prune > 0)
1664            ComputeG_BKZThresh(&c[jj], kk-jj+1);
1665
1666   
1667         cbar = to_double(c[jj]);
1668         utildavec[jj] = uvec[jj] = 1;
1669   
1670         yvec[jj] = vvec[jj] = 0;
1671         Deltavec[jj] = 0;
1672   
1673   
1674         s = t = jj;
1675         deltavec[jj] = 1;
1676   
1677         for (i = jj+1; i <= kk+1; i++) {
1678            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
1679            Deltavec[i] = 0;
1680            vvec[i] = 0;
1681            deltavec[i] = 1;
1682         }
1683
1684         long enum_cnt = 0;
1685   
1686         while (t <= kk) {
1687
1688            ctilda[t] = ctilda[t+1] + 
1689               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*to_double(c[t]);
1690
1691            ForceToMem(&ctilda[t]); // prevents an infinite loop
1692   
1693            if (prune > 0 && t > jj) {
1694               eta = to_double(G_BKZThresh(t-jj));
1695            }
1696            else
1697               eta = 0;
1698
1699            if (ctilda[t] < cbar - eta) {
1700               if (t > jj) {
1701                  double t1;
1702
1703                  t--;
1704                  t1 = 0;
1705                  for (i = t+1; i <= s; i++) {
1706                     t1 += utildavec[i]*to_double(mu[i][t]);
1707                  }
1708
1709
1710                  yvec[t] = t1;
1711                  t1 = -t1;
1712                  if (t1 >= 0)
1713                     t1 = ceil(t1-0.5);
1714                  else
1715                     t1 = floor(t1+0.5);
1716
1717                  utildavec[t] = vvec[t] = t1;
1718                  Deltavec[t] = 0;
1719                  if (utildavec[t] > -yvec[t]) 
1720                     deltavec[t] = -1;
1721                  else
1722                     deltavec[t] = 1;
1723               }
1724               else {
1725                  cbar = ctilda[jj];
1726                  for (i = jj; i <= kk; i++) {
1727                     uvec[i] = utildavec[i];
1728                  }
1729               }
1730            }
1731            else {
1732               t++;
1733               s = max(s, t);
1734               if (t < s) Deltavec[t] = -Deltavec[t];
1735               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1736               utildavec[t] = vvec[t] + Deltavec[t];
1737            }
1738         }
1739
1740         NumIterations++;
1741
1742         h = min(kk+1, m);
1743
1744         quad_float t1;
1745   
1746         if ((delta-8*red_fudge)*c[jj] > cbar*(1+64/NTL_FDOUBLE_PRECISION)) {
1747
1748            clean = 0;
1749
1750            // we treat the case that the new vector is b_s (jj < s <= kk)
1751            // as a special case that appears to occur most of the time.
1752   
1753            s = 0;
1754            for (i = jj+1; i <= kk; i++) {
1755               if (uvec[i] != 0) {
1756                  if (s == 0)
1757                     s = i;
1758                  else
1759                     s = -1;
1760               }
1761            }
1762   
1763            if (s == 0) Error("G_BKZ_QP: internal error");
1764   
1765            if (s > 0) {
1766               // special case
1767
1768               NumTrivial++;
1769   
1770               for (i = s; i > jj; i--) {
1771                  // swap i, i-1
1772                  swap(B(i-1), B(i));
1773                  if (U) swap((*U)(i-1), (*U)(i));
1774                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1775               }
1776   
1777               // cerr << "special case\n";
1778               new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1779                                B1, mu, aux, h, jj, quit, cache);
1780               if (new_m != h) Error("G_BKZ_QP: internal error");
1781               if (quit) break;
1782            }
1783            else {
1784               // the general case
1785
1786               NumNonTrivial++;
1787
1788   
1789               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);
1790
1791               if (U) {
1792                  for (i = 1; i <= m_orig; i++)
1793                     conv((*U)(m+1, i), 0);
1794               }
1795
1796               for (i = jj; i <= kk; i++) {
1797                  if (uvec[i] == 0) continue;
1798                  conv(MU, uvec[i]);
1799                  RowTransform2(B(m+1), B(i), MU);
1800                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
1801               }
1802     
1803               for (i = m+1; i >= jj+1; i--) {
1804                  // swap i, i-1
1805                  swap(B(i-1), B(i));
1806                  if (U) swap((*U)(i-1), (*U)(i));
1807                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1808               }
1809     
1810               for (i = 1; i <= n; i++) {
1811                  conv(B1[jj][i], B(jj, i));
1812                  CheckFinite(&B1[jj][i]);
1813               }
1814     
1815               if (IsZero(B(jj))) Error("G_BKZ_QP: internal error"); 
1816     
1817               // remove linear dependencies
1818   
1819               // cerr << "general case\n";
1820               new_m = ll_G_LLL_QP(B, U, delta, 0, 0, B1, mu, aux,
1821                                  kk+1, jj, quit, cache);
1822             
1823               if (new_m != kk) Error("G_BKZ_QP: internal error"); 
1824
1825               // remove zero vector
1826     
1827               for (i = kk+2; i <= m+1; i++) {
1828                  // swap i, i-1
1829                  swap(B(i-1), B(i));
1830                  if (U) swap((*U)(i-1), (*U)(i));
1831                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1832               }
1833     
1834               quit = 0;
1835               if (check) {
1836                  for (i = 1; i <= kk; i++)
1837                     if ((*check)(B(i))) {
1838                        quit = 1;
1839                        break;
1840                     }
1841               }
1842
1843               if (quit) break;
1844   
1845               if (h > kk) {
1846                  // extend reduced basis
1847   
1848                  new_m = ll_G_LLL_QP(B, U, delta, 0, check,
1849                                   B1, mu, aux, h, h, quit, cache);
1850   
1851                  if (new_m != h) Error("G_BKZ_QP: internal error");
1852                  if (quit) break;
1853               }
1854            }
1855   
1856            z = 0;
1857         }
1858         else {
1859            // G_LLL_QP
1860            // cerr << "progress\n";
1861
1862            NumNoOps++;
1863
1864
1865            if (!clean) {
1866               new_m = ll_G_LLL_QP(B, U, delta, 0, check, B1, mu, aux,
1867                                   h, h, quit, cache);
1868
1869               if (new_m != h) Error("G_BKZ_QP: internal error");
1870               if (quit) break;
1871            }
1872   
1873            z++;
1874         }
1875      }
1876   }
1877
1878
1879   // clean up
1880
1881
1882   if (m_orig > m) {
1883      // for consistency, we move zero vectors to the front
1884
1885      for (i = m+1; i <= m_orig; i++) {
1886         swap(B(i), B(i+1));
1887         if (U) swap((*U)(i), (*U)(i+1));
1888      }
1889
1890      for (i = 0; i < m; i++) {
1891         swap(B(m_orig-i), B(m-i));
1892         if (U) swap((*U)(m_orig-i), (*U)(m-i));
1893      }
1894   }
1895
1896   B.SetDims(m_orig, n);
1897   BB = B;
1898
1899   if (U) {
1900      U->SetDims(m_orig, m_orig);
1901      *UU = *U;
1902   }
1903
1904   for (i = 1; i <= m_orig+1; i++) {
1905      delete [] B1[i];
1906   }
1907
1908   delete [] B1;
1909
1910   for (i = 1; i <= m_orig+1; i++) {
1911      delete [] mu[i];
1912   }
1913
1914   delete [] mu;
1915
1916   for (i = 1; i <= m_orig+1; i++) {
1917      delete [] aux[i];
1918   }
1919
1920   delete [] aux;
1921
1922
1923   delete [] c;
1924   delete [] ctilda;
1925   delete [] vvec;
1926   delete [] yvec;
1927   delete [] uvec;
1928   delete [] utildavec;
1929   delete [] Deltavec;
1930   delete [] deltavec;
1931
1932   return m;
1933}
1934
1935long G_BKZ_QP1(mat_ZZ& BB, mat_ZZ& UU, double delta, 
1936         long beta, long prune, LLLCheckFct check, long verb)
1937{
1938   verbose = verb;
1939   NumSwaps = 0;
1940
1941   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1942   if (beta < 2) Error("G_BKZ_QP: bad block size");
1943
1944   return G_BKZ_QP1(BB, &UU, to_quad_float(delta), beta, prune, check);
1945}
1946
1947long G_BKZ_QP1(mat_ZZ& BB, double delta, 
1948         long beta, long prune, LLLCheckFct check, long verb)
1949{
1950   verbose = verb;
1951   NumSwaps = 0;
1952
1953   if (delta < 0.50 || delta >= 1) Error("G_BKZ_QP: bad delta");
1954   if (beta < 2) Error("G_BKZ_QP: bad block size");
1955
1956   return G_BKZ_QP1(BB, 0, to_quad_float(delta), beta, prune, check);
1957}
1958
1959NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.