source: git/ntl/src/G_LLL_RR.c @ 447abc

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