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