[6ba162] | 1 | #ifndef LLL_CC |
---|
| 2 | #define LLL_CC |
---|
| 3 | |
---|
| 4 | #include "LLL.h" |
---|
| 5 | |
---|
| 6 | // subroutines for the LLL-algorithms |
---|
| 7 | |
---|
| 8 | void 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 | |
---|
| 55 | void 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 | |
---|
| 97 | void 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 | |
---|
| 136 | void 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 | |
---|
[4dcfb4f] | 162 | if(_lambda==BigInt(0)) |
---|
[6ba162] | 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 | |
---|
[36f8d9] | 204 | typedef BigInt* BigIntP; |
---|
[6ba162] | 205 | |
---|
| 206 | |
---|
| 207 | short 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++) |
---|
[4dcfb4f] | 237 | if(b[0][m]!=BigInt(0)) |
---|
[6ba162] | 238 | // nonzero entry detected |
---|
| 239 | r=0; |
---|
| 240 | |
---|
| 241 | if(r==1) |
---|
| 242 | { |
---|
[36f8d9] | 243 | H=new BigIntP[1]; |
---|
[6ba162] | 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 | |
---|
[36f8d9] | 262 | H=new BigIntP[number_of_vectors]; |
---|
[6ba162] | 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 | |
---|
[36f8d9] | 270 | BigInt** lambda=new BigIntP[number_of_vectors]; |
---|
[6ba162] | 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 | |
---|
[4dcfb4f] | 298 | if(t!=BigInt(0)) |
---|
[6ba162] | 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 |
---|
[4dcfb4f] | 345 | if(u!=BigInt(0)) |
---|
[6ba162] | 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) |
---|
[36f8d9] | 413 | H=new BigIntP[r]; |
---|
[6ba162] | 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 | |
---|
| 438 | short 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 | |
---|
[36f8d9] | 477 | BigInt** lambda=new BigIntP[number_of_vectors]; |
---|
[6ba162] | 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 | |
---|
[4dcfb4f] | 534 | if(d[k+1]==BigInt(0)) |
---|
[6ba162] | 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 | |
---|
[a7c8b18] | 549 | //if(4*d[k+1]*d[k-1] < 3*d[k]*d[k] - lambda[k][k-1]*lambda[k][k-1]) |
---|
[3f2e22] | 550 | if((BigInt(4))*d[k+1]*d[k-1] |
---|
| 551 | < (BigInt(3))*d[k]*d[k] - lambda[k][k-1]*lambda[k][k-1]) |
---|
[6ba162] | 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 |
---|