source: git/ntl/src/LLL_FP.c @ 08a955

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