source: git/ntl/src/G_LLL_FP.c @ 54e84ca

fieker-DuValspielwiese
Last change on this file since 54e84ca was de6a29, checked in by Hans Schönemann <hannes@…>, 19 years ago
* hannes: NTL-5.4 git-svn-id: file:///usr/local/Singular/svn/trunk@8693 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 29.1 KB
Line 
1
2#include <NTL/LLL.h>
3#include <NTL/fileio.h>
4#include <NTL/vec_double.h>
5
6
7#include <NTL/new.h>
8
9NTL_START_IMPL
10
11static inline
12void CheckFinite(double *p)
13{
14   if (!IsFinite(p)) Error("G_LLL_FP: numbers too big...use G_LLL_XD");
15}
16
17
18
19static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
20// x = x - y*MU
21{
22   static ZZ T, MU;
23   long k;
24
25   long n = A.length();
26   long i;
27
28   MU = MU1;
29
30   if (MU == 1) {
31      for (i = 1; i <= n; i++)
32         sub(A(i), A(i), B(i));
33
34      return;
35   }
36
37   if (MU == -1) {
38      for (i = 1; i <= n; i++)
39         add(A(i), A(i), B(i));
40
41      return;
42   }
43
44   if (MU == 0) return;
45
46   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
47      k = MakeOdd(MU);
48   else
49      k = 0;
50
51
52   if (MU.WideSinglePrecision()) {
53      long mu1;
54      conv(mu1, MU);
55
56      for (i = 1; i <= n; i++) {
57         mul(T, B(i), mu1);
58         if (k > 0) LeftShift(T, T, k);
59         sub(A(i), A(i), T);
60      }
61   }
62   else {
63      for (i = 1; i <= n; i++) {
64         mul(T, B(i), MU);
65         if (k > 0) LeftShift(T, T, k);
66         sub(A(i), A(i), T);
67      }
68   }
69}
70
71
72#define TR_BND (NTL_FDOUBLE_PRECISION/2.0)
73// Just to be safe!!
74
75static double max_abs(double *v, long n)
76{
77   long i;
78   double res, t;
79
80   res = 0;
81
82   for (i = 1; i <= n; i++) {
83      t = fabs(v[i]);
84      if (t > res) res = t;
85   }
86
87   return res;
88}
89
90
91static void RowTransformStart(double *a, long *in_a, long& in_float, long n)
92{
93   long i;
94   long inf = 1;
95
96   for (i = 1; i <= n; i++) {
97      in_a[i] = (a[i] < TR_BND && a[i] > -TR_BND);
98      inf = inf & in_a[i];
99   }
100
101   in_float = inf;
102}
103
104
105static void RowTransformFinish(vec_ZZ& A, double *a, long *in_a)
106{
107   long n = A.length();
108   long i;
109
110   for (i = 1; i <= n; i++) {
111      if (in_a[i])  {
112         conv(A(i), a[i]);
113      }
114      else {
115         conv(a[i], A(i));
116         CheckFinite(&a[i]);
117      }
118   }
119}
120
121
122static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1, 
123                         double *a, double *b, long *in_a,
124                         double& max_a, double max_b, long& in_float)
125// x = x - y*MU
126{
127   static ZZ T, MU;
128   long k;
129   double mu;
130
131   conv(mu, MU1);
132   CheckFinite(&mu);
133
134   long n = A.length();
135   long i;
136
137   if (in_float) {
138      double mu_abs = fabs(mu);
139      if (mu_abs > 0 && max_b > 0 && (mu_abs >= TR_BND || max_b >= TR_BND)) {
140         in_float = 0;
141      }
142      else {
143         max_a += mu_abs*max_b;
144         if (max_a >= TR_BND)
145            in_float = 0;
146      }
147   }
148
149   if (in_float) {
150      if (mu == 1) {
151         for (i = 1; i <= n; i++)
152            a[i] -= b[i];
153
154         return;
155      }
156
157      if (mu == -1) {
158         for (i = 1; i <= n; i++)
159            a[i] += b[i];
160
161         return;
162      }
163
164      if (mu == 0) return;
165
166      for (i = 1; i <= n; i++)
167         a[i] -= mu*b[i];
168
169
170      return;
171   }
172
173
174   MU = MU1;
175
176   if (MU == 1) {
177      for (i = 1; i <= n; i++) {
178         if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND &&
179             b[i] < TR_BND && b[i] > -TR_BND) {
180
181            a[i] -= b[i];
182         }
183         else {
184            if (in_a[i]) {
185               conv(A(i), a[i]);
186               in_a[i] = 0;
187            }
188         
189            sub(A(i), A(i), B(i));
190         }
191      }
192      return;
193   }
194
195   if (MU == -1) {
196      for (i = 1; i <= n; i++) {
197         if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND &&
198             b[i] < TR_BND && b[i] > -TR_BND) {
199
200            a[i] += b[i];
201         }
202         else {
203            if (in_a[i]) {
204               conv(A(i), a[i]);
205               in_a[i] = 0;
206            }
207         
208            add(A(i), A(i), B(i));
209         }
210      }
211      return;
212   }
213
214   if (MU == 0) return;
215
216   double b_bnd = fabs(TR_BND/mu) - 1;
217   if (b_bnd < 0) b_bnd = 0; 
218
219   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
220      k = MakeOdd(MU);
221   else
222      k = 0;
223
224
225   if (MU.WideSinglePrecision()) {
226      long mu1;
227      conv(mu1, MU);
228
229      if (k > 0) {
230         for (i = 1; i <= n; i++) {
231            if (in_a[i]) {
232               conv(A(i), a[i]);
233               in_a[i] = 0;
234            }
235
236            mul(T, B(i), mu1);
237            LeftShift(T, T, k);
238            sub(A(i), A(i), T);
239         }
240      }
241      else {
242         for (i = 1; i <= n; i++) {
243            if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND &&
244                b[i] < b_bnd && b[i] > -b_bnd) {
245   
246               a[i] -= b[i]*mu;
247            }
248            else {
249               if (in_a[i]) {
250                  conv(A(i), a[i]);
251                  in_a[i] = 0;
252               }
253               mul(T, B(i), mu1);
254               sub(A(i), A(i), T);
255            }
256         }
257      }
258   }
259   else {
260      for (i = 1; i <= n; i++) {
261         if (in_a[i]) {
262            conv(A(i), a[i]);
263            in_a[i] = 0;
264         }
265         mul(T, B(i), MU);
266         if (k > 0) LeftShift(T, T, k);
267         sub(A(i), A(i), T);
268      }
269   }
270}
271
272static void RowTransform2(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1)
273// x = x + y*MU
274
275{
276   static ZZ T, MU;
277   long k;
278
279   long n = A.length();
280   long i;
281
282   MU = MU1;
283
284   if (MU == 1) {
285      for (i = 1; i <= n; i++)
286         add(A(i), A(i), B(i));
287
288      return;
289   }
290
291   if (MU == -1) {
292      for (i = 1; i <= n; i++)
293         sub(A(i), A(i), B(i));
294
295      return;
296   }
297
298   if (MU == 0) return;
299
300   if (NumTwos(MU) >= NTL_ZZ_NBITS) 
301      k = MakeOdd(MU);
302   else
303      k = 0;
304
305   if (MU.WideSinglePrecision()) {
306      long mu1;
307      conv(mu1, MU);
308
309      for (i = 1; i <= n; i++) {
310         mul(T, B(i), mu1);
311         if (k > 0) LeftShift(T, T, k);
312         add(A(i), A(i), T);
313      }
314   }
315   else {
316      for (i = 1; i <= n; i++) {
317         mul(T, B(i), MU);
318         if (k > 0) LeftShift(T, T, k);
319         add(A(i), A(i), T);
320      }
321   }
322}
323
324
325
326class GivensCache_FP {
327public:
328   GivensCache_FP(long m, long n);
329   ~GivensCache_FP();
330
331   void flush();
332   void selective_flush(long l);
333   void swap(long l);
334   void swap();
335   void touch();
336   void incr();
337
338   long sz;
339
340   double **buf;
341   long *bl;
342   long *bv;
343   long bp;
344};
345
346
347GivensCache_FP::GivensCache_FP(long m, long n)
348{
349   sz = min(m, n)/10;
350   if (sz < 2) 
351      sz = 2;
352   else if (sz > 20)
353      sz = 20;
354
355   typedef double *doubleptr;
356
357   long i;
358   buf = NTL_NEW_OP doubleptr[sz]; 
359   if (!buf) Error("out of memory");
360   for (i = 0; i < sz; i++)
361      if (!(buf[i] = NTL_NEW_OP double[n+1])) Error("out of memory");
362
363   bl = NTL_NEW_OP long[sz];
364   if (!bl) Error("out of memory");
365   for (i = 0; i < sz; i++) bl[0] = 0;
366
367   bv = NTL_NEW_OP long[sz];
368   if (!bv) Error("out of memory");
369   for (i = 0; i < sz; i++) bv[0] = 0;
370
371   bp = 0;
372}
373
374GivensCache_FP::~GivensCache_FP()
375{
376   long i;
377
378   for (i = 0; i < sz; i++) delete [] buf[i];
379   delete [] buf;
380   delete [] bl;
381   delete [] bv;
382}
383
384void GivensCache_FP::flush()
385{
386   long i;
387   for (i = 0; i < sz; i++) bl[i] = 0;
388}
389
390void GivensCache_FP::selective_flush(long l)
391{
392   long i;
393
394   for (i = 0; i < sz; i++)
395      if (bl[i] && bv[i] >= l)
396         bl[i] = 0;
397}
398
399void GivensCache_FP::swap(long l)
400{
401   long k = bl[bp];
402   long i;
403
404   i = 0;
405   while (i < sz && bl[i] != l)
406      i++;
407
408   if (i < sz) {
409      bl[bp] = l;
410      bl[i] = k;
411   }
412   else
413      bl[bp] = l;
414
415   selective_flush(l);
416}
417
418void GivensCache_FP::swap()
419{
420   swap(bl[bp] - 1);
421}
422
423void GivensCache_FP::touch()
424{
425   long k = bl[bp];
426   bl[bp] = 0;
427   selective_flush(k);
428}
429
430void GivensCache_FP::incr()
431{
432   long k = bl[bp];
433   long k1 = k+1;
434   long i;
435
436   i = 0;
437   while (i < sz && bl[i] != k1)
438      i++;
439
440   if (i < sz) {
441      bp = i;
442      return;
443   }
444
445   i = 0; 
446   while (i < sz && bl[i] != 0)
447      i++;
448
449   if (i < sz) {
450      bp = i;
451      return;
452   }
453
454   long max_val = 0;
455   long max_index = 0;
456   for (i = 0; i < sz; i++) {
457      long t = labs(bl[i]-k1);
458      if (t > max_val) {
459         max_val = t;
460         max_index = i;
461      }
462   }
463
464   bp = max_index;
465   bl[max_index] = 0;
466}
467
468
469static
470void GivensComputeGS(double **B1, double **mu, double **aux, long k, long n,
471                     GivensCache_FP& cache)
472{
473   long i, j;
474
475   double c, s, a, b, t;
476
477   double *p = mu[k];
478
479   double *pp = cache.buf[cache.bp];
480
481   if (!cache.bl[cache.bp]) {
482      for (j = 1; j <= n; j++)
483         pp[j] = B1[k][j];
484
485      long backoff;
486      backoff = k/4;
487      if (backoff < 2)
488         backoff = 2;
489      else if (backoff > cache.sz + 2)
490         backoff = cache.sz + 2; 
491
492      long ub = k-(backoff-1);
493
494      for (i = 1; i < ub; i++) {
495         double *cptr = mu[i];
496         double *sptr = aux[i];
497   
498         for (j = n; j > i; j--) {
499            c = cptr[j];
500            s = sptr[j];
501   
502            a = c*pp[j-1] - s*pp[j];
503            b = s*pp[j-1] + c*pp[j];
504   
505            pp[j-1] = a;
506            pp[j] = b;
507         }
508   
509         pp[i] = pp[i]/mu[i][i]; 
510      }
511
512      cache.bl[cache.bp] = k;
513      cache.bv[cache.bp] = k-backoff;
514   }
515
516   for (j = 1; j <= n; j++)
517      p[j] = pp[j];
518
519   for (i = max(cache.bv[cache.bp]+1, 1); i < k; i++) {
520      double *cptr = mu[i];
521      double *sptr = aux[i];
522 
523      for (j = n; j > i; j--) {
524         c = cptr[j];
525         s = sptr[j];
526 
527         a = c*p[j-1] - s*p[j];
528         b = s*p[j-1] + c*p[j];
529 
530         p[j-1] = a;
531         p[j] = b;
532      }
533 
534      p[i] = p[i]/mu[i][i];
535   }
536
537   for (j = n; j > k; j--) {
538      a = p[j-1];
539      b = p[j];
540
541      if (b == 0) {
542         c = 1;
543         s = 0;
544      }
545      else if (fabs(b) > fabs(a)) {
546         t = -a/b;
547         s = 1/sqrt(1 + t*t);
548         c = s*t;
549      }
550      else {
551         t = -b/a;
552         c = 1/sqrt(1 + t*t);
553         s = c*t;
554      }
555   
556      p[j-1] = c*a - s*b;
557      p[j] = c;
558      aux[k][j] = s;
559   }
560
561   if (k > n+1) Error("G_LLL_FP: internal error");
562   if (k > n) p[k] = 0;
563
564   for (i = 1; i <= k; i++)
565      CheckFinite(&p[i]);
566}
567
568static double red_fudge = 0;
569static long log_red = 0;
570
571static long verbose = 0;
572
573static unsigned long NumSwaps = 0;
574static double StartTime = 0;
575static double LastTime = 0;
576
577
578
579static void init_red_fudge()
580{
581   long i;
582
583   log_red = long(0.50*NTL_DOUBLE_PRECISION);
584   red_fudge = 1;
585
586   for (i = log_red; i > 0; i--)
587      red_fudge = red_fudge*0.5;
588}
589
590static void inc_red_fudge()
591{
592
593   red_fudge = red_fudge * 2;
594   log_red--;
595
596   
597   if (log_red < 4)
598      Error("G_LLL_FP: too much loss of precision...stop!");
599}
600
601
602
603static
604long ll_G_LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, 
605           LLLCheckFct check, double **B1, double **mu, 
606           double **aux,
607           long m, long init_k, long &quit, GivensCache_FP& cache)
608{
609   long n = B.NumCols();
610
611   long i, j, k, Fc1;
612   ZZ MU;
613   double mu1;
614
615   double t1;
616   ZZ T1;
617   double *tp;
618
619   double half_plus_fudge = 0.5 + red_fudge;
620
621   quit = 0;
622   k = init_k;
623
624   vec_long in_vec_mem;
625   in_vec_mem.SetLength(n+1);
626   long *in_vec = in_vec_mem.elts();
627
628   double *max_b;
629   max_b = NTL_NEW_OP double [m+1];
630   if (!max_b) Error("out of memory in lll_G_LLL_FP");
631
632   for (i = 1; i <= m; i++)
633      max_b[i] = max_abs(B1[i], n);
634
635   long in_float;
636
637   long counter;
638
639   long trigger_index;
640   long small_trigger;
641   long cnt;
642
643
644   long max_k = 0;
645
646   double tt;
647
648   long swap_cnt = 0;
649
650   cache.flush();
651
652   while (k <= m) {
653
654      if (k > max_k) {
655         max_k = k;
656         swap_cnt = 0;
657      }
658
659      GivensComputeGS(B1, mu, aux, k, n, cache);
660
661      counter = 0;
662      trigger_index = k;
663      small_trigger = 0;
664      cnt = 0;
665
666      long sz=0, new_sz;
667
668      do {
669         // size reduction
670
671         counter++;
672         if ((counter & 127) == 0) {
673
674            new_sz = 0;
675            for (j = 1; j <= n; j++)
676               new_sz += NumBits(B(k,j));
677
678            if ((counter >> 7) == 1 || new_sz < sz) {
679               sz = new_sz;
680            }
681         }
682
683         Fc1 = 0;
684   
685         for (j = k-1; j >= 1; j--) {
686            t1 = fabs(mu[k][j]);
687            if (t1 > half_plus_fudge) { 
688
689
690               if (!Fc1) {
691                  if (j > trigger_index || 
692                      (j == trigger_index && small_trigger)) {
693
694                     cnt++;
695
696                     if (cnt > 10) {
697                        inc_red_fudge();
698                        half_plus_fudge = 0.5 + red_fudge;
699                        cnt = 0;
700                     }
701                  }
702
703                  trigger_index = j;
704                  small_trigger = (t1 < 4);
705
706                  Fc1 = 1;
707                  RowTransformStart(B1[k], in_vec, in_float, n);
708               }
709                 
710
711               mu1 = mu[k][j];
712               if (mu1 >= 0)
713                  mu1 = ceil(mu1-0.5);
714               else
715                  mu1 = floor(mu1+0.5);
716   
717               double *mu_k = mu[k];
718               double *mu_j = mu[j];
719   
720               if (mu1 == 1) {
721                  for (i = 1; i <= j-1; i++)
722                     mu_k[i] -= mu_j[i];
723               }
724               else if (mu1 == -1) {
725                  for (i = 1; i <= j-1; i++)
726                     mu_k[i] += mu_j[i];
727               }
728               else {
729                  for (i = 1; i <= j-1; i++)
730                     mu_k[i] -= mu1*mu_j[i];
731               }
732   
733               mu_k[j] -= mu1;
734   
735               conv(MU, mu1);
736
737               RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec,
738                            max_b[k], max_b[j], in_float);
739               if (U) RowTransform((*U)(k), (*U)(j), MU);
740            }
741         }
742
743
744         if (Fc1) {
745            RowTransformFinish(B(k), B1[k], in_vec);
746            max_b[k] = max_abs(B1[k], n);
747            cache.touch();
748            GivensComputeGS(B1, mu, aux, k, n, cache);
749         }
750      } while (Fc1);
751
752      if (check && (*check)(B(k))) 
753         quit = 1;
754
755      if (IsZero(B(k))) {
756         for (i = k; i < m; i++) {
757            // swap i, i+1
758            swap(B(i), B(i+1));
759            tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp;
760            t1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = t1;
761            if (U) swap((*U)(i), (*U)(i+1));
762         }
763
764         cache.flush();
765
766         m--;
767         if (quit) break;
768         continue;
769      }
770
771      if (quit) break;
772
773      if (deep > 0) {
774         // deep insertions
775
776         Error("sorry...deep insertions not implemented");
777      } // end deep insertions
778
779      // test G_LLL reduction condition
780
781      if (k > 1 && 
782         sqrt(delta - mu[k][k-1]*mu[k][k-1])*fabs(mu[k-1][k-1]) > 
783         fabs(mu[k][k])) {
784         // swap rows k, k-1
785
786         swap(B(k), B(k-1));
787         tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp;
788         t1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = t1;
789         if (U) swap((*U)(k), (*U)(k-1));
790
791         cache.swap();
792
793         k--;
794         NumSwaps++;
795         swap_cnt++;
796         // cout << "-\n";
797      }
798      else {
799
800         cache.incr();
801
802         k++;
803         // cout << "+\n";
804      }
805
806   }
807
808   delete [] max_b;
809
810   return m;
811}
812
813
814
815
816
817static
818long G_LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, 
819           LLLCheckFct check)
820{
821   long m = B.NumRows();
822   long n = B.NumCols();
823
824   long i, j;
825   long new_m, dep, quit;
826   ZZ MU;
827
828   ZZ T1;
829
830   init_red_fudge();
831
832   if (U) ident(*U, m);
833
834   double **B1;  // approximates B
835
836   typedef double *doubleptr;
837
838   B1 = NTL_NEW_OP doubleptr[m+1];
839   if (!B1) Error("G_LLL_FP: out of memory");
840
841   for (i = 1; i <= m; i++) {
842      B1[i] = NTL_NEW_OP double[n+1];
843      if (!B1[i]) Error("G_LLL_FP: out of memory");
844   }
845
846   double **mu;
847   mu = NTL_NEW_OP doubleptr[m+1];
848   if (!mu) Error("G_LLL_FP: out of memory");
849
850   for (i = 1; i <= m; i++) {
851      mu[i] = NTL_NEW_OP double[n+2];
852      if (!mu[i]) Error("G_LLL_FP: out of memory");
853   }
854
855   double **aux;
856   aux = NTL_NEW_OP doubleptr[m+1];
857   if (!aux) Error("G_LLL_FP: out of memory");
858
859   for (i = 1; i <= m; i++) {
860      aux[i] = NTL_NEW_OP double[n+1];
861      if (!aux[i]) Error("G_LLL_FP: out of memory");
862   }
863
864   for (i = 1; i <=m; i++)
865      for (j = 1; j <= n; j++) {
866         conv(B1[i][j], B(i, j));
867         CheckFinite(&B1[i][j]);
868      }
869
870         
871   GivensCache_FP cache(m, n);
872
873   new_m = ll_G_LLL_FP(B, U, delta, deep, check, B1, mu, aux, m, 1, quit, cache);
874   dep = m - new_m;
875   m = new_m;
876
877   if (dep > 0) {
878      // for consistency, we move all of the zero rows to the front
879
880      for (i = 0; i < m; i++) {
881         swap(B(m+dep-i), B(m-i));
882         if (U) swap((*U)(m+dep-i), (*U)(m-i));
883      }
884   }
885
886
887   // clean-up
888
889   for (i = 1; i <= m; i++) {
890      delete [] B1[i];
891   }
892
893   delete [] B1;
894
895   for (i = 1; i <= m; i++) {
896      delete [] mu[i];
897   }
898
899   delete [] mu;
900
901   for (i = 1; i <= m; i++) {
902      delete [] aux[i];
903   }
904
905   delete [] aux;
906
907   return m;
908}
909
910         
911
912long G_LLL_FP(mat_ZZ& B, double delta, long deep, LLLCheckFct check, 
913           long verb)
914{
915   verbose = verb;
916   NumSwaps = 0;
917   if (verbose) {
918      StartTime = GetTime();
919      LastTime = StartTime;
920   }
921
922   if (delta < 0.50 || delta >= 1) Error("G_LLL_FP: bad delta");
923   if (deep < 0) Error("G_LLL_FP: bad deep");
924   return G_LLL_FP(B, 0, delta, deep, check);
925}
926
927long G_LLL_FP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, 
928           LLLCheckFct check, long verb)
929{
930   verbose = verb;
931   NumSwaps = 0;
932   if (verbose) {
933      StartTime = GetTime();
934      LastTime = StartTime;
935   }
936
937   if (delta < 0.50 || delta >= 1) Error("G_LLL_FP: bad delta");
938   if (deep < 0) Error("G_LLL_FP: bad deep");
939   return G_LLL_FP(B, &U, delta, deep, check);
940}
941
942
943
944static vec_double G_BKZConstant;
945
946static
947void ComputeG_BKZConstant(long beta, long p)
948{
949   const double c_PI = 3.14159265358979323846264338328;
950   const double LogPI = 1.14472988584940017414342735135;
951
952   G_BKZConstant.SetLength(beta-1);
953
954   vec_double Log;
955   Log.SetLength(beta);
956
957
958   long i, j, k;
959   double x, y;
960
961   for (j = 1; j <= beta; j++)
962      Log(j) = log(double(j));
963
964   for (i = 1; i <= beta-1; i++) {
965      // First, we compute x = gamma(i/2)^{2/i}
966
967      k = i/2;
968
969      if ((i & 1) == 0) { // i even
970         x = 0;
971         for (j = 1; j <= k; j++)
972            x = x + Log(j);
973         
974         x = x * (1/double(k));
975
976         x = exp(x);
977      }
978      else { // i odd
979         x = 0;
980         for (j = k + 2; j <= 2*k + 2; j++)
981            x = x + Log(j);
982
983         x = 0.5*LogPI + x - 2*(k+1)*Log(2);
984
985         x = x * (2.0/double(i));
986
987         x = exp(x);
988      }
989
990      // Second, we compute y = 2^{2*p/i}
991
992      y = -(2*p/double(i))*Log(2);
993      y = exp(y);
994
995      G_BKZConstant(i) = x*y/c_PI;
996   }
997}
998
999static vec_double G_BKZThresh;
1000
1001static 
1002void ComputeG_BKZThresh(double *c, long beta)
1003{
1004   G_BKZThresh.SetLength(beta-1);
1005
1006   long i;
1007   double x;
1008
1009   x = 0;
1010
1011   for (i = 1; i <= beta-1; i++) {
1012      x += log(c[i-1]);
1013      G_BKZThresh(i) = exp(x/double(i))*G_BKZConstant(i);
1014      if (!IsFinite(&G_BKZThresh(i))) G_BKZThresh(i) = 0;
1015   }
1016}
1017
1018static
1019long G_BKZ_FP(mat_ZZ& BB, mat_ZZ* UU, double delta, 
1020         long beta, long prune, LLLCheckFct check)
1021{
1022
1023   
1024
1025   long m = BB.NumRows();
1026   long n = BB.NumCols();
1027   long m_orig = m;
1028   
1029   long i, j;
1030   ZZ MU;
1031
1032   double t1;
1033   ZZ T1;
1034   double *tp;
1035
1036   init_red_fudge();
1037
1038   mat_ZZ B;
1039   B = BB;
1040
1041   B.SetDims(m+1, n);
1042
1043
1044   double **B1;  // approximates B
1045
1046   typedef double *doubleptr;
1047
1048   B1 = NTL_NEW_OP doubleptr[m+2];
1049   if (!B1) Error("G_BKZ_FP: out of memory");
1050
1051   for (i = 1; i <= m+1; i++) {
1052      B1[i] = NTL_NEW_OP double[n+1];
1053      if (!B1[i]) Error("G_BKZ_FP: out of memory");
1054   }
1055
1056   double **mu;
1057   mu = NTL_NEW_OP doubleptr[m+2];
1058   if (!mu) Error("G_LLL_FP: out of memory");
1059
1060   for (i = 1; i <= m+1; i++) {
1061      mu[i] = NTL_NEW_OP double[n+2];
1062      if (!mu[i]) Error("G_BKZ_FP: out of memory");
1063   }
1064
1065   double **aux;
1066   aux = NTL_NEW_OP doubleptr[m+2];
1067   if (!aux) Error("G_LLL_FP: out of memory");
1068
1069   for (i = 1; i <= m+1; i++) {
1070      aux[i] = NTL_NEW_OP double[n+1];
1071      if (!aux[i]) Error("G_BKZ_FP: out of memory");
1072   }
1073
1074
1075   double *c; // squared lengths of Gramm-Schmidt basis vectors
1076
1077   c = NTL_NEW_OP double[m+2];
1078   if (!c) Error("G_BKZ_FP: out of memory");
1079
1080   double cbar;
1081
1082   double *ctilda;
1083   ctilda = NTL_NEW_OP double[m+2];
1084   if (!ctilda) Error("G_BKZ_FP: out of memory");
1085
1086   double *vvec;
1087   vvec = NTL_NEW_OP double[m+2];
1088   if (!vvec) Error("G_BKZ_FP: out of memory");
1089
1090   double *yvec;
1091   yvec = NTL_NEW_OP double[m+2];
1092   if (!yvec) Error("G_BKZ_FP: out of memory");
1093
1094   double *uvec;
1095   uvec = NTL_NEW_OP double[m+2];
1096   if (!uvec) Error("G_BKZ_FP: out of memory");
1097
1098   double *utildavec;
1099   utildavec = NTL_NEW_OP double[m+2];
1100   if (!utildavec) Error("G_BKZ_FP: out of memory");
1101
1102
1103   long *Deltavec;
1104   Deltavec = NTL_NEW_OP long[m+2];
1105   if (!Deltavec) Error("G_BKZ_FP: out of memory");
1106
1107   long *deltavec;
1108   deltavec = NTL_NEW_OP long[m+2];
1109   if (!deltavec) Error("G_BKZ_FP: out of memory");
1110
1111   mat_ZZ Ulocal;
1112   mat_ZZ *U;
1113
1114   if (UU) {
1115      Ulocal.SetDims(m+1, m);
1116      for (i = 1; i <= m; i++)
1117         conv(Ulocal(i, i), 1);
1118      U = &Ulocal;
1119   }
1120   else
1121      U = 0;
1122
1123   long quit;
1124   long new_m;
1125   long z, jj, kk;
1126   long s, t;
1127   long h;
1128   double eta;
1129
1130
1131   for (i = 1; i <=m; i++)
1132      for (j = 1; j <= n; j++) {
1133         conv(B1[i][j], B(i, j));
1134         CheckFinite(&B1[i][j]);
1135      }
1136
1137         
1138   GivensCache_FP cache(m, n);
1139
1140   m = ll_G_LLL_FP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1141
1142   double tt;
1143
1144   double enum_time = 0;
1145   unsigned long NumIterations = 0;
1146   unsigned long NumTrivial = 0;
1147   unsigned long NumNonTrivial = 0;
1148   unsigned long NumNoOps = 0;
1149
1150   long verb = verbose;
1151
1152   verbose = 0;
1153
1154   long clean = 1;
1155
1156   if (m < m_orig) {
1157      for (i = m_orig+1; i >= m+2; i--) {
1158         // swap i, i-1
1159
1160         swap(B(i), B(i-1));
1161         if (U) swap((*U)(i), (*U)(i-1));
1162      }
1163   }
1164
1165   if (!quit && m > 1) {
1166      if (beta > m) beta = m;
1167
1168      if (prune > 0) 
1169         ComputeG_BKZConstant(beta, prune);
1170
1171      z = 0;
1172      jj = 0;
1173   
1174      while (z < m-1) {
1175         jj++;
1176         kk = min(jj+beta-1, m);
1177   
1178         if (jj == m) {
1179            jj = 1;
1180            kk = beta;
1181            clean = 1;
1182         }
1183
1184         // ENUM
1185
1186         double tt1;
1187
1188         if (verb) {
1189            tt1 = GetTime();
1190         }
1191
1192         for (i = jj; i <= kk; i++) {
1193            c[i] = mu[i][i]*mu[i][i];
1194            CheckFinite(&c[i]);
1195         }
1196
1197         if (prune > 0)
1198            ComputeG_BKZThresh(&c[jj], kk-jj+1);
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                  yvec[t] = t1;
1237                  t1 = -t1;
1238                  if (t1 >= 0)
1239                     t1 = ceil(t1-0.5);
1240                  else
1241                     t1 = floor(t1+0.5);
1242                  utildavec[t] = vvec[t] = t1;
1243                  Deltavec[t] = 0;
1244                  if (utildavec[t] > -yvec[t]) 
1245                     deltavec[t] = -1;
1246                  else
1247                     deltavec[t] = 1;
1248               }
1249               else {
1250                  cbar = ctilda[jj];
1251                  for (i = jj; i <= kk; i++) {
1252                     uvec[i] = utildavec[i];
1253                  }
1254               }
1255            }
1256            else {
1257               t++;
1258               s = max(s, t);
1259               if (t < s) Deltavec[t] = -Deltavec[t];
1260               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1261               utildavec[t] = vvec[t] + Deltavec[t];
1262            }
1263         }
1264
1265         if (verb) {
1266            tt1 = GetTime() - tt1;
1267            enum_time += tt1;
1268         }
1269         
1270         NumIterations++;
1271   
1272         h = min(kk+1, m);
1273   
1274         if ((delta - 8*red_fudge)*c[jj] > cbar) {
1275
1276            clean = 0;
1277
1278            // we treat the case that the new vector is b_s (jj < s <= kk)
1279            // as a special case that appears to occur most of the time.
1280   
1281            s = 0;
1282            for (i = jj+1; i <= kk; i++) {
1283               if (uvec[i] != 0) {
1284                  if (s == 0)
1285                     s = i;
1286                  else
1287                     s = -1;
1288               }
1289            }
1290   
1291            if (s == 0) Error("G_BKZ_FP: internal error");
1292   
1293            if (s > 0) {
1294               // special case
1295
1296               NumTrivial++;
1297   
1298               for (i = s; i > jj; i--) {
1299                  // swap i, i-1
1300                  swap(B(i-1), B(i));
1301                  if (U) swap((*U)(i-1), (*U)(i));
1302                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1303               }
1304   
1305               // cerr << "special case\n";
1306               new_m = ll_G_LLL_FP(B, U, delta, 0, check, 
1307                                B1, mu, aux, h, jj, quit, cache);
1308               if (new_m != h) Error("G_BKZ_FP: internal error");
1309               if (quit) break;
1310            }
1311            else {
1312               // the general case
1313
1314               NumNonTrivial++;
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_FP: internal error");
1343     
1344               // remove linear dependencies
1345   
1346               // cerr << "general case\n";
1347               new_m = ll_G_LLL_FP(B, U, delta, 0, 0, B1, mu, aux, 
1348                                  kk+1, jj, quit, cache);
1349             
1350               if (new_m != kk) Error("G_BKZ_FP: 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_FP(B, U, delta, 0, check, 
1376                                   B1, mu, aux, h, h, quit, cache);
1377   
1378                  if (new_m != h) Error("G_BKZ_FP: internal error");
1379                  if (quit) break;
1380               }
1381            }
1382   
1383            z = 0;
1384         }
1385         else {
1386            // G_LLL_FP
1387            // cerr << "progress\n";
1388
1389            NumNoOps++;
1390
1391            if (!clean) {
1392               new_m = ll_G_LLL_FP(B, U, delta, 0, check, B1, mu, aux, 
1393                                   h, h, quit, cache);
1394               if (new_m != h) Error("G_BKZ_FP: internal error");
1395               if (quit) break;
1396            }
1397   
1398            z++;
1399         }
1400      }
1401   }
1402
1403
1404   // clean up
1405
1406
1407   if (m_orig > m) {
1408      // for consistency, we move zero vectors to the front
1409
1410      for (i = m+1; i <= m_orig; i++) {
1411         swap(B(i), B(i+1));
1412         if (U) swap((*U)(i), (*U)(i+1));
1413      }
1414
1415      for (i = 0; i < m; i++) {
1416         swap(B(m_orig-i), B(m-i));
1417         if (U) swap((*U)(m_orig-i), (*U)(m-i));
1418      }
1419   }
1420
1421   B.SetDims(m_orig, n);
1422   BB = B;
1423
1424   if (U) {
1425      U->SetDims(m_orig, m_orig);
1426      *UU = *U;
1427   }
1428
1429   for (i = 1; i <= m+1; i++) {
1430      delete [] B1[i];
1431   }
1432
1433   delete [] B1;
1434
1435   for (i = 1; i <= m+1; i++) {
1436      delete [] mu[i];
1437   }
1438
1439   delete [] mu;
1440
1441   for (i = 1; i <= m+1; i++) {
1442      delete [] aux[i];
1443   }
1444
1445   delete [] aux;
1446
1447   delete [] c;
1448   delete [] ctilda;
1449   delete [] vvec;
1450   delete [] yvec;
1451   delete [] uvec;
1452   delete [] utildavec;
1453   delete [] Deltavec;
1454   delete [] deltavec;
1455
1456   return m;
1457}
1458
1459long G_BKZ_FP(mat_ZZ& BB, mat_ZZ& UU, double delta, 
1460         long beta, long prune, LLLCheckFct check, long verb)
1461{
1462   verbose = verb;
1463   NumSwaps = 0;
1464   if (verbose) {
1465      StartTime = GetTime();
1466      LastTime = StartTime;
1467   }
1468
1469   if (delta < 0.50 || delta >= 1) Error("G_BKZ_FP: bad delta");
1470   if (beta < 2) Error("G_BKZ_FP: bad block size");
1471
1472   return G_BKZ_FP(BB, &UU, delta, beta, prune, check);
1473}
1474
1475long G_BKZ_FP(mat_ZZ& BB, double delta, 
1476         long beta, long prune, LLLCheckFct check, long verb)
1477{
1478   verbose = verb;
1479   NumSwaps = 0;
1480   if (verbose) {
1481      StartTime = GetTime();
1482      LastTime = StartTime;
1483   }
1484
1485   if (delta < 0.50 || delta >= 1) Error("G_BKZ_FP: bad delta");
1486   if (beta < 2) Error("G_BKZ_FP: bad block size");
1487
1488   return G_BKZ_FP(BB, 0, delta, beta, prune, check);
1489}
1490
1491NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.