source: git/ntl/src/G_LLL_XD.c @ 6ce030f

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