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

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