source: git/IntegerProgramming/LLL.cc @ 4dcfb4f

spielwiese
Last change on this file since 4dcfb4f was 4dcfb4f, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: const Bigint&)0 git-svn-id: file:///usr/local/Singular/svn/trunk@7681 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.0 KB
Line 
1#ifndef LLL_CC
2#define LLL_CC
3
4#include "LLL.h"
5
6// subroutines for the LLL-algorithms
7
8void REDI_KB(const short& k, const short& l, BigInt** b,
9             const short& number_of_vectors, const short& vector_dimension,
10             BigInt** H, BigInt* d, BigInt** lambda)
11// the REDI procedure for relations(...) (to compute the Kernel Basis,
12// algorithm 2.7.2 in Cohen's book)
13{
14#ifdef GMP
15  if(abs(BigInt(2)*lambda[k][l])<=d[l+1])
16#else   // GMP
17  if(labs(2*lambda[k][l])<=d[l+1])
18    // labs is the abs-function for long ints
19#endif  // GMP
20    return;
21
22#ifdef GMP
23  BigInt q=(BigInt(2)*lambda[k][l]+d[l+1])/(BigInt(2)*d[l+1]);
24#else   // GMP
25  long q=(long int) floor(((float)(2*lambda[k][l]+d[l+1]))/(2*d[l+1]));
26#endif  // GMP
27
28  // q is the integer quotient of the division
29  // (2*lambda[k][l]+d[l+1])/(2*d[l+1]).
30  // Because of the rounding mode (always towards zero) of GNU C++,
31  // we cannot use the built-in integer division
32  // here; it causes errors when dealing with negative numbers. Therefore
33  // the complicated casts: The divident is first casted to a float which
34  // causes the division result to be a float. This result is explicitly
35  // rounded downwards. As the floor-function returns a double (for range
36  // reasons), this has to be casted to an integer again.
37
38  for(short n=0;n<number_of_vectors;n++)
39    H[k][n]-=q*H[l][n];
40  // H[k]=H[k]-q*H[l]
41
42  for(short m=0;m<vector_dimension;m++)
43    b[k][m]-=q*b[l][m];
44  // b[k]=b[k]-q*b[l]
45
46  lambda[k][l]-=q*d[l+1];
47
48  for(short i=0;i<=l-1;i++)
49    lambda[k][i]-=q*lambda[l][i];
50}
51
52
53
54
55void REDI_IL(const short& k, const short& l, BigInt** b,
56             const short& vector_dimension, BigInt* d, BigInt** lambda)
57// the REDI procedure for the integer LLL algorithm (algorithm 2.6.7 in
58// Cohen's book)
59{
60#ifdef GMP
61  if(abs(BigInt(2)*lambda[k][l])<=d[l+1])
62#else   // GMP
63  if(labs(2*lambda[k][l])<=d[l+1])
64    // labs is the abs-function for long ints
65#endif  // GMP
66    return;
67
68#ifdef GMP
69  BigInt q=(BigInt(2)*lambda[k][l]+d[l+1])/(BigInt(2)*d[l+1]);
70#else   // GMP
71  long q=(long int) floor(((float)(2*lambda[k][l]+d[l+1]))/(2*d[l+1]));
72#endif  // GMP
73
74  // q is the integer quotient of the division
75  // (2*lambda[k][l]+d[l+1])/(2*d[l+1]).
76  // Because of the rounding mode (always towards zero) of GNU C++,
77  // we cannot use the built-in integer division
78  // here; it causes errors when dealing with negative numbers. Therefore
79  // the complicated casts: The divident is first casted to a float which
80  // causes the division result to be a float. This result is explicitly
81  // rounded downwards. As the floor-function returns a double (for range
82  // reasons), this has to be casted to an integer again.
83
84  for(short m=0;m<vector_dimension;m++)
85    b[k][m]-=q*b[l][m];
86  // b[k]=b[k]-q*b[l]
87
88  lambda[k][l]-=q*d[l+1];
89
90  for(short i=0;i<=l-1;i++)
91    lambda[k][i]-=q*lambda[l][i];
92}
93
94
95
96
97void SWAPI(const short& k, const short& k_max, BigInt** b, BigInt* d,
98           BigInt** lambda)
99// the SWAPI procedure of algorithm 2.6.7
100{
101
102  // exchange b[k] and b[k-1]
103  // This can be done efficiently by swapping pointers (not entries).
104  BigInt* swap=b[k];
105  b[k]=b[k-1];
106  b[k-1]=swap;
107
108  if(k>1)
109    for(short j=0;j<=k-2;j++)
110    {
111      // exchange lambda[k][j] and lambda[k-1][j]
112      BigInt swap=lambda[k][j];
113      lambda[k][j]=lambda[k-1][j];
114      lambda[k-1][j]=swap;
115    }
116
117  BigInt _lambda=lambda[k][k-1];
118
119  BigInt B=(d[k-1]*d[k+1] + _lambda*_lambda)/d[k];
120  // It might be better to choose another evaluation order for this formula,
121  // see below.
122
123  for(short i=k+1;i<=k_max;i++)
124  {
125    BigInt t=lambda[i][k];
126    lambda[i][k]=(d[k+1]*lambda[i][k-1] - _lambda*t)/d[k];
127    lambda[i][k-1]=(B*t + _lambda*lambda[i][k])/d[k+1];
128  }
129
130  d[k]=B;
131}
132
133
134
135
136void SWAPK(const short& k, const short& k_max, BigInt** b, BigInt** H,
137           char *f, BigInt* d, BigInt** lambda)
138// the SWAPK procedure of algorithm 2.7.2
139{
140  // exchange H[k] and H[k-1]
141  // This can be done efficiently by swapping pointers (not entries).
142  BigInt *swap=H[k];
143  H[k]=H[k-1];
144  H[k-1]=swap;
145
146  // exchange b[k] and b[k-1] by the same method
147  swap=b[k];
148  b[k]=b[k-1];
149  b[k-1]=swap;
150
151  if(k>1)
152    for(short j=0;j<=k-2;j++)
153    {
154      // exchange lambda[k][j] and lambda[k-1][j]
155      BigInt swap=lambda[k][j];
156      lambda[k][j]=lambda[k-1][j];
157      lambda[k-1][j]=swap;
158    }
159
160  BigInt _lambda=lambda[k][k-1];
161
162  if(_lambda==BigInt(0))
163  {
164    d[k]=d[k-1];
165    f[k-1]=0;
166    f[k]=1;
167    lambda[k][k-1]=0;
168    for(short i=k+1;i<=k_max;i++)
169    {
170      lambda[i][k]=lambda[i][k-1];
171      lambda[i][k-1]=0;
172    }
173  }
174  else
175    // lambda!=0
176  {
177    for(short i=k+1;i<=k_max;i++)
178      lambda[i][k-1]=(_lambda*lambda[i][k-1])/d[k];
179
180    // Multiplie lambda[i][k-1] by _lambda/d[k].
181    // One could also write
182    //   lambda[i][k-1]*=(lambda/d[k]);   (*)
183    // Without a BigInt class, this can prevent overflows when computing
184    // _lambda*lambda[i][k-1].
185    // But examples show that lambda/d[k] is in general not an integer.
186    // So (*) could lead to problems due to the inexact floating point
187    // arithmetic...
188    // Therefore, we chose the secure evaluation order in all such cases.
189
190    BigInt t=d[k+1];
191    d[k]=(_lambda*_lambda)/d[k];
192    d[k+1]=d[k];
193
194    for(short j=k+1;j<=k_max-1;j++)
195      for(short i=j+1;i<=k_max;i++)
196        lambda[i][j]=(lambda[i][j]*d[k])/t;
197
198    for(short j=k+1;j<=k_max;j++)
199      d[j+1]=(d[j+1]*d[k])/t;
200  }
201
202}
203
204typedef BigInt* BigIntP;
205
206
207short relations(BigInt **b, const short& number_of_vectors,
208                const short& vector_dimension, BigInt**& H)
209{
210
211// first check arguments
212
213  if(number_of_vectors<0)
214  {
215    cerr<<"\nWARNING: short relations(BigInt**, const short&, const short&, "
216      "BigInt**):\nargument number_of_vectors out of range"<<endl;
217    return -1;
218  }
219
220  if(vector_dimension<=0)
221  {
222    cerr<<"\nWARNING: short relations(BigInt**, const short&, const short&, "
223      "BigInt**):\nargument vector_dimension out of range"<<endl;
224    return -1;
225  }
226
227
228// consider special case
229
230  if(number_of_vectors==1)
231    // Only one vector which has no relations if it is not zero,
232    // else relation 1.
233  {
234    short r=1;    // Suppose the only column of the matrix is zero.
235
236    for(short m=0;m<vector_dimension;m++)
237      if(b[0][m]!=BigInt(0))
238        // nonzero entry detected
239        r=0;
240
241    if(r==1)
242    {
243      H=new BigIntP[1];
244      H[0]=new BigInt[1];
245      H[0][0]=1;
246      // This is the lattice basis of the relations...
247    }
248
249    return r;
250  }
251
252
253// memory allocation
254
255// The names are chosen (as far as possible) according to Cohen's book.
256// However, for technical reasons, the indices do not run from 1 to
257// (e.g.) number_of_vectors, but from 0 to number_of_vectors-1.
258// Therefore all indices are shifted by -1 in comparison with this book,
259// except from the indices of the array d which has size
260// number_of_vectors+1.
261
262  H=new BigIntP[number_of_vectors];
263  for(short n=0;n<number_of_vectors;n++)
264    H[n]=new BigInt[number_of_vectors];
265
266  char* f=new char[number_of_vectors];
267
268  BigInt* d=new BigInt[number_of_vectors+1];
269
270  BigInt** lambda=new BigIntP[number_of_vectors];
271  for(short n=1;n<number_of_vectors;n++)
272    lambda[n]=new BigInt[n];
273  // We only need lambda[n][k] for n>k.
274
275
276
277// Step 1: Initialization
278
279  short k=1;
280  short k_max=0;
281  // for iteration
282
283  d[0]=1;
284
285  BigInt t=0;
286  for(short m=0;m<vector_dimension;m++)
287    t+=b[0][m]*b[0][m];
288  // Now, t is the scalar product of b[0] with itself.
289
290  for(short n=0;n<number_of_vectors;n++)
291    for(short l=0;l<number_of_vectors;l++)
292      if(n==l)
293        H[n][l]=1;
294      else
295        H[n][l]=0;
296  // Now, H equals the matrix I_(number_of_vectors).
297
298  if(t!=BigInt(0))
299  {
300    d[1]=t;
301    f[0]=1;
302  }
303  else
304  {
305    d[1]=1;
306    f[0]=0;
307  }
308
309
310
311// The other steps are not programmed with "goto" as in Cohen's book.
312// Instead, we enter a do-while-loop which terminates iff
313// k>=number_of_vectors.
314
315  do
316  {
317
318
319// Step 2: Incremental Gram-Schmidt
320
321    if(k>k_max)
322      // else we can immediately go to step 3.
323    {
324      k_max=k;
325
326      for(short j=0;j<=k;j++)
327        if((f[j]==0) && (j<k))
328          lambda[k][j]=0;
329        else
330        {
331          BigInt u=0;
332
333          // compute scalar product of b[k] and b[j]
334          for(short m=0;m<vector_dimension;m++)
335            u+=b[k][m]*b[j][m];
336
337          for(short i=0;i<=j-1;i++)
338            if(f[i]!=0)
339              u=(d[i+1]*u-lambda[k][i]*lambda[j][i])/d[i];
340
341          if(j<k)
342            lambda[k][j]=u;
343          else
344            // j==k
345            if(u!=BigInt(0))
346            {
347              d[k+1]=u;
348              f[k]=1;
349            }
350            else
351              // u==0
352            {
353              d[k+1]=d[k];
354              f[k]=0;
355            }
356        }
357    }
358
359
360// Step 3: Test f[k]==0 and f[k-1]!=0
361
362    do
363    {
364      if(f[k-1]!=0)
365        REDI_KB(k,k-1,b,number_of_vectors,vector_dimension,H,d,lambda);
366
367      if((f[k-1]!=0) && (f[k]==0))
368      {
369        SWAPK(k,k_max,b,H,f,d,lambda);
370
371        if(k>1)
372          k--;
373        else
374          k=1;
375        // k=max(1,k-1)
376      }
377
378      else
379        break;
380    }
381    while(1);
382
383    // Now the conditions above are no longer satisfied.
384
385    for(short l=k-2;l>=0;l--)
386      if(f[l]!=0)
387        REDI_KB(k,l,b,number_of_vectors,vector_dimension,H,d,lambda);
388    k++;
389
390
391// Step 4: Finished?
392
393  }
394  while(k<number_of_vectors);
395
396
397
398// Now we have computed a lattice basis of the relations of the b[i].
399// Prepare the LLL-reduction.
400
401  // Compute the dimension r of the relations.
402  short r=0;
403  for(short n=0;n<number_of_vectors;n++)
404    if(f[n]==0) // n==r!!
405      r++;
406    else
407      break;
408
409  // Delete the part of H that is no longer needed (especially the vectors
410  // H[r],...,H[number_of_vectors-1]).
411  BigInt **aux=H;
412  if(r>0)
413    H=new BigIntP[r];
414  for(short i=0;i<r;i++)
415    H[i]=aux[i];
416
417  for(short i=r;i<number_of_vectors;i++)
418    delete[] aux[i];
419  delete[] aux;
420
421  delete[] f;
422
423  delete[] d;
424
425  for(short i=1;i<number_of_vectors;i++)
426    delete[] lambda[i];
427  delete[] lambda;
428
429  integral_LLL(H,r,number_of_vectors);
430
431  return r;
432
433}
434
435
436
437
438short integral_LLL(BigInt** b, const short& number_of_vectors,
439                  const short& vector_dimension)
440{
441
442// first check arguments
443
444  if(number_of_vectors<0)
445  {
446    cerr<<"\nWARNING: short integral_LL(BigInt**, const short&, const short&):"
447      "\nargument number_of_vectors out of range"<<endl;
448    return -1;
449  }
450
451  if(vector_dimension<=0)
452  {
453    cerr<<"\nWARNING: short integral_LLL(BigInt**, const short&, const "
454      "short&):\nargument vector_dimension out of range"<<endl;
455    return -1;
456  }
457
458
459// consider special case
460
461  if(number_of_vectors<=1)
462    // 0 or 1 input vector, nothing to be done
463    return 0;
464
465
466// memory allocation
467
468// The names are chosen (as far as possible) according to Cohen's book.
469// However, for technical reasons, the indices do not run from 1 to
470// (e.g.) number_of_vectors, but from 0 to number_of_vectors-1.
471// Therefore all indices are shifted by -1 in comparison with this book,
472// except from the indices of the array d which has size
473// number_of_vectors+1.
474
475  BigInt* d=new BigInt[number_of_vectors+1];
476
477  BigInt** lambda=new BigIntP[number_of_vectors];
478  for(short s=1;s<number_of_vectors;s++)
479    lambda[s]=new BigInt[s];
480  // We only need lambda[n][k] for n>k.
481
482
483
484// Step 1: Initialization
485
486  short k=1;
487  short k_max=0;
488  // for iteration
489  d[0]=1;
490
491  d[1]=0;
492  for(short n=0;n<vector_dimension;n++)
493    d[1]+=b[0][n]*b[0][n];
494  // Now, d[1] is the scalar product of b[0] with itself.
495
496
497// The other steps are not programmed with "goto" as in Cohen's book.
498// Instead, we enter a do-while-loop which terminates iff k>r.
499
500  do
501  {
502
503
504// Step 2: Incremental Gram-Schmidt
505
506    if(k>k_max)
507      // else we can immediately go to step 3.
508    {
509      k_max=k;
510
511      for(short j=0;j<=k;j++)
512      {
513        BigInt u=0;
514        // compute scalar product of b[k] and b[j]
515        for(short n=0;n<vector_dimension;n++)
516          u+=b[k][n]*b[j][n];
517
518        for(short i=0;i<=j-1;i++)
519        {
520          u*=d[i+1];
521          u-=lambda[k][i]*lambda[j][i];
522          u/=d[i];
523
524          //u=(d[i+1]*u-lambda[k][i]*lambda[j][i])/d[i];
525        }
526
527        if(j<k)
528          lambda[k][j]=u;
529        else
530          // j==k
531          d[k+1]=u;
532      }
533
534      if(d[k+1]==BigInt(0))
535      {
536        cerr<<"\nERROR: void integral_LLL(BigInt**, const short&, const "
537          "short&):\ninput vectors must be linearly independent"<<endl;
538        return -1;
539      }
540    }
541
542
543// Step 3: Test LLL-condition
544
545    do
546    {
547      REDI_IL(k,k-1,b,vector_dimension,d,lambda);
548
549      //if(4*d[k+1]*d[k-1] < 3*d[k]*d[k] - lambda[k][k-1]*lambda[k][k-1])
550      if((const BigInt&)4*d[k+1]*d[k-1]
551          < (const BigInt&)3*d[k]*d[k] - lambda[k][k-1]*lambda[k][k-1])
552      {
553        SWAPI(k,k_max,b,d,lambda);
554        if(k>1)
555          k--;
556        // k=max(1,k-1)
557      }
558      else
559        break;
560
561    }
562    while(1);
563
564    // Now the condition above is no longer satisfied.
565
566    for(short l=k-2;l>=0;l--)
567      REDI_IL(k,l,b,vector_dimension,d,lambda);
568    k++;
569
570
571
572// Step 4: Finished?
573
574  }
575  while(k<number_of_vectors);
576
577
578// Now, b contains the LLL-reduced lattice basis.
579// Memory cleanup.
580
581  delete[] d;
582
583  for(short i=1;i<number_of_vectors;i++)
584    delete[] lambda[i];
585  delete[] lambda;
586
587  return 0;
588
589}
590#endif  // LLL_CC
Note: See TracBrowser for help on using the repository browser.