[2cfffe] | 1 | |
---|
| 2 | |
---|
| 3 | #include <NTL/GF2EXFactoring.h> |
---|
| 4 | #include <NTL/vec_GF2XVec.h> |
---|
| 5 | #include <NTL/fileio.h> |
---|
| 6 | #include <NTL/FacVec.h> |
---|
| 7 | |
---|
| 8 | #include <stdio.h> |
---|
| 9 | |
---|
| 10 | #include <NTL/new.h> |
---|
| 11 | |
---|
| 12 | NTL_START_IMPL |
---|
| 13 | |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | |
---|
| 17 | static |
---|
| 18 | void IterSqr(GF2E& c, const GF2E& a, long n) |
---|
| 19 | { |
---|
| 20 | GF2E res; |
---|
| 21 | |
---|
| 22 | long i; |
---|
| 23 | |
---|
| 24 | res = a; |
---|
| 25 | |
---|
| 26 | for (i = 0; i < n; i++) |
---|
| 27 | sqr(res, res); |
---|
| 28 | |
---|
| 29 | c = res; |
---|
| 30 | } |
---|
| 31 | |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | void SquareFreeDecomp(vec_pair_GF2EX_long& u, const GF2EX& ff) |
---|
| 35 | { |
---|
| 36 | GF2EX f = ff; |
---|
| 37 | |
---|
| 38 | if (!IsOne(LeadCoeff(f))) |
---|
| 39 | Error("SquareFreeDecomp: bad args"); |
---|
| 40 | |
---|
| 41 | GF2EX r, t, v, tmp1; |
---|
| 42 | long m, j, finished, done; |
---|
| 43 | |
---|
| 44 | u.SetLength(0); |
---|
| 45 | |
---|
| 46 | if (deg(f) == 0) |
---|
| 47 | return; |
---|
| 48 | |
---|
| 49 | m = 1; |
---|
| 50 | finished = 0; |
---|
| 51 | |
---|
| 52 | do { |
---|
| 53 | j = 1; |
---|
| 54 | diff(tmp1, f); |
---|
| 55 | GCD(r, f, tmp1); |
---|
| 56 | div(t, f, r); |
---|
| 57 | |
---|
| 58 | if (deg(t) > 0) { |
---|
| 59 | done = 0; |
---|
| 60 | do { |
---|
| 61 | GCD(v, r, t); |
---|
| 62 | div(tmp1, t, v); |
---|
| 63 | if (deg(tmp1) > 0) append(u, cons(tmp1, j*m)); |
---|
| 64 | if (deg(v) > 0) { |
---|
| 65 | div(r, r, v); |
---|
| 66 | t = v; |
---|
| 67 | j++; |
---|
| 68 | } |
---|
| 69 | else |
---|
| 70 | done = 1; |
---|
| 71 | } while (!done); |
---|
| 72 | if (deg(r) == 0) finished = 1; |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | if (!finished) { |
---|
| 76 | /* r is a square */ |
---|
| 77 | |
---|
| 78 | long k, d; |
---|
| 79 | d = deg(r)/2; |
---|
| 80 | f.rep.SetLength(d+1); |
---|
| 81 | for (k = 0; k <= d; k++) |
---|
| 82 | IterSqr(f.rep[k], r.rep[k*2], GF2E::degree()-1); |
---|
| 83 | m = m*2; |
---|
| 84 | } |
---|
| 85 | } while (!finished); |
---|
| 86 | } |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | static |
---|
| 91 | void NullSpace(long& r, vec_long& D, vec_GF2XVec& M, long verbose) |
---|
| 92 | { |
---|
| 93 | long k, l, n; |
---|
| 94 | long i, j; |
---|
| 95 | long pos; |
---|
| 96 | GF2X t1, t2; |
---|
| 97 | GF2X *x, *y; |
---|
| 98 | |
---|
| 99 | const GF2XModulus& p = GF2E::modulus(); |
---|
| 100 | |
---|
| 101 | n = M.length(); |
---|
| 102 | |
---|
| 103 | D.SetLength(n); |
---|
| 104 | for (j = 0; j < n; j++) D[j] = -1; |
---|
| 105 | |
---|
| 106 | r = 0; |
---|
| 107 | |
---|
| 108 | l = 0; |
---|
| 109 | for (k = 0; k < n; k++) { |
---|
| 110 | |
---|
| 111 | pos = -1; |
---|
| 112 | for (i = l; i < n; i++) { |
---|
| 113 | rem(t1, M[i][k], p); |
---|
| 114 | M[i][k] = t1; |
---|
| 115 | if (pos == -1 && !IsZero(t1)) |
---|
| 116 | pos = i; |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | if (pos != -1) { |
---|
| 120 | swap(M[pos], M[l]); |
---|
| 121 | |
---|
| 122 | // make M[l, k] == -1 mod p, and make row l reduced |
---|
| 123 | |
---|
| 124 | InvMod(t1, M[l][k], p); |
---|
| 125 | for (j = k+1; j < n; j++) { |
---|
| 126 | rem(t2, M[l][j], p); |
---|
| 127 | MulMod(M[l][j], t2, t1, p); |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | for (i = l+1; i < n; i++) { |
---|
| 131 | // M[i] = M[i] + M[l]*M[i,k] |
---|
| 132 | |
---|
| 133 | t1 = M[i][k]; // this is already reduced |
---|
| 134 | |
---|
| 135 | x = M[i].elts() + (k+1); |
---|
| 136 | y = M[l].elts() + (k+1); |
---|
| 137 | |
---|
| 138 | for (j = k+1; j < n; j++, x++, y++) { |
---|
| 139 | // *x = *x + (*y)*t1 |
---|
| 140 | |
---|
| 141 | mul(t2, *y, t1); |
---|
| 142 | add(*x, *x, t2); |
---|
| 143 | } |
---|
| 144 | } |
---|
| 145 | |
---|
| 146 | D[k] = l; // variable k is defined by row l |
---|
| 147 | l++; |
---|
| 148 | |
---|
| 149 | } |
---|
| 150 | else { |
---|
| 151 | r++; |
---|
| 152 | } |
---|
| 153 | } |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | static |
---|
| 159 | void BuildMatrix(vec_GF2XVec& M, long n, const GF2EX& g, const GF2EXModulus& F, |
---|
| 160 | long verbose) |
---|
| 161 | { |
---|
| 162 | long i, j, m; |
---|
| 163 | GF2EX h; |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | M.SetLength(n); |
---|
| 167 | for (i = 0; i < n; i++) |
---|
| 168 | M[i].SetSize(n, 2*GF2E::WordLength()); |
---|
| 169 | |
---|
| 170 | set(h); |
---|
| 171 | for (j = 0; j < n; j++) { |
---|
[eb5966] | 172 | |
---|
[2cfffe] | 173 | m = deg(h); |
---|
| 174 | for (i = 0; i < n; i++) { |
---|
| 175 | if (i <= m) |
---|
| 176 | M[i][j] = rep(h.rep[i]); |
---|
| 177 | else |
---|
| 178 | clear(M[i][j]); |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | if (j < n-1) |
---|
| 182 | MulMod(h, h, g, F); |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | for (i = 0; i < n; i++) |
---|
| 186 | add(M[i][i], M[i][i], 1); |
---|
| 187 | |
---|
| 188 | } |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | static |
---|
| 192 | void TraceMap(GF2EX& h, const GF2EX& a, const GF2EXModulus& F) |
---|
| 193 | |
---|
| 194 | // one could consider making a version based on modular composition, |
---|
| 195 | // as in ComposeFrobeniusMap... |
---|
| 196 | |
---|
| 197 | { |
---|
| 198 | GF2EX res, tmp; |
---|
| 199 | |
---|
| 200 | res = a; |
---|
| 201 | tmp = a; |
---|
| 202 | |
---|
| 203 | long i; |
---|
| 204 | for (i = 0; i < GF2E::degree()-1; i++) { |
---|
| 205 | SqrMod(tmp, tmp, F); |
---|
| 206 | add(res, res, tmp); |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | h = res; |
---|
| 210 | } |
---|
| 211 | |
---|
| 212 | void PlainFrobeniusMap(GF2EX& h, const GF2EXModulus& F) |
---|
| 213 | { |
---|
| 214 | GF2EX res; |
---|
| 215 | |
---|
| 216 | SetX(res); |
---|
| 217 | long i; |
---|
| 218 | for (i = 0; i < GF2E::degree(); i++) |
---|
| 219 | SqrMod(res, res, F); |
---|
| 220 | |
---|
| 221 | h = res; |
---|
| 222 | } |
---|
| 223 | |
---|
| 224 | long UseComposeFrobenius(long d, long n) |
---|
| 225 | { |
---|
| 226 | long i; |
---|
| 227 | i = 1; |
---|
| 228 | while (i <= d) i = i << 1; |
---|
| 229 | i = i >> 1; |
---|
| 230 | |
---|
| 231 | i = i >> 1; |
---|
| 232 | long m = 1; |
---|
| 233 | |
---|
| 234 | long dz; |
---|
| 235 | |
---|
| 236 | if (n == 2) { |
---|
| 237 | dz = 1; |
---|
| 238 | } |
---|
| 239 | else { |
---|
| 240 | while (i) { |
---|
| 241 | long m1 = 2*m; |
---|
| 242 | if (i & d) m1++; |
---|
| 243 | |
---|
| 244 | if (m1 >= NTL_BITS_PER_LONG-1 || (1L << m1) >= n) break; |
---|
| 245 | |
---|
| 246 | m = m1; |
---|
| 247 | i = i >> 1; |
---|
| 248 | } |
---|
| 249 | |
---|
| 250 | dz = 1L << m; |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | long rootn = SqrRoot(n); |
---|
| 254 | long cnt = 0; |
---|
| 255 | |
---|
| 256 | if (i) { |
---|
| 257 | cnt += SqrRoot(dz+1); |
---|
| 258 | i = i >> 1; |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | while (i) { |
---|
| 262 | cnt += rootn; |
---|
| 263 | i = i >> 1; |
---|
| 264 | } |
---|
| 265 | |
---|
| 266 | return 4*cnt <= d; |
---|
| 267 | } |
---|
| 268 | |
---|
| 269 | void ComposeFrobeniusMap(GF2EX& y, const GF2EXModulus& F) |
---|
| 270 | { |
---|
| 271 | long d = GF2E::degree(); |
---|
| 272 | long n = deg(F); |
---|
| 273 | |
---|
| 274 | long i; |
---|
| 275 | i = 1; |
---|
| 276 | while (i <= d) i = i << 1; |
---|
| 277 | i = i >> 1; |
---|
| 278 | |
---|
| 279 | GF2EX z(INIT_SIZE, n), z1(INIT_SIZE, n); |
---|
| 280 | |
---|
| 281 | i = i >> 1; |
---|
| 282 | long m = 1; |
---|
| 283 | |
---|
| 284 | if (n == 2) { |
---|
| 285 | SetX(z); |
---|
| 286 | SqrMod(z, z, F); |
---|
| 287 | } |
---|
| 288 | else { |
---|
| 289 | while (i) { |
---|
| 290 | long m1 = 2*m; |
---|
| 291 | if (i & d) m1++; |
---|
| 292 | |
---|
| 293 | if (m1 >= NTL_BITS_PER_LONG-1 || (1L << m1) >= n) break; |
---|
| 294 | |
---|
| 295 | m = m1; |
---|
| 296 | i = i >> 1; |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | clear(z); |
---|
| 300 | SetCoeff(z, 1L << m); |
---|
| 301 | } |
---|
| 302 | |
---|
| 303 | |
---|
| 304 | while (i) { |
---|
| 305 | z1 = z; |
---|
| 306 | |
---|
| 307 | long j, k, dz; |
---|
| 308 | dz = deg(z); |
---|
| 309 | |
---|
| 310 | for (j = 0; j <= dz; j++) |
---|
| 311 | for (k = 0; k < m; k++) |
---|
| 312 | sqr(z1.rep[j], z1.rep[j]); |
---|
| 313 | |
---|
| 314 | CompMod(z, z1, z, F); |
---|
| 315 | m = 2*m; |
---|
| 316 | |
---|
| 317 | if (d & i) { |
---|
| 318 | SqrMod(z, z, F); |
---|
| 319 | m++; |
---|
| 320 | } |
---|
| 321 | |
---|
| 322 | i = i >> 1; |
---|
| 323 | } |
---|
| 324 | |
---|
| 325 | y = z; |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | void FrobeniusMap(GF2EX& h, const GF2EXModulus& F) |
---|
| 329 | { |
---|
| 330 | long n = deg(F); |
---|
| 331 | long d = GF2E::degree(); |
---|
| 332 | |
---|
| 333 | if (n == 1) { |
---|
| 334 | h = ConstTerm(F); |
---|
| 335 | return; |
---|
| 336 | } |
---|
| 337 | |
---|
| 338 | if (UseComposeFrobenius(d, n)) |
---|
| 339 | ComposeFrobeniusMap(h, F); |
---|
| 340 | else |
---|
| 341 | PlainFrobeniusMap(h, F); |
---|
| 342 | } |
---|
| 343 | |
---|
| 344 | |
---|
| 345 | |
---|
| 346 | |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | |
---|
| 350 | static |
---|
| 351 | void RecFindRoots(vec_GF2E& x, const GF2EX& f) |
---|
| 352 | { |
---|
| 353 | if (deg(f) == 0) return; |
---|
| 354 | |
---|
| 355 | if (deg(f) == 1) { |
---|
| 356 | long k = x.length(); |
---|
| 357 | x.SetLength(k+1); |
---|
| 358 | x[k] = ConstTerm(f); |
---|
| 359 | return; |
---|
| 360 | } |
---|
| 361 | |
---|
| 362 | GF2EX h; |
---|
| 363 | |
---|
| 364 | GF2E r; |
---|
| 365 | |
---|
| 366 | |
---|
| 367 | { |
---|
| 368 | GF2EXModulus F; |
---|
| 369 | build(F, f); |
---|
| 370 | |
---|
| 371 | do { |
---|
| 372 | random(r); |
---|
| 373 | clear(h); |
---|
| 374 | SetCoeff(h, 1, r); |
---|
| 375 | TraceMap(h, h, F); |
---|
| 376 | GCD(h, h, f); |
---|
| 377 | } while (deg(h) <= 0 || deg(h) == deg(f)); |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | RecFindRoots(x, h); |
---|
| 381 | div(h, f, h); |
---|
| 382 | RecFindRoots(x, h); |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | void FindRoots(vec_GF2E& x, const GF2EX& ff) |
---|
| 386 | { |
---|
| 387 | GF2EX f = ff; |
---|
| 388 | |
---|
| 389 | if (!IsOne(LeadCoeff(f))) |
---|
| 390 | Error("FindRoots: bad args"); |
---|
| 391 | |
---|
| 392 | x.SetMaxLength(deg(f)); |
---|
| 393 | x.SetLength(0); |
---|
| 394 | RecFindRoots(x, f); |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | |
---|
| 398 | static |
---|
| 399 | void RandomBasisElt(GF2EX& g, const vec_long& D, const vec_GF2XVec& M) |
---|
| 400 | { |
---|
| 401 | static GF2X t1, t2; |
---|
| 402 | |
---|
| 403 | long n = D.length(); |
---|
| 404 | |
---|
| 405 | long i, j, s; |
---|
| 406 | |
---|
| 407 | g.rep.SetLength(n); |
---|
| 408 | |
---|
| 409 | vec_GF2E& v = g.rep; |
---|
| 410 | |
---|
| 411 | for (j = n-1; j >= 0; j--) { |
---|
| 412 | if (D[j] == -1) |
---|
| 413 | random(v[j]); |
---|
| 414 | else { |
---|
| 415 | i = D[j]; |
---|
| 416 | |
---|
| 417 | // v[j] = sum_{s=j+1}^{n-1} v[s]*M[i,s] |
---|
| 418 | |
---|
| 419 | clear(t1); |
---|
| 420 | |
---|
| 421 | for (s = j+1; s < n; s++) { |
---|
| 422 | mul(t2, rep(v[s]), M[i][s]); |
---|
| 423 | add(t1, t1, t2); |
---|
| 424 | } |
---|
| 425 | |
---|
| 426 | conv(v[j], t1); |
---|
| 427 | } |
---|
| 428 | } |
---|
[eb5966] | 429 | |
---|
| 430 | g.normalize(); |
---|
[2cfffe] | 431 | } |
---|
| 432 | |
---|
| 433 | |
---|
| 434 | |
---|
| 435 | static |
---|
| 436 | void split(GF2EX& f1, GF2EX& g1, GF2EX& f2, GF2EX& g2, |
---|
| 437 | const GF2EX& f, const GF2EX& g, |
---|
| 438 | const vec_GF2E& roots, long lo, long mid) |
---|
| 439 | { |
---|
| 440 | long r = mid-lo+1; |
---|
| 441 | |
---|
| 442 | GF2EXModulus F; |
---|
| 443 | build(F, f); |
---|
| 444 | |
---|
| 445 | vec_GF2E lroots(INIT_SIZE, r); |
---|
| 446 | long i; |
---|
| 447 | |
---|
| 448 | for (i = 0; i < r; i++) |
---|
| 449 | lroots[i] = roots[lo+i]; |
---|
| 450 | |
---|
| 451 | |
---|
| 452 | GF2EX h, a, d; |
---|
| 453 | BuildFromRoots(h, lroots); |
---|
| 454 | CompMod(a, h, g, F); |
---|
| 455 | |
---|
| 456 | |
---|
| 457 | GCD(f1, a, f); |
---|
| 458 | |
---|
| 459 | div(f2, f, f1); |
---|
| 460 | |
---|
| 461 | rem(g1, g, f1); |
---|
| 462 | rem(g2, g, f2); |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | static |
---|
| 466 | void RecFindFactors(vec_GF2EX& factors, const GF2EX& f, const GF2EX& g, |
---|
| 467 | const vec_GF2E& roots, long lo, long hi) |
---|
| 468 | { |
---|
| 469 | long r = hi-lo+1; |
---|
| 470 | |
---|
| 471 | if (r == 0) return; |
---|
| 472 | |
---|
| 473 | if (r == 1) { |
---|
| 474 | append(factors, f); |
---|
| 475 | return; |
---|
| 476 | } |
---|
| 477 | |
---|
| 478 | GF2EX f1, g1, f2, g2; |
---|
| 479 | |
---|
| 480 | long mid = (lo+hi)/2; |
---|
| 481 | |
---|
| 482 | split(f1, g1, f2, g2, f, g, roots, lo, mid); |
---|
| 483 | |
---|
| 484 | RecFindFactors(factors, f1, g1, roots, lo, mid); |
---|
| 485 | RecFindFactors(factors, f2, g2, roots, mid+1, hi); |
---|
| 486 | } |
---|
| 487 | |
---|
| 488 | |
---|
| 489 | static |
---|
| 490 | void FindFactors(vec_GF2EX& factors, const GF2EX& f, const GF2EX& g, |
---|
| 491 | const vec_GF2E& roots) |
---|
| 492 | { |
---|
| 493 | long r = roots.length(); |
---|
| 494 | |
---|
| 495 | factors.SetMaxLength(r); |
---|
| 496 | factors.SetLength(0); |
---|
| 497 | |
---|
| 498 | RecFindFactors(factors, f, g, roots, 0, r-1); |
---|
| 499 | } |
---|
| 500 | |
---|
| 501 | #if 0 |
---|
| 502 | |
---|
| 503 | static |
---|
| 504 | void IterFindFactors(vec_GF2EX& factors, const GF2EX& f, |
---|
| 505 | const GF2EX& g, const vec_GF2E& roots) |
---|
| 506 | { |
---|
| 507 | long r = roots.length(); |
---|
| 508 | long i; |
---|
| 509 | GF2EX h; |
---|
| 510 | |
---|
| 511 | factors.SetLength(r); |
---|
| 512 | |
---|
| 513 | for (i = 0; i < r; i++) { |
---|
| 514 | add(h, g, roots[i]); |
---|
| 515 | GCD(factors[i], f, h); |
---|
| 516 | } |
---|
| 517 | } |
---|
| 518 | |
---|
| 519 | #endif |
---|
| 520 | |
---|
| 521 | |
---|
| 522 | |
---|
| 523 | |
---|
| 524 | void SFBerlekamp(vec_GF2EX& factors, const GF2EX& ff, long verbose) |
---|
| 525 | { |
---|
| 526 | GF2EX f = ff; |
---|
| 527 | |
---|
| 528 | if (!IsOne(LeadCoeff(f))) |
---|
| 529 | Error("SFBerlekamp: bad args"); |
---|
| 530 | |
---|
| 531 | if (deg(f) == 0) { |
---|
| 532 | factors.SetLength(0); |
---|
| 533 | return; |
---|
| 534 | } |
---|
| 535 | |
---|
| 536 | if (deg(f) == 1) { |
---|
| 537 | factors.SetLength(1); |
---|
| 538 | factors[0] = f; |
---|
| 539 | return; |
---|
| 540 | } |
---|
| 541 | |
---|
| 542 | double t; |
---|
| 543 | |
---|
| 544 | long n = deg(f); |
---|
| 545 | |
---|
| 546 | GF2EXModulus F; |
---|
| 547 | |
---|
| 548 | build(F, f); |
---|
| 549 | |
---|
| 550 | GF2EX g, h; |
---|
| 551 | |
---|
| 552 | FrobeniusMap(g, F); |
---|
| 553 | |
---|
| 554 | vec_long D; |
---|
| 555 | long r; |
---|
| 556 | |
---|
| 557 | vec_GF2XVec M; |
---|
| 558 | |
---|
| 559 | BuildMatrix(M, n, g, F, verbose); |
---|
| 560 | |
---|
| 561 | NullSpace(r, D, M, verbose); |
---|
| 562 | |
---|
[09da99] | 563 | |
---|
[de6a29] | 564 | |
---|
[2cfffe] | 565 | if (r == 1) { |
---|
| 566 | factors.SetLength(1); |
---|
| 567 | factors[0] = f; |
---|
| 568 | return; |
---|
| 569 | } |
---|
| 570 | |
---|
[de6a29] | 571 | |
---|
[2cfffe] | 572 | vec_GF2E roots; |
---|
| 573 | |
---|
| 574 | RandomBasisElt(g, D, M); |
---|
| 575 | MinPolyMod(h, g, F, r); |
---|
| 576 | if (deg(h) == r) M.kill(); |
---|
| 577 | FindRoots(roots, h); |
---|
| 578 | FindFactors(factors, f, g, roots); |
---|
| 579 | |
---|
| 580 | GF2EX g1; |
---|
| 581 | vec_GF2EX S, S1; |
---|
| 582 | long i; |
---|
| 583 | |
---|
| 584 | while (factors.length() < r) { |
---|
| 585 | RandomBasisElt(g, D, M); |
---|
| 586 | S.kill(); |
---|
| 587 | for (i = 0; i < factors.length(); i++) { |
---|
| 588 | const GF2EX& f = factors[i]; |
---|
| 589 | if (deg(f) == 1) { |
---|
| 590 | append(S, f); |
---|
| 591 | continue; |
---|
| 592 | } |
---|
| 593 | build(F, f); |
---|
| 594 | rem(g1, g, F); |
---|
| 595 | if (deg(g1) <= 0) { |
---|
| 596 | append(S, f); |
---|
| 597 | continue; |
---|
| 598 | } |
---|
| 599 | MinPolyMod(h, g1, F, min(deg(f), r-factors.length()+1)); |
---|
| 600 | FindRoots(roots, h); |
---|
| 601 | S1.kill(); |
---|
| 602 | FindFactors(S1, f, g1, roots); |
---|
| 603 | append(S, S1); |
---|
| 604 | } |
---|
| 605 | swap(factors, S); |
---|
| 606 | } |
---|
[eb5966] | 607 | |
---|
[de6a29] | 608 | |
---|
[2cfffe] | 609 | } |
---|
| 610 | |
---|
| 611 | |
---|
| 612 | void berlekamp(vec_pair_GF2EX_long& factors, const GF2EX& f, long verbose) |
---|
| 613 | { |
---|
| 614 | double t; |
---|
| 615 | vec_pair_GF2EX_long sfd; |
---|
| 616 | vec_GF2EX x; |
---|
| 617 | |
---|
| 618 | if (!IsOne(LeadCoeff(f))) |
---|
| 619 | Error("berlekamp: bad args"); |
---|
| 620 | |
---|
| 621 | |
---|
| 622 | SquareFreeDecomp(sfd, f); |
---|
| 623 | |
---|
| 624 | factors.SetLength(0); |
---|
| 625 | |
---|
| 626 | long i, j; |
---|
| 627 | |
---|
| 628 | for (i = 0; i < sfd.length(); i++) { |
---|
| 629 | |
---|
| 630 | SFBerlekamp(x, sfd[i].a, verbose); |
---|
| 631 | |
---|
| 632 | for (j = 0; j < x.length(); j++) |
---|
| 633 | append(factors, cons(x[j], sfd[i].b)); |
---|
| 634 | } |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | |
---|
| 638 | |
---|
| 639 | static |
---|
| 640 | void AddFactor(vec_pair_GF2EX_long& factors, const GF2EX& g, long d, long verbose) |
---|
| 641 | { |
---|
| 642 | append(factors, cons(g, d)); |
---|
| 643 | } |
---|
| 644 | |
---|
| 645 | static |
---|
| 646 | void ProcessTable(GF2EX& f, vec_pair_GF2EX_long& factors, |
---|
| 647 | const GF2EXModulus& F, long limit, const vec_GF2EX& tbl, |
---|
| 648 | long d, long verbose) |
---|
| 649 | |
---|
| 650 | { |
---|
| 651 | if (limit == 0) return; |
---|
| 652 | |
---|
| 653 | GF2EX t1; |
---|
| 654 | |
---|
| 655 | if (limit == 1) { |
---|
| 656 | GCD(t1, f, tbl[0]); |
---|
| 657 | if (deg(t1) > 0) { |
---|
| 658 | AddFactor(factors, t1, d, verbose); |
---|
| 659 | div(f, f, t1); |
---|
| 660 | } |
---|
| 661 | |
---|
| 662 | return; |
---|
| 663 | } |
---|
| 664 | |
---|
| 665 | long i; |
---|
| 666 | |
---|
| 667 | t1 = tbl[0]; |
---|
| 668 | for (i = 1; i < limit; i++) |
---|
| 669 | MulMod(t1, t1, tbl[i], F); |
---|
| 670 | |
---|
| 671 | GCD(t1, f, t1); |
---|
| 672 | |
---|
| 673 | if (deg(t1) == 0) return; |
---|
| 674 | |
---|
| 675 | div(f, f, t1); |
---|
| 676 | |
---|
| 677 | GF2EX t2; |
---|
| 678 | |
---|
| 679 | i = 0; |
---|
| 680 | d = d - limit + 1; |
---|
| 681 | |
---|
| 682 | while (2*d <= deg(t1)) { |
---|
| 683 | GCD(t2, tbl[i], t1); |
---|
| 684 | if (deg(t2) > 0) { |
---|
| 685 | AddFactor(factors, t2, d, verbose); |
---|
| 686 | div(t1, t1, t2); |
---|
| 687 | } |
---|
| 688 | |
---|
| 689 | i++; |
---|
| 690 | d++; |
---|
| 691 | } |
---|
| 692 | |
---|
| 693 | if (deg(t1) > 0) |
---|
| 694 | AddFactor(factors, t1, deg(t1), verbose); |
---|
| 695 | } |
---|
| 696 | |
---|
| 697 | |
---|
| 698 | void TraceMap(GF2EX& w, const GF2EX& a, long d, const GF2EXModulus& F, |
---|
| 699 | const GF2EX& b) |
---|
| 700 | |
---|
| 701 | { |
---|
| 702 | if (d < 0) Error("TraceMap: bad args"); |
---|
| 703 | |
---|
| 704 | GF2EX y, z, t; |
---|
| 705 | |
---|
| 706 | z = b; |
---|
| 707 | y = a; |
---|
| 708 | clear(w); |
---|
| 709 | |
---|
| 710 | while (d) { |
---|
| 711 | if (d == 1) { |
---|
| 712 | if (IsZero(w)) |
---|
| 713 | w = y; |
---|
| 714 | else { |
---|
| 715 | CompMod(w, w, z, F); |
---|
| 716 | add(w, w, y); |
---|
| 717 | } |
---|
| 718 | } |
---|
| 719 | else if ((d & 1) == 0) { |
---|
| 720 | Comp2Mod(z, t, z, y, z, F); |
---|
| 721 | add(y, t, y); |
---|
| 722 | } |
---|
| 723 | else if (IsZero(w)) { |
---|
| 724 | w = y; |
---|
| 725 | Comp2Mod(z, t, z, y, z, F); |
---|
| 726 | add(y, t, y); |
---|
| 727 | } |
---|
| 728 | else { |
---|
| 729 | Comp3Mod(z, t, w, z, y, w, z, F); |
---|
| 730 | add(w, w, y); |
---|
| 731 | add(y, t, y); |
---|
| 732 | } |
---|
| 733 | |
---|
| 734 | d = d >> 1; |
---|
| 735 | } |
---|
| 736 | } |
---|
| 737 | |
---|
| 738 | |
---|
| 739 | void PowerCompose(GF2EX& y, const GF2EX& h, long q, const GF2EXModulus& F) |
---|
| 740 | { |
---|
| 741 | if (q < 0) Error("powerCompose: bad args"); |
---|
| 742 | |
---|
| 743 | GF2EX z(INIT_SIZE, F.n); |
---|
| 744 | long sw; |
---|
| 745 | |
---|
| 746 | z = h; |
---|
| 747 | SetX(y); |
---|
| 748 | |
---|
| 749 | while (q) { |
---|
| 750 | sw = 0; |
---|
| 751 | |
---|
| 752 | if (q > 1) sw = 2; |
---|
| 753 | if (q & 1) { |
---|
| 754 | if (IsX(y)) |
---|
| 755 | y = z; |
---|
| 756 | else |
---|
| 757 | sw = sw | 1; |
---|
| 758 | } |
---|
| 759 | |
---|
| 760 | switch (sw) { |
---|
| 761 | case 0: |
---|
| 762 | break; |
---|
| 763 | |
---|
| 764 | case 1: |
---|
| 765 | CompMod(y, y, z, F); |
---|
| 766 | break; |
---|
| 767 | |
---|
| 768 | case 2: |
---|
| 769 | CompMod(z, z, z, F); |
---|
| 770 | break; |
---|
| 771 | |
---|
| 772 | case 3: |
---|
| 773 | Comp2Mod(y, z, y, z, z, F); |
---|
| 774 | break; |
---|
| 775 | } |
---|
| 776 | |
---|
| 777 | q = q >> 1; |
---|
| 778 | } |
---|
| 779 | } |
---|
| 780 | |
---|
| 781 | |
---|
| 782 | long ProbIrredTest(const GF2EX& f, long iter) |
---|
| 783 | { |
---|
| 784 | long n = deg(f); |
---|
| 785 | |
---|
| 786 | if (n <= 0) return 0; |
---|
| 787 | if (n == 1) return 1; |
---|
| 788 | |
---|
| 789 | GF2EXModulus F; |
---|
| 790 | |
---|
| 791 | build(F, f); |
---|
| 792 | |
---|
| 793 | GF2EX b, r, s; |
---|
| 794 | |
---|
| 795 | FrobeniusMap(b, F); |
---|
| 796 | |
---|
| 797 | long all_zero = 1; |
---|
| 798 | |
---|
| 799 | long i; |
---|
| 800 | |
---|
| 801 | for (i = 0; i < iter; i++) { |
---|
| 802 | random(r, n); |
---|
| 803 | TraceMap(s, r, n, F, b); |
---|
| 804 | |
---|
| 805 | all_zero = all_zero && IsZero(s); |
---|
| 806 | |
---|
| 807 | if (deg(s) > 0) return 0; |
---|
| 808 | } |
---|
| 809 | |
---|
| 810 | if (!all_zero || (n & 1)) return 1; |
---|
| 811 | |
---|
| 812 | PowerCompose(s, b, n/2, F); |
---|
| 813 | return !IsX(s); |
---|
| 814 | } |
---|
| 815 | |
---|
| 816 | |
---|
| 817 | long GF2EX_BlockingFactor = 10; |
---|
| 818 | |
---|
| 819 | void DDF(vec_pair_GF2EX_long& factors, const GF2EX& ff, const GF2EX& hh, |
---|
| 820 | long verbose) |
---|
| 821 | { |
---|
| 822 | GF2EX f = ff; |
---|
| 823 | GF2EX h = hh; |
---|
| 824 | |
---|
| 825 | if (!IsOne(LeadCoeff(f))) |
---|
| 826 | Error("DDF: bad args"); |
---|
| 827 | |
---|
| 828 | factors.SetLength(0); |
---|
| 829 | |
---|
| 830 | if (deg(f) == 0) |
---|
| 831 | return; |
---|
| 832 | |
---|
| 833 | if (deg(f) == 1) { |
---|
| 834 | AddFactor(factors, f, 1, verbose); |
---|
| 835 | return; |
---|
| 836 | } |
---|
| 837 | |
---|
| 838 | long CompTableSize = 2*SqrRoot(deg(f)); |
---|
| 839 | |
---|
| 840 | long GCDTableSize = GF2EX_BlockingFactor; |
---|
| 841 | |
---|
| 842 | GF2EXModulus F; |
---|
| 843 | build(F, f); |
---|
| 844 | |
---|
| 845 | GF2EXArgument H; |
---|
| 846 | |
---|
| 847 | build(H, h, F, min(CompTableSize, deg(f))); |
---|
| 848 | |
---|
| 849 | long i, d, limit, old_n; |
---|
| 850 | GF2EX g, X; |
---|
| 851 | |
---|
| 852 | |
---|
| 853 | vec_GF2EX tbl(INIT_SIZE, GCDTableSize); |
---|
| 854 | |
---|
| 855 | SetX(X); |
---|
| 856 | |
---|
| 857 | i = 0; |
---|
| 858 | g = h; |
---|
| 859 | d = 1; |
---|
| 860 | limit = GCDTableSize; |
---|
| 861 | |
---|
| 862 | |
---|
| 863 | while (2*d <= deg(f)) { |
---|
| 864 | |
---|
| 865 | old_n = deg(f); |
---|
| 866 | add(tbl[i], g, X); |
---|
| 867 | i++; |
---|
| 868 | if (i == limit) { |
---|
| 869 | ProcessTable(f, factors, F, i, tbl, d, verbose); |
---|
| 870 | i = 0; |
---|
| 871 | } |
---|
| 872 | |
---|
| 873 | d = d + 1; |
---|
| 874 | if (2*d <= deg(f)) { |
---|
| 875 | // we need to go further |
---|
| 876 | |
---|
| 877 | if (deg(f) < old_n) { |
---|
| 878 | // f has changed |
---|
| 879 | |
---|
| 880 | build(F, f); |
---|
| 881 | rem(h, h, f); |
---|
| 882 | rem(g, g, f); |
---|
| 883 | build(H, h, F, min(CompTableSize, deg(f))); |
---|
| 884 | } |
---|
| 885 | |
---|
| 886 | CompMod(g, g, H, F); |
---|
| 887 | } |
---|
| 888 | } |
---|
| 889 | |
---|
| 890 | ProcessTable(f, factors, F, i, tbl, d-1, verbose); |
---|
| 891 | |
---|
| 892 | if (!IsOne(f)) AddFactor(factors, f, deg(f), verbose); |
---|
| 893 | } |
---|
| 894 | |
---|
| 895 | |
---|
| 896 | |
---|
| 897 | void RootEDF(vec_GF2EX& factors, const GF2EX& f, long verbose) |
---|
| 898 | { |
---|
| 899 | vec_GF2E roots; |
---|
| 900 | double t; |
---|
| 901 | |
---|
| 902 | FindRoots(roots, f); |
---|
| 903 | |
---|
| 904 | long r = roots.length(); |
---|
| 905 | factors.SetLength(r); |
---|
| 906 | for (long j = 0; j < r; j++) { |
---|
| 907 | SetX(factors[j]); |
---|
| 908 | add(factors[j], factors[j], roots[j]); |
---|
| 909 | } |
---|
| 910 | } |
---|
| 911 | |
---|
| 912 | static |
---|
| 913 | void EDFSplit(vec_GF2EX& v, const GF2EX& f, const GF2EX& b, long d) |
---|
| 914 | { |
---|
| 915 | GF2EX a, g, h; |
---|
| 916 | GF2EXModulus F; |
---|
| 917 | vec_GF2E roots; |
---|
| 918 | |
---|
| 919 | build(F, f); |
---|
| 920 | long n = F.n; |
---|
| 921 | long r = n/d; |
---|
| 922 | random(a, n); |
---|
| 923 | TraceMap(g, a, d, F, b); |
---|
| 924 | MinPolyMod(h, g, F, r); |
---|
| 925 | FindRoots(roots, h); |
---|
| 926 | FindFactors(v, f, g, roots); |
---|
| 927 | } |
---|
| 928 | |
---|
| 929 | static |
---|
| 930 | void RecEDF(vec_GF2EX& factors, const GF2EX& f, const GF2EX& b, long d, |
---|
| 931 | long verbose) |
---|
| 932 | { |
---|
| 933 | vec_GF2EX v; |
---|
| 934 | long i; |
---|
| 935 | GF2EX bb; |
---|
| 936 | |
---|
| 937 | EDFSplit(v, f, b, d); |
---|
| 938 | for (i = 0; i < v.length(); i++) { |
---|
| 939 | if (deg(v[i]) == d) { |
---|
| 940 | append(factors, v[i]); |
---|
| 941 | } |
---|
| 942 | else { |
---|
| 943 | GF2EX bb; |
---|
| 944 | rem(bb, b, v[i]); |
---|
| 945 | RecEDF(factors, v[i], bb, d, verbose); |
---|
| 946 | } |
---|
| 947 | } |
---|
| 948 | } |
---|
| 949 | |
---|
| 950 | |
---|
| 951 | void EDF(vec_GF2EX& factors, const GF2EX& ff, const GF2EX& bb, |
---|
| 952 | long d, long verbose) |
---|
| 953 | |
---|
| 954 | { |
---|
| 955 | GF2EX f = ff; |
---|
| 956 | GF2EX b = bb; |
---|
| 957 | |
---|
| 958 | if (!IsOne(LeadCoeff(f))) |
---|
| 959 | Error("EDF: bad args"); |
---|
| 960 | |
---|
| 961 | long n = deg(f); |
---|
| 962 | long r = n/d; |
---|
| 963 | |
---|
| 964 | if (r == 0) { |
---|
| 965 | factors.SetLength(0); |
---|
| 966 | return; |
---|
| 967 | } |
---|
| 968 | |
---|
| 969 | if (r == 1) { |
---|
| 970 | factors.SetLength(1); |
---|
| 971 | factors[0] = f; |
---|
| 972 | return; |
---|
| 973 | } |
---|
| 974 | |
---|
| 975 | if (d == 1) { |
---|
| 976 | RootEDF(factors, f, verbose); |
---|
| 977 | return; |
---|
| 978 | } |
---|
| 979 | |
---|
| 980 | |
---|
| 981 | factors.SetLength(0); |
---|
| 982 | |
---|
| 983 | RecEDF(factors, f, b, d, verbose); |
---|
| 984 | |
---|
| 985 | } |
---|
| 986 | |
---|
| 987 | |
---|
| 988 | void SFCanZass(vec_GF2EX& factors, const GF2EX& ff, long verbose) |
---|
| 989 | { |
---|
| 990 | GF2EX f = ff; |
---|
| 991 | |
---|
| 992 | if (!IsOne(LeadCoeff(f))) |
---|
| 993 | Error("SFCanZass: bad args"); |
---|
| 994 | |
---|
| 995 | if (deg(f) == 0) { |
---|
| 996 | factors.SetLength(0); |
---|
| 997 | return; |
---|
| 998 | } |
---|
| 999 | |
---|
| 1000 | if (deg(f) == 1) { |
---|
| 1001 | factors.SetLength(1); |
---|
| 1002 | factors[0] = f; |
---|
| 1003 | return; |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | factors.SetLength(0); |
---|
| 1007 | |
---|
| 1008 | double t; |
---|
| 1009 | |
---|
| 1010 | |
---|
| 1011 | GF2EXModulus F; |
---|
| 1012 | build(F, f); |
---|
| 1013 | |
---|
| 1014 | GF2EX h; |
---|
| 1015 | |
---|
| 1016 | FrobeniusMap(h, F); |
---|
| 1017 | |
---|
| 1018 | vec_pair_GF2EX_long u; |
---|
| 1019 | NewDDF(u, f, h, verbose); |
---|
| 1020 | |
---|
| 1021 | GF2EX hh; |
---|
| 1022 | vec_GF2EX v; |
---|
| 1023 | |
---|
| 1024 | long i; |
---|
| 1025 | for (i = 0; i < u.length(); i++) { |
---|
| 1026 | const GF2EX& g = u[i].a; |
---|
| 1027 | long d = u[i].b; |
---|
| 1028 | long r = deg(g)/d; |
---|
| 1029 | |
---|
| 1030 | if (r == 1) { |
---|
| 1031 | // g is already irreducible |
---|
| 1032 | |
---|
| 1033 | append(factors, g); |
---|
| 1034 | } |
---|
| 1035 | else { |
---|
| 1036 | // must perform EDF |
---|
| 1037 | |
---|
| 1038 | if (d == 1) { |
---|
| 1039 | // root finding |
---|
| 1040 | RootEDF(v, g, verbose); |
---|
| 1041 | append(factors, v); |
---|
| 1042 | } |
---|
| 1043 | else { |
---|
| 1044 | // general case |
---|
| 1045 | rem(hh, h, g); |
---|
| 1046 | EDF(v, g, hh, d, verbose); |
---|
| 1047 | append(factors, v); |
---|
| 1048 | } |
---|
| 1049 | } |
---|
| 1050 | } |
---|
| 1051 | } |
---|
| 1052 | |
---|
| 1053 | void CanZass(vec_pair_GF2EX_long& factors, const GF2EX& f, long verbose) |
---|
| 1054 | { |
---|
| 1055 | if (!IsOne(LeadCoeff(f))) |
---|
| 1056 | Error("CanZass: bad args"); |
---|
| 1057 | |
---|
| 1058 | double t; |
---|
| 1059 | vec_pair_GF2EX_long sfd; |
---|
| 1060 | vec_GF2EX x; |
---|
| 1061 | |
---|
| 1062 | |
---|
| 1063 | SquareFreeDecomp(sfd, f); |
---|
| 1064 | |
---|
| 1065 | factors.SetLength(0); |
---|
| 1066 | |
---|
| 1067 | long i, j; |
---|
| 1068 | |
---|
| 1069 | for (i = 0; i < sfd.length(); i++) { |
---|
| 1070 | |
---|
| 1071 | SFCanZass(x, sfd[i].a, verbose); |
---|
| 1072 | |
---|
| 1073 | for (j = 0; j < x.length(); j++) |
---|
| 1074 | append(factors, cons(x[j], sfd[i].b)); |
---|
| 1075 | } |
---|
| 1076 | } |
---|
| 1077 | |
---|
| 1078 | void mul(GF2EX& f, const vec_pair_GF2EX_long& v) |
---|
| 1079 | { |
---|
| 1080 | long i, j, n; |
---|
| 1081 | |
---|
| 1082 | n = 0; |
---|
| 1083 | for (i = 0; i < v.length(); i++) |
---|
| 1084 | n += v[i].b*deg(v[i].a); |
---|
| 1085 | |
---|
| 1086 | GF2EX g(INIT_SIZE, n+1); |
---|
| 1087 | |
---|
| 1088 | set(g); |
---|
| 1089 | for (i = 0; i < v.length(); i++) |
---|
| 1090 | for (j = 0; j < v[i].b; j++) { |
---|
| 1091 | mul(g, g, v[i].a); |
---|
| 1092 | } |
---|
| 1093 | |
---|
| 1094 | f = g; |
---|
| 1095 | } |
---|
| 1096 | |
---|
| 1097 | |
---|
| 1098 | |
---|
| 1099 | |
---|
| 1100 | static |
---|
| 1101 | long BaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F) |
---|
| 1102 | { |
---|
| 1103 | long b, e; |
---|
| 1104 | GF2EX lh(INIT_SIZE, F.n); |
---|
| 1105 | |
---|
| 1106 | lh = h; |
---|
| 1107 | b = 1; |
---|
| 1108 | e = 0; |
---|
| 1109 | while (e < a-1 && !IsX(lh)) { |
---|
| 1110 | e++; |
---|
| 1111 | b *= q; |
---|
| 1112 | PowerCompose(lh, lh, q, F); |
---|
| 1113 | } |
---|
| 1114 | |
---|
| 1115 | if (!IsX(lh)) b *= q; |
---|
| 1116 | |
---|
| 1117 | return b; |
---|
| 1118 | } |
---|
| 1119 | |
---|
| 1120 | |
---|
| 1121 | |
---|
| 1122 | static |
---|
| 1123 | void TandemPowerCompose(GF2EX& y1, GF2EX& y2, const GF2EX& h, |
---|
| 1124 | long q1, long q2, const GF2EXModulus& F) |
---|
| 1125 | { |
---|
| 1126 | GF2EX z(INIT_SIZE, F.n); |
---|
| 1127 | long sw; |
---|
| 1128 | |
---|
| 1129 | z = h; |
---|
| 1130 | SetX(y1); |
---|
| 1131 | SetX(y2); |
---|
| 1132 | |
---|
| 1133 | while (q1 || q2) { |
---|
| 1134 | sw = 0; |
---|
| 1135 | |
---|
| 1136 | if (q1 > 1 || q2 > 1) sw = 4; |
---|
| 1137 | |
---|
| 1138 | if (q1 & 1) { |
---|
| 1139 | if (IsX(y1)) |
---|
| 1140 | y1 = z; |
---|
| 1141 | else |
---|
| 1142 | sw = sw | 2; |
---|
| 1143 | } |
---|
| 1144 | |
---|
| 1145 | if (q2 & 1) { |
---|
| 1146 | if (IsX(y2)) |
---|
| 1147 | y2 = z; |
---|
| 1148 | else |
---|
| 1149 | sw = sw | 1; |
---|
| 1150 | } |
---|
| 1151 | |
---|
| 1152 | switch (sw) { |
---|
| 1153 | case 0: |
---|
| 1154 | break; |
---|
| 1155 | |
---|
| 1156 | case 1: |
---|
| 1157 | CompMod(y2, y2, z, F); |
---|
| 1158 | break; |
---|
| 1159 | |
---|
| 1160 | case 2: |
---|
| 1161 | CompMod(y1, y1, z, F); |
---|
| 1162 | break; |
---|
| 1163 | |
---|
| 1164 | case 3: |
---|
| 1165 | Comp2Mod(y1, y2, y1, y2, z, F); |
---|
| 1166 | break; |
---|
| 1167 | |
---|
| 1168 | case 4: |
---|
| 1169 | CompMod(z, z, z, F); |
---|
| 1170 | break; |
---|
| 1171 | |
---|
| 1172 | case 5: |
---|
| 1173 | Comp2Mod(z, y2, z, y2, z, F); |
---|
| 1174 | break; |
---|
| 1175 | |
---|
| 1176 | case 6: |
---|
| 1177 | Comp2Mod(z, y1, z, y1, z, F); |
---|
| 1178 | break; |
---|
| 1179 | |
---|
| 1180 | case 7: |
---|
| 1181 | Comp3Mod(z, y1, y2, z, y1, y2, z, F); |
---|
| 1182 | break; |
---|
| 1183 | } |
---|
| 1184 | |
---|
| 1185 | q1 = q1 >> 1; |
---|
| 1186 | q2 = q2 >> 1; |
---|
| 1187 | } |
---|
| 1188 | } |
---|
| 1189 | |
---|
| 1190 | |
---|
| 1191 | |
---|
| 1192 | static |
---|
| 1193 | long RecComputeDegree(long u, const GF2EX& h, const GF2EXModulus& F, |
---|
| 1194 | FacVec& fvec) |
---|
| 1195 | { |
---|
| 1196 | if (IsX(h)) return 1; |
---|
| 1197 | |
---|
| 1198 | if (fvec[u].link == -1) return BaseCase(h, fvec[u].q, fvec[u].a, F); |
---|
| 1199 | |
---|
| 1200 | GF2EX h1, h2; |
---|
| 1201 | long q1, q2, r1, r2; |
---|
| 1202 | |
---|
| 1203 | q1 = fvec[fvec[u].link].val; |
---|
| 1204 | q2 = fvec[fvec[u].link+1].val; |
---|
| 1205 | |
---|
| 1206 | TandemPowerCompose(h1, h2, h, q1, q2, F); |
---|
| 1207 | r1 = RecComputeDegree(fvec[u].link, h2, F, fvec); |
---|
| 1208 | r2 = RecComputeDegree(fvec[u].link+1, h1, F, fvec); |
---|
| 1209 | return r1*r2; |
---|
| 1210 | } |
---|
| 1211 | |
---|
| 1212 | |
---|
| 1213 | |
---|
| 1214 | |
---|
| 1215 | long RecComputeDegree(const GF2EX& h, const GF2EXModulus& F) |
---|
| 1216 | // f = F.f is assumed to be an "equal degree" polynomial |
---|
| 1217 | // h = X^p mod f |
---|
| 1218 | // the common degree of the irreducible factors of f is computed |
---|
| 1219 | { |
---|
| 1220 | if (F.n == 1 || IsX(h)) |
---|
| 1221 | return 1; |
---|
| 1222 | |
---|
| 1223 | FacVec fvec; |
---|
| 1224 | |
---|
| 1225 | FactorInt(fvec, F.n); |
---|
| 1226 | |
---|
| 1227 | return RecComputeDegree(fvec.length()-1, h, F, fvec); |
---|
| 1228 | } |
---|
| 1229 | |
---|
| 1230 | |
---|
| 1231 | void FindRoot(GF2E& root, const GF2EX& ff) |
---|
| 1232 | // finds a root of ff. |
---|
| 1233 | // assumes that ff is monic and splits into distinct linear factors |
---|
| 1234 | |
---|
| 1235 | { |
---|
| 1236 | GF2EXModulus F; |
---|
| 1237 | GF2EX h, h1, f; |
---|
| 1238 | GF2E r; |
---|
| 1239 | |
---|
| 1240 | f = ff; |
---|
| 1241 | |
---|
| 1242 | if (!IsOne(LeadCoeff(f))) |
---|
| 1243 | Error("FindRoot: bad args"); |
---|
| 1244 | |
---|
| 1245 | if (deg(f) == 0) |
---|
| 1246 | Error("FindRoot: bad args"); |
---|
| 1247 | |
---|
| 1248 | |
---|
| 1249 | while (deg(f) > 1) { |
---|
| 1250 | build(F, f); |
---|
| 1251 | random(r); |
---|
| 1252 | clear(h); |
---|
| 1253 | SetCoeff(h, 1, r); |
---|
| 1254 | TraceMap(h, h, F); |
---|
| 1255 | GCD(h, h, f); |
---|
| 1256 | if (deg(h) > 0 && deg(h) < deg(f)) { |
---|
| 1257 | if (deg(h) > deg(f)/2) |
---|
| 1258 | div(f, f, h); |
---|
| 1259 | else |
---|
| 1260 | f = h; |
---|
| 1261 | } |
---|
| 1262 | } |
---|
| 1263 | |
---|
| 1264 | root = ConstTerm(f); |
---|
| 1265 | } |
---|
| 1266 | |
---|
| 1267 | |
---|
| 1268 | static |
---|
| 1269 | long power(long a, long e) |
---|
| 1270 | { |
---|
| 1271 | long i, res; |
---|
| 1272 | |
---|
| 1273 | res = 1; |
---|
| 1274 | for (i = 1; i <= e; i++) |
---|
| 1275 | res = res * a; |
---|
| 1276 | |
---|
| 1277 | return res; |
---|
| 1278 | } |
---|
| 1279 | |
---|
| 1280 | |
---|
| 1281 | static |
---|
| 1282 | long IrredBaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F) |
---|
| 1283 | { |
---|
| 1284 | long e; |
---|
| 1285 | GF2EX X, s, d; |
---|
| 1286 | |
---|
| 1287 | e = power(q, a-1); |
---|
| 1288 | PowerCompose(s, h, e, F); |
---|
| 1289 | SetX(X); |
---|
| 1290 | add(s, s, X); |
---|
| 1291 | GCD(d, F.f, s); |
---|
| 1292 | return IsOne(d); |
---|
| 1293 | } |
---|
| 1294 | |
---|
| 1295 | |
---|
| 1296 | static |
---|
| 1297 | long RecIrredTest(long u, const GF2EX& h, const GF2EXModulus& F, |
---|
| 1298 | const FacVec& fvec) |
---|
| 1299 | { |
---|
| 1300 | long q1, q2; |
---|
| 1301 | GF2EX h1, h2; |
---|
| 1302 | |
---|
| 1303 | if (IsX(h)) return 0; |
---|
| 1304 | |
---|
| 1305 | if (fvec[u].link == -1) { |
---|
| 1306 | return IrredBaseCase(h, fvec[u].q, fvec[u].a, F); |
---|
| 1307 | } |
---|
| 1308 | |
---|
| 1309 | |
---|
| 1310 | q1 = fvec[fvec[u].link].val; |
---|
| 1311 | q2 = fvec[fvec[u].link+1].val; |
---|
| 1312 | |
---|
| 1313 | TandemPowerCompose(h1, h2, h, q1, q2, F); |
---|
| 1314 | return RecIrredTest(fvec[u].link, h2, F, fvec) |
---|
| 1315 | && RecIrredTest(fvec[u].link+1, h1, F, fvec); |
---|
| 1316 | } |
---|
| 1317 | |
---|
| 1318 | long DetIrredTest(const GF2EX& f) |
---|
| 1319 | { |
---|
| 1320 | if (deg(f) <= 0) return 0; |
---|
| 1321 | if (deg(f) == 1) return 1; |
---|
| 1322 | |
---|
| 1323 | GF2EXModulus F; |
---|
| 1324 | |
---|
| 1325 | build(F, f); |
---|
| 1326 | |
---|
| 1327 | GF2EX h; |
---|
| 1328 | |
---|
| 1329 | FrobeniusMap(h, F); |
---|
| 1330 | |
---|
| 1331 | GF2EX s; |
---|
| 1332 | PowerCompose(s, h, F.n, F); |
---|
| 1333 | if (!IsX(s)) return 0; |
---|
| 1334 | |
---|
| 1335 | FacVec fvec; |
---|
| 1336 | |
---|
| 1337 | FactorInt(fvec, F.n); |
---|
| 1338 | |
---|
| 1339 | return RecIrredTest(fvec.length()-1, h, F, fvec); |
---|
| 1340 | } |
---|
| 1341 | |
---|
| 1342 | |
---|
| 1343 | |
---|
| 1344 | long IterIrredTest(const GF2EX& f) |
---|
| 1345 | { |
---|
| 1346 | if (deg(f) <= 0) return 0; |
---|
| 1347 | if (deg(f) == 1) return 1; |
---|
| 1348 | |
---|
| 1349 | GF2EXModulus F; |
---|
| 1350 | |
---|
| 1351 | build(F, f); |
---|
| 1352 | |
---|
| 1353 | GF2EX h; |
---|
| 1354 | |
---|
| 1355 | FrobeniusMap(h, F); |
---|
| 1356 | |
---|
| 1357 | long CompTableSize = 2*SqrRoot(deg(f)); |
---|
| 1358 | |
---|
| 1359 | GF2EXArgument H; |
---|
| 1360 | |
---|
| 1361 | build(H, h, F, CompTableSize); |
---|
| 1362 | |
---|
| 1363 | long i, d, limit, limit_sqr; |
---|
| 1364 | GF2EX g, X, t, prod; |
---|
| 1365 | |
---|
| 1366 | |
---|
| 1367 | SetX(X); |
---|
| 1368 | |
---|
| 1369 | i = 0; |
---|
| 1370 | g = h; |
---|
| 1371 | d = 1; |
---|
| 1372 | limit = 2; |
---|
| 1373 | limit_sqr = limit*limit; |
---|
| 1374 | |
---|
| 1375 | set(prod); |
---|
| 1376 | |
---|
| 1377 | |
---|
| 1378 | while (2*d <= deg(f)) { |
---|
| 1379 | add(t, g, X); |
---|
| 1380 | MulMod(prod, prod, t, F); |
---|
| 1381 | i++; |
---|
| 1382 | if (i == limit_sqr) { |
---|
| 1383 | GCD(t, f, prod); |
---|
| 1384 | if (!IsOne(t)) return 0; |
---|
| 1385 | |
---|
| 1386 | set(prod); |
---|
| 1387 | limit++; |
---|
| 1388 | limit_sqr = limit*limit; |
---|
| 1389 | i = 0; |
---|
| 1390 | } |
---|
| 1391 | |
---|
| 1392 | d = d + 1; |
---|
| 1393 | if (2*d <= deg(f)) { |
---|
| 1394 | CompMod(g, g, H, F); |
---|
| 1395 | } |
---|
| 1396 | } |
---|
| 1397 | |
---|
| 1398 | if (i > 0) { |
---|
| 1399 | GCD(t, f, prod); |
---|
| 1400 | if (!IsOne(t)) return 0; |
---|
| 1401 | } |
---|
| 1402 | |
---|
| 1403 | return 1; |
---|
| 1404 | } |
---|
| 1405 | |
---|
| 1406 | static |
---|
| 1407 | void MulByXPlusY(vec_GF2EX& h, const GF2EX& f, const GF2EX& g) |
---|
| 1408 | // h represents the bivariate polynomial h[0] + h[1]*Y + ... + h[n-1]*Y^k, |
---|
| 1409 | // where the h[i]'s are polynomials in X, each of degree < deg(f), |
---|
| 1410 | // and k < deg(g). |
---|
| 1411 | // h is replaced by the bivariate polynomial h*(X+Y) (mod f(X), g(Y)). |
---|
| 1412 | |
---|
| 1413 | { |
---|
| 1414 | long n = deg(g); |
---|
| 1415 | long k = h.length()-1; |
---|
| 1416 | |
---|
| 1417 | if (k < 0) return; |
---|
| 1418 | |
---|
| 1419 | if (k < n-1) { |
---|
| 1420 | h.SetLength(k+2); |
---|
| 1421 | h[k+1] = h[k]; |
---|
| 1422 | for (long i = k; i >= 1; i--) { |
---|
| 1423 | MulByXMod(h[i], h[i], f); |
---|
| 1424 | add(h[i], h[i], h[i-1]); |
---|
| 1425 | } |
---|
| 1426 | MulByXMod(h[0], h[0], f); |
---|
| 1427 | } |
---|
| 1428 | else { |
---|
| 1429 | GF2EX b, t; |
---|
| 1430 | |
---|
| 1431 | b = h[n-1]; |
---|
| 1432 | for (long i = n-1; i >= 1; i--) { |
---|
| 1433 | mul(t, b, g.rep[i]); |
---|
| 1434 | MulByXMod(h[i], h[i], f); |
---|
| 1435 | add(h[i], h[i], h[i-1]); |
---|
| 1436 | add(h[i], h[i], t); |
---|
| 1437 | } |
---|
| 1438 | mul(t, b, g.rep[0]); |
---|
| 1439 | MulByXMod(h[0], h[0], f); |
---|
| 1440 | add(h[0], h[0], t); |
---|
| 1441 | } |
---|
| 1442 | |
---|
| 1443 | // normalize |
---|
| 1444 | |
---|
| 1445 | k = h.length()-1; |
---|
| 1446 | while (k >= 0 && IsZero(h[k])) k--; |
---|
| 1447 | h.SetLength(k+1); |
---|
| 1448 | } |
---|
| 1449 | |
---|
| 1450 | |
---|
| 1451 | static |
---|
| 1452 | void IrredCombine(GF2EX& x, const GF2EX& f, const GF2EX& g) |
---|
| 1453 | { |
---|
| 1454 | if (deg(f) < deg(g)) { |
---|
| 1455 | IrredCombine(x, g, f); |
---|
| 1456 | return; |
---|
| 1457 | } |
---|
| 1458 | |
---|
| 1459 | // deg(f) >= deg(g)...not necessary, but maybe a little more |
---|
| 1460 | // time & space efficient |
---|
| 1461 | |
---|
| 1462 | long df = deg(f); |
---|
| 1463 | long dg = deg(g); |
---|
| 1464 | long m = df*dg; |
---|
| 1465 | |
---|
| 1466 | vec_GF2EX h(INIT_SIZE, dg); |
---|
| 1467 | |
---|
| 1468 | long i; |
---|
| 1469 | for (i = 0; i < dg; i++) h[i].SetMaxLength(df); |
---|
| 1470 | |
---|
| 1471 | h.SetLength(1); |
---|
| 1472 | set(h[0]); |
---|
| 1473 | |
---|
| 1474 | vec_GF2E a; |
---|
| 1475 | |
---|
| 1476 | a.SetLength(2*m); |
---|
| 1477 | |
---|
| 1478 | for (i = 0; i < 2*m; i++) { |
---|
| 1479 | a[i] = ConstTerm(h[0]); |
---|
| 1480 | if (i < 2*m-1) |
---|
| 1481 | MulByXPlusY(h, f, g); |
---|
| 1482 | } |
---|
| 1483 | |
---|
| 1484 | MinPolySeq(x, a, m); |
---|
| 1485 | } |
---|
| 1486 | |
---|
| 1487 | |
---|
| 1488 | static |
---|
| 1489 | void BuildPrimePowerIrred(GF2EX& f, long q, long e) |
---|
| 1490 | { |
---|
| 1491 | long n = power(q, e); |
---|
| 1492 | |
---|
| 1493 | do { |
---|
| 1494 | random(f, n); |
---|
| 1495 | SetCoeff(f, n); |
---|
| 1496 | } while (!IterIrredTest(f)); |
---|
| 1497 | } |
---|
| 1498 | |
---|
| 1499 | static |
---|
| 1500 | void RecBuildIrred(GF2EX& f, long u, const FacVec& fvec) |
---|
| 1501 | { |
---|
| 1502 | if (fvec[u].link == -1) |
---|
| 1503 | BuildPrimePowerIrred(f, fvec[u].q, fvec[u].a); |
---|
| 1504 | else { |
---|
| 1505 | GF2EX g, h; |
---|
| 1506 | RecBuildIrred(g, fvec[u].link, fvec); |
---|
| 1507 | RecBuildIrred(h, fvec[u].link+1, fvec); |
---|
| 1508 | IrredCombine(f, g, h); |
---|
| 1509 | } |
---|
| 1510 | } |
---|
| 1511 | |
---|
| 1512 | |
---|
| 1513 | void BuildIrred(GF2EX& f, long n) |
---|
| 1514 | { |
---|
| 1515 | if (n <= 0) |
---|
| 1516 | Error("BuildIrred: n must be positive"); |
---|
| 1517 | |
---|
[09da99] | 1518 | if (NTL_OVERFLOW(n, 1, 0)) |
---|
[2cfffe] | 1519 | Error("overflow in BuildIrred"); |
---|
| 1520 | |
---|
| 1521 | if (n == 1) { |
---|
| 1522 | SetX(f); |
---|
| 1523 | return; |
---|
| 1524 | } |
---|
| 1525 | |
---|
| 1526 | FacVec fvec; |
---|
| 1527 | |
---|
| 1528 | FactorInt(fvec, n); |
---|
| 1529 | |
---|
| 1530 | RecBuildIrred(f, fvec.length()-1, fvec); |
---|
| 1531 | } |
---|
| 1532 | |
---|
| 1533 | |
---|
| 1534 | |
---|
| 1535 | #if 0 |
---|
| 1536 | void BuildIrred(GF2EX& f, long n) |
---|
| 1537 | { |
---|
| 1538 | if (n <= 0) |
---|
| 1539 | Error("BuildIrred: n must be positive"); |
---|
| 1540 | |
---|
[09da99] | 1541 | if (NTL_OVERFLOW(n, 1, 0)) |
---|
[2cfffe] | 1542 | Error("overflow in BuildIrred"); |
---|
| 1543 | |
---|
| 1544 | if (n == 1) { |
---|
| 1545 | SetX(f); |
---|
| 1546 | return; |
---|
| 1547 | } |
---|
| 1548 | |
---|
| 1549 | GF2EX g; |
---|
| 1550 | |
---|
| 1551 | do { |
---|
| 1552 | random(g, n); |
---|
| 1553 | SetCoeff(g, n); |
---|
| 1554 | } while (!IterIrredTest(g)); |
---|
| 1555 | |
---|
| 1556 | f = g; |
---|
| 1557 | |
---|
| 1558 | } |
---|
| 1559 | #endif |
---|
| 1560 | |
---|
| 1561 | |
---|
| 1562 | |
---|
| 1563 | void BuildRandomIrred(GF2EX& f, const GF2EX& g) |
---|
| 1564 | { |
---|
| 1565 | GF2EXModulus G; |
---|
| 1566 | GF2EX h, ff; |
---|
| 1567 | |
---|
| 1568 | build(G, g); |
---|
| 1569 | do { |
---|
| 1570 | random(h, deg(g)); |
---|
| 1571 | IrredPolyMod(ff, h, G); |
---|
| 1572 | } while (deg(ff) < deg(g)); |
---|
| 1573 | |
---|
| 1574 | f = ff; |
---|
| 1575 | } |
---|
| 1576 | |
---|
| 1577 | |
---|
| 1578 | /************* NEW DDF ****************/ |
---|
| 1579 | |
---|
| 1580 | long GF2EX_GCDTableSize = 4; |
---|
| 1581 | char GF2EX_stem[256] = ""; |
---|
| 1582 | |
---|
| 1583 | double GF2EXFileThresh = 256; |
---|
| 1584 | |
---|
| 1585 | static vec_GF2EX BabyStepFile; |
---|
| 1586 | static vec_GF2EX GiantStepFile; |
---|
| 1587 | |
---|
| 1588 | static |
---|
| 1589 | double CalcTableSize(long n, long k) |
---|
| 1590 | { |
---|
[09da99] | 1591 | double sz = GF2E::storage(); |
---|
[2cfffe] | 1592 | sz = sz * n; |
---|
| 1593 | sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_GF2E); |
---|
| 1594 | sz = sz * k; |
---|
| 1595 | sz = sz/1024; |
---|
| 1596 | return sz; |
---|
| 1597 | } |
---|
| 1598 | |
---|
| 1599 | |
---|
| 1600 | |
---|
| 1601 | static |
---|
| 1602 | void GenerateBabySteps(GF2EX& h1, const GF2EX& f, const GF2EX& h, long k, |
---|
| 1603 | long verbose) |
---|
| 1604 | |
---|
| 1605 | { |
---|
[de6a29] | 1606 | double t; |
---|
[eb5966] | 1607 | |
---|
[2cfffe] | 1608 | GF2EXModulus F; |
---|
| 1609 | build(F, f); |
---|
| 1610 | |
---|
| 1611 | GF2EXArgument H; |
---|
| 1612 | |
---|
| 1613 | #if 0 |
---|
| 1614 | double n2 = sqrt(double(F.n)); |
---|
| 1615 | double n4 = sqrt(n2); |
---|
| 1616 | double n34 = n2*n4; |
---|
| 1617 | long sz = long(ceil(n34/sqrt(sqrt(2.0)))); |
---|
| 1618 | #else |
---|
| 1619 | long sz = 2*SqrRoot(F.n); |
---|
| 1620 | #endif |
---|
| 1621 | |
---|
| 1622 | build(H, h, F, sz); |
---|
| 1623 | |
---|
| 1624 | |
---|
| 1625 | h1 = h; |
---|
| 1626 | |
---|
| 1627 | long i; |
---|
| 1628 | |
---|
| 1629 | long HexOutput = GF2X::HexOutput; |
---|
| 1630 | GF2X::HexOutput = 1; |
---|
| 1631 | |
---|
| 1632 | BabyStepFile.kill(); |
---|
| 1633 | BabyStepFile.SetLength(k-1); |
---|
| 1634 | |
---|
| 1635 | for (i = 1; i <= k-1; i++) { |
---|
[eb5966] | 1636 | BabyStepFile(i) = h1; |
---|
[2cfffe] | 1637 | |
---|
| 1638 | CompMod(h1, h1, H, F); |
---|
| 1639 | } |
---|
| 1640 | |
---|
| 1641 | GF2X::HexOutput = HexOutput; |
---|
| 1642 | } |
---|
| 1643 | |
---|
| 1644 | |
---|
| 1645 | static |
---|
| 1646 | void GenerateGiantSteps(const GF2EX& f, const GF2EX& h, long l, long verbose) |
---|
| 1647 | { |
---|
| 1648 | |
---|
[de6a29] | 1649 | double t; |
---|
| 1650 | |
---|
| 1651 | |
---|
[2cfffe] | 1652 | GF2EXModulus F; |
---|
| 1653 | build(F, f); |
---|
| 1654 | |
---|
| 1655 | GF2EXArgument H; |
---|
| 1656 | |
---|
| 1657 | #if 0 |
---|
| 1658 | double n2 = sqrt(double(F.n)); |
---|
| 1659 | double n4 = sqrt(n2); |
---|
| 1660 | double n34 = n2*n4; |
---|
| 1661 | long sz = long(ceil(n34/sqrt(sqrt(2.0)))); |
---|
| 1662 | #else |
---|
| 1663 | long sz = 2*SqrRoot(F.n); |
---|
| 1664 | #endif |
---|
| 1665 | |
---|
| 1666 | build(H, h, F, sz); |
---|
| 1667 | |
---|
| 1668 | GF2EX h1; |
---|
| 1669 | |
---|
| 1670 | h1 = h; |
---|
| 1671 | |
---|
| 1672 | long i; |
---|
| 1673 | |
---|
| 1674 | long HexOutput = GF2X::HexOutput; |
---|
| 1675 | GF2X::HexOutput = 1; |
---|
| 1676 | |
---|
| 1677 | GiantStepFile.kill(); |
---|
| 1678 | GiantStepFile.SetLength(l); |
---|
| 1679 | |
---|
| 1680 | for (i = 1; i <= l-1; i++) { |
---|
[eb5966] | 1681 | GiantStepFile(i) = h1; |
---|
[2cfffe] | 1682 | |
---|
| 1683 | CompMod(h1, h1, H, F); |
---|
[eb5966] | 1684 | } |
---|
| 1685 | |
---|
| 1686 | GiantStepFile(i) = h1; |
---|
[2cfffe] | 1687 | |
---|
| 1688 | GF2X::HexOutput = HexOutput; |
---|
| 1689 | } |
---|
| 1690 | |
---|
| 1691 | static |
---|
| 1692 | void FileCleanup(long k, long l) |
---|
| 1693 | { |
---|
| 1694 | BabyStepFile.kill(); |
---|
| 1695 | GiantStepFile.kill(); |
---|
| 1696 | } |
---|
| 1697 | |
---|
| 1698 | |
---|
| 1699 | static |
---|
| 1700 | void NewAddFactor(vec_pair_GF2EX_long& u, const GF2EX& g, long m, long verbose) |
---|
| 1701 | { |
---|
| 1702 | long len = u.length(); |
---|
| 1703 | |
---|
| 1704 | u.SetLength(len+1); |
---|
| 1705 | u[len].a = g; |
---|
| 1706 | u[len].b = m; |
---|
| 1707 | |
---|
| 1708 | } |
---|
| 1709 | |
---|
| 1710 | |
---|
| 1711 | |
---|
| 1712 | |
---|
| 1713 | static |
---|
| 1714 | void NewProcessTable(vec_pair_GF2EX_long& u, GF2EX& f, const GF2EXModulus& F, |
---|
| 1715 | vec_GF2EX& buf, long size, long StartInterval, |
---|
| 1716 | long IntervalLength, long verbose) |
---|
| 1717 | |
---|
| 1718 | { |
---|
| 1719 | if (size == 0) return; |
---|
| 1720 | |
---|
| 1721 | GF2EX& g = buf[size-1]; |
---|
| 1722 | |
---|
| 1723 | long i; |
---|
| 1724 | |
---|
| 1725 | for (i = 0; i < size-1; i++) |
---|
| 1726 | MulMod(g, g, buf[i], F); |
---|
| 1727 | |
---|
| 1728 | GCD(g, f, g); |
---|
| 1729 | |
---|
| 1730 | if (deg(g) == 0) return; |
---|
| 1731 | |
---|
| 1732 | div(f, f, g); |
---|
| 1733 | |
---|
| 1734 | long d = (StartInterval-1)*IntervalLength + 1; |
---|
| 1735 | i = 0; |
---|
| 1736 | long interval = StartInterval; |
---|
| 1737 | |
---|
| 1738 | while (i < size-1 && 2*d <= deg(g)) { |
---|
| 1739 | GCD(buf[i], buf[i], g); |
---|
| 1740 | if (deg(buf[i]) > 0) { |
---|
| 1741 | NewAddFactor(u, buf[i], interval, verbose); |
---|
| 1742 | div(g, g, buf[i]); |
---|
| 1743 | } |
---|
| 1744 | |
---|
| 1745 | i++; |
---|
| 1746 | interval++; |
---|
| 1747 | d += IntervalLength; |
---|
| 1748 | } |
---|
| 1749 | |
---|
| 1750 | if (deg(g) > 0) { |
---|
| 1751 | if (i == size-1) |
---|
| 1752 | NewAddFactor(u, g, interval, verbose); |
---|
| 1753 | else |
---|
| 1754 | NewAddFactor(u, g, (deg(g)+IntervalLength-1)/IntervalLength, verbose); |
---|
| 1755 | } |
---|
| 1756 | } |
---|
| 1757 | |
---|
| 1758 | |
---|
| 1759 | static |
---|
| 1760 | void FetchGiantStep(GF2EX& g, long gs, const GF2EXModulus& F) |
---|
| 1761 | { |
---|
[de6a29] | 1762 | g = GiantStepFile(gs); |
---|
[2cfffe] | 1763 | |
---|
| 1764 | rem(g, g, F); |
---|
| 1765 | } |
---|
| 1766 | |
---|
| 1767 | |
---|
| 1768 | static |
---|
| 1769 | void FetchBabySteps(vec_GF2EX& v, long k) |
---|
| 1770 | { |
---|
| 1771 | v.SetLength(k); |
---|
| 1772 | |
---|
| 1773 | SetX(v[0]); |
---|
| 1774 | |
---|
| 1775 | long i; |
---|
| 1776 | for (i = 1; i <= k-1; i++) { |
---|
[eb5966] | 1777 | v[i] = BabyStepFile(i); |
---|
[2cfffe] | 1778 | } |
---|
| 1779 | } |
---|
| 1780 | |
---|
| 1781 | |
---|
| 1782 | |
---|
| 1783 | static |
---|
| 1784 | void GiantRefine(vec_pair_GF2EX_long& u, const GF2EX& ff, long k, long l, |
---|
| 1785 | long verbose) |
---|
| 1786 | |
---|
| 1787 | { |
---|
| 1788 | u.SetLength(0); |
---|
| 1789 | |
---|
| 1790 | vec_GF2EX BabyStep; |
---|
| 1791 | |
---|
| 1792 | FetchBabySteps(BabyStep, k); |
---|
| 1793 | |
---|
| 1794 | vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize); |
---|
| 1795 | |
---|
| 1796 | GF2EX f; |
---|
| 1797 | f = ff; |
---|
| 1798 | |
---|
| 1799 | GF2EXModulus F; |
---|
| 1800 | build(F, f); |
---|
| 1801 | |
---|
| 1802 | GF2EX g; |
---|
| 1803 | GF2EX h; |
---|
| 1804 | |
---|
| 1805 | long size = 0; |
---|
| 1806 | |
---|
| 1807 | long first_gs; |
---|
| 1808 | |
---|
| 1809 | long d = 1; |
---|
| 1810 | |
---|
| 1811 | while (2*d <= deg(f)) { |
---|
| 1812 | |
---|
| 1813 | long old_n = deg(f); |
---|
| 1814 | |
---|
| 1815 | long gs = (d+k-1)/k; |
---|
| 1816 | long bs = gs*k - d; |
---|
| 1817 | |
---|
| 1818 | if (bs == k-1) { |
---|
| 1819 | size++; |
---|
| 1820 | if (size == 1) first_gs = gs; |
---|
| 1821 | FetchGiantStep(g, gs, F); |
---|
| 1822 | add(buf[size-1], g, BabyStep[bs]); |
---|
| 1823 | } |
---|
| 1824 | else { |
---|
| 1825 | add(h, g, BabyStep[bs]); |
---|
| 1826 | MulMod(buf[size-1], buf[size-1], h, F); |
---|
| 1827 | } |
---|
| 1828 | |
---|
| 1829 | if (size == GF2EX_GCDTableSize && bs == 0) { |
---|
| 1830 | NewProcessTable(u, f, F, buf, size, first_gs, k, verbose); |
---|
| 1831 | size = 0; |
---|
| 1832 | } |
---|
| 1833 | |
---|
| 1834 | d++; |
---|
| 1835 | |
---|
| 1836 | if (2*d <= deg(f) && deg(f) < old_n) { |
---|
| 1837 | build(F, f); |
---|
| 1838 | |
---|
| 1839 | long i; |
---|
| 1840 | for (i = 1; i <= k-1; i++) |
---|
| 1841 | rem(BabyStep[i], BabyStep[i], F); |
---|
| 1842 | } |
---|
| 1843 | } |
---|
| 1844 | |
---|
| 1845 | if (size > 0) { |
---|
| 1846 | NewProcessTable(u, f, F, buf, size, first_gs, k, verbose); |
---|
| 1847 | } |
---|
| 1848 | |
---|
| 1849 | if (deg(f) > 0) |
---|
| 1850 | NewAddFactor(u, f, 0, verbose); |
---|
| 1851 | |
---|
| 1852 | } |
---|
| 1853 | |
---|
| 1854 | |
---|
| 1855 | static |
---|
| 1856 | void IntervalRefine(vec_pair_GF2EX_long& factors, const GF2EX& ff, |
---|
| 1857 | long k, long gs, const vec_GF2EX& BabyStep, long verbose) |
---|
| 1858 | |
---|
| 1859 | { |
---|
| 1860 | vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize); |
---|
| 1861 | |
---|
| 1862 | GF2EX f; |
---|
| 1863 | f = ff; |
---|
| 1864 | |
---|
| 1865 | GF2EXModulus F; |
---|
| 1866 | build(F, f); |
---|
| 1867 | |
---|
| 1868 | GF2EX g; |
---|
| 1869 | |
---|
| 1870 | FetchGiantStep(g, gs, F); |
---|
| 1871 | |
---|
| 1872 | long size = 0; |
---|
| 1873 | |
---|
| 1874 | long first_d; |
---|
| 1875 | |
---|
| 1876 | long d = (gs-1)*k + 1; |
---|
| 1877 | long bs = k-1; |
---|
| 1878 | |
---|
| 1879 | while (bs >= 0 && 2*d <= deg(f)) { |
---|
| 1880 | |
---|
| 1881 | long old_n = deg(f); |
---|
| 1882 | |
---|
| 1883 | if (size == 0) first_d = d; |
---|
| 1884 | rem(buf[size], BabyStep[bs], F); |
---|
| 1885 | add(buf[size], buf[size], g); |
---|
| 1886 | size++; |
---|
| 1887 | |
---|
| 1888 | if (size == GF2EX_GCDTableSize) { |
---|
| 1889 | NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose); |
---|
| 1890 | size = 0; |
---|
| 1891 | } |
---|
| 1892 | |
---|
| 1893 | d++; |
---|
| 1894 | bs--; |
---|
| 1895 | |
---|
| 1896 | if (bs >= 0 && 2*d <= deg(f) && deg(f) < old_n) { |
---|
| 1897 | build(F, f); |
---|
| 1898 | rem(g, g, F); |
---|
| 1899 | } |
---|
| 1900 | } |
---|
| 1901 | |
---|
| 1902 | NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose); |
---|
| 1903 | |
---|
| 1904 | if (deg(f) > 0) |
---|
| 1905 | NewAddFactor(factors, f, deg(f), verbose); |
---|
| 1906 | } |
---|
| 1907 | |
---|
| 1908 | |
---|
| 1909 | |
---|
| 1910 | |
---|
| 1911 | static |
---|
| 1912 | void BabyRefine(vec_pair_GF2EX_long& factors, const vec_pair_GF2EX_long& u, |
---|
| 1913 | long k, long l, long verbose) |
---|
| 1914 | |
---|
| 1915 | { |
---|
[eb5966] | 1916 | |
---|
[2cfffe] | 1917 | factors.SetLength(0); |
---|
| 1918 | |
---|
| 1919 | vec_GF2EX BabyStep; |
---|
| 1920 | |
---|
| 1921 | long i; |
---|
| 1922 | for (i = 0; i < u.length(); i++) { |
---|
| 1923 | const GF2EX& g = u[i].a; |
---|
| 1924 | long gs = u[i].b; |
---|
| 1925 | |
---|
| 1926 | if (gs == 0 || 2*((gs-1)*k+1) > deg(g)) |
---|
| 1927 | NewAddFactor(factors, g, deg(g), verbose); |
---|
| 1928 | else { |
---|
| 1929 | if (BabyStep.length() == 0) |
---|
| 1930 | FetchBabySteps(BabyStep, k); |
---|
| 1931 | IntervalRefine(factors, g, k, gs, BabyStep, verbose); |
---|
| 1932 | } |
---|
| 1933 | } |
---|
| 1934 | |
---|
| 1935 | } |
---|
| 1936 | |
---|
| 1937 | |
---|
| 1938 | |
---|
| 1939 | |
---|
| 1940 | |
---|
| 1941 | |
---|
| 1942 | void NewDDF(vec_pair_GF2EX_long& factors, |
---|
| 1943 | const GF2EX& f, |
---|
| 1944 | const GF2EX& h, |
---|
| 1945 | long verbose) |
---|
| 1946 | |
---|
| 1947 | { |
---|
| 1948 | if (!IsOne(LeadCoeff(f))) |
---|
| 1949 | Error("NewDDF: bad args"); |
---|
| 1950 | |
---|
| 1951 | if (deg(f) == 0) { |
---|
| 1952 | factors.SetLength(0); |
---|
| 1953 | return; |
---|
| 1954 | } |
---|
| 1955 | |
---|
| 1956 | if (deg(f) == 1) { |
---|
| 1957 | factors.SetLength(0); |
---|
| 1958 | append(factors, cons(f, 1)); |
---|
| 1959 | return; |
---|
| 1960 | } |
---|
| 1961 | |
---|
| 1962 | if (!GF2EX_stem[0]) |
---|
| 1963 | sprintf(GF2EX_stem, "ddf-%ld", RandomBnd(10000)); |
---|
| 1964 | |
---|
| 1965 | long B = deg(f)/2; |
---|
| 1966 | long k = SqrRoot(B); |
---|
| 1967 | long l = (B+k-1)/k; |
---|
| 1968 | |
---|
| 1969 | GF2EX h1; |
---|
| 1970 | |
---|
| 1971 | GenerateBabySteps(h1, f, h, k, verbose); |
---|
| 1972 | |
---|
| 1973 | GenerateGiantSteps(f, h1, l, verbose); |
---|
| 1974 | |
---|
| 1975 | vec_pair_GF2EX_long u; |
---|
| 1976 | GiantRefine(u, f, k, l, verbose); |
---|
| 1977 | BabyRefine(factors, u, k, l, verbose); |
---|
| 1978 | |
---|
| 1979 | FileCleanup(k, l); |
---|
| 1980 | } |
---|
| 1981 | |
---|
| 1982 | long IterComputeDegree(const GF2EX& h, const GF2EXModulus& F) |
---|
| 1983 | { |
---|
| 1984 | long n = deg(F); |
---|
| 1985 | |
---|
| 1986 | if (n == 1 || IsX(h)) return 1; |
---|
| 1987 | |
---|
| 1988 | long B = n/2; |
---|
| 1989 | long k = SqrRoot(B); |
---|
| 1990 | long l = (B+k-1)/k; |
---|
| 1991 | |
---|
| 1992 | |
---|
| 1993 | GF2EXArgument H; |
---|
| 1994 | |
---|
| 1995 | #if 0 |
---|
| 1996 | double n2 = sqrt(double(n)); |
---|
| 1997 | double n4 = sqrt(n2); |
---|
| 1998 | double n34 = n2*n4; |
---|
| 1999 | long sz = long(ceil(n34/sqrt(sqrt(2.0)))); |
---|
| 2000 | #else |
---|
| 2001 | long sz = 2*SqrRoot(F.n); |
---|
| 2002 | #endif |
---|
| 2003 | |
---|
| 2004 | build(H, h, F, sz); |
---|
| 2005 | |
---|
| 2006 | GF2EX h1; |
---|
| 2007 | h1 = h; |
---|
| 2008 | |
---|
| 2009 | vec_GF2EX baby; |
---|
| 2010 | baby.SetLength(k); |
---|
| 2011 | |
---|
| 2012 | SetX(baby[0]); |
---|
| 2013 | |
---|
| 2014 | long i; |
---|
| 2015 | |
---|
| 2016 | for (i = 1; i <= k-1; i++) { |
---|
| 2017 | baby[i] = h1; |
---|
| 2018 | CompMod(h1, h1, H, F); |
---|
| 2019 | if (IsX(h1)) return i+1; |
---|
| 2020 | } |
---|
| 2021 | |
---|
| 2022 | build(H, h1, F, sz); |
---|
| 2023 | |
---|
| 2024 | long j; |
---|
| 2025 | |
---|
| 2026 | for (j = 2; j <= l; j++) { |
---|
| 2027 | CompMod(h1, h1, H, F); |
---|
| 2028 | |
---|
| 2029 | for (i = k-1; i >= 0; i--) { |
---|
| 2030 | if (h1 == baby[i]) |
---|
| 2031 | return j*k-i; |
---|
| 2032 | } |
---|
| 2033 | } |
---|
| 2034 | |
---|
| 2035 | return n; |
---|
| 2036 | } |
---|
| 2037 | |
---|
| 2038 | NTL_END_IMPL |
---|