source: git/ntl/src/G_LLL_QP.c @ 2cfffe

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