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

spielwiese
Last change on this file since 33a041 was 311902, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: 5.3.2-C++-fix git-svn-id: file:///usr/local/Singular/svn/trunk@7474 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 31.2 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 G_LLLStatus(long max_k, double t, long m, const mat_ZZ& B)
580{
581   ZZ t1;
582   long i;
583   double prodlen = 0;
584
585   for (i = 1; i <= m; i++) {
586      InnerProduct(t1, B(i), B(i));
587      if (!IsZero(t1))
588         prodlen += log(t1);
589   }
590
591   LastTime = t;
592   
593}
594
595static void init_red_fudge()
596{
597   long i;
598
599   log_red = long(0.50*NTL_DOUBLE_PRECISION);
600   red_fudge = 1;
601
602   for (i = log_red; i > 0; i--)
603      red_fudge = red_fudge*0.5;
604}
605
606static void inc_red_fudge()
607{
608
609   red_fudge = red_fudge * 2;
610   log_red--;
611
612   
613   if (log_red < 4)
614      Error("G_LLL_FP: too much loss of precision...stop!");
615}
616
617
618#if 0
619
620static void print_mus(double **mu, long k)
621{
622   long i;
623
624   for (i = k-1; i >= 1; i--)
625      cerr << mu[k][i] << " ";
626   cerr << "\n";
627}
628
629#endif
630
631
632
633static
634long ll_G_LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, 
635           LLLCheckFct check, double **B1, double **mu, 
636           double **aux,
637           long m, long init_k, long &quit, GivensCache_FP& cache)
638{
639   long n = B.NumCols();
640
641   long i, j, k, Fc1;
642   ZZ MU;
643   double mu1;
644
645   double t1;
646   ZZ T1;
647   double *tp;
648
649   double half_plus_fudge = 0.5 + red_fudge;
650
651   quit = 0;
652   k = init_k;
653
654   vec_long in_vec_mem;
655   in_vec_mem.SetLength(n+1);
656   long *in_vec = in_vec_mem.elts();
657
658   double *max_b;
659   max_b = NTL_NEW_OP double [m+1];
660   if (!max_b) Error("out of memory in lll_G_LLL_FP");
661
662   for (i = 1; i <= m; i++)
663      max_b[i] = max_abs(B1[i], n);
664
665   long in_float;
666
667   long counter;
668
669   long trigger_index;
670   long small_trigger;
671   long cnt;
672
673
674   long max_k = 0;
675
676   double tt;
677
678   long swap_cnt = 0;
679
680   cache.flush();
681
682   while (k <= m) {
683
684      if (k > max_k) {
685         max_k = k;
686         swap_cnt = 0;
687      }
688
689      if (verbose) {
690         tt = GetTime();
691
692         if (tt > LastTime + LLLStatusInterval)
693            G_LLLStatus(max_k, tt, m, B);
694      }
695
696      GivensComputeGS(B1, mu, aux, k, n, cache);
697
698      if (swap_cnt > 200000) {
699         //cerr << "G_LLL_FP: swap loop?\n";
700         swap_cnt = 0;
701      }
702
703      counter = 0;
704      trigger_index = k;
705      small_trigger = 0;
706      cnt = 0;
707
708      long sz=0, new_sz;
709
710      do {
711         // size reduction
712
713         counter++;
714         if ((counter & 127) == 0) {
715
716            new_sz = 0;
717            for (j = 1; j <= n; j++)
718               new_sz += NumBits(B(k,j));
719
720            if ((counter >> 7) == 1 || new_sz < sz) {
721               sz = new_sz;
722            }
723            //else {
724            //   cerr << "G_LLL_FP: warning--infinite loop? (" << k << ")\n";
725            //}
726         }
727
728         Fc1 = 0;
729   
730         for (j = k-1; j >= 1; j--) {
731            t1 = fabs(mu[k][j]);
732            if (t1 > half_plus_fudge) { 
733
734
735               if (!Fc1) {
736                  if (j > trigger_index || 
737                      (j == trigger_index && small_trigger)) {
738
739                     cnt++;
740
741                     if (cnt > 10) {
742                        inc_red_fudge();
743                        half_plus_fudge = 0.5 + red_fudge;
744                        cnt = 0;
745                     }
746                  }
747
748                  trigger_index = j;
749                  small_trigger = (t1 < 4);
750
751                  Fc1 = 1;
752                  RowTransformStart(B1[k], in_vec, in_float, n);
753               }
754                 
755
756               mu1 = mu[k][j];
757               if (mu1 >= 0)
758                  mu1 = ceil(mu1-0.5);
759               else
760                  mu1 = floor(mu1+0.5);
761   
762               double *mu_k = mu[k];
763               double *mu_j = mu[j];
764   
765               if (mu1 == 1) {
766                  for (i = 1; i <= j-1; i++)
767                     mu_k[i] -= mu_j[i];
768               }
769               else if (mu1 == -1) {
770                  for (i = 1; i <= j-1; i++)
771                     mu_k[i] += mu_j[i];
772               }
773               else {
774                  for (i = 1; i <= j-1; i++)
775                     mu_k[i] -= mu1*mu_j[i];
776               }
777   
778               mu_k[j] -= mu1;
779   
780               conv(MU, mu1);
781
782               RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec,
783                            max_b[k], max_b[j], in_float);
784               if (U) RowTransform((*U)(k), (*U)(j), MU);
785            }
786         }
787
788
789         if (Fc1) {
790            RowTransformFinish(B(k), B1[k], in_vec);
791            max_b[k] = max_abs(B1[k], n);
792            cache.touch();
793            GivensComputeGS(B1, mu, aux, k, n, cache);
794         }
795      } while (Fc1);
796
797      if (check && (*check)(B(k))) 
798         quit = 1;
799
800      if (IsZero(B(k))) {
801         for (i = k; i < m; i++) {
802            // swap i, i+1
803            swap(B(i), B(i+1));
804            tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp;
805            t1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = t1;
806            if (U) swap((*U)(i), (*U)(i+1));
807         }
808
809         cache.flush();
810
811         m--;
812         if (quit) break;
813         continue;
814      }
815
816      if (quit) break;
817
818      if (deep > 0) {
819         // deep insertions
820
821         Error("sorry...deep insertions not implemented");
822      } // end deep insertions
823
824      // test G_LLL reduction condition
825
826      if (k > 1 && 
827         sqrt(delta - mu[k][k-1]*mu[k][k-1])*fabs(mu[k-1][k-1]) > 
828         fabs(mu[k][k])) {
829         // swap rows k, k-1
830
831         swap(B(k), B(k-1));
832         tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp;
833         t1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = t1;
834         if (U) swap((*U)(k), (*U)(k-1));
835
836         cache.swap();
837
838         k--;
839         NumSwaps++;
840         swap_cnt++;
841         // cout << "-\n";
842      }
843      else {
844
845         cache.incr();
846
847         k++;
848         // cout << "+\n";
849      }
850
851   }
852
853   if (verbose) {
854      G_LLLStatus(m+1, GetTime(), m, B);
855   }
856
857
858   delete [] max_b;
859
860   return m;
861}
862
863
864
865
866
867static
868long G_LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, 
869           LLLCheckFct check)
870{
871   long m = B.NumRows();
872   long n = B.NumCols();
873
874   long i, j;
875   long new_m, dep, quit;
876   ZZ MU;
877
878   ZZ T1;
879
880   init_red_fudge();
881
882   if (U) ident(*U, m);
883
884   double **B1;  // approximates B
885
886   typedef double *doubleptr;
887
888   B1 = NTL_NEW_OP doubleptr[m+1];
889   if (!B1) Error("G_LLL_FP: out of memory");
890
891   for (i = 1; i <= m; i++) {
892      B1[i] = NTL_NEW_OP double[n+1];
893      if (!B1[i]) Error("G_LLL_FP: out of memory");
894   }
895
896   double **mu;
897   mu = NTL_NEW_OP doubleptr[m+1];
898   if (!mu) Error("G_LLL_FP: out of memory");
899
900   for (i = 1; i <= m; i++) {
901      mu[i] = NTL_NEW_OP double[n+2];
902      if (!mu[i]) Error("G_LLL_FP: out of memory");
903   }
904
905   double **aux;
906   aux = NTL_NEW_OP doubleptr[m+1];
907   if (!aux) Error("G_LLL_FP: out of memory");
908
909   for (i = 1; i <= m; i++) {
910      aux[i] = NTL_NEW_OP double[n+1];
911      if (!aux[i]) Error("G_LLL_FP: out of memory");
912   }
913
914   for (i = 1; i <=m; i++)
915      for (j = 1; j <= n; j++) {
916         conv(B1[i][j], B(i, j));
917         CheckFinite(&B1[i][j]);
918      }
919
920         
921   GivensCache_FP cache(m, n);
922
923   new_m = ll_G_LLL_FP(B, U, delta, deep, check, B1, mu, aux, m, 1, quit, cache);
924   dep = m - new_m;
925   m = new_m;
926
927   if (dep > 0) {
928      // for consistency, we move all of the zero rows to the front
929
930      for (i = 0; i < m; i++) {
931         swap(B(m+dep-i), B(m-i));
932         if (U) swap((*U)(m+dep-i), (*U)(m-i));
933      }
934   }
935
936
937   // clean-up
938
939   for (i = 1; i <= m; i++) {
940      delete [] B1[i];
941   }
942
943   delete [] B1;
944
945   for (i = 1; i <= m; i++) {
946      delete [] mu[i];
947   }
948
949   delete [] mu;
950
951   for (i = 1; i <= m; i++) {
952      delete [] aux[i];
953   }
954
955   delete [] aux;
956
957   return m;
958}
959
960         
961
962long G_LLL_FP(mat_ZZ& B, double delta, long deep, LLLCheckFct check, 
963           long verb)
964{
965   verbose = verb;
966   NumSwaps = 0;
967   if (verbose) {
968      StartTime = GetTime();
969      LastTime = StartTime;
970   }
971
972   if (delta < 0.50 || delta >= 1) Error("G_LLL_FP: bad delta");
973   if (deep < 0) Error("G_LLL_FP: bad deep");
974   return G_LLL_FP(B, 0, delta, deep, check);
975}
976
977long G_LLL_FP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, 
978           LLLCheckFct check, long verb)
979{
980   verbose = verb;
981   NumSwaps = 0;
982   if (verbose) {
983      StartTime = GetTime();
984      LastTime = StartTime;
985   }
986
987   if (delta < 0.50 || delta >= 1) Error("G_LLL_FP: bad delta");
988   if (deep < 0) Error("G_LLL_FP: bad deep");
989   return G_LLL_FP(B, &U, delta, deep, check);
990}
991
992
993
994static vec_double G_BKZConstant;
995
996static
997void ComputeG_BKZConstant(long beta, long p)
998{
999   const double c_PI = 3.14159265358979323846264338328;
1000   const double LogPI = 1.14472988584940017414342735135;
1001
1002   G_BKZConstant.SetLength(beta-1);
1003
1004   vec_double Log;
1005   Log.SetLength(beta);
1006
1007
1008   long i, j, k;
1009   double x, y;
1010
1011   for (j = 1; j <= beta; j++)
1012      Log(j) = log(double(j));
1013
1014   for (i = 1; i <= beta-1; i++) {
1015      // First, we compute x = gamma(i/2)^{2/i}
1016
1017      k = i/2;
1018
1019      if ((i & 1) == 0) { // i even
1020         x = 0;
1021         for (j = 1; j <= k; j++)
1022            x = x + Log(j);
1023         
1024         x = x * (1/double(k));
1025
1026         x = exp(x);
1027      }
1028      else { // i odd
1029         x = 0;
1030         for (j = k + 2; j <= 2*k + 2; j++)
1031            x = x + Log(j);
1032
1033         x = 0.5*LogPI + x - 2*(k+1)*Log(2);
1034
1035         x = x * (2.0/double(i));
1036
1037         x = exp(x);
1038      }
1039
1040      // Second, we compute y = 2^{2*p/i}
1041
1042      y = -(2*p/double(i))*Log(2);
1043      y = exp(y);
1044
1045      G_BKZConstant(i) = x*y/c_PI;
1046   }
1047}
1048
1049static vec_double G_BKZThresh;
1050
1051static 
1052void ComputeG_BKZThresh(double *c, long beta)
1053{
1054   G_BKZThresh.SetLength(beta-1);
1055
1056   long i;
1057   double x;
1058
1059   x = 0;
1060
1061   for (i = 1; i <= beta-1; i++) {
1062      x += log(c[i-1]);
1063      G_BKZThresh(i) = exp(x/double(i))*G_BKZConstant(i);
1064      if (!IsFinite(&G_BKZThresh(i))) G_BKZThresh(i) = 0;
1065   }
1066}
1067
1068static 
1069void G_BKZStatus(double tt, double enum_time, unsigned long NumIterations, 
1070               unsigned long NumTrivial, unsigned long NumNonTrivial, 
1071               unsigned long NumNoOps, long m, 
1072               const mat_ZZ& B)
1073{
1074
1075   ZZ t1;
1076   long i;
1077   double prodlen = 0;
1078
1079   for (i = 1; i <= m; i++) {
1080      InnerProduct(t1, B(i), B(i));
1081      if (!IsZero(t1))
1082         prodlen += log(t1);
1083   }
1084
1085   LastTime = tt;
1086   
1087}
1088
1089
1090
1091static
1092long G_BKZ_FP(mat_ZZ& BB, mat_ZZ* UU, double delta, 
1093         long beta, long prune, LLLCheckFct check)
1094{
1095
1096   
1097
1098   long m = BB.NumRows();
1099   long n = BB.NumCols();
1100   long m_orig = m;
1101   
1102   long i, j;
1103   ZZ MU;
1104
1105   double t1;
1106   ZZ T1;
1107   double *tp;
1108
1109   init_red_fudge();
1110
1111   mat_ZZ B;
1112   B = BB;
1113
1114   B.SetDims(m+1, n);
1115
1116
1117   double **B1;  // approximates B
1118
1119   typedef double *doubleptr;
1120
1121   B1 = NTL_NEW_OP doubleptr[m+2];
1122   if (!B1) Error("G_BKZ_FP: out of memory");
1123
1124   for (i = 1; i <= m+1; i++) {
1125      B1[i] = NTL_NEW_OP double[n+1];
1126      if (!B1[i]) Error("G_BKZ_FP: out of memory");
1127   }
1128
1129   double **mu;
1130   mu = NTL_NEW_OP doubleptr[m+2];
1131   if (!mu) Error("G_LLL_FP: out of memory");
1132
1133   for (i = 1; i <= m+1; i++) {
1134      mu[i] = NTL_NEW_OP double[n+2];
1135      if (!mu[i]) Error("G_BKZ_FP: out of memory");
1136   }
1137
1138   double **aux;
1139   aux = NTL_NEW_OP doubleptr[m+2];
1140   if (!aux) Error("G_LLL_FP: out of memory");
1141
1142   for (i = 1; i <= m+1; i++) {
1143      aux[i] = NTL_NEW_OP double[n+1];
1144      if (!aux[i]) Error("G_BKZ_FP: out of memory");
1145   }
1146
1147
1148   double *c; // squared lengths of Gramm-Schmidt basis vectors
1149
1150   c = NTL_NEW_OP double[m+2];
1151   if (!c) Error("G_BKZ_FP: out of memory");
1152
1153   double cbar;
1154
1155   double *ctilda;
1156   ctilda = NTL_NEW_OP double[m+2];
1157   if (!ctilda) Error("G_BKZ_FP: out of memory");
1158
1159   double *vvec;
1160   vvec = NTL_NEW_OP double[m+2];
1161   if (!vvec) Error("G_BKZ_FP: out of memory");
1162
1163   double *yvec;
1164   yvec = NTL_NEW_OP double[m+2];
1165   if (!yvec) Error("G_BKZ_FP: out of memory");
1166
1167   double *uvec;
1168   uvec = NTL_NEW_OP double[m+2];
1169   if (!uvec) Error("G_BKZ_FP: out of memory");
1170
1171   double *utildavec;
1172   utildavec = NTL_NEW_OP double[m+2];
1173   if (!utildavec) Error("G_BKZ_FP: out of memory");
1174
1175
1176   long *Deltavec;
1177   Deltavec = NTL_NEW_OP long[m+2];
1178   if (!Deltavec) Error("G_BKZ_FP: out of memory");
1179
1180   long *deltavec;
1181   deltavec = NTL_NEW_OP long[m+2];
1182   if (!deltavec) Error("G_BKZ_FP: out of memory");
1183
1184   mat_ZZ Ulocal;
1185   mat_ZZ *U;
1186
1187   if (UU) {
1188      Ulocal.SetDims(m+1, m);
1189      for (i = 1; i <= m; i++)
1190         conv(Ulocal(i, i), 1);
1191      U = &Ulocal;
1192   }
1193   else
1194      U = 0;
1195
1196   long quit;
1197   long new_m;
1198   long z, jj, kk;
1199   long s, t;
1200   long h;
1201   double eta;
1202
1203
1204   for (i = 1; i <=m; i++)
1205      for (j = 1; j <= n; j++) {
1206         conv(B1[i][j], B(i, j));
1207         CheckFinite(&B1[i][j]);
1208      }
1209
1210         
1211   GivensCache_FP cache(m, n);
1212
1213   m = ll_G_LLL_FP(B, U, delta, 0, check, B1, mu, aux, m, 1, quit, cache);
1214
1215   double tt;
1216
1217   double enum_time = 0;
1218   unsigned long NumIterations = 0;
1219   unsigned long NumTrivial = 0;
1220   unsigned long NumNonTrivial = 0;
1221   unsigned long NumNoOps = 0;
1222
1223   long verb = verbose;
1224
1225   verbose = 0;
1226
1227   long clean = 1;
1228
1229   if (m < m_orig) {
1230      for (i = m_orig+1; i >= m+2; i--) {
1231         // swap i, i-1
1232
1233         swap(B(i), B(i-1));
1234         if (U) swap((*U)(i), (*U)(i-1));
1235      }
1236   }
1237
1238   if (!quit && m > 1) {
1239      if (beta > m) beta = m;
1240
1241      if (prune > 0) 
1242         ComputeG_BKZConstant(beta, prune);
1243
1244      z = 0;
1245      jj = 0;
1246   
1247      while (z < m-1) {
1248         jj++;
1249         kk = min(jj+beta-1, m);
1250   
1251         if (jj == m) {
1252            jj = 1;
1253            kk = beta;
1254            clean = 1;
1255         }
1256
1257         if (verb) {
1258            tt = GetTime();
1259            if (tt > LastTime + LLLStatusInterval)
1260               G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1261                         NumNonTrivial, NumNoOps, m, B);
1262         }
1263
1264   
1265         // ENUM
1266
1267         double tt1;
1268
1269         if (verb) {
1270            tt1 = GetTime();
1271         }
1272
1273         for (i = jj; i <= kk; i++) {
1274            c[i] = mu[i][i]*mu[i][i];
1275            CheckFinite(&c[i]);
1276         }
1277
1278         if (prune > 0)
1279            ComputeG_BKZThresh(&c[jj], kk-jj+1);
1280   
1281         cbar = c[jj];
1282         utildavec[jj] = uvec[jj] = 1;
1283   
1284         yvec[jj] = vvec[jj] = 0;
1285         Deltavec[jj] = 0;
1286   
1287   
1288         s = t = jj;
1289         deltavec[jj] = 1;
1290   
1291         for (i = jj+1; i <= kk+1; i++) {
1292            ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0;
1293            Deltavec[i] = 0;
1294            vvec[i] = 0;
1295            deltavec[i] = 1;
1296         }
1297
1298         long enum_cnt = 0;
1299   
1300         while (t <= kk) {
1301            if (verb) {
1302               enum_cnt++;
1303               if (enum_cnt > 100000) {
1304                  enum_cnt = 0;
1305                  tt = GetTime();
1306                  if (tt > LastTime + LLLStatusInterval) {
1307                     enum_time += tt - tt1;
1308                     tt1 = tt;
1309                     G_BKZStatus(tt, enum_time, NumIterations, NumTrivial,
1310                               NumNonTrivial, NumNoOps, m, B);
1311                  }
1312               }
1313            }
1314
1315            ctilda[t] = ctilda[t+1] + 
1316               (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*c[t];
1317   
1318            if (prune > 0 && t > jj) {
1319               eta = G_BKZThresh(t-jj);
1320            }
1321            else
1322               eta = 0;
1323   
1324            if (ctilda[t] < cbar - eta) {
1325               if (t > jj) {
1326                  t--;
1327                  t1 = 0;
1328                  for (i = t+1; i <= s; i++)
1329                     t1 += utildavec[i]*mu[i][t];
1330                  yvec[t] = t1;
1331                  t1 = -t1;
1332                  if (t1 >= 0)
1333                     t1 = ceil(t1-0.5);
1334                  else
1335                     t1 = floor(t1+0.5);
1336                  utildavec[t] = vvec[t] = t1;
1337                  Deltavec[t] = 0;
1338                  if (utildavec[t] > -yvec[t]) 
1339                     deltavec[t] = -1;
1340                  else
1341                     deltavec[t] = 1;
1342               }
1343               else {
1344                  cbar = ctilda[jj];
1345                  for (i = jj; i <= kk; i++) {
1346                     uvec[i] = utildavec[i];
1347                  }
1348               }
1349            }
1350            else {
1351               t++;
1352               s = max(s, t);
1353               if (t < s) Deltavec[t] = -Deltavec[t];
1354               if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t];
1355               utildavec[t] = vvec[t] + Deltavec[t];
1356            }
1357         }
1358
1359         if (verb) {
1360            tt1 = GetTime() - tt1;
1361            enum_time += tt1;
1362         }
1363         
1364         NumIterations++;
1365   
1366         h = min(kk+1, m);
1367   
1368         if ((delta - 8*red_fudge)*c[jj] > cbar) {
1369
1370            clean = 0;
1371
1372            // we treat the case that the new vector is b_s (jj < s <= kk)
1373            // as a special case that appears to occur most of the time.
1374   
1375            s = 0;
1376            for (i = jj+1; i <= kk; i++) {
1377               if (uvec[i] != 0) {
1378                  if (s == 0)
1379                     s = i;
1380                  else
1381                     s = -1;
1382               }
1383            }
1384   
1385            if (s == 0) Error("G_BKZ_FP: internal error");
1386   
1387            if (s > 0) {
1388               // special case
1389
1390               NumTrivial++;
1391   
1392               for (i = s; i > jj; i--) {
1393                  // swap i, i-1
1394                  swap(B(i-1), B(i));
1395                  if (U) swap((*U)(i-1), (*U)(i));
1396                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1397               }
1398   
1399               // cerr << "special case\n";
1400               new_m = ll_G_LLL_FP(B, U, delta, 0, check, 
1401                                B1, mu, aux, h, jj, quit, cache);
1402               if (new_m != h) Error("G_BKZ_FP: internal error");
1403               if (quit) break;
1404            }
1405            else {
1406               // the general case
1407
1408               NumNonTrivial++;
1409   
1410               for (i = 1; i <= n; i++) conv(B(m+1, i), 0);
1411
1412               if (U) {
1413                  for (i = 1; i <= m_orig; i++)
1414                     conv((*U)(m+1, i), 0);
1415               }
1416
1417               for (i = jj; i <= kk; i++) {
1418                  if (uvec[i] == 0) continue;
1419                  conv(MU, uvec[i]);
1420                  RowTransform2(B(m+1), B(i), MU);
1421                  if (U) RowTransform2((*U)(m+1), (*U)(i), MU);
1422               }
1423     
1424               for (i = m+1; i >= jj+1; i--) {
1425                  // swap i, i-1
1426                  swap(B(i-1), B(i));
1427                  if (U) swap((*U)(i-1), (*U)(i));
1428                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1429               }
1430     
1431               for (i = 1; i <= n; i++) {
1432                  conv(B1[jj][i], B(jj, i));
1433                  CheckFinite(&B1[jj][i]);
1434               }
1435
1436               if (IsZero(B(jj))) Error("G_BKZ_FP: internal error");
1437     
1438               // remove linear dependencies
1439   
1440               // cerr << "general case\n";
1441               new_m = ll_G_LLL_FP(B, U, delta, 0, 0, B1, mu, aux, 
1442                                  kk+1, jj, quit, cache);
1443             
1444               if (new_m != kk) Error("G_BKZ_FP: internal error"); 
1445
1446               // remove zero vector
1447     
1448               for (i = kk+2; i <= m+1; i++) {
1449                  // swap i, i-1
1450                  swap(B(i-1), B(i));
1451                  if (U) swap((*U)(i-1), (*U)(i));
1452                  tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp;
1453               }
1454     
1455               quit = 0;
1456               if (check) {
1457                  for (i = 1; i <= kk; i++)
1458                     if ((*check)(B(i))) {
1459                        quit = 1;
1460                        break;
1461                     }
1462               }
1463
1464               if (quit) break;
1465   
1466               if (h > kk) {
1467                  // extend reduced basis
1468   
1469                  new_m = ll_G_LLL_FP(B, U, delta, 0, check, 
1470                                   B1, mu, aux, h, h, quit, cache);
1471   
1472                  if (new_m != h) Error("G_BKZ_FP: internal error");
1473                  if (quit) break;
1474               }
1475            }
1476   
1477            z = 0;
1478         }
1479         else {
1480            // G_LLL_FP
1481            // cerr << "progress\n";
1482
1483            NumNoOps++;
1484
1485            if (!clean) {
1486               new_m = ll_G_LLL_FP(B, U, delta, 0, check, B1, mu, aux, 
1487                                   h, h, quit, cache);
1488               if (new_m != h) Error("G_BKZ_FP: internal error");
1489               if (quit) break;
1490            }
1491   
1492            z++;
1493         }
1494      }
1495   }
1496
1497
1498   if (verb) {
1499      G_BKZStatus(GetTime(), enum_time, NumIterations, NumTrivial, NumNonTrivial, 
1500                NumNoOps, m, B);
1501   }
1502
1503   // clean up
1504
1505
1506   if (m_orig > m) {
1507      // for consistency, we move zero vectors to the front
1508
1509      for (i = m+1; i <= m_orig; i++) {
1510         swap(B(i), B(i+1));
1511         if (U) swap((*U)(i), (*U)(i+1));
1512      }
1513
1514      for (i = 0; i < m; i++) {
1515         swap(B(m_orig-i), B(m-i));
1516         if (U) swap((*U)(m_orig-i), (*U)(m-i));
1517      }
1518   }
1519
1520   B.SetDims(m_orig, n);
1521   BB = B;
1522
1523   if (U) {
1524      U->SetDims(m_orig, m_orig);
1525      *UU = *U;
1526   }
1527
1528   for (i = 1; i <= m+1; i++) {
1529      delete [] B1[i];
1530   }
1531
1532   delete [] B1;
1533
1534   for (i = 1; i <= m+1; i++) {
1535      delete [] mu[i];
1536   }
1537
1538   delete [] mu;
1539
1540   for (i = 1; i <= m+1; i++) {
1541      delete [] aux[i];
1542   }
1543
1544   delete [] aux;
1545
1546   delete [] c;
1547   delete [] ctilda;
1548   delete [] vvec;
1549   delete [] yvec;
1550   delete [] uvec;
1551   delete [] utildavec;
1552   delete [] Deltavec;
1553   delete [] deltavec;
1554
1555   return m;
1556}
1557
1558long G_BKZ_FP(mat_ZZ& BB, mat_ZZ& UU, double delta, 
1559         long beta, long prune, LLLCheckFct check, long verb)
1560{
1561   verbose = verb;
1562   NumSwaps = 0;
1563   if (verbose) {
1564      StartTime = GetTime();
1565      LastTime = StartTime;
1566   }
1567
1568   if (delta < 0.50 || delta >= 1) Error("G_BKZ_FP: bad delta");
1569   if (beta < 2) Error("G_BKZ_FP: bad block size");
1570
1571   return G_BKZ_FP(BB, &UU, delta, beta, prune, check);
1572}
1573
1574long G_BKZ_FP(mat_ZZ& BB, double delta, 
1575         long beta, long prune, LLLCheckFct check, long verb)
1576{
1577   verbose = verb;
1578   NumSwaps = 0;
1579   if (verbose) {
1580      StartTime = GetTime();
1581      LastTime = StartTime;
1582   }
1583
1584   if (delta < 0.50 || delta >= 1) Error("G_BKZ_FP: bad delta");
1585   if (beta < 2) Error("G_BKZ_FP: bad block size");
1586
1587   return G_BKZ_FP(BB, 0, delta, beta, prune, check);
1588}
1589
1590NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.