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