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

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