[2cfffe] | 1 | |
---|
| 2 | #include <NTL/LLL.h> |
---|
| 3 | #include <NTL/fileio.h> |
---|
| 4 | #include <NTL/vec_double.h> |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | #include <NTL/new.h> |
---|
| 8 | |
---|
| 9 | NTL_START_IMPL |
---|
| 10 | |
---|
[09da99] | 11 | static inline |
---|
| 12 | void CheckFinite(double *p) |
---|
| 13 | { |
---|
| 14 | if (!IsFinite(p)) Error("LLL_FP: numbers too big...use LLL_XD"); |
---|
| 15 | } |
---|
[2cfffe] | 16 | |
---|
| 17 | static double InnerProduct(double *a, double *b, long n) |
---|
| 18 | { |
---|
| 19 | double s; |
---|
| 20 | long i; |
---|
| 21 | |
---|
| 22 | s = 0; |
---|
| 23 | for (i = 1; i <= n; i++) |
---|
| 24 | s += a[i]*b[i]; |
---|
| 25 | |
---|
| 26 | return s; |
---|
| 27 | } |
---|
| 28 | |
---|
| 29 | static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1) |
---|
| 30 | // x = x - y*MU |
---|
| 31 | { |
---|
| 32 | static ZZ T, MU; |
---|
| 33 | long k; |
---|
| 34 | |
---|
| 35 | long n = A.length(); |
---|
| 36 | long i; |
---|
| 37 | |
---|
| 38 | MU = MU1; |
---|
| 39 | |
---|
| 40 | if (MU == 1) { |
---|
| 41 | for (i = 1; i <= n; i++) |
---|
| 42 | sub(A(i), A(i), B(i)); |
---|
| 43 | |
---|
| 44 | return; |
---|
| 45 | } |
---|
| 46 | |
---|
| 47 | if (MU == -1) { |
---|
| 48 | for (i = 1; i <= n; i++) |
---|
| 49 | add(A(i), A(i), B(i)); |
---|
| 50 | |
---|
| 51 | return; |
---|
| 52 | } |
---|
| 53 | |
---|
| 54 | if (MU == 0) return; |
---|
| 55 | |
---|
| 56 | if (NumTwos(MU) >= NTL_ZZ_NBITS) |
---|
| 57 | k = MakeOdd(MU); |
---|
| 58 | else |
---|
| 59 | k = 0; |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | if (MU.WideSinglePrecision()) { |
---|
| 63 | long mu1; |
---|
| 64 | conv(mu1, MU); |
---|
| 65 | |
---|
| 66 | for (i = 1; i <= n; i++) { |
---|
| 67 | mul(T, B(i), mu1); |
---|
| 68 | if (k > 0) LeftShift(T, T, k); |
---|
| 69 | sub(A(i), A(i), T); |
---|
| 70 | } |
---|
| 71 | } |
---|
| 72 | else { |
---|
| 73 | for (i = 1; i <= n; i++) { |
---|
| 74 | mul(T, B(i), MU); |
---|
| 75 | if (k > 0) LeftShift(T, T, k); |
---|
| 76 | sub(A(i), A(i), T); |
---|
| 77 | } |
---|
| 78 | } |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | #define TR_BND (NTL_FDOUBLE_PRECISION/2.0) |
---|
| 83 | // Just to be safe!! |
---|
| 84 | |
---|
| 85 | static double max_abs(double *v, long n) |
---|
| 86 | { |
---|
| 87 | long i; |
---|
| 88 | double res, t; |
---|
| 89 | |
---|
| 90 | res = 0; |
---|
| 91 | |
---|
| 92 | for (i = 1; i <= n; i++) { |
---|
| 93 | t = fabs(v[i]); |
---|
| 94 | if (t > res) res = t; |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | return res; |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | static void RowTransformStart(double *a, long *in_a, long& in_float, long n) |
---|
| 102 | { |
---|
| 103 | long i; |
---|
| 104 | long inf = 1; |
---|
| 105 | |
---|
| 106 | for (i = 1; i <= n; i++) { |
---|
| 107 | in_a[i] = (a[i] < TR_BND && a[i] > -TR_BND); |
---|
| 108 | inf = inf & in_a[i]; |
---|
| 109 | } |
---|
| 110 | |
---|
| 111 | in_float = inf; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | static void RowTransformFinish(vec_ZZ& A, double *a, long *in_a) |
---|
| 116 | { |
---|
| 117 | long n = A.length(); |
---|
| 118 | long i; |
---|
| 119 | |
---|
| 120 | for (i = 1; i <= n; i++) { |
---|
| 121 | if (in_a[i]) { |
---|
| 122 | conv(A(i), a[i]); |
---|
| 123 | } |
---|
| 124 | else { |
---|
| 125 | conv(a[i], A(i)); |
---|
[09da99] | 126 | CheckFinite(&a[i]); |
---|
[2cfffe] | 127 | } |
---|
| 128 | } |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | static void RowTransform(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1, |
---|
[09da99] | 133 | double *a, double *b, long *in_a, |
---|
[2cfffe] | 134 | double& max_a, double max_b, long& in_float) |
---|
| 135 | // x = x - y*MU |
---|
| 136 | { |
---|
| 137 | static ZZ T, MU; |
---|
| 138 | long k; |
---|
[09da99] | 139 | double mu; |
---|
| 140 | |
---|
| 141 | conv(mu, MU1); |
---|
| 142 | CheckFinite(&mu); |
---|
[2cfffe] | 143 | |
---|
| 144 | long n = A.length(); |
---|
| 145 | long i; |
---|
| 146 | |
---|
| 147 | if (in_float) { |
---|
[09da99] | 148 | double mu_abs = fabs(mu); |
---|
| 149 | if (mu_abs > 0 && max_b > 0 && (mu_abs >= TR_BND || max_b >= TR_BND)) { |
---|
[2cfffe] | 150 | in_float = 0; |
---|
| 151 | } |
---|
[09da99] | 152 | else { |
---|
| 153 | max_a += mu_abs*max_b; |
---|
| 154 | if (max_a >= TR_BND) |
---|
| 155 | in_float = 0; |
---|
| 156 | } |
---|
[2cfffe] | 157 | } |
---|
| 158 | |
---|
| 159 | if (in_float) { |
---|
| 160 | if (mu == 1) { |
---|
| 161 | for (i = 1; i <= n; i++) |
---|
| 162 | a[i] -= b[i]; |
---|
| 163 | |
---|
| 164 | return; |
---|
| 165 | } |
---|
| 166 | |
---|
| 167 | if (mu == -1) { |
---|
| 168 | for (i = 1; i <= n; i++) |
---|
| 169 | a[i] += b[i]; |
---|
| 170 | |
---|
| 171 | return; |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | if (mu == 0) return; |
---|
| 175 | |
---|
| 176 | for (i = 1; i <= n; i++) |
---|
| 177 | a[i] -= mu*b[i]; |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | return; |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | |
---|
| 184 | MU = MU1; |
---|
| 185 | |
---|
| 186 | if (MU == 1) { |
---|
| 187 | for (i = 1; i <= n; i++) { |
---|
| 188 | if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND && |
---|
| 189 | b[i] < TR_BND && b[i] > -TR_BND) { |
---|
| 190 | |
---|
| 191 | a[i] -= b[i]; |
---|
| 192 | } |
---|
| 193 | else { |
---|
| 194 | if (in_a[i]) { |
---|
| 195 | conv(A(i), a[i]); |
---|
| 196 | in_a[i] = 0; |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | sub(A(i), A(i), B(i)); |
---|
| 200 | } |
---|
| 201 | } |
---|
| 202 | return; |
---|
| 203 | } |
---|
| 204 | |
---|
| 205 | if (MU == -1) { |
---|
| 206 | for (i = 1; i <= n; i++) { |
---|
| 207 | if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND && |
---|
| 208 | b[i] < TR_BND && b[i] > -TR_BND) { |
---|
| 209 | |
---|
| 210 | a[i] += b[i]; |
---|
| 211 | } |
---|
| 212 | else { |
---|
| 213 | if (in_a[i]) { |
---|
| 214 | conv(A(i), a[i]); |
---|
| 215 | in_a[i] = 0; |
---|
| 216 | } |
---|
| 217 | |
---|
| 218 | add(A(i), A(i), B(i)); |
---|
| 219 | } |
---|
| 220 | } |
---|
| 221 | return; |
---|
| 222 | } |
---|
| 223 | |
---|
| 224 | if (MU == 0) return; |
---|
| 225 | |
---|
| 226 | double b_bnd = fabs(TR_BND/mu) - 1; |
---|
| 227 | if (b_bnd < 0) b_bnd = 0; |
---|
| 228 | |
---|
| 229 | if (NumTwos(MU) >= NTL_ZZ_NBITS) |
---|
| 230 | k = MakeOdd(MU); |
---|
| 231 | else |
---|
| 232 | k = 0; |
---|
| 233 | |
---|
| 234 | |
---|
| 235 | if (MU.WideSinglePrecision()) { |
---|
| 236 | long mu1; |
---|
| 237 | conv(mu1, MU); |
---|
| 238 | |
---|
| 239 | if (k > 0) { |
---|
| 240 | for (i = 1; i <= n; i++) { |
---|
| 241 | if (in_a[i]) { |
---|
| 242 | conv(A(i), a[i]); |
---|
| 243 | in_a[i] = 0; |
---|
| 244 | } |
---|
| 245 | |
---|
| 246 | mul(T, B(i), mu1); |
---|
| 247 | LeftShift(T, T, k); |
---|
| 248 | sub(A(i), A(i), T); |
---|
| 249 | } |
---|
| 250 | } |
---|
| 251 | else { |
---|
| 252 | for (i = 1; i <= n; i++) { |
---|
| 253 | if (in_a[i] && a[i] < TR_BND && a[i] > -TR_BND && |
---|
| 254 | b[i] < b_bnd && b[i] > -b_bnd) { |
---|
| 255 | |
---|
| 256 | a[i] -= b[i]*mu; |
---|
| 257 | } |
---|
| 258 | else { |
---|
| 259 | if (in_a[i]) { |
---|
| 260 | conv(A(i), a[i]); |
---|
| 261 | in_a[i] = 0; |
---|
| 262 | } |
---|
| 263 | mul(T, B(i), mu1); |
---|
| 264 | sub(A(i), A(i), T); |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | } |
---|
| 268 | } |
---|
| 269 | else { |
---|
| 270 | for (i = 1; i <= n; i++) { |
---|
| 271 | if (in_a[i]) { |
---|
| 272 | conv(A(i), a[i]); |
---|
| 273 | in_a[i] = 0; |
---|
| 274 | } |
---|
| 275 | mul(T, B(i), MU); |
---|
| 276 | if (k > 0) LeftShift(T, T, k); |
---|
| 277 | sub(A(i), A(i), T); |
---|
| 278 | } |
---|
| 279 | } |
---|
| 280 | } |
---|
| 281 | |
---|
| 282 | static void RowTransform2(vec_ZZ& A, vec_ZZ& B, const ZZ& MU1) |
---|
| 283 | // x = x + y*MU |
---|
| 284 | |
---|
| 285 | { |
---|
| 286 | static ZZ T, MU; |
---|
| 287 | long k; |
---|
| 288 | |
---|
| 289 | long n = A.length(); |
---|
| 290 | long i; |
---|
| 291 | |
---|
| 292 | MU = MU1; |
---|
| 293 | |
---|
| 294 | if (MU == 1) { |
---|
| 295 | for (i = 1; i <= n; i++) |
---|
| 296 | add(A(i), A(i), B(i)); |
---|
| 297 | |
---|
| 298 | return; |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | if (MU == -1) { |
---|
| 302 | for (i = 1; i <= n; i++) |
---|
| 303 | sub(A(i), A(i), B(i)); |
---|
| 304 | |
---|
| 305 | return; |
---|
| 306 | } |
---|
| 307 | |
---|
| 308 | if (MU == 0) return; |
---|
| 309 | |
---|
| 310 | if (NumTwos(MU) >= NTL_ZZ_NBITS) |
---|
| 311 | k = MakeOdd(MU); |
---|
| 312 | else |
---|
| 313 | k = 0; |
---|
| 314 | |
---|
| 315 | if (MU.WideSinglePrecision()) { |
---|
| 316 | long mu1; |
---|
| 317 | conv(mu1, MU); |
---|
| 318 | |
---|
| 319 | for (i = 1; i <= n; i++) { |
---|
| 320 | mul(T, B(i), mu1); |
---|
| 321 | if (k > 0) LeftShift(T, T, k); |
---|
| 322 | add(A(i), A(i), T); |
---|
| 323 | } |
---|
| 324 | } |
---|
| 325 | else { |
---|
| 326 | for (i = 1; i <= n; i++) { |
---|
| 327 | mul(T, B(i), MU); |
---|
| 328 | if (k > 0) LeftShift(T, T, k); |
---|
| 329 | add(A(i), A(i), T); |
---|
| 330 | } |
---|
| 331 | } |
---|
| 332 | } |
---|
| 333 | |
---|
| 334 | static |
---|
| 335 | void ComputeGS(mat_ZZ& B, double **B1, double **mu, double *b, |
---|
| 336 | double *c, long k, double bound, long st, double *buf) |
---|
| 337 | |
---|
| 338 | { |
---|
| 339 | long n = B.NumCols(); |
---|
| 340 | long i, j; |
---|
| 341 | double s, t1, y, t; |
---|
| 342 | |
---|
| 343 | ZZ T1; |
---|
| 344 | long test; |
---|
| 345 | |
---|
| 346 | double *mu_k = mu[k]; |
---|
| 347 | |
---|
| 348 | if (st < k) { |
---|
| 349 | for (i = 1; i < st; i++) |
---|
| 350 | buf[i] = mu_k[i]*c[i]; |
---|
| 351 | } |
---|
| 352 | |
---|
| 353 | for (j = st; j <= k-1; j++) { |
---|
| 354 | s = InnerProduct(B1[k], B1[j], n); |
---|
| 355 | |
---|
| 356 | // test = b[k]*b[j] >= NTL_FDOUBLE_PRECISION^2 |
---|
| 357 | |
---|
| 358 | test = (b[k]/NTL_FDOUBLE_PRECISION >= NTL_FDOUBLE_PRECISION/b[j]); |
---|
| 359 | |
---|
| 360 | // test = test && s^2 <= b[k]*b[j]/bound, |
---|
| 361 | // but we compute it in a strange way to avoid overflow |
---|
| 362 | |
---|
| 363 | if (test && (y = fabs(s)) != 0) { |
---|
| 364 | t = y/b[j]; |
---|
| 365 | t1 = b[k]/y; |
---|
| 366 | if (t <= 1) |
---|
| 367 | test = (t*bound <= t1); |
---|
| 368 | else if (t1 >= 1) |
---|
| 369 | test = (t <= t1/bound); |
---|
| 370 | else |
---|
| 371 | test = 0; |
---|
| 372 | } |
---|
| 373 | |
---|
| 374 | if (test) { |
---|
| 375 | InnerProduct(T1, B(k), B(j)); |
---|
| 376 | conv(s, T1); |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | double *mu_j = mu[j]; |
---|
| 380 | |
---|
| 381 | t1 = 0; |
---|
| 382 | for (i = 1; i <= j-1; i++) { |
---|
| 383 | t1 += mu_j[i]*buf[i]; |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | mu_k[j] = (buf[j] = (s - t1))/c[j]; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | #if (!NTL_EXT_DOUBLE) |
---|
| 390 | |
---|
| 391 | // Kahan summation |
---|
| 392 | |
---|
| 393 | double c1; |
---|
| 394 | |
---|
| 395 | s = c1 = 0; |
---|
| 396 | for (j = 1; j <= k-1; j++) { |
---|
| 397 | y = mu_k[j]*buf[j] - c1; |
---|
| 398 | t = s+y; |
---|
| 399 | c1 = t-s; |
---|
| 400 | c1 = c1-y; |
---|
| 401 | s = t; |
---|
| 402 | } |
---|
| 403 | |
---|
| 404 | |
---|
| 405 | #else |
---|
| 406 | |
---|
| 407 | s = 0; |
---|
| 408 | for (j = 1; j <= k-1; j++) |
---|
| 409 | s += mu_k[j]*buf[j]; |
---|
| 410 | |
---|
| 411 | #endif |
---|
| 412 | |
---|
| 413 | c[k] = b[k] - s; |
---|
| 414 | } |
---|
| 415 | |
---|
| 416 | static double red_fudge = 0; |
---|
| 417 | static long log_red = 0; |
---|
| 418 | |
---|
| 419 | static long verbose = 0; |
---|
| 420 | |
---|
| 421 | char *LLLDumpFile = 0; |
---|
| 422 | |
---|
| 423 | static unsigned long NumSwaps = 0; |
---|
| 424 | static double RR_GS_time = 0; |
---|
| 425 | static double StartTime = 0; |
---|
| 426 | static double LastTime = 0; |
---|
| 427 | |
---|
| 428 | |
---|
| 429 | |
---|
| 430 | static void init_red_fudge() |
---|
| 431 | { |
---|
| 432 | long i; |
---|
| 433 | |
---|
| 434 | log_red = long(0.50*NTL_DOUBLE_PRECISION); |
---|
| 435 | red_fudge = 1; |
---|
| 436 | |
---|
| 437 | for (i = log_red; i > 0; i--) |
---|
| 438 | red_fudge = red_fudge*0.5; |
---|
| 439 | } |
---|
| 440 | |
---|
| 441 | static void inc_red_fudge() |
---|
| 442 | { |
---|
| 443 | |
---|
| 444 | red_fudge = red_fudge * 2; |
---|
| 445 | log_red--; |
---|
| 446 | |
---|
| 447 | |
---|
| 448 | if (log_red < 4) |
---|
| 449 | Error("LLL_FP: too much loss of precision...stop!"); |
---|
| 450 | } |
---|
| 451 | |
---|
| 452 | |
---|
| 453 | void ComputeGS(const mat_ZZ& B, mat_RR& B1, |
---|
| 454 | mat_RR& mu, vec_RR& b, |
---|
| 455 | vec_RR& c, long k, const RR& bound, long st, |
---|
| 456 | vec_RR& buf, const RR& bound2); |
---|
| 457 | |
---|
| 458 | |
---|
| 459 | |
---|
| 460 | static void RR_GS(mat_ZZ& B, double **B1, double **mu, |
---|
| 461 | double *b, double *c, double *buf, long prec, |
---|
| 462 | long rr_st, long k, long m_orig, |
---|
| 463 | mat_RR& rr_B1, mat_RR& rr_mu, |
---|
| 464 | vec_RR& rr_b, vec_RR& rr_c) |
---|
| 465 | { |
---|
| 466 | double tt; |
---|
| 467 | |
---|
| 468 | tt = GetTime(); |
---|
| 469 | |
---|
| 470 | if (rr_st > k) Error("LLL_FP: can not continue!!!"); |
---|
| 471 | |
---|
| 472 | long old_p = RR::precision(); |
---|
| 473 | RR::SetPrecision(prec); |
---|
| 474 | |
---|
| 475 | long n = B.NumCols(); |
---|
| 476 | |
---|
| 477 | rr_B1.SetDims(k, n); |
---|
| 478 | rr_mu.SetDims(k, m_orig); |
---|
| 479 | rr_b.SetLength(k); |
---|
| 480 | rr_c.SetLength(k); |
---|
| 481 | |
---|
| 482 | vec_RR rr_buf; |
---|
| 483 | rr_buf.SetLength(k); |
---|
| 484 | |
---|
| 485 | long i, j; |
---|
| 486 | |
---|
| 487 | for (i = rr_st; i <= k; i++) |
---|
| 488 | for (j = 1; j <= n; j++) |
---|
| 489 | conv(rr_B1(i, j), B(i, j)); |
---|
| 490 | |
---|
| 491 | for (i = rr_st; i <= k; i++) |
---|
| 492 | InnerProduct(rr_b(i), rr_B1(i), rr_B1(i)); |
---|
| 493 | |
---|
| 494 | |
---|
| 495 | |
---|
| 496 | RR bound; |
---|
| 497 | power2(bound, 2*long(0.15*RR::precision())); |
---|
| 498 | |
---|
| 499 | RR bound2; |
---|
| 500 | power2(bound2, 2*RR::precision()); |
---|
| 501 | |
---|
| 502 | for (i = rr_st; i <= k; i++) |
---|
| 503 | ComputeGS(B, rr_B1, rr_mu, rr_b, rr_c, i, bound, 1, rr_buf, bound2); |
---|
| 504 | |
---|
| 505 | for (i = rr_st; i <= k; i++) |
---|
[09da99] | 506 | for (j = 1; j <= n; j++) { |
---|
[2cfffe] | 507 | conv(B1[i][j], rr_B1(i,j)); |
---|
[09da99] | 508 | CheckFinite(&B1[i][j]); |
---|
| 509 | } |
---|
[2cfffe] | 510 | |
---|
| 511 | for (i = rr_st; i <= k; i++) |
---|
| 512 | for (j = 1; j <= i-1; j++) { |
---|
| 513 | conv(mu[i][j], rr_mu(i,j)); |
---|
| 514 | } |
---|
| 515 | |
---|
| 516 | for (i = rr_st; i <= k; i++) { |
---|
| 517 | conv(b[i], rr_b(i)); |
---|
[09da99] | 518 | CheckFinite(&b[i]); |
---|
[2cfffe] | 519 | } |
---|
| 520 | |
---|
| 521 | |
---|
| 522 | for (i = rr_st; i <= k; i++) { |
---|
| 523 | conv(c[i], rr_c(i)); |
---|
[09da99] | 524 | CheckFinite(&c[i]); |
---|
[2cfffe] | 525 | } |
---|
| 526 | |
---|
| 527 | for (i = 1; i <= k-1; i++) { |
---|
| 528 | conv(buf[i], rr_buf[i]); |
---|
| 529 | } |
---|
| 530 | |
---|
| 531 | |
---|
| 532 | RR::SetPrecision(old_p); |
---|
| 533 | |
---|
| 534 | tt = GetTime()-tt; |
---|
| 535 | RR_GS_time += tt; |
---|
| 536 | } |
---|
| 537 | |
---|
| 538 | void ComputeGS(const mat_ZZ& B, mat_RR& mu, vec_RR& c) |
---|
| 539 | { |
---|
| 540 | long n = B.NumCols(); |
---|
| 541 | long k = B.NumRows(); |
---|
| 542 | |
---|
| 543 | mat_RR B1; |
---|
| 544 | vec_RR b; |
---|
| 545 | |
---|
| 546 | B1.SetDims(k, n); |
---|
| 547 | mu.SetDims(k, k); |
---|
| 548 | b.SetLength(k); |
---|
| 549 | c.SetLength(k); |
---|
| 550 | |
---|
| 551 | vec_RR buf; |
---|
| 552 | buf.SetLength(k); |
---|
| 553 | |
---|
| 554 | long i, j; |
---|
| 555 | |
---|
| 556 | for (i = 1; i <= k; i++) |
---|
| 557 | for (j = 1; j <= n; j++) |
---|
| 558 | conv(B1(i, j), B(i, j)); |
---|
| 559 | |
---|
| 560 | for (i = 1; i <= k; i++) |
---|
| 561 | InnerProduct(b(i), B1(i), B1(i)); |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | |
---|
| 565 | RR bound; |
---|
| 566 | power2(bound, 2*long(0.15*RR::precision())); |
---|
| 567 | |
---|
| 568 | RR bound2; |
---|
| 569 | power2(bound2, 2*RR::precision()); |
---|
| 570 | |
---|
| 571 | |
---|
| 572 | for (i = 1; i <= k; i++) |
---|
| 573 | ComputeGS(B, B1, mu, b, c, i, bound, 1, buf, bound2); |
---|
| 574 | |
---|
| 575 | } |
---|
| 576 | |
---|
| 577 | |
---|
| 578 | |
---|
| 579 | |
---|
| 580 | |
---|
| 581 | static |
---|
| 582 | long ll_LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, |
---|
| 583 | LLLCheckFct check, double **B1, double **mu, |
---|
| 584 | double *b, double *c, |
---|
| 585 | long m, long init_k, long &quit) |
---|
| 586 | { |
---|
| 587 | long n = B.NumCols(); |
---|
| 588 | |
---|
| 589 | long i, j, k, Fc1; |
---|
| 590 | ZZ MU; |
---|
| 591 | double mu1; |
---|
| 592 | |
---|
| 593 | double t1; |
---|
| 594 | ZZ T1; |
---|
| 595 | double *tp; |
---|
| 596 | |
---|
| 597 | |
---|
| 598 | static double bound = 0; |
---|
| 599 | |
---|
| 600 | if (bound == 0) { |
---|
| 601 | // we tolerate a 15% loss of precision in computing |
---|
| 602 | // inner products in ComputeGS. |
---|
| 603 | |
---|
| 604 | bound = 1; |
---|
| 605 | for (i = 2*long(0.15*NTL_DOUBLE_PRECISION); i > 0; i--) |
---|
| 606 | bound = bound * 2; |
---|
| 607 | } |
---|
| 608 | |
---|
| 609 | double half_plus_fudge = 0.5 + red_fudge; |
---|
| 610 | |
---|
| 611 | quit = 0; |
---|
| 612 | k = init_k; |
---|
| 613 | |
---|
| 614 | |
---|
| 615 | vec_long st_mem; |
---|
| 616 | st_mem.SetLength(m+2); |
---|
| 617 | long *st = st_mem.elts(); |
---|
| 618 | |
---|
| 619 | for (i = 1; i < k; i++) |
---|
| 620 | st[i] = i; |
---|
| 621 | |
---|
| 622 | for (i = k; i <= m+1; i++) |
---|
| 623 | st[i] = 1; |
---|
| 624 | |
---|
| 625 | double *buf; |
---|
| 626 | buf = NTL_NEW_OP double [m+1]; |
---|
| 627 | if (!buf) Error("out of memory in lll_LLL_FP"); |
---|
| 628 | |
---|
| 629 | vec_long in_vec_mem; |
---|
| 630 | in_vec_mem.SetLength(n+1); |
---|
| 631 | long *in_vec = in_vec_mem.elts(); |
---|
| 632 | |
---|
| 633 | double *max_b; |
---|
| 634 | max_b = NTL_NEW_OP double [m+1]; |
---|
| 635 | if (!max_b) Error("out of memory in lll_LLL_FP"); |
---|
| 636 | |
---|
| 637 | for (i = 1; i <= m; i++) |
---|
| 638 | max_b[i] = max_abs(B1[i], n); |
---|
| 639 | |
---|
| 640 | long in_float; |
---|
| 641 | |
---|
| 642 | long rst; |
---|
| 643 | long counter; |
---|
| 644 | long start_over; |
---|
| 645 | |
---|
| 646 | long trigger_index; |
---|
| 647 | long small_trigger; |
---|
| 648 | long cnt; |
---|
| 649 | |
---|
| 650 | mat_RR rr_B1; |
---|
| 651 | mat_RR rr_mu; |
---|
| 652 | vec_RR rr_c; |
---|
| 653 | vec_RR rr_b; |
---|
| 654 | |
---|
| 655 | long m_orig = m; |
---|
| 656 | |
---|
| 657 | long rr_st = 1; |
---|
| 658 | |
---|
| 659 | long max_k = 0; |
---|
| 660 | |
---|
| 661 | long prec = RR::precision(); |
---|
| 662 | |
---|
| 663 | double tt; |
---|
| 664 | |
---|
| 665 | long swap_cnt = 0; |
---|
| 666 | |
---|
| 667 | |
---|
| 668 | while (k <= m) { |
---|
| 669 | |
---|
| 670 | if (k > max_k) { |
---|
| 671 | max_k = k; |
---|
| 672 | swap_cnt = 0; |
---|
| 673 | } |
---|
| 674 | |
---|
| 675 | if (k < rr_st) rr_st = k; |
---|
| 676 | |
---|
| 677 | if (st[k] == k) |
---|
| 678 | rst = 1; |
---|
| 679 | else |
---|
| 680 | rst = k; |
---|
| 681 | |
---|
| 682 | if (st[k] < st[k+1]) st[k+1] = st[k]; |
---|
| 683 | ComputeGS(B, B1, mu, b, c, k, bound, st[k], buf); |
---|
[09da99] | 684 | CheckFinite(&c[k]); |
---|
[2cfffe] | 685 | st[k] = k; |
---|
| 686 | |
---|
| 687 | if (swap_cnt > 200000) { |
---|
| 688 | RR_GS(B, B1, mu, b, c, buf, prec, |
---|
| 689 | rr_st, k, m_orig, rr_B1, rr_mu, rr_b, rr_c); |
---|
| 690 | if (rr_st < st[k+1]) st[k+1] = rr_st; |
---|
| 691 | rr_st = k+1; |
---|
| 692 | rst = k; |
---|
| 693 | swap_cnt = 0; |
---|
| 694 | } |
---|
| 695 | |
---|
| 696 | counter = 0; |
---|
| 697 | trigger_index = k; |
---|
| 698 | small_trigger = 0; |
---|
| 699 | cnt = 0; |
---|
| 700 | |
---|
| 701 | long thresh = 10; |
---|
| 702 | long sz=0, new_sz; |
---|
| 703 | |
---|
| 704 | long did_rr_gs = 0; |
---|
| 705 | |
---|
| 706 | |
---|
| 707 | do { |
---|
| 708 | // size reduction |
---|
| 709 | |
---|
| 710 | counter++; |
---|
| 711 | if ((counter & 127) == 0) { |
---|
| 712 | |
---|
| 713 | new_sz = 0; |
---|
| 714 | for (j = 1; j <= n; j++) |
---|
| 715 | new_sz += NumBits(B(k,j)); |
---|
| 716 | |
---|
| 717 | if ((counter >> 7) == 1 || new_sz < sz) { |
---|
| 718 | sz = new_sz; |
---|
| 719 | } |
---|
| 720 | } |
---|
| 721 | |
---|
| 722 | Fc1 = 0; |
---|
| 723 | start_over = 0; |
---|
| 724 | |
---|
| 725 | for (j = rst-1; j >= 1; j--) { |
---|
| 726 | t1 = fabs(mu[k][j]); |
---|
| 727 | if (t1 > half_plus_fudge) { |
---|
| 728 | |
---|
| 729 | |
---|
| 730 | if (!Fc1) { |
---|
| 731 | if (j > trigger_index || |
---|
| 732 | (j == trigger_index && small_trigger)) { |
---|
| 733 | |
---|
| 734 | cnt++; |
---|
| 735 | |
---|
| 736 | if (cnt > thresh) { |
---|
| 737 | if (log_red <= 15) { |
---|
| 738 | |
---|
| 739 | while (log_red > 10) |
---|
| 740 | inc_red_fudge(); |
---|
| 741 | |
---|
| 742 | half_plus_fudge = 0.5 + red_fudge; |
---|
| 743 | |
---|
| 744 | if (!did_rr_gs) { |
---|
| 745 | RR_GS(B, B1, mu, b, c, buf, prec, |
---|
| 746 | rr_st, k, m_orig, rr_B1, rr_mu, rr_b, rr_c); |
---|
| 747 | if (rr_st < st[k+1]) st[k+1] = rr_st; |
---|
| 748 | rr_st = k+1; |
---|
| 749 | did_rr_gs = 1; |
---|
| 750 | rst = k; |
---|
| 751 | trigger_index = k; |
---|
| 752 | small_trigger = 0; |
---|
| 753 | start_over = 1; |
---|
| 754 | break; |
---|
| 755 | } |
---|
| 756 | } |
---|
| 757 | else { |
---|
| 758 | inc_red_fudge(); |
---|
| 759 | half_plus_fudge = 0.5 + red_fudge; |
---|
| 760 | cnt = 0; |
---|
| 761 | } |
---|
| 762 | } |
---|
| 763 | } |
---|
| 764 | |
---|
| 765 | trigger_index = j; |
---|
| 766 | small_trigger = (t1 < 4); |
---|
| 767 | |
---|
| 768 | Fc1 = 1; |
---|
| 769 | if (k < rr_st) rr_st = k; |
---|
| 770 | RowTransformStart(B1[k], in_vec, in_float, n); |
---|
| 771 | } |
---|
| 772 | |
---|
| 773 | |
---|
| 774 | mu1 = mu[k][j]; |
---|
| 775 | if (mu1 >= 0) |
---|
| 776 | mu1 = ceil(mu1-0.5); |
---|
| 777 | else |
---|
| 778 | mu1 = floor(mu1+0.5); |
---|
| 779 | |
---|
| 780 | double *mu_k = mu[k]; |
---|
| 781 | double *mu_j = mu[j]; |
---|
| 782 | |
---|
| 783 | if (mu1 == 1) { |
---|
| 784 | for (i = 1; i <= j-1; i++) |
---|
| 785 | mu_k[i] -= mu_j[i]; |
---|
| 786 | } |
---|
| 787 | else if (mu1 == -1) { |
---|
| 788 | for (i = 1; i <= j-1; i++) |
---|
| 789 | mu_k[i] += mu_j[i]; |
---|
| 790 | } |
---|
| 791 | else { |
---|
| 792 | for (i = 1; i <= j-1; i++) |
---|
| 793 | mu_k[i] -= mu1*mu_j[i]; |
---|
| 794 | } |
---|
| 795 | |
---|
| 796 | mu_k[j] -= mu1; |
---|
| 797 | |
---|
| 798 | conv(MU, mu1); |
---|
| 799 | |
---|
[09da99] | 800 | RowTransform(B(k), B(j), MU, B1[k], B1[j], in_vec, |
---|
[2cfffe] | 801 | max_b[k], max_b[j], in_float); |
---|
| 802 | if (U) RowTransform((*U)(k), (*U)(j), MU); |
---|
| 803 | } |
---|
| 804 | } |
---|
| 805 | |
---|
| 806 | |
---|
| 807 | if (Fc1) { |
---|
| 808 | RowTransformFinish(B(k), B1[k], in_vec); |
---|
| 809 | max_b[k] = max_abs(B1[k], n); |
---|
| 810 | |
---|
| 811 | if (!did_rr_gs) { |
---|
| 812 | b[k] = InnerProduct(B1[k], B1[k], n); |
---|
[09da99] | 813 | CheckFinite(&b[k]); |
---|
| 814 | |
---|
[2cfffe] | 815 | ComputeGS(B, B1, mu, b, c, k, bound, 1, buf); |
---|
[09da99] | 816 | CheckFinite(&c[k]); |
---|
[2cfffe] | 817 | } |
---|
| 818 | else { |
---|
| 819 | RR_GS(B, B1, mu, b, c, buf, prec, |
---|
| 820 | rr_st, k, m_orig, rr_B1, rr_mu, rr_b, rr_c); |
---|
| 821 | rr_st = k+1; |
---|
| 822 | } |
---|
| 823 | |
---|
| 824 | rst = k; |
---|
| 825 | } |
---|
| 826 | } while (Fc1 || start_over); |
---|
| 827 | |
---|
| 828 | if (check && (*check)(B(k))) |
---|
| 829 | quit = 1; |
---|
| 830 | |
---|
| 831 | if (b[k] == 0) { |
---|
| 832 | for (i = k; i < m; i++) { |
---|
| 833 | // swap i, i+1 |
---|
| 834 | swap(B(i), B(i+1)); |
---|
| 835 | tp = B1[i]; B1[i] = B1[i+1]; B1[i+1] = tp; |
---|
| 836 | t1 = b[i]; b[i] = b[i+1]; b[i+1] = t1; |
---|
| 837 | t1 = max_b[i]; max_b[i] = max_b[i+1]; max_b[i+1] = t1; |
---|
| 838 | if (U) swap((*U)(i), (*U)(i+1)); |
---|
| 839 | } |
---|
| 840 | |
---|
| 841 | for (i = k; i <= m+1; i++) st[i] = 1; |
---|
| 842 | if (k < rr_st) rr_st = k; |
---|
| 843 | |
---|
| 844 | m--; |
---|
| 845 | if (quit) break; |
---|
| 846 | continue; |
---|
| 847 | } |
---|
| 848 | |
---|
| 849 | if (quit) break; |
---|
| 850 | |
---|
| 851 | if (deep > 0) { |
---|
| 852 | // deep insertions |
---|
| 853 | |
---|
| 854 | double cc = b[k]; |
---|
| 855 | long l = 1; |
---|
| 856 | while (l <= k-1 && delta*c[l] <= cc) { |
---|
| 857 | cc = cc - mu[k][l]*mu[k][l]*c[l]; |
---|
| 858 | l++; |
---|
| 859 | } |
---|
| 860 | |
---|
| 861 | if (l <= k-1 && (l <= deep || k-l <= deep)) { |
---|
| 862 | // deep insertion at position l |
---|
| 863 | |
---|
| 864 | for (i = k; i > l; i--) { |
---|
| 865 | // swap rows i, i-1 |
---|
| 866 | swap(B(i), B(i-1)); |
---|
| 867 | tp = B1[i]; B1[i] = B1[i-1]; B1[i-1] = tp; |
---|
| 868 | tp = mu[i]; mu[i] = mu[i-1]; mu[i-1] = tp; |
---|
| 869 | t1 = b[i]; b[i] = b[i-1]; b[i-1] = t1; |
---|
| 870 | t1 = max_b[i]; max_b[i] = max_b[i-1]; max_b[i-1] = t1; |
---|
| 871 | if (U) swap((*U)(i), (*U)(i-1)); |
---|
| 872 | } |
---|
| 873 | |
---|
| 874 | k = l; |
---|
| 875 | NumSwaps++; |
---|
| 876 | swap_cnt++; |
---|
| 877 | continue; |
---|
| 878 | } |
---|
| 879 | } // end deep insertions |
---|
| 880 | |
---|
| 881 | // test LLL reduction condition |
---|
| 882 | |
---|
| 883 | if (k > 1 && delta*c[k-1] > c[k] + mu[k][k-1]*mu[k][k-1]*c[k-1]) { |
---|
| 884 | // swap rows k, k-1 |
---|
| 885 | swap(B(k), B(k-1)); |
---|
| 886 | tp = B1[k]; B1[k] = B1[k-1]; B1[k-1] = tp; |
---|
| 887 | tp = mu[k]; mu[k] = mu[k-1]; mu[k-1] = tp; |
---|
| 888 | t1 = b[k]; b[k] = b[k-1]; b[k-1] = t1; |
---|
| 889 | t1 = max_b[k]; max_b[k] = max_b[k-1]; max_b[k-1] = t1; |
---|
| 890 | if (U) swap((*U)(k), (*U)(k-1)); |
---|
| 891 | |
---|
| 892 | k--; |
---|
| 893 | NumSwaps++; |
---|
| 894 | swap_cnt++; |
---|
| 895 | // cout << "-\n"; |
---|
| 896 | } |
---|
| 897 | else { |
---|
| 898 | |
---|
| 899 | k++; |
---|
| 900 | // cout << "+\n"; |
---|
| 901 | } |
---|
| 902 | |
---|
| 903 | } |
---|
| 904 | |
---|
| 905 | delete [] buf; |
---|
| 906 | delete [] max_b; |
---|
| 907 | |
---|
| 908 | return m; |
---|
| 909 | } |
---|
| 910 | |
---|
| 911 | |
---|
| 912 | |
---|
| 913 | |
---|
| 914 | |
---|
| 915 | static |
---|
| 916 | long LLL_FP(mat_ZZ& B, mat_ZZ* U, double delta, long deep, |
---|
| 917 | LLLCheckFct check) |
---|
| 918 | { |
---|
| 919 | long m = B.NumRows(); |
---|
| 920 | long n = B.NumCols(); |
---|
| 921 | |
---|
| 922 | long i, j; |
---|
| 923 | long new_m, dep, quit; |
---|
| 924 | ZZ MU; |
---|
| 925 | |
---|
| 926 | ZZ T1; |
---|
| 927 | |
---|
| 928 | init_red_fudge(); |
---|
| 929 | |
---|
| 930 | if (U) ident(*U, m); |
---|
| 931 | |
---|
| 932 | double **B1; // approximates B |
---|
| 933 | |
---|
| 934 | typedef double *doubleptr; |
---|
| 935 | |
---|
| 936 | B1 = NTL_NEW_OP doubleptr[m+1]; |
---|
| 937 | if (!B1) Error("LLL_FP: out of memory"); |
---|
| 938 | |
---|
| 939 | for (i = 1; i <= m; i++) { |
---|
| 940 | B1[i] = NTL_NEW_OP double[n+1]; |
---|
| 941 | if (!B1[i]) Error("LLL_FP: out of memory"); |
---|
| 942 | } |
---|
| 943 | |
---|
| 944 | double **mu; |
---|
| 945 | mu = NTL_NEW_OP doubleptr[m+1]; |
---|
| 946 | if (!mu) Error("LLL_FP: out of memory"); |
---|
| 947 | |
---|
| 948 | for (i = 1; i <= m; i++) { |
---|
| 949 | mu[i] = NTL_NEW_OP double[m+1]; |
---|
| 950 | if (!mu[i]) Error("LLL_FP: out of memory"); |
---|
| 951 | } |
---|
| 952 | |
---|
| 953 | double *c; // squared lengths of Gramm-Schmidt basis vectors |
---|
| 954 | |
---|
| 955 | c = NTL_NEW_OP double[m+1]; |
---|
| 956 | if (!c) Error("LLL_FP: out of memory"); |
---|
| 957 | |
---|
| 958 | double *b; // squared lengths of basis vectors |
---|
| 959 | |
---|
| 960 | b = NTL_NEW_OP double[m+1]; |
---|
| 961 | if (!b) Error("LLL_FP: out of memory"); |
---|
| 962 | |
---|
| 963 | |
---|
| 964 | for (i = 1; i <=m; i++) |
---|
[09da99] | 965 | for (j = 1; j <= n; j++) { |
---|
[2cfffe] | 966 | conv(B1[i][j], B(i, j)); |
---|
[09da99] | 967 | CheckFinite(&B1[i][j]); |
---|
| 968 | } |
---|
[2cfffe] | 969 | |
---|
| 970 | |
---|
| 971 | for (i = 1; i <= m; i++) { |
---|
| 972 | b[i] = InnerProduct(B1[i], B1[i], n); |
---|
[09da99] | 973 | CheckFinite(&b[i]); |
---|
[2cfffe] | 974 | } |
---|
| 975 | |
---|
| 976 | new_m = ll_LLL_FP(B, U, delta, deep, check, B1, mu, b, c, m, 1, quit); |
---|
| 977 | dep = m - new_m; |
---|
| 978 | m = new_m; |
---|
| 979 | |
---|
| 980 | if (dep > 0) { |
---|
| 981 | // for consistency, we move all of the zero rows to the front |
---|
| 982 | |
---|
| 983 | for (i = 0; i < m; i++) { |
---|
| 984 | swap(B(m+dep-i), B(m-i)); |
---|
| 985 | if (U) swap((*U)(m+dep-i), (*U)(m-i)); |
---|
| 986 | } |
---|
| 987 | } |
---|
| 988 | |
---|
| 989 | |
---|
| 990 | // clean-up |
---|
| 991 | |
---|
| 992 | for (i = 1; i <= m; i++) { |
---|
| 993 | delete [] B1[i]; |
---|
| 994 | } |
---|
| 995 | |
---|
| 996 | delete [] B1; |
---|
| 997 | |
---|
| 998 | for (i = 1; i <= m; i++) { |
---|
| 999 | delete [] mu[i]; |
---|
| 1000 | } |
---|
| 1001 | |
---|
| 1002 | delete [] mu; |
---|
| 1003 | |
---|
| 1004 | delete [] c; |
---|
| 1005 | |
---|
| 1006 | delete [] b; |
---|
| 1007 | |
---|
| 1008 | return m; |
---|
| 1009 | } |
---|
| 1010 | |
---|
| 1011 | |
---|
| 1012 | |
---|
| 1013 | long LLL_FP(mat_ZZ& B, double delta, long deep, LLLCheckFct check, |
---|
| 1014 | long verb) |
---|
| 1015 | { |
---|
| 1016 | verbose = verb; |
---|
| 1017 | RR_GS_time = 0; |
---|
| 1018 | NumSwaps = 0; |
---|
| 1019 | if (verbose) { |
---|
| 1020 | StartTime = GetTime(); |
---|
| 1021 | LastTime = StartTime; |
---|
| 1022 | } |
---|
| 1023 | |
---|
| 1024 | if (delta < 0.50 || delta >= 1) Error("LLL_FP: bad delta"); |
---|
| 1025 | if (deep < 0) Error("LLL_FP: bad deep"); |
---|
| 1026 | return LLL_FP(B, 0, delta, deep, check); |
---|
| 1027 | } |
---|
| 1028 | |
---|
| 1029 | long LLL_FP(mat_ZZ& B, mat_ZZ& U, double delta, long deep, |
---|
| 1030 | LLLCheckFct check, long verb) |
---|
| 1031 | { |
---|
| 1032 | verbose = verb; |
---|
| 1033 | RR_GS_time = 0; |
---|
| 1034 | NumSwaps = 0; |
---|
| 1035 | if (verbose) { |
---|
| 1036 | StartTime = GetTime(); |
---|
| 1037 | LastTime = StartTime; |
---|
| 1038 | } |
---|
| 1039 | |
---|
| 1040 | if (delta < 0.50 || delta >= 1) Error("LLL_FP: bad delta"); |
---|
| 1041 | if (deep < 0) Error("LLL_FP: bad deep"); |
---|
| 1042 | return LLL_FP(B, &U, delta, deep, check); |
---|
| 1043 | } |
---|
| 1044 | |
---|
| 1045 | |
---|
| 1046 | |
---|
| 1047 | static vec_double BKZConstant; |
---|
| 1048 | |
---|
| 1049 | static |
---|
| 1050 | void ComputeBKZConstant(long beta, long p) |
---|
| 1051 | { |
---|
| 1052 | const double c_PI = 3.14159265358979323846264338328; |
---|
| 1053 | const double LogPI = 1.14472988584940017414342735135; |
---|
| 1054 | |
---|
| 1055 | BKZConstant.SetLength(beta-1); |
---|
| 1056 | |
---|
| 1057 | vec_double Log; |
---|
| 1058 | Log.SetLength(beta); |
---|
| 1059 | |
---|
| 1060 | |
---|
| 1061 | long i, j, k; |
---|
| 1062 | double x, y; |
---|
| 1063 | |
---|
| 1064 | for (j = 1; j <= beta; j++) |
---|
| 1065 | Log(j) = log(double(j)); |
---|
| 1066 | |
---|
| 1067 | for (i = 1; i <= beta-1; i++) { |
---|
| 1068 | // First, we compute x = gamma(i/2)^{2/i} |
---|
| 1069 | |
---|
| 1070 | k = i/2; |
---|
| 1071 | |
---|
| 1072 | if ((i & 1) == 0) { // i even |
---|
| 1073 | x = 0; |
---|
| 1074 | for (j = 1; j <= k; j++) |
---|
| 1075 | x = x + Log(j); |
---|
| 1076 | |
---|
| 1077 | x = x * (1/double(k)); |
---|
| 1078 | |
---|
| 1079 | x = exp(x); |
---|
| 1080 | } |
---|
| 1081 | else { // i odd |
---|
| 1082 | x = 0; |
---|
| 1083 | for (j = k + 2; j <= 2*k + 2; j++) |
---|
| 1084 | x = x + Log(j); |
---|
| 1085 | |
---|
| 1086 | x = 0.5*LogPI + x - 2*(k+1)*Log(2); |
---|
| 1087 | |
---|
| 1088 | x = x * (2.0/double(i)); |
---|
| 1089 | |
---|
| 1090 | x = exp(x); |
---|
| 1091 | } |
---|
| 1092 | |
---|
| 1093 | // Second, we compute y = 2^{2*p/i} |
---|
| 1094 | |
---|
| 1095 | y = -(2*p/double(i))*Log(2); |
---|
| 1096 | y = exp(y); |
---|
| 1097 | |
---|
| 1098 | BKZConstant(i) = x*y/c_PI; |
---|
| 1099 | } |
---|
| 1100 | } |
---|
| 1101 | |
---|
| 1102 | static vec_double BKZThresh; |
---|
| 1103 | |
---|
| 1104 | static |
---|
| 1105 | void ComputeBKZThresh(double *c, long beta) |
---|
| 1106 | { |
---|
| 1107 | BKZThresh.SetLength(beta-1); |
---|
| 1108 | |
---|
| 1109 | long i; |
---|
| 1110 | double x; |
---|
| 1111 | |
---|
| 1112 | x = 0; |
---|
| 1113 | |
---|
| 1114 | for (i = 1; i <= beta-1; i++) { |
---|
| 1115 | x += log(c[i-1]); |
---|
| 1116 | BKZThresh(i) = exp(x/double(i))*BKZConstant(i); |
---|
| 1117 | if (!IsFinite(&BKZThresh(i))) BKZThresh(i) = 0; |
---|
| 1118 | } |
---|
| 1119 | } |
---|
| 1120 | |
---|
| 1121 | static |
---|
| 1122 | long BKZ_FP(mat_ZZ& BB, mat_ZZ* UU, double delta, |
---|
| 1123 | long beta, long prune, LLLCheckFct check) |
---|
| 1124 | { |
---|
| 1125 | |
---|
| 1126 | |
---|
| 1127 | |
---|
| 1128 | long m = BB.NumRows(); |
---|
| 1129 | long n = BB.NumCols(); |
---|
| 1130 | long m_orig = m; |
---|
| 1131 | |
---|
| 1132 | long i, j; |
---|
| 1133 | ZZ MU; |
---|
| 1134 | |
---|
| 1135 | double t1; |
---|
| 1136 | ZZ T1; |
---|
| 1137 | double *tp; |
---|
| 1138 | |
---|
| 1139 | init_red_fudge(); |
---|
| 1140 | |
---|
| 1141 | mat_ZZ B; |
---|
| 1142 | B = BB; |
---|
| 1143 | |
---|
| 1144 | B.SetDims(m+1, n); |
---|
| 1145 | |
---|
| 1146 | |
---|
| 1147 | double **B1; // approximates B |
---|
| 1148 | |
---|
| 1149 | typedef double *doubleptr; |
---|
| 1150 | |
---|
| 1151 | B1 = NTL_NEW_OP doubleptr[m+2]; |
---|
| 1152 | if (!B1) Error("BKZ_FP: out of memory"); |
---|
| 1153 | |
---|
| 1154 | for (i = 1; i <= m+1; i++) { |
---|
| 1155 | B1[i] = NTL_NEW_OP double[n+1]; |
---|
| 1156 | if (!B1[i]) Error("BKZ_FP: out of memory"); |
---|
| 1157 | } |
---|
| 1158 | |
---|
| 1159 | double **mu; |
---|
| 1160 | mu = NTL_NEW_OP doubleptr[m+2]; |
---|
| 1161 | if (!mu) Error("LLL_FP: out of memory"); |
---|
| 1162 | |
---|
| 1163 | for (i = 1; i <= m+1; i++) { |
---|
| 1164 | mu[i] = NTL_NEW_OP double[m+1]; |
---|
| 1165 | if (!mu[i]) Error("BKZ_FP: out of memory"); |
---|
| 1166 | } |
---|
| 1167 | |
---|
| 1168 | |
---|
| 1169 | double *c; // squared lengths of Gramm-Schmidt basis vectors |
---|
| 1170 | |
---|
| 1171 | c = NTL_NEW_OP double[m+2]; |
---|
| 1172 | if (!c) Error("BKZ_FP: out of memory"); |
---|
| 1173 | |
---|
| 1174 | double *b; // squared lengths of basis vectors |
---|
| 1175 | |
---|
| 1176 | b = NTL_NEW_OP double[m+2]; |
---|
| 1177 | if (!b) Error("BKZ_FP: out of memory"); |
---|
| 1178 | |
---|
| 1179 | double cbar; |
---|
| 1180 | |
---|
| 1181 | double *ctilda; |
---|
| 1182 | ctilda = NTL_NEW_OP double[m+2]; |
---|
| 1183 | if (!ctilda) Error("BKZ_FP: out of memory"); |
---|
| 1184 | |
---|
| 1185 | double *vvec; |
---|
| 1186 | vvec = NTL_NEW_OP double[m+2]; |
---|
| 1187 | if (!vvec) Error("BKZ_FP: out of memory"); |
---|
| 1188 | |
---|
| 1189 | double *yvec; |
---|
| 1190 | yvec = NTL_NEW_OP double[m+2]; |
---|
| 1191 | if (!yvec) Error("BKZ_FP: out of memory"); |
---|
| 1192 | |
---|
| 1193 | double *uvec; |
---|
| 1194 | uvec = NTL_NEW_OP double[m+2]; |
---|
| 1195 | if (!uvec) Error("BKZ_FP: out of memory"); |
---|
| 1196 | |
---|
| 1197 | double *utildavec; |
---|
| 1198 | utildavec = NTL_NEW_OP double[m+2]; |
---|
| 1199 | if (!utildavec) Error("BKZ_FP: out of memory"); |
---|
| 1200 | |
---|
| 1201 | |
---|
| 1202 | long *Deltavec; |
---|
| 1203 | Deltavec = NTL_NEW_OP long[m+2]; |
---|
| 1204 | if (!Deltavec) Error("BKZ_FP: out of memory"); |
---|
| 1205 | |
---|
| 1206 | long *deltavec; |
---|
| 1207 | deltavec = NTL_NEW_OP long[m+2]; |
---|
| 1208 | if (!deltavec) Error("BKZ_FP: out of memory"); |
---|
| 1209 | |
---|
| 1210 | mat_ZZ Ulocal; |
---|
| 1211 | mat_ZZ *U; |
---|
| 1212 | |
---|
| 1213 | if (UU) { |
---|
| 1214 | Ulocal.SetDims(m+1, m); |
---|
| 1215 | for (i = 1; i <= m; i++) |
---|
| 1216 | conv(Ulocal(i, i), 1); |
---|
| 1217 | U = &Ulocal; |
---|
| 1218 | } |
---|
| 1219 | else |
---|
| 1220 | U = 0; |
---|
| 1221 | |
---|
| 1222 | long quit; |
---|
| 1223 | long new_m; |
---|
| 1224 | long z, jj, kk; |
---|
| 1225 | long s, t; |
---|
| 1226 | long h; |
---|
| 1227 | double eta; |
---|
| 1228 | |
---|
| 1229 | |
---|
| 1230 | for (i = 1; i <=m; i++) |
---|
[09da99] | 1231 | for (j = 1; j <= n; j++) { |
---|
[2cfffe] | 1232 | conv(B1[i][j], B(i, j)); |
---|
[09da99] | 1233 | CheckFinite(&B1[i][j]); |
---|
| 1234 | } |
---|
[2cfffe] | 1235 | |
---|
| 1236 | |
---|
| 1237 | for (i = 1; i <= m; i++) { |
---|
| 1238 | b[i] = InnerProduct(B1[i], B1[i], n); |
---|
[09da99] | 1239 | CheckFinite(&b[i]); |
---|
[2cfffe] | 1240 | } |
---|
| 1241 | |
---|
| 1242 | |
---|
| 1243 | |
---|
| 1244 | m = ll_LLL_FP(B, U, delta, 0, check, B1, mu, b, c, m, 1, quit); |
---|
| 1245 | |
---|
| 1246 | double tt; |
---|
| 1247 | |
---|
| 1248 | double enum_time = 0; |
---|
[09da99] | 1249 | unsigned long NumIterations = 0; |
---|
| 1250 | unsigned long NumTrivial = 0; |
---|
| 1251 | unsigned long NumNonTrivial = 0; |
---|
| 1252 | unsigned long NumNoOps = 0; |
---|
[2cfffe] | 1253 | |
---|
| 1254 | long verb = verbose; |
---|
| 1255 | |
---|
| 1256 | verbose = 0; |
---|
| 1257 | |
---|
| 1258 | long clean = 1; |
---|
| 1259 | |
---|
| 1260 | if (m < m_orig) { |
---|
| 1261 | for (i = m_orig+1; i >= m+2; i--) { |
---|
| 1262 | // swap i, i-1 |
---|
| 1263 | |
---|
| 1264 | swap(B(i), B(i-1)); |
---|
| 1265 | if (U) swap((*U)(i), (*U)(i-1)); |
---|
| 1266 | } |
---|
| 1267 | } |
---|
| 1268 | |
---|
| 1269 | if (!quit && m > 1) { |
---|
| 1270 | if (beta > m) beta = m; |
---|
| 1271 | |
---|
| 1272 | if (prune > 0) |
---|
| 1273 | ComputeBKZConstant(beta, prune); |
---|
| 1274 | |
---|
| 1275 | z = 0; |
---|
| 1276 | jj = 0; |
---|
| 1277 | |
---|
| 1278 | while (z < m-1) { |
---|
| 1279 | jj++; |
---|
| 1280 | kk = min(jj+beta-1, m); |
---|
| 1281 | |
---|
| 1282 | if (jj == m) { |
---|
| 1283 | jj = 1; |
---|
| 1284 | kk = beta; |
---|
| 1285 | clean = 1; |
---|
| 1286 | } |
---|
| 1287 | |
---|
| 1288 | // ENUM |
---|
| 1289 | |
---|
| 1290 | double tt1; |
---|
| 1291 | |
---|
| 1292 | if (verb) { |
---|
| 1293 | tt1 = GetTime(); |
---|
| 1294 | } |
---|
| 1295 | |
---|
| 1296 | |
---|
| 1297 | if (prune > 0) |
---|
| 1298 | ComputeBKZThresh(&c[jj], kk-jj+1); |
---|
| 1299 | |
---|
| 1300 | |
---|
| 1301 | cbar = c[jj]; |
---|
| 1302 | utildavec[jj] = uvec[jj] = 1; |
---|
| 1303 | |
---|
| 1304 | yvec[jj] = vvec[jj] = 0; |
---|
| 1305 | Deltavec[jj] = 0; |
---|
| 1306 | |
---|
| 1307 | |
---|
| 1308 | s = t = jj; |
---|
| 1309 | deltavec[jj] = 1; |
---|
| 1310 | |
---|
| 1311 | for (i = jj+1; i <= kk+1; i++) { |
---|
| 1312 | ctilda[i] = uvec[i] = utildavec[i] = yvec[i] = 0; |
---|
| 1313 | Deltavec[i] = 0; |
---|
| 1314 | vvec[i] = 0; |
---|
| 1315 | deltavec[i] = 1; |
---|
| 1316 | } |
---|
| 1317 | |
---|
| 1318 | long enum_cnt = 0; |
---|
| 1319 | |
---|
| 1320 | while (t <= kk) { |
---|
| 1321 | |
---|
| 1322 | ctilda[t] = ctilda[t+1] + |
---|
| 1323 | (yvec[t]+utildavec[t])*(yvec[t]+utildavec[t])*c[t]; |
---|
| 1324 | |
---|
| 1325 | if (prune > 0 && t > jj) { |
---|
| 1326 | eta = BKZThresh(t-jj); |
---|
| 1327 | } |
---|
| 1328 | else |
---|
| 1329 | eta = 0; |
---|
| 1330 | |
---|
| 1331 | if (ctilda[t] < cbar - eta) { |
---|
| 1332 | if (t > jj) { |
---|
| 1333 | t--; |
---|
| 1334 | t1 = 0; |
---|
| 1335 | for (i = t+1; i <= s; i++) |
---|
| 1336 | t1 += utildavec[i]*mu[i][t]; |
---|
| 1337 | yvec[t] = t1; |
---|
| 1338 | t1 = -t1; |
---|
| 1339 | if (t1 >= 0) |
---|
| 1340 | t1 = ceil(t1-0.5); |
---|
| 1341 | else |
---|
| 1342 | t1 = floor(t1+0.5); |
---|
| 1343 | utildavec[t] = vvec[t] = t1; |
---|
| 1344 | Deltavec[t] = 0; |
---|
| 1345 | if (utildavec[t] > -yvec[t]) |
---|
| 1346 | deltavec[t] = -1; |
---|
| 1347 | else |
---|
| 1348 | deltavec[t] = 1; |
---|
| 1349 | } |
---|
| 1350 | else { |
---|
| 1351 | cbar = ctilda[jj]; |
---|
| 1352 | for (i = jj; i <= kk; i++) { |
---|
| 1353 | uvec[i] = utildavec[i]; |
---|
| 1354 | } |
---|
| 1355 | } |
---|
| 1356 | } |
---|
| 1357 | else { |
---|
| 1358 | t++; |
---|
| 1359 | s = max(s, t); |
---|
| 1360 | if (t < s) Deltavec[t] = -Deltavec[t]; |
---|
| 1361 | if (Deltavec[t]*deltavec[t] >= 0) Deltavec[t] += deltavec[t]; |
---|
| 1362 | utildavec[t] = vvec[t] + Deltavec[t]; |
---|
| 1363 | } |
---|
| 1364 | } |
---|
| 1365 | |
---|
| 1366 | if (verb) { |
---|
| 1367 | tt1 = GetTime() - tt1; |
---|
| 1368 | enum_time += tt1; |
---|
| 1369 | } |
---|
| 1370 | |
---|
| 1371 | NumIterations++; |
---|
| 1372 | |
---|
| 1373 | h = min(kk+1, m); |
---|
| 1374 | |
---|
| 1375 | if ((delta - 8*red_fudge)*c[jj] > cbar) { |
---|
| 1376 | |
---|
| 1377 | clean = 0; |
---|
| 1378 | |
---|
| 1379 | // we treat the case that the new vector is b_s (jj < s <= kk) |
---|
| 1380 | // as a special case that appears to occur most of the time. |
---|
| 1381 | |
---|
| 1382 | s = 0; |
---|
| 1383 | for (i = jj+1; i <= kk; i++) { |
---|
| 1384 | if (uvec[i] != 0) { |
---|
| 1385 | if (s == 0) |
---|
| 1386 | s = i; |
---|
| 1387 | else |
---|
| 1388 | s = -1; |
---|
| 1389 | } |
---|
| 1390 | } |
---|
| 1391 | |
---|
| 1392 | if (s == 0) Error("BKZ_FP: internal error"); |
---|
| 1393 | |
---|
| 1394 | if (s > 0) { |
---|
| 1395 | // special case |
---|
| 1396 | |
---|
| 1397 | NumTrivial++; |
---|
| 1398 | |
---|
| 1399 | for (i = s; i > jj; i--) { |
---|
| 1400 | // swap i, i-1 |
---|
| 1401 | swap(B(i-1), B(i)); |
---|
| 1402 | if (U) swap((*U)(i-1), (*U)(i)); |
---|
| 1403 | tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp; |
---|
| 1404 | t1 = b[i-1]; b[i-1] = b[i]; b[i] = t1; |
---|
| 1405 | } |
---|
| 1406 | |
---|
| 1407 | // cerr << "special case\n"; |
---|
| 1408 | new_m = ll_LLL_FP(B, U, delta, 0, check, |
---|
| 1409 | B1, mu, b, c, h, jj, quit); |
---|
| 1410 | if (new_m != h) Error("BKZ_FP: internal error"); |
---|
| 1411 | if (quit) break; |
---|
| 1412 | } |
---|
| 1413 | else { |
---|
| 1414 | // the general case |
---|
| 1415 | |
---|
| 1416 | NumNonTrivial++; |
---|
| 1417 | |
---|
| 1418 | for (i = 1; i <= n; i++) conv(B(m+1, i), 0); |
---|
| 1419 | |
---|
| 1420 | if (U) { |
---|
| 1421 | for (i = 1; i <= m_orig; i++) |
---|
| 1422 | conv((*U)(m+1, i), 0); |
---|
| 1423 | } |
---|
| 1424 | |
---|
| 1425 | for (i = jj; i <= kk; i++) { |
---|
| 1426 | if (uvec[i] == 0) continue; |
---|
| 1427 | conv(MU, uvec[i]); |
---|
| 1428 | RowTransform2(B(m+1), B(i), MU); |
---|
| 1429 | if (U) RowTransform2((*U)(m+1), (*U)(i), MU); |
---|
| 1430 | } |
---|
| 1431 | |
---|
| 1432 | for (i = m+1; i >= jj+1; i--) { |
---|
| 1433 | // swap i, i-1 |
---|
| 1434 | swap(B(i-1), B(i)); |
---|
| 1435 | if (U) swap((*U)(i-1), (*U)(i)); |
---|
| 1436 | tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp; |
---|
| 1437 | t1 = b[i-1]; b[i-1] = b[i]; b[i] = t1; |
---|
| 1438 | } |
---|
| 1439 | |
---|
[09da99] | 1440 | for (i = 1; i <= n; i++) { |
---|
[2cfffe] | 1441 | conv(B1[jj][i], B(jj, i)); |
---|
[09da99] | 1442 | CheckFinite(&B1[jj][i]); |
---|
| 1443 | } |
---|
[2cfffe] | 1444 | |
---|
| 1445 | b[jj] = InnerProduct(B1[jj], B1[jj], n); |
---|
[09da99] | 1446 | CheckFinite(&b[jj]); |
---|
[2cfffe] | 1447 | |
---|
| 1448 | if (b[jj] == 0) Error("BKZ_FP: internal error"); |
---|
| 1449 | |
---|
| 1450 | // remove linear dependencies |
---|
| 1451 | |
---|
| 1452 | // cerr << "general case\n"; |
---|
| 1453 | new_m = ll_LLL_FP(B, U, delta, 0, 0, B1, mu, b, c, kk+1, jj, quit); |
---|
| 1454 | |
---|
| 1455 | if (new_m != kk) Error("BKZ_FP: internal error"); |
---|
| 1456 | |
---|
| 1457 | // remove zero vector |
---|
| 1458 | |
---|
| 1459 | for (i = kk+2; i <= m+1; i++) { |
---|
| 1460 | // swap i, i-1 |
---|
| 1461 | swap(B(i-1), B(i)); |
---|
| 1462 | if (U) swap((*U)(i-1), (*U)(i)); |
---|
| 1463 | tp = B1[i-1]; B1[i-1] = B1[i]; B1[i] = tp; |
---|
| 1464 | t1 = b[i-1]; b[i-1] = b[i]; b[i] = t1; |
---|
| 1465 | } |
---|
| 1466 | |
---|
| 1467 | quit = 0; |
---|
| 1468 | if (check) { |
---|
| 1469 | for (i = 1; i <= kk; i++) |
---|
| 1470 | if ((*check)(B(i))) { |
---|
| 1471 | quit = 1; |
---|
| 1472 | break; |
---|
| 1473 | } |
---|
| 1474 | } |
---|
| 1475 | |
---|
| 1476 | if (quit) break; |
---|
| 1477 | |
---|
| 1478 | if (h > kk) { |
---|
| 1479 | // extend reduced basis |
---|
| 1480 | |
---|
| 1481 | new_m = ll_LLL_FP(B, U, delta, 0, check, |
---|
| 1482 | B1, mu, b, c, h, h, quit); |
---|
| 1483 | |
---|
| 1484 | if (new_m != h) Error("BKZ_FP: internal error"); |
---|
| 1485 | if (quit) break; |
---|
| 1486 | } |
---|
| 1487 | } |
---|
| 1488 | |
---|
| 1489 | z = 0; |
---|
| 1490 | } |
---|
| 1491 | else { |
---|
| 1492 | // LLL_FP |
---|
| 1493 | // cerr << "progress\n"; |
---|
| 1494 | |
---|
| 1495 | NumNoOps++; |
---|
| 1496 | |
---|
| 1497 | if (!clean) { |
---|
| 1498 | new_m = |
---|
| 1499 | ll_LLL_FP(B, U, delta, 0, check, B1, mu, b, c, h, h, quit); |
---|
| 1500 | if (new_m != h) Error("BKZ_FP: internal error"); |
---|
| 1501 | if (quit) break; |
---|
| 1502 | } |
---|
| 1503 | |
---|
| 1504 | z++; |
---|
| 1505 | } |
---|
| 1506 | } |
---|
| 1507 | } |
---|
| 1508 | |
---|
| 1509 | |
---|
| 1510 | // clean up |
---|
| 1511 | |
---|
| 1512 | |
---|
| 1513 | if (m_orig > m) { |
---|
| 1514 | // for consistency, we move zero vectors to the front |
---|
| 1515 | |
---|
| 1516 | for (i = m+1; i <= m_orig; i++) { |
---|
| 1517 | swap(B(i), B(i+1)); |
---|
| 1518 | if (U) swap((*U)(i), (*U)(i+1)); |
---|
| 1519 | } |
---|
| 1520 | |
---|
| 1521 | for (i = 0; i < m; i++) { |
---|
| 1522 | swap(B(m_orig-i), B(m-i)); |
---|
| 1523 | if (U) swap((*U)(m_orig-i), (*U)(m-i)); |
---|
| 1524 | } |
---|
| 1525 | } |
---|
| 1526 | |
---|
| 1527 | B.SetDims(m_orig, n); |
---|
| 1528 | BB = B; |
---|
| 1529 | |
---|
| 1530 | if (U) { |
---|
| 1531 | U->SetDims(m_orig, m_orig); |
---|
| 1532 | *UU = *U; |
---|
| 1533 | } |
---|
| 1534 | |
---|
| 1535 | for (i = 1; i <= m+1; i++) { |
---|
| 1536 | delete [] B1[i]; |
---|
| 1537 | } |
---|
| 1538 | |
---|
| 1539 | delete [] B1; |
---|
| 1540 | |
---|
| 1541 | for (i = 1; i <= m+1; i++) { |
---|
| 1542 | delete [] mu[i]; |
---|
| 1543 | } |
---|
| 1544 | |
---|
| 1545 | delete [] mu; |
---|
| 1546 | |
---|
| 1547 | delete [] c; |
---|
| 1548 | delete [] b; |
---|
| 1549 | delete [] ctilda; |
---|
| 1550 | delete [] vvec; |
---|
| 1551 | delete [] yvec; |
---|
| 1552 | delete [] uvec; |
---|
| 1553 | delete [] utildavec; |
---|
| 1554 | delete [] Deltavec; |
---|
| 1555 | delete [] deltavec; |
---|
| 1556 | |
---|
| 1557 | return m; |
---|
| 1558 | } |
---|
| 1559 | |
---|
| 1560 | long BKZ_FP(mat_ZZ& BB, mat_ZZ& UU, double delta, |
---|
| 1561 | long beta, long prune, LLLCheckFct check, long verb) |
---|
| 1562 | { |
---|
| 1563 | verbose = verb; |
---|
| 1564 | RR_GS_time = 0; |
---|
| 1565 | NumSwaps = 0; |
---|
| 1566 | if (verbose) { |
---|
| 1567 | StartTime = GetTime(); |
---|
| 1568 | LastTime = StartTime; |
---|
| 1569 | } |
---|
| 1570 | |
---|
| 1571 | if (delta < 0.50 || delta >= 1) Error("BKZ_FP: bad delta"); |
---|
| 1572 | if (beta < 2) Error("BKZ_FP: bad block size"); |
---|
| 1573 | |
---|
| 1574 | return BKZ_FP(BB, &UU, delta, beta, prune, check); |
---|
| 1575 | } |
---|
| 1576 | |
---|
| 1577 | long BKZ_FP(mat_ZZ& BB, double delta, |
---|
| 1578 | long beta, long prune, LLLCheckFct check, long verb) |
---|
| 1579 | { |
---|
| 1580 | verbose = verb; |
---|
| 1581 | RR_GS_time = 0; |
---|
| 1582 | NumSwaps = 0; |
---|
| 1583 | if (verbose) { |
---|
| 1584 | StartTime = GetTime(); |
---|
| 1585 | LastTime = StartTime; |
---|
| 1586 | } |
---|
| 1587 | |
---|
| 1588 | if (delta < 0.50 || delta >= 1) Error("BKZ_FP: bad delta"); |
---|
| 1589 | if (beta < 2) Error("BKZ_FP: bad block size"); |
---|
| 1590 | |
---|
| 1591 | return BKZ_FP(BB, 0, delta, beta, prune, check); |
---|
| 1592 | } |
---|
| 1593 | |
---|
| 1594 | |
---|
| 1595 | NTL_END_IMPL |
---|