[2cfffe] | 1 | |
---|
| 2 | #include <NTL/ZZ_pX.h> |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | // The mul & sqr routines use routines from ZZX, |
---|
| 6 | // which is faster for small degree polynomials. |
---|
| 7 | // Define this macro to revert to old strategy. |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | #ifndef NTL_OLD_ZZ_pX_MUL |
---|
| 11 | |
---|
| 12 | #include <NTL/ZZX.h> |
---|
| 13 | |
---|
| 14 | #endif |
---|
| 15 | |
---|
| 16 | #include <NTL/new.h> |
---|
| 17 | |
---|
| 18 | |
---|
[de6a29] | 19 | #if (defined(NTL_GMP_LIP) || defined(NTL_GMP_HACK)) |
---|
| 20 | #define KARX 200 |
---|
| 21 | #else |
---|
| 22 | #define KARX 80 |
---|
| 23 | #endif |
---|
| 24 | |
---|
[2cfffe] | 25 | |
---|
| 26 | NTL_START_IMPL |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | const ZZ_pX& ZZ_pX::zero() |
---|
| 32 | { |
---|
| 33 | static ZZ_pX z; |
---|
| 34 | return z; |
---|
| 35 | } |
---|
| 36 | |
---|
| 37 | |
---|
| 38 | ZZ_pX& ZZ_pX::operator=(long a) |
---|
| 39 | { |
---|
| 40 | conv(*this, a); |
---|
| 41 | return *this; |
---|
| 42 | } |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | ZZ_pX& ZZ_pX::operator=(const ZZ_p& a) |
---|
| 46 | { |
---|
| 47 | conv(*this, a); |
---|
| 48 | return *this; |
---|
| 49 | } |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | void ZZ_pX::normalize() |
---|
| 53 | { |
---|
| 54 | long n; |
---|
| 55 | const ZZ_p* p; |
---|
| 56 | |
---|
| 57 | n = rep.length(); |
---|
| 58 | if (n == 0) return; |
---|
[09da99] | 59 | p = rep.elts() + n; |
---|
| 60 | while (n > 0 && IsZero(*--p)) { |
---|
[2cfffe] | 61 | n--; |
---|
| 62 | } |
---|
| 63 | rep.SetLength(n); |
---|
| 64 | } |
---|
| 65 | |
---|
| 66 | |
---|
| 67 | long IsZero(const ZZ_pX& a) |
---|
| 68 | { |
---|
| 69 | return a.rep.length() == 0; |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | |
---|
| 73 | long IsOne(const ZZ_pX& a) |
---|
| 74 | { |
---|
| 75 | return a.rep.length() == 1 && IsOne(a.rep[0]); |
---|
| 76 | } |
---|
| 77 | |
---|
| 78 | void GetCoeff(ZZ_p& x, const ZZ_pX& a, long i) |
---|
| 79 | { |
---|
| 80 | if (i < 0 || i > deg(a)) |
---|
| 81 | clear(x); |
---|
| 82 | else |
---|
| 83 | x = a.rep[i]; |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | void SetCoeff(ZZ_pX& x, long i, const ZZ_p& a) |
---|
| 87 | { |
---|
| 88 | long j, m; |
---|
| 89 | |
---|
| 90 | if (i < 0) |
---|
| 91 | Error("SetCoeff: negative index"); |
---|
| 92 | |
---|
[09da99] | 93 | if (NTL_OVERFLOW(i, 1, 0)) |
---|
[2cfffe] | 94 | Error("overflow in SetCoeff"); |
---|
| 95 | |
---|
| 96 | m = deg(x); |
---|
| 97 | |
---|
| 98 | if (i > m) { |
---|
[09da99] | 99 | /* careful: a may alias a coefficient of x */ |
---|
[2cfffe] | 100 | |
---|
[09da99] | 101 | long alloc = x.rep.allocated(); |
---|
| 102 | |
---|
| 103 | if (alloc > 0 && i >= alloc) { |
---|
| 104 | ZZ_pTemp aa_tmp; ZZ_p& aa = aa_tmp.val(); |
---|
| 105 | aa = a; |
---|
| 106 | x.rep.SetLength(i+1); |
---|
| 107 | x.rep[i] = aa; |
---|
| 108 | } |
---|
| 109 | else { |
---|
| 110 | x.rep.SetLength(i+1); |
---|
[2cfffe] | 111 | x.rep[i] = a; |
---|
[09da99] | 112 | } |
---|
[2cfffe] | 113 | |
---|
| 114 | for (j = m+1; j < i; j++) |
---|
| 115 | clear(x.rep[j]); |
---|
| 116 | } |
---|
| 117 | else |
---|
| 118 | x.rep[i] = a; |
---|
| 119 | |
---|
| 120 | x.normalize(); |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | void SetCoeff(ZZ_pX& x, long i, long a) |
---|
| 124 | { |
---|
| 125 | if (a == 1) |
---|
| 126 | SetCoeff(x, i); |
---|
| 127 | else { |
---|
| 128 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 129 | conv(T, a); |
---|
| 130 | SetCoeff(x, i, T); |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | void SetCoeff(ZZ_pX& x, long i) |
---|
| 135 | { |
---|
| 136 | long j, m; |
---|
| 137 | |
---|
| 138 | if (i < 0) |
---|
| 139 | Error("coefficient index out of range"); |
---|
| 140 | |
---|
[09da99] | 141 | if (NTL_OVERFLOW(i, 1, 0)) |
---|
[2cfffe] | 142 | Error("overflow in SetCoeff"); |
---|
| 143 | |
---|
| 144 | m = deg(x); |
---|
| 145 | |
---|
| 146 | if (i > m) { |
---|
| 147 | x.rep.SetLength(i+1); |
---|
| 148 | for (j = m+1; j < i; j++) |
---|
| 149 | clear(x.rep[j]); |
---|
| 150 | } |
---|
| 151 | set(x.rep[i]); |
---|
| 152 | x.normalize(); |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | void SetX(ZZ_pX& x) |
---|
| 157 | { |
---|
| 158 | clear(x); |
---|
| 159 | SetCoeff(x, 1); |
---|
| 160 | } |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | long IsX(const ZZ_pX& a) |
---|
| 164 | { |
---|
| 165 | return deg(a) == 1 && IsOne(LeadCoeff(a)) && IsZero(ConstTerm(a)); |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | const ZZ_p& coeff(const ZZ_pX& a, long i) |
---|
| 171 | { |
---|
| 172 | if (i < 0 || i > deg(a)) |
---|
| 173 | return ZZ_p::zero(); |
---|
| 174 | else |
---|
| 175 | return a.rep[i]; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | const ZZ_p& LeadCoeff(const ZZ_pX& a) |
---|
| 180 | { |
---|
| 181 | if (IsZero(a)) |
---|
| 182 | return ZZ_p::zero(); |
---|
| 183 | else |
---|
| 184 | return a.rep[deg(a)]; |
---|
| 185 | } |
---|
| 186 | |
---|
| 187 | const ZZ_p& ConstTerm(const ZZ_pX& a) |
---|
| 188 | { |
---|
| 189 | if (IsZero(a)) |
---|
| 190 | return ZZ_p::zero(); |
---|
| 191 | else |
---|
| 192 | return a.rep[0]; |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | |
---|
| 197 | void conv(ZZ_pX& x, const ZZ_p& a) |
---|
| 198 | { |
---|
| 199 | if (IsZero(a)) |
---|
| 200 | x.rep.SetLength(0); |
---|
| 201 | else { |
---|
| 202 | x.rep.SetLength(1); |
---|
| 203 | x.rep[0] = a; |
---|
| 204 | |
---|
| 205 | // note: if a aliases x.rep[i], i > 0, this code |
---|
| 206 | // will still work, since is is assumed that |
---|
| 207 | // SetLength(1) will not relocate or destroy x.rep[i] |
---|
| 208 | } |
---|
| 209 | } |
---|
| 210 | |
---|
| 211 | void conv(ZZ_pX& x, long a) |
---|
| 212 | { |
---|
| 213 | if (a == 0) |
---|
| 214 | clear(x); |
---|
| 215 | else if (a == 1) |
---|
| 216 | set(x); |
---|
| 217 | else { |
---|
| 218 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 219 | conv(T, a); |
---|
| 220 | conv(x, T); |
---|
| 221 | } |
---|
| 222 | } |
---|
| 223 | |
---|
| 224 | void conv(ZZ_pX& x, const ZZ& a) |
---|
| 225 | { |
---|
| 226 | if (IsZero(a)) |
---|
| 227 | clear(x); |
---|
| 228 | else { |
---|
| 229 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 230 | conv(T, a); |
---|
| 231 | conv(x, T); |
---|
| 232 | } |
---|
| 233 | } |
---|
| 234 | |
---|
| 235 | void conv(ZZ_pX& x, const vec_ZZ_p& a) |
---|
| 236 | { |
---|
| 237 | x.rep = a; |
---|
| 238 | x.normalize(); |
---|
| 239 | } |
---|
| 240 | |
---|
| 241 | |
---|
| 242 | void add(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 243 | { |
---|
| 244 | long da = deg(a); |
---|
| 245 | long db = deg(b); |
---|
| 246 | long minab = min(da, db); |
---|
| 247 | long maxab = max(da, db); |
---|
| 248 | x.rep.SetLength(maxab+1); |
---|
| 249 | |
---|
| 250 | long i; |
---|
| 251 | const ZZ_p *ap, *bp; |
---|
| 252 | ZZ_p* xp; |
---|
| 253 | |
---|
| 254 | for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts(); |
---|
| 255 | i; i--, ap++, bp++, xp++) |
---|
| 256 | add(*xp, (*ap), (*bp)); |
---|
| 257 | |
---|
| 258 | if (da > minab && &x != &a) |
---|
| 259 | for (i = da-minab; i; i--, xp++, ap++) |
---|
| 260 | *xp = *ap; |
---|
| 261 | else if (db > minab && &x != &b) |
---|
| 262 | for (i = db-minab; i; i--, xp++, bp++) |
---|
| 263 | *xp = *bp; |
---|
| 264 | else |
---|
| 265 | x.normalize(); |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | |
---|
| 269 | void add(ZZ_pX& x, const ZZ_pX& a, const ZZ_p& b) |
---|
| 270 | { |
---|
| 271 | long n = a.rep.length(); |
---|
| 272 | if (n == 0) { |
---|
| 273 | conv(x, b); |
---|
| 274 | } |
---|
| 275 | else if (&x == &a) { |
---|
| 276 | add(x.rep[0], a.rep[0], b); |
---|
| 277 | x.normalize(); |
---|
| 278 | } |
---|
| 279 | else if (x.rep.MaxLength() == 0) { |
---|
| 280 | x = a; |
---|
| 281 | add(x.rep[0], a.rep[0], b); |
---|
| 282 | x.normalize(); |
---|
| 283 | } |
---|
| 284 | else { |
---|
| 285 | // ugly...b could alias a coeff of x |
---|
| 286 | |
---|
| 287 | ZZ_p *xp = x.rep.elts(); |
---|
| 288 | add(xp[0], a.rep[0], b); |
---|
| 289 | x.rep.SetLength(n); |
---|
| 290 | xp = x.rep.elts(); |
---|
| 291 | const ZZ_p *ap = a.rep.elts(); |
---|
| 292 | long i; |
---|
| 293 | for (i = 1; i < n; i++) |
---|
| 294 | xp[i] = ap[i]; |
---|
| 295 | x.normalize(); |
---|
| 296 | } |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | void add(ZZ_pX& x, const ZZ_pX& a, long b) |
---|
| 300 | { |
---|
| 301 | if (a.rep.length() == 0) { |
---|
| 302 | conv(x, b); |
---|
| 303 | } |
---|
| 304 | else { |
---|
| 305 | if (&x != &a) x = a; |
---|
| 306 | add(x.rep[0], x.rep[0], b); |
---|
| 307 | x.normalize(); |
---|
| 308 | } |
---|
| 309 | } |
---|
| 310 | |
---|
| 311 | void sub(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 312 | { |
---|
| 313 | long da = deg(a); |
---|
| 314 | long db = deg(b); |
---|
| 315 | long minab = min(da, db); |
---|
| 316 | long maxab = max(da, db); |
---|
| 317 | x.rep.SetLength(maxab+1); |
---|
| 318 | |
---|
| 319 | long i; |
---|
| 320 | const ZZ_p *ap, *bp; |
---|
| 321 | ZZ_p* xp; |
---|
| 322 | |
---|
| 323 | for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts(); |
---|
| 324 | i; i--, ap++, bp++, xp++) |
---|
| 325 | sub(*xp, (*ap), (*bp)); |
---|
| 326 | |
---|
| 327 | if (da > minab && &x != &a) |
---|
| 328 | for (i = da-minab; i; i--, xp++, ap++) |
---|
| 329 | *xp = *ap; |
---|
| 330 | else if (db > minab) |
---|
| 331 | for (i = db-minab; i; i--, xp++, bp++) |
---|
| 332 | negate(*xp, *bp); |
---|
| 333 | else |
---|
| 334 | x.normalize(); |
---|
| 335 | |
---|
| 336 | } |
---|
| 337 | |
---|
| 338 | void sub(ZZ_pX& x, const ZZ_pX& a, const ZZ_p& b) |
---|
| 339 | { |
---|
| 340 | long n = a.rep.length(); |
---|
| 341 | if (n == 0) { |
---|
| 342 | conv(x, b); |
---|
| 343 | negate(x, x); |
---|
| 344 | } |
---|
| 345 | else if (&x == &a) { |
---|
| 346 | sub(x.rep[0], a.rep[0], b); |
---|
| 347 | x.normalize(); |
---|
| 348 | } |
---|
| 349 | else if (x.rep.MaxLength() == 0) { |
---|
| 350 | x = a; |
---|
| 351 | sub(x.rep[0], a.rep[0], b); |
---|
| 352 | x.normalize(); |
---|
| 353 | } |
---|
| 354 | else { |
---|
| 355 | // ugly...b could alias a coeff of x |
---|
| 356 | |
---|
| 357 | ZZ_p *xp = x.rep.elts(); |
---|
| 358 | sub(xp[0], a.rep[0], b); |
---|
| 359 | x.rep.SetLength(n); |
---|
| 360 | xp = x.rep.elts(); |
---|
| 361 | const ZZ_p *ap = a.rep.elts(); |
---|
| 362 | long i; |
---|
| 363 | for (i = 1; i < n; i++) |
---|
| 364 | xp[i] = ap[i]; |
---|
| 365 | x.normalize(); |
---|
| 366 | } |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | void sub(ZZ_pX& x, const ZZ_pX& a, long b) |
---|
| 370 | { |
---|
| 371 | if (b == 0) { |
---|
| 372 | x = a; |
---|
| 373 | return; |
---|
| 374 | } |
---|
| 375 | |
---|
| 376 | if (a.rep.length() == 0) { |
---|
| 377 | x.rep.SetLength(1); |
---|
| 378 | x.rep[0] = b; |
---|
| 379 | negate(x.rep[0], x.rep[0]); |
---|
| 380 | } |
---|
| 381 | else { |
---|
| 382 | if (&x != &a) x = a; |
---|
| 383 | sub(x.rep[0], x.rep[0], b); |
---|
| 384 | } |
---|
| 385 | x.normalize(); |
---|
| 386 | } |
---|
| 387 | |
---|
| 388 | void sub(ZZ_pX& x, const ZZ_p& a, const ZZ_pX& b) |
---|
| 389 | { |
---|
| 390 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 391 | T = a; |
---|
| 392 | |
---|
| 393 | negate(x, b); |
---|
| 394 | add(x, x, T); |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | void sub(ZZ_pX& x, long a, const ZZ_pX& b) |
---|
| 398 | { |
---|
| 399 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 400 | T = a; |
---|
| 401 | |
---|
| 402 | negate(x, b); |
---|
| 403 | add(x, x, T); |
---|
| 404 | } |
---|
| 405 | |
---|
| 406 | void negate(ZZ_pX& x, const ZZ_pX& a) |
---|
| 407 | { |
---|
| 408 | long n = a.rep.length(); |
---|
| 409 | x.rep.SetLength(n); |
---|
| 410 | |
---|
| 411 | const ZZ_p* ap = a.rep.elts(); |
---|
| 412 | ZZ_p* xp = x.rep.elts(); |
---|
| 413 | long i; |
---|
| 414 | |
---|
| 415 | for (i = n; i; i--, ap++, xp++) |
---|
| 416 | negate((*xp), (*ap)); |
---|
| 417 | |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | #ifndef NTL_OLD_ZZ_pX_MUL |
---|
| 422 | |
---|
| 423 | // These crossovers are tuned for a Pentium, but hopefully |
---|
| 424 | // they should be OK on other machines as well. |
---|
| 425 | |
---|
| 426 | |
---|
| 427 | const long SS_kbound = 40; |
---|
| 428 | const double SS_rbound = 1.25; |
---|
| 429 | |
---|
| 430 | |
---|
| 431 | void mul(ZZ_pX& c, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 432 | { |
---|
| 433 | if (IsZero(a) || IsZero(b)) { |
---|
| 434 | clear(c); |
---|
| 435 | return; |
---|
| 436 | } |
---|
| 437 | |
---|
| 438 | if (&a == &b) { |
---|
| 439 | sqr(c, a); |
---|
| 440 | return; |
---|
| 441 | } |
---|
| 442 | |
---|
| 443 | long k = ZZ_p::ModulusSize(); |
---|
| 444 | long s = min(deg(a), deg(b)) + 1; |
---|
| 445 | |
---|
| 446 | if (s == 1 || (k == 1 && s < 40) || (k == 2 && s < 20) || |
---|
| 447 | (k == 3 && s < 12) || (k <= 5 && s < 8) || |
---|
| 448 | (k <= 12 && s < 4) ) { |
---|
| 449 | PlainMul(c, a, b); |
---|
| 450 | } |
---|
[de6a29] | 451 | else if (s < KARX) { |
---|
[2cfffe] | 452 | ZZX A, B, C; |
---|
| 453 | conv(A, a); |
---|
| 454 | conv(B, b); |
---|
| 455 | KarMul(C, A, B); |
---|
| 456 | conv(c, C); |
---|
| 457 | } |
---|
| 458 | else { |
---|
| 459 | long mbits; |
---|
| 460 | mbits = NumBits(ZZ_p::modulus()); |
---|
| 461 | if (k >= SS_kbound && |
---|
| 462 | SSRatio(deg(a), mbits, deg(b), mbits) < SS_rbound) { |
---|
| 463 | ZZX A, B, C; |
---|
| 464 | conv(A, a); |
---|
| 465 | conv(B, b); |
---|
| 466 | SSMul(C, A, B); |
---|
| 467 | conv(c, C); |
---|
| 468 | } |
---|
| 469 | else { |
---|
| 470 | FFTMul(c, a, b); |
---|
| 471 | } |
---|
| 472 | } |
---|
| 473 | } |
---|
| 474 | |
---|
| 475 | void sqr(ZZ_pX& c, const ZZ_pX& a) |
---|
| 476 | { |
---|
| 477 | if (IsZero(a)) { |
---|
| 478 | clear(c); |
---|
| 479 | return; |
---|
| 480 | } |
---|
| 481 | |
---|
| 482 | long k = ZZ_p::ModulusSize(); |
---|
| 483 | long s = deg(a) + 1; |
---|
| 484 | |
---|
| 485 | if (s == 1 || (k == 1 && s < 50) || (k == 2 && s < 25) || |
---|
| 486 | (k == 3 && s < 25) || (k <= 6 && s < 12) || |
---|
| 487 | (k <= 8 && s < 8) || (k == 9 && s < 6) || |
---|
| 488 | (k <= 30 && s < 4) ) { |
---|
| 489 | |
---|
| 490 | PlainSqr(c, a); |
---|
| 491 | } |
---|
| 492 | else if (s < 80) { |
---|
| 493 | ZZX C, A; |
---|
| 494 | conv(A, a); |
---|
| 495 | KarSqr(C, A); |
---|
| 496 | conv(c, C); |
---|
| 497 | } |
---|
| 498 | else { |
---|
| 499 | long mbits; |
---|
| 500 | mbits = NumBits(ZZ_p::modulus()); |
---|
| 501 | if (k >= SS_kbound && |
---|
| 502 | SSRatio(deg(a), mbits, deg(a), mbits) < SS_rbound) { |
---|
| 503 | ZZX A, C; |
---|
| 504 | conv(A, a); |
---|
| 505 | SSSqr(C, A); |
---|
| 506 | conv(c, C); |
---|
| 507 | } |
---|
| 508 | else { |
---|
| 509 | FFTSqr(c, a); |
---|
| 510 | } |
---|
| 511 | } |
---|
| 512 | } |
---|
| 513 | |
---|
| 514 | #else |
---|
| 515 | |
---|
| 516 | void mul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 517 | { |
---|
| 518 | if (&a == &b) { |
---|
| 519 | sqr(x, a); |
---|
| 520 | return; |
---|
| 521 | } |
---|
| 522 | |
---|
| 523 | if (deg(a) > NTL_ZZ_pX_FFT_CROSSOVER && deg(b) > NTL_ZZ_pX_FFT_CROSSOVER) |
---|
| 524 | FFTMul(x, a, b); |
---|
| 525 | else |
---|
| 526 | PlainMul(x, a, b); |
---|
| 527 | } |
---|
| 528 | |
---|
| 529 | void sqr(ZZ_pX& x, const ZZ_pX& a) |
---|
| 530 | { |
---|
| 531 | if (deg(a) > NTL_ZZ_pX_FFT_CROSSOVER) |
---|
| 532 | FFTSqr(x, a); |
---|
| 533 | else |
---|
| 534 | PlainSqr(x, a); |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | |
---|
| 538 | #endif |
---|
| 539 | |
---|
| 540 | |
---|
| 541 | void PlainMul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 542 | { |
---|
| 543 | long da = deg(a); |
---|
| 544 | long db = deg(b); |
---|
| 545 | |
---|
| 546 | if (da < 0 || db < 0) { |
---|
| 547 | clear(x); |
---|
| 548 | return; |
---|
| 549 | } |
---|
| 550 | |
---|
| 551 | if (da == 0) { |
---|
| 552 | mul(x, b, a.rep[0]); |
---|
| 553 | return; |
---|
| 554 | } |
---|
| 555 | |
---|
| 556 | if (db == 0) { |
---|
| 557 | mul(x, a, b.rep[0]); |
---|
| 558 | return; |
---|
| 559 | } |
---|
| 560 | |
---|
| 561 | long d = da+db; |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | |
---|
| 565 | const ZZ_p *ap, *bp; |
---|
| 566 | ZZ_p *xp; |
---|
| 567 | |
---|
| 568 | ZZ_pX la, lb; |
---|
| 569 | |
---|
| 570 | if (&x == &a) { |
---|
| 571 | la = a; |
---|
| 572 | ap = la.rep.elts(); |
---|
| 573 | } |
---|
| 574 | else |
---|
| 575 | ap = a.rep.elts(); |
---|
| 576 | |
---|
| 577 | if (&x == &b) { |
---|
| 578 | lb = b; |
---|
| 579 | bp = lb.rep.elts(); |
---|
| 580 | } |
---|
| 581 | else |
---|
| 582 | bp = b.rep.elts(); |
---|
| 583 | |
---|
| 584 | x.rep.SetLength(d+1); |
---|
| 585 | |
---|
| 586 | xp = x.rep.elts(); |
---|
| 587 | |
---|
| 588 | long i, j, jmin, jmax; |
---|
| 589 | static ZZ t, accum; |
---|
| 590 | |
---|
| 591 | for (i = 0; i <= d; i++) { |
---|
| 592 | jmin = max(0, i-db); |
---|
| 593 | jmax = min(da, i); |
---|
| 594 | clear(accum); |
---|
| 595 | for (j = jmin; j <= jmax; j++) { |
---|
| 596 | mul(t, rep(ap[j]), rep(bp[i-j])); |
---|
| 597 | add(accum, accum, t); |
---|
| 598 | } |
---|
| 599 | conv(xp[i], accum); |
---|
| 600 | } |
---|
| 601 | x.normalize(); |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | void PlainSqr(ZZ_pX& x, const ZZ_pX& a) |
---|
| 605 | { |
---|
| 606 | long da = deg(a); |
---|
| 607 | |
---|
| 608 | if (da < 0) { |
---|
| 609 | clear(x); |
---|
| 610 | return; |
---|
| 611 | } |
---|
| 612 | |
---|
| 613 | long d = 2*da; |
---|
| 614 | |
---|
| 615 | const ZZ_p *ap; |
---|
| 616 | ZZ_p *xp; |
---|
| 617 | |
---|
| 618 | ZZ_pX la; |
---|
| 619 | |
---|
| 620 | if (&x == &a) { |
---|
| 621 | la = a; |
---|
| 622 | ap = la.rep.elts(); |
---|
| 623 | } |
---|
| 624 | else |
---|
| 625 | ap = a.rep.elts(); |
---|
| 626 | |
---|
| 627 | |
---|
| 628 | x.rep.SetLength(d+1); |
---|
| 629 | |
---|
| 630 | xp = x.rep.elts(); |
---|
| 631 | |
---|
| 632 | long i, j, jmin, jmax; |
---|
| 633 | long m, m2; |
---|
| 634 | static ZZ t, accum; |
---|
| 635 | |
---|
| 636 | for (i = 0; i <= d; i++) { |
---|
| 637 | jmin = max(0, i-da); |
---|
| 638 | jmax = min(da, i); |
---|
| 639 | m = jmax - jmin + 1; |
---|
| 640 | m2 = m >> 1; |
---|
| 641 | jmax = jmin + m2 - 1; |
---|
| 642 | clear(accum); |
---|
| 643 | for (j = jmin; j <= jmax; j++) { |
---|
| 644 | mul(t, rep(ap[j]), rep(ap[i-j])); |
---|
| 645 | add(accum, accum, t); |
---|
| 646 | } |
---|
| 647 | add(accum, accum, accum); |
---|
| 648 | if (m & 1) { |
---|
| 649 | sqr(t, rep(ap[jmax + 1])); |
---|
| 650 | add(accum, accum, t); |
---|
| 651 | } |
---|
| 652 | |
---|
| 653 | conv(xp[i], accum); |
---|
| 654 | } |
---|
| 655 | |
---|
| 656 | x.normalize(); |
---|
| 657 | } |
---|
| 658 | |
---|
| 659 | void PlainDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 660 | { |
---|
| 661 | long da, db, dq, i, j, LCIsOne; |
---|
| 662 | const ZZ_p *bp; |
---|
| 663 | ZZ_p *qp; |
---|
| 664 | ZZ *xp; |
---|
| 665 | |
---|
| 666 | |
---|
| 667 | ZZ_p LCInv, t; |
---|
| 668 | static ZZ s; |
---|
| 669 | |
---|
| 670 | da = deg(a); |
---|
| 671 | db = deg(b); |
---|
| 672 | |
---|
| 673 | if (db < 0) Error("ZZ_pX: division by zero"); |
---|
| 674 | |
---|
| 675 | if (da < db) { |
---|
| 676 | r = a; |
---|
| 677 | clear(q); |
---|
| 678 | return; |
---|
| 679 | } |
---|
| 680 | |
---|
| 681 | ZZ_pX lb; |
---|
| 682 | |
---|
| 683 | if (&q == &b) { |
---|
| 684 | lb = b; |
---|
| 685 | bp = lb.rep.elts(); |
---|
| 686 | } |
---|
| 687 | else |
---|
| 688 | bp = b.rep.elts(); |
---|
| 689 | |
---|
| 690 | if (IsOne(bp[db])) |
---|
| 691 | LCIsOne = 1; |
---|
| 692 | else { |
---|
| 693 | LCIsOne = 0; |
---|
| 694 | inv(LCInv, bp[db]); |
---|
| 695 | } |
---|
| 696 | |
---|
| 697 | ZZVec x(da + 1, ZZ_pInfo->ExtendedModulusSize); |
---|
| 698 | |
---|
| 699 | for (i = 0; i <= da; i++) |
---|
| 700 | x[i] = rep(a.rep[i]); |
---|
| 701 | |
---|
| 702 | xp = x.elts(); |
---|
| 703 | |
---|
| 704 | dq = da - db; |
---|
| 705 | q.rep.SetLength(dq+1); |
---|
| 706 | qp = q.rep.elts(); |
---|
| 707 | |
---|
| 708 | for (i = dq; i >= 0; i--) { |
---|
| 709 | conv(t, xp[i+db]); |
---|
| 710 | if (!LCIsOne) |
---|
| 711 | mul(t, t, LCInv); |
---|
| 712 | qp[i] = t; |
---|
| 713 | negate(t, t); |
---|
| 714 | |
---|
| 715 | for (j = db-1; j >= 0; j--) { |
---|
| 716 | mul(s, rep(t), rep(bp[j])); |
---|
| 717 | add(xp[i+j], xp[i+j], s); |
---|
| 718 | } |
---|
| 719 | } |
---|
| 720 | |
---|
| 721 | r.rep.SetLength(db); |
---|
| 722 | for (i = 0; i < db; i++) |
---|
| 723 | conv(r.rep[i], xp[i]); |
---|
| 724 | r.normalize(); |
---|
| 725 | } |
---|
| 726 | |
---|
| 727 | |
---|
| 728 | void PlainRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b, ZZVec& x) |
---|
| 729 | { |
---|
| 730 | long da, db, dq, i, j, LCIsOne; |
---|
| 731 | const ZZ_p *bp; |
---|
| 732 | ZZ *xp; |
---|
| 733 | |
---|
| 734 | |
---|
| 735 | ZZ_p LCInv, t; |
---|
| 736 | static ZZ s; |
---|
| 737 | |
---|
| 738 | da = deg(a); |
---|
| 739 | db = deg(b); |
---|
| 740 | |
---|
| 741 | if (db < 0) Error("ZZ_pX: division by zero"); |
---|
| 742 | |
---|
| 743 | if (da < db) { |
---|
| 744 | r = a; |
---|
| 745 | return; |
---|
| 746 | } |
---|
| 747 | |
---|
| 748 | bp = b.rep.elts(); |
---|
| 749 | |
---|
| 750 | if (IsOne(bp[db])) |
---|
| 751 | LCIsOne = 1; |
---|
| 752 | else { |
---|
| 753 | LCIsOne = 0; |
---|
| 754 | inv(LCInv, bp[db]); |
---|
| 755 | } |
---|
| 756 | |
---|
| 757 | for (i = 0; i <= da; i++) |
---|
| 758 | x[i] = rep(a.rep[i]); |
---|
| 759 | |
---|
| 760 | xp = x.elts(); |
---|
| 761 | |
---|
| 762 | dq = da - db; |
---|
| 763 | |
---|
| 764 | for (i = dq; i >= 0; i--) { |
---|
| 765 | conv(t, xp[i+db]); |
---|
| 766 | if (!LCIsOne) |
---|
| 767 | mul(t, t, LCInv); |
---|
| 768 | negate(t, t); |
---|
| 769 | |
---|
| 770 | for (j = db-1; j >= 0; j--) { |
---|
| 771 | mul(s, rep(t), rep(bp[j])); |
---|
| 772 | add(xp[i+j], xp[i+j], s); |
---|
| 773 | } |
---|
| 774 | } |
---|
| 775 | |
---|
| 776 | r.rep.SetLength(db); |
---|
| 777 | for (i = 0; i < db; i++) |
---|
| 778 | conv(r.rep[i], xp[i]); |
---|
| 779 | r.normalize(); |
---|
| 780 | } |
---|
| 781 | |
---|
| 782 | |
---|
| 783 | void PlainDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b, ZZVec& x) |
---|
| 784 | { |
---|
| 785 | long da, db, dq, i, j, LCIsOne; |
---|
| 786 | const ZZ_p *bp; |
---|
| 787 | ZZ_p *qp; |
---|
| 788 | ZZ *xp; |
---|
| 789 | |
---|
| 790 | |
---|
| 791 | ZZ_p LCInv, t; |
---|
| 792 | static ZZ s; |
---|
| 793 | |
---|
| 794 | da = deg(a); |
---|
| 795 | db = deg(b); |
---|
| 796 | |
---|
| 797 | if (db < 0) Error("ZZ_pX: division by zero"); |
---|
| 798 | |
---|
| 799 | if (da < db) { |
---|
| 800 | r = a; |
---|
| 801 | clear(q); |
---|
| 802 | return; |
---|
| 803 | } |
---|
| 804 | |
---|
| 805 | ZZ_pX lb; |
---|
| 806 | |
---|
| 807 | if (&q == &b) { |
---|
| 808 | lb = b; |
---|
| 809 | bp = lb.rep.elts(); |
---|
| 810 | } |
---|
| 811 | else |
---|
| 812 | bp = b.rep.elts(); |
---|
| 813 | |
---|
| 814 | if (IsOne(bp[db])) |
---|
| 815 | LCIsOne = 1; |
---|
| 816 | else { |
---|
| 817 | LCIsOne = 0; |
---|
| 818 | inv(LCInv, bp[db]); |
---|
| 819 | } |
---|
| 820 | |
---|
| 821 | for (i = 0; i <= da; i++) |
---|
| 822 | x[i] = rep(a.rep[i]); |
---|
| 823 | |
---|
| 824 | xp = x.elts(); |
---|
| 825 | |
---|
| 826 | dq = da - db; |
---|
| 827 | q.rep.SetLength(dq+1); |
---|
| 828 | qp = q.rep.elts(); |
---|
| 829 | |
---|
| 830 | for (i = dq; i >= 0; i--) { |
---|
| 831 | conv(t, xp[i+db]); |
---|
| 832 | if (!LCIsOne) |
---|
| 833 | mul(t, t, LCInv); |
---|
| 834 | qp[i] = t; |
---|
| 835 | negate(t, t); |
---|
| 836 | |
---|
| 837 | for (j = db-1; j >= 0; j--) { |
---|
| 838 | mul(s, rep(t), rep(bp[j])); |
---|
| 839 | add(xp[i+j], xp[i+j], s); |
---|
| 840 | } |
---|
| 841 | } |
---|
| 842 | |
---|
| 843 | r.rep.SetLength(db); |
---|
| 844 | for (i = 0; i < db; i++) |
---|
| 845 | conv(r.rep[i], xp[i]); |
---|
| 846 | r.normalize(); |
---|
| 847 | } |
---|
| 848 | |
---|
| 849 | |
---|
| 850 | void PlainDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 851 | { |
---|
| 852 | long da, db, dq, i, j, LCIsOne; |
---|
| 853 | const ZZ_p *bp; |
---|
| 854 | ZZ_p *qp; |
---|
| 855 | ZZ *xp; |
---|
| 856 | |
---|
| 857 | |
---|
| 858 | ZZ_p LCInv, t; |
---|
| 859 | static ZZ s; |
---|
| 860 | |
---|
| 861 | da = deg(a); |
---|
| 862 | db = deg(b); |
---|
| 863 | |
---|
| 864 | if (db < 0) Error("ZZ_pX: division by zero"); |
---|
| 865 | |
---|
| 866 | if (da < db) { |
---|
| 867 | clear(q); |
---|
| 868 | return; |
---|
| 869 | } |
---|
| 870 | |
---|
| 871 | ZZ_pX lb; |
---|
| 872 | |
---|
| 873 | if (&q == &b) { |
---|
| 874 | lb = b; |
---|
| 875 | bp = lb.rep.elts(); |
---|
| 876 | } |
---|
| 877 | else |
---|
| 878 | bp = b.rep.elts(); |
---|
| 879 | |
---|
| 880 | if (IsOne(bp[db])) |
---|
| 881 | LCIsOne = 1; |
---|
| 882 | else { |
---|
| 883 | LCIsOne = 0; |
---|
| 884 | inv(LCInv, bp[db]); |
---|
| 885 | } |
---|
| 886 | |
---|
| 887 | ZZVec x(da + 1 - db, ZZ_pInfo->ExtendedModulusSize); |
---|
| 888 | |
---|
| 889 | for (i = db; i <= da; i++) |
---|
| 890 | x[i-db] = rep(a.rep[i]); |
---|
| 891 | |
---|
| 892 | xp = x.elts(); |
---|
| 893 | |
---|
| 894 | dq = da - db; |
---|
| 895 | q.rep.SetLength(dq+1); |
---|
| 896 | qp = q.rep.elts(); |
---|
| 897 | |
---|
| 898 | for (i = dq; i >= 0; i--) { |
---|
| 899 | conv(t, xp[i]); |
---|
| 900 | if (!LCIsOne) |
---|
| 901 | mul(t, t, LCInv); |
---|
| 902 | qp[i] = t; |
---|
| 903 | negate(t, t); |
---|
| 904 | |
---|
| 905 | long lastj = max(0, db-i); |
---|
| 906 | |
---|
| 907 | for (j = db-1; j >= lastj; j--) { |
---|
| 908 | mul(s, rep(t), rep(bp[j])); |
---|
| 909 | add(xp[i+j-db], xp[i+j-db], s); |
---|
| 910 | } |
---|
| 911 | } |
---|
| 912 | } |
---|
| 913 | |
---|
| 914 | void PlainRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 915 | { |
---|
| 916 | long da, db, dq, i, j, LCIsOne; |
---|
| 917 | const ZZ_p *bp; |
---|
| 918 | ZZ *xp; |
---|
| 919 | |
---|
| 920 | |
---|
| 921 | ZZ_p LCInv, t; |
---|
| 922 | static ZZ s; |
---|
| 923 | |
---|
| 924 | da = deg(a); |
---|
| 925 | db = deg(b); |
---|
| 926 | |
---|
| 927 | if (db < 0) Error("ZZ_pX: division by zero"); |
---|
| 928 | |
---|
| 929 | if (da < db) { |
---|
| 930 | r = a; |
---|
| 931 | return; |
---|
| 932 | } |
---|
| 933 | |
---|
| 934 | bp = b.rep.elts(); |
---|
| 935 | |
---|
| 936 | if (IsOne(bp[db])) |
---|
| 937 | LCIsOne = 1; |
---|
| 938 | else { |
---|
| 939 | LCIsOne = 0; |
---|
| 940 | inv(LCInv, bp[db]); |
---|
| 941 | } |
---|
| 942 | |
---|
| 943 | ZZVec x(da + 1, ZZ_pInfo->ExtendedModulusSize); |
---|
| 944 | |
---|
| 945 | for (i = 0; i <= da; i++) |
---|
| 946 | x[i] = rep(a.rep[i]); |
---|
| 947 | |
---|
| 948 | xp = x.elts(); |
---|
| 949 | |
---|
| 950 | dq = da - db; |
---|
| 951 | |
---|
| 952 | for (i = dq; i >= 0; i--) { |
---|
| 953 | conv(t, xp[i+db]); |
---|
| 954 | if (!LCIsOne) |
---|
| 955 | mul(t, t, LCInv); |
---|
| 956 | negate(t, t); |
---|
| 957 | |
---|
| 958 | for (j = db-1; j >= 0; j--) { |
---|
| 959 | mul(s, rep(t), rep(bp[j])); |
---|
| 960 | add(xp[i+j], xp[i+j], s); |
---|
| 961 | } |
---|
| 962 | } |
---|
| 963 | |
---|
| 964 | r.rep.SetLength(db); |
---|
| 965 | for (i = 0; i < db; i++) |
---|
| 966 | conv(r.rep[i], xp[i]); |
---|
| 967 | r.normalize(); |
---|
| 968 | } |
---|
| 969 | |
---|
| 970 | void mul(ZZ_pX& x, const ZZ_pX& a, const ZZ_p& b) |
---|
| 971 | { |
---|
| 972 | if (IsZero(b)) { |
---|
| 973 | clear(x); |
---|
| 974 | return; |
---|
| 975 | } |
---|
| 976 | |
---|
| 977 | if (IsOne(b)) { |
---|
| 978 | x = a; |
---|
| 979 | return; |
---|
| 980 | } |
---|
| 981 | |
---|
| 982 | ZZ_pTemp TT; ZZ_p& t = TT.val(); |
---|
| 983 | |
---|
| 984 | long i, da; |
---|
| 985 | |
---|
| 986 | const ZZ_p *ap; |
---|
| 987 | ZZ_p* xp; |
---|
| 988 | |
---|
| 989 | |
---|
| 990 | t = b; |
---|
| 991 | |
---|
| 992 | da = deg(a); |
---|
| 993 | x.rep.SetLength(da+1); |
---|
| 994 | ap = a.rep.elts(); |
---|
| 995 | xp = x.rep.elts(); |
---|
| 996 | |
---|
| 997 | for (i = 0; i <= da; i++) |
---|
| 998 | mul(xp[i], ap[i], t); |
---|
| 999 | |
---|
| 1000 | x.normalize(); |
---|
| 1001 | } |
---|
| 1002 | |
---|
| 1003 | void mul(ZZ_pX& x, const ZZ_pX& a, long b) |
---|
| 1004 | { |
---|
| 1005 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 1006 | conv(T, b); |
---|
| 1007 | mul(x, a, T); |
---|
| 1008 | } |
---|
| 1009 | |
---|
| 1010 | |
---|
| 1011 | void PlainGCD(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 1012 | { |
---|
| 1013 | ZZ_p t; |
---|
| 1014 | |
---|
| 1015 | if (IsZero(b)) |
---|
| 1016 | x = a; |
---|
| 1017 | else if (IsZero(a)) |
---|
| 1018 | x = b; |
---|
| 1019 | else { |
---|
| 1020 | long n = max(deg(a),deg(b)) + 1; |
---|
| 1021 | ZZ_pX u(INIT_SIZE, n), v(INIT_SIZE, n); |
---|
| 1022 | ZZVec tmp(n, ZZ_pInfo->ExtendedModulusSize); |
---|
| 1023 | |
---|
| 1024 | u = a; |
---|
| 1025 | v = b; |
---|
| 1026 | do { |
---|
| 1027 | PlainRem(u, u, v, tmp); |
---|
| 1028 | swap(u, v); |
---|
| 1029 | } while (!IsZero(v)); |
---|
| 1030 | |
---|
| 1031 | x = u; |
---|
| 1032 | } |
---|
| 1033 | |
---|
| 1034 | if (IsZero(x)) return; |
---|
| 1035 | if (IsOne(LeadCoeff(x))) return; |
---|
| 1036 | |
---|
| 1037 | /* make gcd monic */ |
---|
| 1038 | |
---|
| 1039 | |
---|
| 1040 | inv(t, LeadCoeff(x)); |
---|
| 1041 | mul(x, x, t); |
---|
| 1042 | } |
---|
| 1043 | |
---|
| 1044 | |
---|
| 1045 | |
---|
| 1046 | |
---|
| 1047 | |
---|
| 1048 | void PlainXGCD(ZZ_pX& d, ZZ_pX& s, ZZ_pX& t, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 1049 | { |
---|
| 1050 | ZZ_p z; |
---|
| 1051 | |
---|
| 1052 | |
---|
| 1053 | if (IsZero(b)) { |
---|
| 1054 | set(s); |
---|
| 1055 | clear(t); |
---|
| 1056 | d = a; |
---|
| 1057 | } |
---|
| 1058 | else if (IsZero(a)) { |
---|
| 1059 | clear(s); |
---|
| 1060 | set(t); |
---|
| 1061 | d = b; |
---|
| 1062 | } |
---|
| 1063 | else { |
---|
| 1064 | long e = max(deg(a), deg(b)) + 1; |
---|
| 1065 | |
---|
| 1066 | ZZ_pX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e), |
---|
| 1067 | u0(INIT_SIZE, e), v0(INIT_SIZE, e), |
---|
| 1068 | u1(INIT_SIZE, e), v1(INIT_SIZE, e), |
---|
| 1069 | u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e); |
---|
| 1070 | |
---|
| 1071 | |
---|
| 1072 | set(u1); clear(v1); |
---|
| 1073 | clear(u2); set(v2); |
---|
| 1074 | u = a; v = b; |
---|
| 1075 | |
---|
| 1076 | do { |
---|
| 1077 | DivRem(q, u, u, v); |
---|
| 1078 | swap(u, v); |
---|
| 1079 | u0 = u2; |
---|
| 1080 | v0 = v2; |
---|
| 1081 | mul(temp, q, u2); |
---|
| 1082 | sub(u2, u1, temp); |
---|
| 1083 | mul(temp, q, v2); |
---|
| 1084 | sub(v2, v1, temp); |
---|
| 1085 | u1 = u0; |
---|
| 1086 | v1 = v0; |
---|
| 1087 | } while (!IsZero(v)); |
---|
| 1088 | |
---|
| 1089 | d = u; |
---|
| 1090 | s = u1; |
---|
| 1091 | t = v1; |
---|
| 1092 | } |
---|
| 1093 | |
---|
| 1094 | if (IsZero(d)) return; |
---|
| 1095 | if (IsOne(LeadCoeff(d))) return; |
---|
| 1096 | |
---|
| 1097 | /* make gcd monic */ |
---|
| 1098 | |
---|
| 1099 | inv(z, LeadCoeff(d)); |
---|
| 1100 | mul(d, d, z); |
---|
| 1101 | mul(s, s, z); |
---|
| 1102 | mul(t, t, z); |
---|
| 1103 | } |
---|
| 1104 | |
---|
| 1105 | |
---|
| 1106 | void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, const ZZ_pX& f) |
---|
| 1107 | { |
---|
| 1108 | if (deg(a) >= deg(f) || deg(b) >= deg(f) || deg(f) == 0) |
---|
| 1109 | Error("MulMod: bad args"); |
---|
| 1110 | |
---|
| 1111 | ZZ_pX t; |
---|
| 1112 | |
---|
| 1113 | mul(t, a, b); |
---|
| 1114 | rem(x, t, f); |
---|
| 1115 | } |
---|
| 1116 | |
---|
| 1117 | void SqrMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f) |
---|
| 1118 | { |
---|
| 1119 | if (deg(a) >= deg(f) || deg(f) == 0) Error("SqrMod: bad args"); |
---|
| 1120 | |
---|
| 1121 | ZZ_pX t; |
---|
| 1122 | |
---|
| 1123 | sqr(t, a); |
---|
| 1124 | rem(x, t, f); |
---|
| 1125 | } |
---|
| 1126 | |
---|
| 1127 | |
---|
| 1128 | void InvMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f) |
---|
| 1129 | { |
---|
| 1130 | if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args"); |
---|
| 1131 | |
---|
| 1132 | ZZ_pX d, t; |
---|
| 1133 | |
---|
| 1134 | XGCD(d, x, t, a, f); |
---|
| 1135 | if (!IsOne(d)) |
---|
| 1136 | Error("ZZ_pX InvMod: can't compute multiplicative inverse"); |
---|
| 1137 | } |
---|
| 1138 | |
---|
| 1139 | long InvModStatus(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f) |
---|
| 1140 | { |
---|
| 1141 | if (deg(a) >= deg(f) || deg(f) == 0) Error("InvModStatus: bad args"); |
---|
| 1142 | ZZ_pX d, t; |
---|
| 1143 | |
---|
| 1144 | XGCD(d, x, t, a, f); |
---|
| 1145 | if (!IsOne(d)) { |
---|
| 1146 | x = d; |
---|
| 1147 | return 1; |
---|
| 1148 | } |
---|
| 1149 | else |
---|
| 1150 | return 0; |
---|
| 1151 | } |
---|
| 1152 | |
---|
| 1153 | |
---|
| 1154 | |
---|
| 1155 | |
---|
| 1156 | static |
---|
| 1157 | void MulByXModAux(ZZ_pX& h, const ZZ_pX& a, const ZZ_pX& f) |
---|
| 1158 | { |
---|
| 1159 | long i, n, m; |
---|
| 1160 | ZZ_p* hh; |
---|
| 1161 | const ZZ_p *aa, *ff; |
---|
| 1162 | |
---|
| 1163 | ZZ_p t, z; |
---|
| 1164 | |
---|
| 1165 | n = deg(f); |
---|
| 1166 | m = deg(a); |
---|
| 1167 | |
---|
| 1168 | if (m >= n || n == 0) Error("MulByXMod: bad args"); |
---|
| 1169 | |
---|
| 1170 | if (m < 0) { |
---|
| 1171 | clear(h); |
---|
| 1172 | return; |
---|
| 1173 | } |
---|
| 1174 | |
---|
| 1175 | if (m < n-1) { |
---|
| 1176 | h.rep.SetLength(m+2); |
---|
| 1177 | hh = h.rep.elts(); |
---|
| 1178 | aa = a.rep.elts(); |
---|
| 1179 | for (i = m+1; i >= 1; i--) |
---|
| 1180 | hh[i] = aa[i-1]; |
---|
| 1181 | clear(hh[0]); |
---|
| 1182 | } |
---|
| 1183 | else { |
---|
| 1184 | h.rep.SetLength(n); |
---|
| 1185 | hh = h.rep.elts(); |
---|
| 1186 | aa = a.rep.elts(); |
---|
| 1187 | ff = f.rep.elts(); |
---|
| 1188 | negate(z, aa[n-1]); |
---|
| 1189 | if (!IsOne(ff[n])) |
---|
| 1190 | div(z, z, ff[n]); |
---|
| 1191 | for (i = n-1; i >= 1; i--) { |
---|
| 1192 | mul(t, z, ff[i]); |
---|
| 1193 | add(hh[i], aa[i-1], t); |
---|
| 1194 | } |
---|
| 1195 | mul(hh[0], z, ff[0]); |
---|
| 1196 | h.normalize(); |
---|
| 1197 | } |
---|
| 1198 | } |
---|
| 1199 | |
---|
| 1200 | |
---|
| 1201 | void MulByXMod(ZZ_pX& h, const ZZ_pX& a, const ZZ_pX& f) |
---|
| 1202 | { |
---|
| 1203 | if (&h == &f) { |
---|
| 1204 | ZZ_pX hh; |
---|
| 1205 | MulByXModAux(hh, a, f); |
---|
| 1206 | h = hh; |
---|
| 1207 | } |
---|
| 1208 | else |
---|
| 1209 | MulByXModAux(h, a, f); |
---|
| 1210 | } |
---|
| 1211 | |
---|
| 1212 | |
---|
| 1213 | |
---|
| 1214 | void random(ZZ_pX& x, long n) |
---|
| 1215 | { |
---|
| 1216 | long i; |
---|
| 1217 | |
---|
| 1218 | x.rep.SetLength(n); |
---|
| 1219 | |
---|
| 1220 | for (i = 0; i < n; i++) |
---|
| 1221 | random(x.rep[i]); |
---|
| 1222 | |
---|
| 1223 | x.normalize(); |
---|
| 1224 | } |
---|
| 1225 | |
---|
| 1226 | |
---|
| 1227 | void FFTRep::SetSize(long NewK) |
---|
| 1228 | { |
---|
| 1229 | |
---|
| 1230 | if (NewK < -1 || NewK >= NTL_BITS_PER_LONG-1) |
---|
| 1231 | Error("bad arg to FFTRep::SetSize()"); |
---|
| 1232 | |
---|
| 1233 | if (NewK <= MaxK) { |
---|
| 1234 | k = NewK; |
---|
| 1235 | return; |
---|
| 1236 | } |
---|
| 1237 | |
---|
| 1238 | ZZ_pInfo->check(); |
---|
| 1239 | |
---|
| 1240 | if (MaxK == -1) |
---|
| 1241 | NumPrimes = ZZ_pInfo->NumPrimes; |
---|
| 1242 | else { |
---|
| 1243 | if (NumPrimes != ZZ_pInfo->NumPrimes) |
---|
| 1244 | Error("FFTRep: inconsistent use"); |
---|
| 1245 | } |
---|
| 1246 | |
---|
| 1247 | long i, n; |
---|
| 1248 | |
---|
| 1249 | if (MaxK == -1) { |
---|
[09da99] | 1250 | tbl = (long **) NTL_MALLOC(NumPrimes, sizeof(long *), 0); |
---|
[2cfffe] | 1251 | if (!tbl) |
---|
| 1252 | Error("out of space in FFTRep::SetSize()"); |
---|
| 1253 | } |
---|
| 1254 | else { |
---|
| 1255 | for (i = 0; i < NumPrimes; i++) |
---|
| 1256 | free(tbl[i]); |
---|
| 1257 | } |
---|
| 1258 | |
---|
| 1259 | n = 1L << NewK; |
---|
| 1260 | |
---|
| 1261 | for (i = 0; i < NumPrimes; i++) { |
---|
[09da99] | 1262 | if ( !(tbl[i] = (long *) NTL_MALLOC(n, sizeof(long), 0)) ) |
---|
[2cfffe] | 1263 | Error("out of space in FFTRep::SetSize()"); |
---|
| 1264 | } |
---|
| 1265 | |
---|
| 1266 | k = MaxK = NewK; |
---|
| 1267 | } |
---|
| 1268 | |
---|
| 1269 | FFTRep::FFTRep(const FFTRep& R) |
---|
| 1270 | { |
---|
| 1271 | k = MaxK = R.k; |
---|
| 1272 | tbl = 0; |
---|
| 1273 | NumPrimes = 0; |
---|
| 1274 | |
---|
| 1275 | if (k < 0) return; |
---|
| 1276 | |
---|
| 1277 | NumPrimes = R.NumPrimes; |
---|
| 1278 | |
---|
| 1279 | long i, j, n; |
---|
| 1280 | |
---|
[09da99] | 1281 | tbl = (long **) NTL_MALLOC(NumPrimes, sizeof(long *), 0); |
---|
[2cfffe] | 1282 | if (!tbl) |
---|
| 1283 | Error("out of space in FFTRep"); |
---|
| 1284 | |
---|
| 1285 | n = 1L << k; |
---|
| 1286 | |
---|
| 1287 | for (i = 0; i < NumPrimes; i++) { |
---|
[09da99] | 1288 | if ( !(tbl[i] = (long *) NTL_MALLOC(n, sizeof(long), 0)) ) |
---|
[2cfffe] | 1289 | Error("out of space in FFTRep"); |
---|
| 1290 | |
---|
| 1291 | for (j = 0; j < n; j++) |
---|
| 1292 | tbl[i][j] = R.tbl[i][j]; |
---|
| 1293 | } |
---|
| 1294 | } |
---|
| 1295 | |
---|
| 1296 | FFTRep& FFTRep::operator=(const FFTRep& R) |
---|
| 1297 | { |
---|
| 1298 | if (this == &R) return *this; |
---|
| 1299 | |
---|
| 1300 | if (MaxK >= 0 && R.MaxK >= 0 && NumPrimes != R.NumPrimes) |
---|
| 1301 | Error("FFTRep: inconsistent use"); |
---|
| 1302 | |
---|
| 1303 | if (R.k < 0) { |
---|
| 1304 | k = -1; |
---|
| 1305 | return *this; |
---|
| 1306 | } |
---|
| 1307 | |
---|
| 1308 | NumPrimes = R.NumPrimes; |
---|
| 1309 | |
---|
| 1310 | if (R.k > MaxK) { |
---|
| 1311 | long i, n; |
---|
| 1312 | |
---|
| 1313 | if (MaxK == -1) { |
---|
[09da99] | 1314 | tbl = (long **) NTL_MALLOC(NumPrimes, sizeof(long *), 0); |
---|
[2cfffe] | 1315 | if (!tbl) |
---|
| 1316 | Error("out of space in FFTRep"); |
---|
| 1317 | } |
---|
| 1318 | else { |
---|
| 1319 | for (i = 0; i < NumPrimes; i++) |
---|
| 1320 | free(tbl[i]); |
---|
| 1321 | } |
---|
| 1322 | |
---|
| 1323 | n = 1L << R.k; |
---|
| 1324 | |
---|
| 1325 | for (i = 0; i < NumPrimes; i++) { |
---|
[09da99] | 1326 | if ( !(tbl[i] = (long *) NTL_MALLOC(n, sizeof(long), 0)) ) |
---|
[2cfffe] | 1327 | Error("out of space in FFTRep"); |
---|
| 1328 | } |
---|
| 1329 | |
---|
| 1330 | k = MaxK = R.k; |
---|
| 1331 | } |
---|
| 1332 | else { |
---|
| 1333 | k = R.k; |
---|
| 1334 | } |
---|
| 1335 | |
---|
| 1336 | long i, j, n; |
---|
| 1337 | |
---|
| 1338 | n = 1L << k; |
---|
| 1339 | |
---|
| 1340 | for (i = 0; i < NumPrimes; i++) |
---|
| 1341 | for (j = 0; j < n; j++) |
---|
| 1342 | tbl[i][j] = R.tbl[i][j]; |
---|
| 1343 | |
---|
| 1344 | return *this; |
---|
| 1345 | } |
---|
| 1346 | |
---|
| 1347 | FFTRep::~FFTRep() |
---|
| 1348 | { |
---|
| 1349 | if (MaxK == -1) |
---|
| 1350 | return; |
---|
| 1351 | |
---|
| 1352 | for (long i = 0; i < NumPrimes; i++) |
---|
| 1353 | free(tbl[i]); |
---|
| 1354 | |
---|
| 1355 | free(tbl); |
---|
| 1356 | } |
---|
| 1357 | |
---|
| 1358 | |
---|
| 1359 | |
---|
| 1360 | void ZZ_pXModRep::SetSize(long NewN) |
---|
| 1361 | { |
---|
| 1362 | ZZ_pInfo->check(); |
---|
| 1363 | |
---|
| 1364 | NumPrimes = ZZ_pInfo->NumPrimes; |
---|
| 1365 | |
---|
| 1366 | if (NewN < 0) |
---|
| 1367 | Error("bad arg to ZZ_pXModRep::SetSize()"); |
---|
| 1368 | |
---|
| 1369 | if (NewN <= MaxN) { |
---|
| 1370 | n = NewN; |
---|
| 1371 | return; |
---|
| 1372 | } |
---|
| 1373 | |
---|
| 1374 | long i; |
---|
| 1375 | |
---|
| 1376 | |
---|
| 1377 | if (MaxN == 0) { |
---|
[09da99] | 1378 | tbl = (long **) NTL_MALLOC(ZZ_pInfo->NumPrimes, sizeof(long *), 0); |
---|
[2cfffe] | 1379 | if (!tbl) |
---|
| 1380 | Error("out of space in ZZ_pXModRep::SetSize()"); |
---|
| 1381 | } |
---|
| 1382 | else { |
---|
| 1383 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1384 | free(tbl[i]); |
---|
| 1385 | } |
---|
| 1386 | |
---|
| 1387 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
[09da99] | 1388 | if ( !(tbl[i] = (long *) NTL_MALLOC(NewN, sizeof(long), 0)) ) |
---|
[2cfffe] | 1389 | Error("out of space in ZZ_pXModRep::SetSize()"); |
---|
| 1390 | } |
---|
| 1391 | |
---|
| 1392 | n = MaxN = NewN; |
---|
| 1393 | } |
---|
| 1394 | |
---|
| 1395 | ZZ_pXModRep::~ZZ_pXModRep() |
---|
| 1396 | { |
---|
| 1397 | if (MaxN == 0) |
---|
| 1398 | return; |
---|
| 1399 | |
---|
| 1400 | long i; |
---|
| 1401 | for (i = 0; i < NumPrimes; i++) |
---|
| 1402 | free(tbl[i]); |
---|
| 1403 | |
---|
| 1404 | free(tbl); |
---|
| 1405 | } |
---|
| 1406 | |
---|
| 1407 | |
---|
| 1408 | static vec_long ModularRepBuf; |
---|
| 1409 | static vec_long FFTBuf; |
---|
| 1410 | |
---|
| 1411 | |
---|
| 1412 | |
---|
| 1413 | |
---|
| 1414 | void ToModularRep(vec_long& x, const ZZ_p& a) |
---|
| 1415 | { |
---|
| 1416 | ZZ_pInfo->check(); |
---|
| 1417 | ZZ_p_rem_struct_eval(ZZ_pInfo->rem_struct, &x[0], rep(a)); |
---|
| 1418 | } |
---|
| 1419 | |
---|
| 1420 | |
---|
| 1421 | // NOTE: earlier versions used Kahan summation... |
---|
| 1422 | // we no longer do this, as it is less portable than I thought. |
---|
| 1423 | |
---|
| 1424 | |
---|
| 1425 | |
---|
| 1426 | void FromModularRep(ZZ_p& x, const vec_long& a) |
---|
| 1427 | { |
---|
| 1428 | ZZ_pInfo->check(); |
---|
| 1429 | |
---|
| 1430 | long n = ZZ_pInfo->NumPrimes; |
---|
| 1431 | static ZZ q, s, t; |
---|
| 1432 | long i; |
---|
| 1433 | double y; |
---|
| 1434 | |
---|
| 1435 | if (ZZ_p_crt_struct_special(ZZ_pInfo->crt_struct)) { |
---|
| 1436 | ZZ_p_crt_struct_eval(ZZ_pInfo->crt_struct, t, &a[0]); |
---|
| 1437 | x.LoopHole() = t; |
---|
| 1438 | return; |
---|
| 1439 | } |
---|
| 1440 | |
---|
| 1441 | |
---|
| 1442 | if (ZZ_pInfo->QuickCRT) { |
---|
| 1443 | y = 0; |
---|
| 1444 | for (i = 0; i < n; i++) |
---|
| 1445 | y += ((double) a[i])*ZZ_pInfo->x[i]; |
---|
| 1446 | |
---|
| 1447 | conv(q, (y + 0.5)); |
---|
| 1448 | } else { |
---|
| 1449 | long Q, r; |
---|
| 1450 | static ZZ qq; |
---|
| 1451 | |
---|
| 1452 | y = 0; |
---|
| 1453 | |
---|
| 1454 | clear(q); |
---|
| 1455 | |
---|
| 1456 | for (i = 0; i < n; i++) { |
---|
| 1457 | r = MulDivRem(Q, a[i], ZZ_pInfo->u[i], FFTPrime[i], ZZ_pInfo->x[i]); |
---|
| 1458 | add(q, q, Q); |
---|
| 1459 | y += r*FFTPrimeInv[i]; |
---|
| 1460 | } |
---|
| 1461 | |
---|
| 1462 | conv(qq, (y + 0.5)); |
---|
| 1463 | add(q, q, qq); |
---|
| 1464 | } |
---|
| 1465 | |
---|
| 1466 | ZZ_p_crt_struct_eval(ZZ_pInfo->crt_struct, t, &a[0]); |
---|
| 1467 | |
---|
| 1468 | mul(s, q, ZZ_pInfo->MinusMModP); |
---|
| 1469 | add(t, t, s); |
---|
| 1470 | |
---|
| 1471 | conv(x, t); |
---|
| 1472 | } |
---|
| 1473 | |
---|
| 1474 | |
---|
| 1475 | |
---|
| 1476 | |
---|
| 1477 | void ToFFTRep(FFTRep& y, const ZZ_pX& x, long k, long lo, long hi) |
---|
| 1478 | // computes an n = 2^k point convolution. |
---|
| 1479 | // if deg(x) >= 2^k, then x is first reduced modulo X^n-1. |
---|
| 1480 | { |
---|
| 1481 | ZZ_pInfo->check(); |
---|
| 1482 | |
---|
| 1483 | long n, i, j, m, j1; |
---|
| 1484 | vec_long& t = ModularRepBuf; |
---|
[09da99] | 1485 | vec_long& s = FFTBuf; |
---|
[2cfffe] | 1486 | ZZ_p accum; |
---|
| 1487 | |
---|
| 1488 | |
---|
| 1489 | if (k > ZZ_pInfo->MaxRoot) |
---|
| 1490 | Error("Polynomial too big for FFT"); |
---|
| 1491 | |
---|
| 1492 | if (lo < 0) |
---|
| 1493 | Error("bad arg to ToFFTRep"); |
---|
| 1494 | |
---|
| 1495 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1496 | |
---|
| 1497 | hi = min(hi, deg(x)); |
---|
| 1498 | |
---|
| 1499 | y.SetSize(k); |
---|
| 1500 | |
---|
| 1501 | n = 1L << k; |
---|
| 1502 | |
---|
| 1503 | m = max(hi-lo + 1, 0); |
---|
| 1504 | |
---|
| 1505 | const ZZ_p *xx = x.rep.elts(); |
---|
| 1506 | |
---|
| 1507 | for (j = 0; j < n; j++) { |
---|
| 1508 | if (j >= m) { |
---|
| 1509 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1510 | y.tbl[i][j] = 0; |
---|
| 1511 | } |
---|
| 1512 | else { |
---|
| 1513 | accum = xx[j+lo]; |
---|
| 1514 | for (j1 = j + n; j1 < m; j1 += n) |
---|
| 1515 | add(accum, accum, xx[j1+lo]); |
---|
| 1516 | ToModularRep(t, accum); |
---|
| 1517 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1518 | y.tbl[i][j] = t[i]; |
---|
| 1519 | } |
---|
| 1520 | } |
---|
| 1521 | } |
---|
| 1522 | |
---|
| 1523 | |
---|
| 1524 | s.SetLength(n); |
---|
| 1525 | long *sp = s.elts(); |
---|
| 1526 | |
---|
| 1527 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1528 | long *Root = &RootTable[i][0]; |
---|
| 1529 | long *yp = &y.tbl[i][0]; |
---|
| 1530 | FFT(sp, yp, y.k, FFTPrime[i], Root); |
---|
| 1531 | for (j = 0; j < n; j++) |
---|
| 1532 | yp[j] = sp[j]; |
---|
| 1533 | } |
---|
| 1534 | } |
---|
| 1535 | |
---|
| 1536 | |
---|
| 1537 | |
---|
| 1538 | void RevToFFTRep(FFTRep& y, const vec_ZZ_p& x, |
---|
| 1539 | long k, long lo, long hi, long offset) |
---|
| 1540 | // computes an n = 2^k point convolution of X^offset*x[lo..hi] mod X^n-1 |
---|
| 1541 | // using "inverted" evaluation points. |
---|
| 1542 | |
---|
| 1543 | { |
---|
| 1544 | ZZ_pInfo->check(); |
---|
| 1545 | |
---|
| 1546 | long n, i, j, m, j1; |
---|
| 1547 | vec_long& t = ModularRepBuf; |
---|
| 1548 | vec_long& s = FFTBuf; |
---|
| 1549 | ZZ_p accum; |
---|
| 1550 | |
---|
| 1551 | if (k > ZZ_pInfo->MaxRoot) |
---|
| 1552 | Error("Polynomial too big for FFT"); |
---|
| 1553 | |
---|
| 1554 | if (lo < 0) |
---|
| 1555 | Error("bad arg to ToFFTRep"); |
---|
| 1556 | |
---|
| 1557 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1558 | |
---|
| 1559 | hi = min(hi, x.length()-1); |
---|
| 1560 | |
---|
| 1561 | y.SetSize(k); |
---|
| 1562 | |
---|
| 1563 | n = 1L << k; |
---|
| 1564 | |
---|
| 1565 | m = max(hi-lo + 1, 0); |
---|
| 1566 | |
---|
| 1567 | const ZZ_p *xx = x.elts(); |
---|
| 1568 | |
---|
| 1569 | offset = offset & (n-1); |
---|
| 1570 | |
---|
| 1571 | for (j = 0; j < n; j++) { |
---|
| 1572 | if (j >= m) { |
---|
| 1573 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1574 | y.tbl[i][offset] = 0; |
---|
| 1575 | } |
---|
| 1576 | else { |
---|
| 1577 | accum = xx[j+lo]; |
---|
| 1578 | for (j1 = j + n; j1 < m; j1 += n) |
---|
| 1579 | add(accum, accum, xx[j1+lo]); |
---|
| 1580 | ToModularRep(t, accum); |
---|
| 1581 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1582 | y.tbl[i][offset] = t[i]; |
---|
| 1583 | |
---|
| 1584 | } |
---|
| 1585 | } |
---|
| 1586 | |
---|
| 1587 | offset = (offset + 1) & (n-1); |
---|
| 1588 | } |
---|
| 1589 | |
---|
| 1590 | |
---|
| 1591 | s.SetLength(n); |
---|
| 1592 | long *sp = s.elts(); |
---|
| 1593 | |
---|
| 1594 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1595 | long *Root = &RootInvTable[i][0]; |
---|
| 1596 | long *yp = &y.tbl[i][0]; |
---|
| 1597 | long w = TwoInvTable[i][k]; |
---|
| 1598 | long q = FFTPrime[i]; |
---|
| 1599 | double qinv = ((double) 1)/((double) q); |
---|
| 1600 | FFT(sp, yp, y.k, q, Root); |
---|
| 1601 | for (j = 0; j < n; j++) |
---|
| 1602 | yp[j] = MulMod(sp[j], w, q, qinv); |
---|
| 1603 | } |
---|
| 1604 | |
---|
| 1605 | } |
---|
| 1606 | |
---|
| 1607 | void FromFFTRep(ZZ_pX& x, FFTRep& y, long lo, long hi) |
---|
| 1608 | |
---|
| 1609 | // converts from FFT-representation to coefficient representation |
---|
| 1610 | // only the coefficients lo..hi are computed |
---|
| 1611 | |
---|
| 1612 | |
---|
| 1613 | { |
---|
| 1614 | ZZ_pInfo->check(); |
---|
| 1615 | |
---|
| 1616 | long k, n, i, j, l; |
---|
| 1617 | |
---|
| 1618 | vec_long& t = ModularRepBuf; |
---|
| 1619 | vec_long& s = FFTBuf;; |
---|
| 1620 | |
---|
| 1621 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1622 | |
---|
| 1623 | k = y.k; |
---|
| 1624 | n = (1L << k); |
---|
| 1625 | |
---|
| 1626 | s.SetLength(n); |
---|
| 1627 | long *sp = s.elts(); |
---|
| 1628 | |
---|
| 1629 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1630 | long *yp = &y.tbl[i][0]; |
---|
| 1631 | long q = FFTPrime[i]; |
---|
| 1632 | double qinv = ((double) 1)/((double) q); |
---|
| 1633 | long w = TwoInvTable[i][k]; |
---|
| 1634 | long *Root = &RootInvTable[i][0]; |
---|
| 1635 | |
---|
| 1636 | FFT(sp, yp, k, q, Root); |
---|
| 1637 | |
---|
| 1638 | for (j = 0; j < n; j++) yp[j] = MulMod(sp[j], w, q, qinv); |
---|
| 1639 | } |
---|
| 1640 | |
---|
| 1641 | hi = min(hi, n-1); |
---|
| 1642 | l = hi-lo+1; |
---|
| 1643 | l = max(l, 0); |
---|
| 1644 | x.rep.SetLength(l); |
---|
| 1645 | |
---|
| 1646 | for (j = 0; j < l; j++) { |
---|
| 1647 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1648 | t[i] = y.tbl[i][j+lo]; |
---|
| 1649 | |
---|
| 1650 | FromModularRep(x.rep[j], t); |
---|
| 1651 | } |
---|
| 1652 | |
---|
| 1653 | x.normalize(); |
---|
| 1654 | } |
---|
| 1655 | |
---|
| 1656 | void RevFromFFTRep(vec_ZZ_p& x, FFTRep& y, long lo, long hi) |
---|
| 1657 | |
---|
| 1658 | // converts from FFT-representation to coefficient representation |
---|
| 1659 | // using "inverted" evaluation points. |
---|
| 1660 | // only the coefficients lo..hi are computed |
---|
| 1661 | |
---|
| 1662 | |
---|
| 1663 | { |
---|
| 1664 | ZZ_pInfo->check(); |
---|
| 1665 | |
---|
| 1666 | long k, n, i, j, l; |
---|
| 1667 | |
---|
| 1668 | vec_long& t = ModularRepBuf; |
---|
| 1669 | vec_long& s = FFTBuf; |
---|
| 1670 | |
---|
| 1671 | k = y.k; |
---|
| 1672 | n = (1L << k); |
---|
| 1673 | |
---|
| 1674 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1675 | s.SetLength(n); |
---|
| 1676 | long *sp = s.elts(); |
---|
| 1677 | |
---|
| 1678 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1679 | long *yp = &y.tbl[i][0]; |
---|
| 1680 | long q = FFTPrime[i]; |
---|
| 1681 | long *Root = &RootTable[i][0]; |
---|
| 1682 | |
---|
| 1683 | FFT(sp, yp, k, q, Root); |
---|
| 1684 | for (j = 0; j < n; j++) |
---|
| 1685 | yp[j] = sp[j]; |
---|
| 1686 | } |
---|
| 1687 | |
---|
| 1688 | hi = min(hi, n-1); |
---|
| 1689 | l = hi-lo+1; |
---|
| 1690 | l = max(l, 0); |
---|
| 1691 | x.SetLength(l); |
---|
| 1692 | |
---|
| 1693 | for (j = 0; j < l; j++) { |
---|
| 1694 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1695 | t[i] = y.tbl[i][j+lo]; |
---|
| 1696 | |
---|
| 1697 | FromModularRep(x[j], t); |
---|
| 1698 | } |
---|
| 1699 | } |
---|
| 1700 | |
---|
| 1701 | void NDFromFFTRep(ZZ_pX& x, const FFTRep& y, long lo, long hi, FFTRep& z) |
---|
| 1702 | { |
---|
| 1703 | ZZ_pInfo->check(); |
---|
| 1704 | |
---|
| 1705 | long k, n, i, j, l; |
---|
| 1706 | |
---|
| 1707 | vec_long& t = ModularRepBuf; |
---|
| 1708 | |
---|
| 1709 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1710 | k = y.k; |
---|
| 1711 | n = (1L << k); |
---|
| 1712 | |
---|
| 1713 | z.SetSize(k); |
---|
| 1714 | |
---|
| 1715 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1716 | long *zp = &z.tbl[i][0]; |
---|
| 1717 | long q = FFTPrime[i]; |
---|
| 1718 | double qinv = ((double) 1)/((double) q); |
---|
| 1719 | long w = TwoInvTable[i][k]; |
---|
| 1720 | long *Root = &RootInvTable[i][0]; |
---|
| 1721 | |
---|
| 1722 | FFT(zp, &y.tbl[i][0], k, q, Root); |
---|
| 1723 | |
---|
| 1724 | for (j = 0; j < n; j++) zp[j] = MulMod(zp[j], w, q, qinv); |
---|
| 1725 | } |
---|
| 1726 | |
---|
| 1727 | hi = min(hi, n-1); |
---|
| 1728 | l = hi-lo+1; |
---|
| 1729 | l = max(l, 0); |
---|
| 1730 | x.rep.SetLength(l); |
---|
| 1731 | |
---|
| 1732 | for (j = 0; j < l; j++) { |
---|
| 1733 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1734 | t[i] = z.tbl[i][j+lo]; |
---|
| 1735 | |
---|
| 1736 | FromModularRep(x.rep[j], t); |
---|
| 1737 | } |
---|
| 1738 | |
---|
| 1739 | x.normalize(); |
---|
| 1740 | } |
---|
| 1741 | |
---|
| 1742 | void NDFromFFTRep(ZZ_pX& x, FFTRep& y, long lo, long hi) |
---|
| 1743 | { |
---|
| 1744 | FFTRep z; |
---|
| 1745 | NDFromFFTRep(x, y, lo, hi, z); |
---|
| 1746 | } |
---|
| 1747 | |
---|
| 1748 | void FromFFTRep(ZZ_p* x, FFTRep& y, long lo, long hi) |
---|
| 1749 | |
---|
| 1750 | // converts from FFT-representation to coefficient representation |
---|
| 1751 | // only the coefficients lo..hi are computed |
---|
| 1752 | |
---|
| 1753 | |
---|
| 1754 | { |
---|
| 1755 | ZZ_pInfo->check(); |
---|
| 1756 | |
---|
| 1757 | long k, n, i, j; |
---|
| 1758 | |
---|
| 1759 | vec_long& t = ModularRepBuf; |
---|
| 1760 | vec_long& s = FFTBuf; |
---|
| 1761 | |
---|
| 1762 | k = y.k; |
---|
| 1763 | n = (1L << k); |
---|
| 1764 | |
---|
| 1765 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1766 | s.SetLength(n); |
---|
| 1767 | long *sp = s.elts(); |
---|
| 1768 | |
---|
| 1769 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1770 | long *yp = &y.tbl[i][0]; |
---|
| 1771 | long q = FFTPrime[i]; |
---|
| 1772 | double qinv = ((double) 1)/((double) q); |
---|
| 1773 | long w = TwoInvTable[i][k]; |
---|
| 1774 | long *Root = &RootInvTable[i][0]; |
---|
| 1775 | |
---|
| 1776 | FFT(sp, yp, k, q, Root); |
---|
| 1777 | |
---|
| 1778 | for (j = 0; j < n; j++) yp[j] = MulMod(sp[j], w, q, qinv); |
---|
| 1779 | } |
---|
| 1780 | |
---|
| 1781 | for (j = lo; j <= hi; j++) { |
---|
| 1782 | if (j >= n) |
---|
| 1783 | clear(x[j-lo]); |
---|
| 1784 | else { |
---|
| 1785 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1786 | t[i] = y.tbl[i][j]; |
---|
| 1787 | |
---|
| 1788 | FromModularRep(x[j-lo], t); |
---|
| 1789 | } |
---|
| 1790 | } |
---|
| 1791 | } |
---|
| 1792 | |
---|
| 1793 | |
---|
| 1794 | void mul(FFTRep& z, const FFTRep& x, const FFTRep& y) |
---|
| 1795 | { |
---|
| 1796 | ZZ_pInfo->check(); |
---|
| 1797 | |
---|
| 1798 | long k, n, i, j; |
---|
| 1799 | |
---|
| 1800 | if (x.k != y.k) Error("FFT rep mismatch"); |
---|
| 1801 | |
---|
| 1802 | k = x.k; |
---|
| 1803 | n = 1L << k; |
---|
| 1804 | |
---|
| 1805 | z.SetSize(k); |
---|
| 1806 | |
---|
| 1807 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1808 | long *zp = &z.tbl[i][0]; |
---|
| 1809 | const long *xp = &x.tbl[i][0]; |
---|
| 1810 | const long *yp = &y.tbl[i][0]; |
---|
| 1811 | long q = FFTPrime[i]; |
---|
| 1812 | double qinv = ((double) 1)/((double) q); |
---|
| 1813 | |
---|
| 1814 | for (j = 0; j < n; j++) |
---|
| 1815 | zp[j] = MulMod(xp[j], yp[j], q, qinv); |
---|
| 1816 | } |
---|
| 1817 | |
---|
| 1818 | } |
---|
| 1819 | |
---|
| 1820 | void sub(FFTRep& z, const FFTRep& x, const FFTRep& y) |
---|
| 1821 | { |
---|
| 1822 | ZZ_pInfo->check(); |
---|
| 1823 | |
---|
| 1824 | long k, n, i, j; |
---|
| 1825 | |
---|
| 1826 | if (x.k != y.k) Error("FFT rep mismatch"); |
---|
| 1827 | |
---|
| 1828 | k = x.k; |
---|
| 1829 | n = 1L << k; |
---|
| 1830 | |
---|
| 1831 | z.SetSize(k); |
---|
| 1832 | |
---|
| 1833 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1834 | long *zp = &z.tbl[i][0]; |
---|
| 1835 | const long *xp = &x.tbl[i][0]; |
---|
| 1836 | const long *yp = &y.tbl[i][0]; |
---|
| 1837 | long q = FFTPrime[i]; |
---|
| 1838 | |
---|
| 1839 | for (j = 0; j < n; j++) |
---|
| 1840 | zp[j] = SubMod(xp[j], yp[j], q); |
---|
| 1841 | } |
---|
| 1842 | } |
---|
| 1843 | |
---|
| 1844 | void add(FFTRep& z, const FFTRep& x, const FFTRep& y) |
---|
| 1845 | { |
---|
| 1846 | ZZ_pInfo->check(); |
---|
| 1847 | |
---|
| 1848 | long k, n, i, j; |
---|
| 1849 | |
---|
| 1850 | if (x.k != y.k) Error("FFT rep mismatch"); |
---|
| 1851 | |
---|
| 1852 | k = x.k; |
---|
| 1853 | n = 1L << k; |
---|
| 1854 | |
---|
| 1855 | z.SetSize(k); |
---|
| 1856 | |
---|
| 1857 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1858 | long *zp = &z.tbl[i][0]; |
---|
| 1859 | const long *xp = &x.tbl[i][0]; |
---|
| 1860 | const long *yp = &y.tbl[i][0]; |
---|
| 1861 | long q = FFTPrime[i]; |
---|
| 1862 | |
---|
| 1863 | for (j = 0; j < n; j++) |
---|
| 1864 | zp[j] = AddMod(xp[j], yp[j], q); |
---|
| 1865 | } |
---|
| 1866 | } |
---|
| 1867 | |
---|
| 1868 | |
---|
| 1869 | void reduce(FFTRep& x, const FFTRep& a, long k) |
---|
| 1870 | // reduces a 2^l point FFT-rep to a 2^k point FFT-rep |
---|
| 1871 | // input may alias output |
---|
| 1872 | { |
---|
| 1873 | ZZ_pInfo->check(); |
---|
| 1874 | |
---|
| 1875 | long i, j, l, n; |
---|
| 1876 | long* xp; |
---|
| 1877 | const long* ap; |
---|
| 1878 | |
---|
| 1879 | l = a.k; |
---|
| 1880 | n = 1L << k; |
---|
| 1881 | |
---|
| 1882 | if (l < k) Error("reduce: bad operands"); |
---|
| 1883 | |
---|
| 1884 | x.SetSize(k); |
---|
| 1885 | |
---|
| 1886 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1887 | ap = &a.tbl[i][0]; |
---|
| 1888 | xp = &x.tbl[i][0]; |
---|
| 1889 | for (j = 0; j < n; j++) |
---|
| 1890 | xp[j] = ap[j << (l-k)]; |
---|
| 1891 | } |
---|
| 1892 | } |
---|
| 1893 | |
---|
| 1894 | void AddExpand(FFTRep& x, const FFTRep& a) |
---|
| 1895 | // x = x + (an "expanded" version of a) |
---|
| 1896 | { |
---|
| 1897 | ZZ_pInfo->check(); |
---|
| 1898 | |
---|
| 1899 | long i, j, l, k, n; |
---|
| 1900 | |
---|
| 1901 | l = x.k; |
---|
| 1902 | k = a.k; |
---|
| 1903 | n = 1L << k; |
---|
| 1904 | |
---|
| 1905 | if (l < k) Error("AddExpand: bad args"); |
---|
| 1906 | |
---|
| 1907 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) { |
---|
| 1908 | long q = FFTPrime[i]; |
---|
| 1909 | const long *ap = &a.tbl[i][0]; |
---|
| 1910 | long *xp = &x.tbl[i][0]; |
---|
| 1911 | for (j = 0; j < n; j++) { |
---|
| 1912 | long j1 = j << (l-k); |
---|
| 1913 | xp[j1] = AddMod(xp[j1], ap[j], q); |
---|
| 1914 | } |
---|
| 1915 | } |
---|
| 1916 | } |
---|
| 1917 | |
---|
| 1918 | |
---|
| 1919 | |
---|
| 1920 | void ToZZ_pXModRep(ZZ_pXModRep& y, const ZZ_pX& x, long lo, long hi) |
---|
| 1921 | { |
---|
| 1922 | ZZ_pInfo->check(); |
---|
| 1923 | |
---|
| 1924 | long n, i, j; |
---|
| 1925 | vec_long& t = ModularRepBuf; |
---|
| 1926 | |
---|
| 1927 | t.SetLength(ZZ_pInfo->NumPrimes); |
---|
| 1928 | |
---|
| 1929 | if (lo < 0) |
---|
| 1930 | Error("bad arg to ToZZ_pXModRep"); |
---|
| 1931 | hi = min(hi, deg(x)); |
---|
| 1932 | n = max(hi-lo+1, 0); |
---|
| 1933 | |
---|
| 1934 | y.SetSize(n); |
---|
| 1935 | |
---|
| 1936 | const ZZ_p *xx = x.rep.elts(); |
---|
| 1937 | |
---|
| 1938 | for (j = 0; j < n; j++) { |
---|
| 1939 | ToModularRep(t, xx[j+lo]); |
---|
| 1940 | for (i = 0; i < ZZ_pInfo->NumPrimes; i++) |
---|
| 1941 | y.tbl[i][j] = t[i]; |
---|
| 1942 | } |
---|
| 1943 | } |
---|
| 1944 | |
---|
| 1945 | |
---|
| 1946 | void ToFFTRep(FFTRep& x, const ZZ_pXModRep& a, long k, long lo, long hi) |
---|
| 1947 | { |
---|
| 1948 | ZZ_pInfo->check(); |
---|
| 1949 | |
---|
| 1950 | vec_long s; |
---|
| 1951 | long n, m, i, j; |
---|
| 1952 | |
---|
| 1953 | if (k < 0 || lo < 0) |
---|
| 1954 | Error("bad args to ToFFTRep"); |
---|
| 1955 | |
---|
| 1956 | if (hi > a.n-1) hi = a.n-1; |
---|
| 1957 | |
---|
| 1958 | n = 1L << k; |
---|
| 1959 | m = max(hi-lo+1, 0); |
---|
| 1960 | |
---|
| 1961 | if (m > n) |
---|
| 1962 | Error("bad args to ToFFTRep"); |
---|
| 1963 | |
---|
| 1964 | s.SetLength(n); |
---|
| 1965 | long *sp = s.elts(); |
---|
| 1966 | |
---|
| 1967 | x.SetSize(k); |
---|
| 1968 | |
---|
| 1969 | long NumPrimes = ZZ_pInfo->NumPrimes; |
---|
| 1970 | |
---|
| 1971 | for (i = 0; i < NumPrimes; i++) { |
---|
| 1972 | long *Root = &RootTable[i][0]; |
---|
| 1973 | long *xp = &x.tbl[i][0]; |
---|
| 1974 | long *ap = (m == 0 ? 0 : &a.tbl[i][0]); |
---|
| 1975 | for (j = 0; j < m; j++) |
---|
| 1976 | sp[j] = ap[lo+j]; |
---|
| 1977 | for (j = m; j < n; j++) |
---|
| 1978 | sp[j] = 0; |
---|
| 1979 | |
---|
| 1980 | FFT(xp, sp, k, FFTPrime[i], Root); |
---|
| 1981 | } |
---|
| 1982 | } |
---|
| 1983 | |
---|
| 1984 | |
---|
| 1985 | |
---|
| 1986 | |
---|
| 1987 | |
---|
| 1988 | |
---|
| 1989 | void FFTMul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 1990 | { |
---|
| 1991 | long k, d; |
---|
| 1992 | |
---|
| 1993 | if (IsZero(a) || IsZero(b)) { |
---|
| 1994 | clear(x); |
---|
| 1995 | return; |
---|
| 1996 | } |
---|
| 1997 | |
---|
| 1998 | d = deg(a) + deg(b); |
---|
| 1999 | k = NextPowerOfTwo(d+1); |
---|
| 2000 | |
---|
| 2001 | FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k); |
---|
| 2002 | |
---|
| 2003 | ToFFTRep(R1, a, k); |
---|
| 2004 | ToFFTRep(R2, b, k); |
---|
| 2005 | mul(R1, R1, R2); |
---|
| 2006 | FromFFTRep(x, R1, 0, d); |
---|
| 2007 | } |
---|
| 2008 | |
---|
| 2009 | void FFTSqr(ZZ_pX& x, const ZZ_pX& a) |
---|
| 2010 | { |
---|
| 2011 | long k, d; |
---|
| 2012 | |
---|
| 2013 | if (IsZero(a)) { |
---|
| 2014 | clear(x); |
---|
| 2015 | return; |
---|
| 2016 | } |
---|
| 2017 | |
---|
| 2018 | d = 2*deg(a); |
---|
| 2019 | k = NextPowerOfTwo(d+1); |
---|
| 2020 | |
---|
| 2021 | FFTRep R1(INIT_SIZE, k); |
---|
| 2022 | |
---|
| 2023 | ToFFTRep(R1, a, k); |
---|
| 2024 | mul(R1, R1, R1); |
---|
| 2025 | FromFFTRep(x, R1, 0, d); |
---|
| 2026 | } |
---|
| 2027 | |
---|
| 2028 | |
---|
| 2029 | void CopyReverse(ZZ_pX& x, const ZZ_pX& a, long lo, long hi) |
---|
| 2030 | |
---|
| 2031 | // x[0..hi-lo] = reverse(a[lo..hi]), with zero fill |
---|
| 2032 | // input may not alias output |
---|
| 2033 | |
---|
| 2034 | { |
---|
| 2035 | long i, j, n, m; |
---|
| 2036 | |
---|
| 2037 | n = hi-lo+1; |
---|
| 2038 | m = a.rep.length(); |
---|
| 2039 | |
---|
| 2040 | x.rep.SetLength(n); |
---|
| 2041 | |
---|
| 2042 | const ZZ_p* ap = a.rep.elts(); |
---|
| 2043 | ZZ_p* xp = x.rep.elts(); |
---|
| 2044 | |
---|
| 2045 | for (i = 0; i < n; i++) { |
---|
| 2046 | j = hi-i; |
---|
| 2047 | if (j < 0 || j >= m) |
---|
| 2048 | clear(xp[i]); |
---|
| 2049 | else |
---|
| 2050 | xp[i] = ap[j]; |
---|
| 2051 | } |
---|
| 2052 | |
---|
| 2053 | x.normalize(); |
---|
| 2054 | } |
---|
| 2055 | |
---|
| 2056 | void copy(ZZ_pX& x, const ZZ_pX& a, long lo, long hi) |
---|
| 2057 | |
---|
| 2058 | // x[0..hi-lo] = a[lo..hi], with zero fill |
---|
| 2059 | // input may not alias output |
---|
| 2060 | |
---|
| 2061 | { |
---|
| 2062 | long i, j, n, m; |
---|
| 2063 | |
---|
| 2064 | n = hi-lo+1; |
---|
| 2065 | m = a.rep.length(); |
---|
| 2066 | |
---|
| 2067 | x.rep.SetLength(n); |
---|
| 2068 | |
---|
| 2069 | const ZZ_p* ap = a.rep.elts(); |
---|
| 2070 | ZZ_p* xp = x.rep.elts(); |
---|
| 2071 | |
---|
| 2072 | for (i = 0; i < n; i++) { |
---|
| 2073 | j = lo + i; |
---|
| 2074 | if (j < 0 || j >= m) |
---|
| 2075 | clear(xp[i]); |
---|
| 2076 | else |
---|
| 2077 | xp[i] = ap[j]; |
---|
| 2078 | } |
---|
| 2079 | |
---|
| 2080 | x.normalize(); |
---|
| 2081 | } |
---|
| 2082 | |
---|
| 2083 | |
---|
| 2084 | void rem21(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2085 | { |
---|
| 2086 | long i, da, ds, n, kk; |
---|
| 2087 | |
---|
| 2088 | da = deg(a); |
---|
| 2089 | n = F.n; |
---|
| 2090 | |
---|
| 2091 | if (da > 2*n-2) |
---|
| 2092 | Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)"); |
---|
| 2093 | |
---|
| 2094 | |
---|
| 2095 | if (da < n) { |
---|
| 2096 | x = a; |
---|
| 2097 | return; |
---|
| 2098 | } |
---|
| 2099 | |
---|
| 2100 | if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2101 | PlainRem(x, a, F.f); |
---|
| 2102 | return; |
---|
| 2103 | } |
---|
| 2104 | |
---|
| 2105 | FFTRep R1(INIT_SIZE, F.l); |
---|
| 2106 | ZZ_pX P1(INIT_SIZE, n); |
---|
| 2107 | |
---|
| 2108 | ToFFTRep(R1, a, F.l, n, 2*(n-1)); |
---|
| 2109 | mul(R1, R1, F.HRep); |
---|
| 2110 | FromFFTRep(P1, R1, n-2, 2*n-4); |
---|
| 2111 | |
---|
| 2112 | ToFFTRep(R1, P1, F.k); |
---|
| 2113 | mul(R1, R1, F.FRep); |
---|
| 2114 | FromFFTRep(P1, R1, 0, n-1); |
---|
| 2115 | |
---|
| 2116 | ds = deg(P1); |
---|
| 2117 | |
---|
| 2118 | kk = 1L << F.k; |
---|
| 2119 | |
---|
| 2120 | x.rep.SetLength(n); |
---|
| 2121 | const ZZ_p* aa = a.rep.elts(); |
---|
| 2122 | const ZZ_p* ss = P1.rep.elts(); |
---|
| 2123 | ZZ_p* xx = x.rep.elts(); |
---|
| 2124 | |
---|
| 2125 | for (i = 0; i < n; i++) { |
---|
| 2126 | if (i <= ds) |
---|
| 2127 | sub(xx[i], aa[i], ss[i]); |
---|
| 2128 | else |
---|
| 2129 | xx[i] = aa[i]; |
---|
| 2130 | |
---|
| 2131 | if (i + kk <= da) |
---|
| 2132 | add(xx[i], xx[i], aa[i+kk]); |
---|
| 2133 | } |
---|
| 2134 | |
---|
| 2135 | x.normalize(); |
---|
| 2136 | } |
---|
| 2137 | |
---|
| 2138 | void DivRem21(ZZ_pX& q, ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2139 | { |
---|
| 2140 | long i, da, ds, n, kk; |
---|
| 2141 | |
---|
| 2142 | da = deg(a); |
---|
| 2143 | n = F.n; |
---|
| 2144 | |
---|
| 2145 | if (da > 2*n-2) |
---|
| 2146 | Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)"); |
---|
| 2147 | |
---|
| 2148 | |
---|
| 2149 | if (da < n) { |
---|
| 2150 | x = a; |
---|
| 2151 | clear(q); |
---|
| 2152 | return; |
---|
| 2153 | } |
---|
| 2154 | |
---|
| 2155 | if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2156 | PlainDivRem(q, x, a, F.f); |
---|
| 2157 | return; |
---|
| 2158 | } |
---|
| 2159 | |
---|
| 2160 | FFTRep R1(INIT_SIZE, F.l); |
---|
| 2161 | ZZ_pX P1(INIT_SIZE, n), qq; |
---|
| 2162 | |
---|
| 2163 | ToFFTRep(R1, a, F.l, n, 2*(n-1)); |
---|
| 2164 | mul(R1, R1, F.HRep); |
---|
| 2165 | FromFFTRep(P1, R1, n-2, 2*n-4); |
---|
| 2166 | qq = P1; |
---|
| 2167 | |
---|
| 2168 | ToFFTRep(R1, P1, F.k); |
---|
| 2169 | mul(R1, R1, F.FRep); |
---|
| 2170 | FromFFTRep(P1, R1, 0, n-1); |
---|
| 2171 | |
---|
| 2172 | ds = deg(P1); |
---|
| 2173 | |
---|
| 2174 | kk = 1L << F.k; |
---|
| 2175 | |
---|
| 2176 | x.rep.SetLength(n); |
---|
| 2177 | const ZZ_p* aa = a.rep.elts(); |
---|
| 2178 | const ZZ_p* ss = P1.rep.elts(); |
---|
| 2179 | ZZ_p* xx = x.rep.elts(); |
---|
| 2180 | |
---|
| 2181 | for (i = 0; i < n; i++) { |
---|
| 2182 | if (i <= ds) |
---|
| 2183 | sub(xx[i], aa[i], ss[i]); |
---|
| 2184 | else |
---|
| 2185 | xx[i] = aa[i]; |
---|
| 2186 | |
---|
| 2187 | if (i + kk <= da) |
---|
| 2188 | add(xx[i], xx[i], aa[i+kk]); |
---|
| 2189 | } |
---|
| 2190 | |
---|
| 2191 | x.normalize(); |
---|
| 2192 | q = qq; |
---|
| 2193 | } |
---|
| 2194 | |
---|
| 2195 | void div21(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2196 | { |
---|
| 2197 | long da, n; |
---|
| 2198 | |
---|
| 2199 | da = deg(a); |
---|
| 2200 | n = F.n; |
---|
| 2201 | |
---|
| 2202 | if (da > 2*n-2) |
---|
| 2203 | Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)"); |
---|
| 2204 | |
---|
| 2205 | |
---|
| 2206 | if (da < n) { |
---|
| 2207 | clear(x); |
---|
| 2208 | return; |
---|
| 2209 | } |
---|
| 2210 | |
---|
| 2211 | if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2212 | PlainDiv(x, a, F.f); |
---|
| 2213 | return; |
---|
| 2214 | } |
---|
| 2215 | |
---|
| 2216 | FFTRep R1(INIT_SIZE, F.l); |
---|
| 2217 | ZZ_pX P1(INIT_SIZE, n); |
---|
| 2218 | |
---|
| 2219 | ToFFTRep(R1, a, F.l, n, 2*(n-1)); |
---|
| 2220 | mul(R1, R1, F.HRep); |
---|
| 2221 | FromFFTRep(x, R1, n-2, 2*n-4); |
---|
| 2222 | } |
---|
| 2223 | |
---|
| 2224 | |
---|
| 2225 | void rem(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2226 | { |
---|
| 2227 | long da = deg(a); |
---|
| 2228 | long n = F.n; |
---|
| 2229 | |
---|
| 2230 | if (n < 0) Error("rem: unitialized modulus"); |
---|
| 2231 | |
---|
| 2232 | if (da <= 2*n-2) { |
---|
| 2233 | rem21(x, a, F); |
---|
| 2234 | return; |
---|
| 2235 | } |
---|
| 2236 | else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2237 | PlainRem(x, a, F.f); |
---|
| 2238 | return; |
---|
| 2239 | } |
---|
| 2240 | |
---|
| 2241 | ZZ_pX buf(INIT_SIZE, 2*n-1); |
---|
| 2242 | |
---|
| 2243 | long a_len = da+1; |
---|
| 2244 | |
---|
| 2245 | while (a_len > 0) { |
---|
| 2246 | long old_buf_len = buf.rep.length(); |
---|
| 2247 | long amt = min(2*n-1-old_buf_len, a_len); |
---|
| 2248 | |
---|
| 2249 | buf.rep.SetLength(old_buf_len+amt); |
---|
| 2250 | |
---|
| 2251 | long i; |
---|
| 2252 | |
---|
| 2253 | for (i = old_buf_len+amt-1; i >= amt; i--) |
---|
| 2254 | buf.rep[i] = buf.rep[i-amt]; |
---|
| 2255 | |
---|
| 2256 | for (i = amt-1; i >= 0; i--) |
---|
| 2257 | buf.rep[i] = a.rep[a_len-amt+i]; |
---|
| 2258 | |
---|
| 2259 | buf.normalize(); |
---|
| 2260 | |
---|
| 2261 | rem21(buf, buf, F); |
---|
| 2262 | |
---|
| 2263 | a_len -= amt; |
---|
| 2264 | } |
---|
| 2265 | |
---|
| 2266 | x = buf; |
---|
| 2267 | } |
---|
| 2268 | |
---|
| 2269 | void DivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2270 | { |
---|
| 2271 | long da = deg(a); |
---|
| 2272 | long n = F.n; |
---|
| 2273 | |
---|
| 2274 | if (n < 0) Error("uninitialized modulus"); |
---|
| 2275 | |
---|
| 2276 | if (da <= 2*n-2) { |
---|
| 2277 | DivRem21(q, r, a, F); |
---|
| 2278 | return; |
---|
| 2279 | } |
---|
| 2280 | else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2281 | PlainDivRem(q, r, a, F.f); |
---|
| 2282 | return; |
---|
| 2283 | } |
---|
| 2284 | |
---|
| 2285 | ZZ_pX buf(INIT_SIZE, 2*n-1); |
---|
| 2286 | ZZ_pX qbuf(INIT_SIZE, n-1); |
---|
| 2287 | |
---|
| 2288 | ZZ_pX qq; |
---|
| 2289 | qq.rep.SetLength(da-n+1); |
---|
| 2290 | |
---|
| 2291 | long a_len = da+1; |
---|
| 2292 | long q_hi = da-n+1; |
---|
| 2293 | |
---|
| 2294 | while (a_len > 0) { |
---|
| 2295 | long old_buf_len = buf.rep.length(); |
---|
| 2296 | long amt = min(2*n-1-old_buf_len, a_len); |
---|
| 2297 | |
---|
| 2298 | buf.rep.SetLength(old_buf_len+amt); |
---|
| 2299 | |
---|
| 2300 | long i; |
---|
| 2301 | |
---|
| 2302 | for (i = old_buf_len+amt-1; i >= amt; i--) |
---|
| 2303 | buf.rep[i] = buf.rep[i-amt]; |
---|
| 2304 | |
---|
| 2305 | for (i = amt-1; i >= 0; i--) |
---|
| 2306 | buf.rep[i] = a.rep[a_len-amt+i]; |
---|
| 2307 | |
---|
| 2308 | buf.normalize(); |
---|
| 2309 | |
---|
| 2310 | DivRem21(qbuf, buf, buf, F); |
---|
| 2311 | long dl = qbuf.rep.length(); |
---|
| 2312 | a_len = a_len - amt; |
---|
| 2313 | for(i = 0; i < dl; i++) |
---|
| 2314 | qq.rep[a_len+i] = qbuf.rep[i]; |
---|
| 2315 | for(i = dl+a_len; i < q_hi; i++) |
---|
| 2316 | clear(qq.rep[i]); |
---|
| 2317 | q_hi = a_len; |
---|
| 2318 | } |
---|
| 2319 | |
---|
| 2320 | r = buf; |
---|
| 2321 | |
---|
| 2322 | qq.normalize(); |
---|
| 2323 | q = qq; |
---|
| 2324 | } |
---|
| 2325 | |
---|
| 2326 | void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2327 | { |
---|
| 2328 | long da = deg(a); |
---|
| 2329 | long n = F.n; |
---|
| 2330 | |
---|
| 2331 | if (n < 0) Error("uninitialized modulus"); |
---|
| 2332 | |
---|
| 2333 | if (da <= 2*n-2) { |
---|
| 2334 | div21(q, a, F); |
---|
| 2335 | return; |
---|
| 2336 | } |
---|
| 2337 | else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2338 | PlainDiv(q, a, F.f); |
---|
| 2339 | return; |
---|
| 2340 | } |
---|
| 2341 | |
---|
| 2342 | ZZ_pX buf(INIT_SIZE, 2*n-1); |
---|
| 2343 | ZZ_pX qbuf(INIT_SIZE, n-1); |
---|
| 2344 | |
---|
| 2345 | ZZ_pX qq; |
---|
| 2346 | qq.rep.SetLength(da-n+1); |
---|
| 2347 | |
---|
| 2348 | long a_len = da+1; |
---|
| 2349 | long q_hi = da-n+1; |
---|
| 2350 | |
---|
| 2351 | while (a_len > 0) { |
---|
| 2352 | long old_buf_len = buf.rep.length(); |
---|
| 2353 | long amt = min(2*n-1-old_buf_len, a_len); |
---|
| 2354 | |
---|
| 2355 | buf.rep.SetLength(old_buf_len+amt); |
---|
| 2356 | |
---|
| 2357 | long i; |
---|
| 2358 | |
---|
| 2359 | for (i = old_buf_len+amt-1; i >= amt; i--) |
---|
| 2360 | buf.rep[i] = buf.rep[i-amt]; |
---|
| 2361 | |
---|
| 2362 | for (i = amt-1; i >= 0; i--) |
---|
| 2363 | buf.rep[i] = a.rep[a_len-amt+i]; |
---|
| 2364 | |
---|
| 2365 | buf.normalize(); |
---|
| 2366 | |
---|
| 2367 | a_len = a_len - amt; |
---|
| 2368 | if (a_len > 0) |
---|
| 2369 | DivRem21(qbuf, buf, buf, F); |
---|
| 2370 | else |
---|
| 2371 | div21(qbuf, buf, F); |
---|
| 2372 | |
---|
| 2373 | long dl = qbuf.rep.length(); |
---|
| 2374 | for(i = 0; i < dl; i++) |
---|
| 2375 | qq.rep[a_len+i] = qbuf.rep[i]; |
---|
| 2376 | for(i = dl+a_len; i < q_hi; i++) |
---|
| 2377 | clear(qq.rep[i]); |
---|
| 2378 | q_hi = a_len; |
---|
| 2379 | } |
---|
| 2380 | |
---|
| 2381 | qq.normalize(); |
---|
| 2382 | q = qq; |
---|
| 2383 | } |
---|
| 2384 | |
---|
| 2385 | |
---|
| 2386 | |
---|
| 2387 | void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, const ZZ_pXModulus& F) |
---|
| 2388 | { |
---|
| 2389 | long da, db, d, n, k; |
---|
| 2390 | |
---|
| 2391 | da = deg(a); |
---|
| 2392 | db = deg(b); |
---|
| 2393 | n = F.n; |
---|
| 2394 | |
---|
| 2395 | if (n < 0) Error("MulMod: uninitialized modulus"); |
---|
| 2396 | |
---|
| 2397 | if (da >= n || db >= n) |
---|
| 2398 | Error("bad args to MulMod(ZZ_pX,ZZ_pX,ZZ_pX,ZZ_pXModulus)"); |
---|
| 2399 | |
---|
| 2400 | if (da < 0 || db < 0) { |
---|
| 2401 | clear(x); |
---|
| 2402 | return; |
---|
| 2403 | } |
---|
| 2404 | |
---|
| 2405 | if (!F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER || db <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2406 | ZZ_pX P1; |
---|
| 2407 | mul(P1, a, b); |
---|
| 2408 | rem(x, P1, F); |
---|
| 2409 | return; |
---|
| 2410 | } |
---|
| 2411 | |
---|
| 2412 | d = da + db + 1; |
---|
| 2413 | |
---|
| 2414 | k = NextPowerOfTwo(d); |
---|
| 2415 | k = max(k, F.k); |
---|
| 2416 | |
---|
| 2417 | FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, F.l); |
---|
| 2418 | ZZ_pX P1(INIT_SIZE, n); |
---|
| 2419 | |
---|
| 2420 | ToFFTRep(R1, a, k); |
---|
| 2421 | ToFFTRep(R2, b, k); |
---|
| 2422 | mul(R1, R1, R2); |
---|
| 2423 | NDFromFFTRep(P1, R1, n, d-1, R2); // save R1 for future use |
---|
| 2424 | |
---|
| 2425 | ToFFTRep(R2, P1, F.l); |
---|
| 2426 | mul(R2, R2, F.HRep); |
---|
| 2427 | FromFFTRep(P1, R2, n-2, 2*n-4); |
---|
| 2428 | |
---|
| 2429 | ToFFTRep(R2, P1, F.k); |
---|
| 2430 | mul(R2, R2, F.FRep); |
---|
| 2431 | reduce(R1, R1, F.k); |
---|
| 2432 | sub(R1, R1, R2); |
---|
| 2433 | FromFFTRep(x, R1, 0, n-1); |
---|
| 2434 | } |
---|
| 2435 | |
---|
| 2436 | void SqrMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F) |
---|
| 2437 | { |
---|
| 2438 | long da, d, n, k; |
---|
| 2439 | |
---|
| 2440 | da = deg(a); |
---|
| 2441 | n = F.n; |
---|
| 2442 | |
---|
| 2443 | if (n < 0) Error("SqrMod: uninitailized modulus"); |
---|
| 2444 | |
---|
| 2445 | if (da >= n) |
---|
| 2446 | Error("bad args to SqrMod(ZZ_pX,ZZ_pX,ZZ_pXModulus)"); |
---|
| 2447 | |
---|
| 2448 | if (!F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2449 | ZZ_pX P1; |
---|
| 2450 | sqr(P1, a); |
---|
| 2451 | rem(x, P1, F); |
---|
| 2452 | return; |
---|
| 2453 | } |
---|
| 2454 | |
---|
| 2455 | |
---|
| 2456 | d = 2*da + 1; |
---|
| 2457 | |
---|
| 2458 | k = NextPowerOfTwo(d); |
---|
| 2459 | k = max(k, F.k); |
---|
| 2460 | |
---|
| 2461 | FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, F.l); |
---|
| 2462 | ZZ_pX P1(INIT_SIZE, n); |
---|
| 2463 | |
---|
| 2464 | ToFFTRep(R1, a, k); |
---|
| 2465 | mul(R1, R1, R1); |
---|
| 2466 | NDFromFFTRep(P1, R1, n, d-1, R2); // save R1 for future use |
---|
| 2467 | |
---|
| 2468 | ToFFTRep(R2, P1, F.l); |
---|
| 2469 | mul(R2, R2, F.HRep); |
---|
| 2470 | FromFFTRep(P1, R2, n-2, 2*n-4); |
---|
| 2471 | |
---|
| 2472 | ToFFTRep(R2, P1, F.k); |
---|
| 2473 | mul(R2, R2, F.FRep); |
---|
| 2474 | reduce(R1, R1, F.k); |
---|
| 2475 | sub(R1, R1, R2); |
---|
| 2476 | FromFFTRep(x, R1, 0, n-1); |
---|
| 2477 | } |
---|
| 2478 | |
---|
| 2479 | void PlainInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m) |
---|
| 2480 | |
---|
| 2481 | /* x = (1/a) % X^m, input not output, constant term a is nonzero */ |
---|
| 2482 | |
---|
| 2483 | { |
---|
| 2484 | long i, k, n, lb; |
---|
| 2485 | static ZZ v, t; |
---|
| 2486 | ZZ_p s; |
---|
| 2487 | const ZZ_p* ap; |
---|
| 2488 | ZZ_p* xp; |
---|
| 2489 | |
---|
| 2490 | |
---|
| 2491 | n = deg(a); |
---|
| 2492 | |
---|
| 2493 | if (n < 0) Error("division by zero"); |
---|
| 2494 | |
---|
| 2495 | inv(s, ConstTerm(a)); |
---|
| 2496 | |
---|
| 2497 | if (n == 0) { |
---|
| 2498 | conv(x, s); |
---|
| 2499 | return; |
---|
| 2500 | } |
---|
| 2501 | |
---|
| 2502 | ap = a.rep.elts(); |
---|
| 2503 | x.rep.SetLength(m); |
---|
| 2504 | xp = x.rep.elts(); |
---|
| 2505 | |
---|
| 2506 | xp[0] = s; |
---|
| 2507 | |
---|
| 2508 | long is_one = IsOne(s); |
---|
| 2509 | |
---|
| 2510 | for (k = 1; k < m; k++) { |
---|
| 2511 | clear(v); |
---|
| 2512 | lb = max(k-n, 0); |
---|
| 2513 | for (i = lb; i <= k-1; i++) { |
---|
| 2514 | mul(t, rep(xp[i]), rep(ap[k-i])); |
---|
| 2515 | add(v, v, t); |
---|
| 2516 | } |
---|
| 2517 | conv(xp[k], v); |
---|
| 2518 | negate(xp[k], xp[k]); |
---|
| 2519 | if (!is_one) mul(xp[k], xp[k], s); |
---|
| 2520 | } |
---|
| 2521 | |
---|
| 2522 | x.normalize(); |
---|
| 2523 | } |
---|
| 2524 | |
---|
| 2525 | |
---|
| 2526 | void trunc(ZZ_pX& x, const ZZ_pX& a, long m) |
---|
| 2527 | |
---|
| 2528 | // x = a % X^m, output may alias input |
---|
| 2529 | |
---|
| 2530 | { |
---|
| 2531 | if (m < 0) Error("trunc: bad args"); |
---|
| 2532 | |
---|
| 2533 | if (&x == &a) { |
---|
| 2534 | if (x.rep.length() > m) { |
---|
| 2535 | x.rep.SetLength(m); |
---|
| 2536 | x.normalize(); |
---|
| 2537 | } |
---|
| 2538 | } |
---|
| 2539 | else { |
---|
| 2540 | long n; |
---|
| 2541 | long i; |
---|
| 2542 | ZZ_p* xp; |
---|
| 2543 | const ZZ_p* ap; |
---|
| 2544 | |
---|
| 2545 | n = min(a.rep.length(), m); |
---|
| 2546 | x.rep.SetLength(n); |
---|
| 2547 | |
---|
| 2548 | xp = x.rep.elts(); |
---|
| 2549 | ap = a.rep.elts(); |
---|
| 2550 | |
---|
| 2551 | for (i = 0; i < n; i++) xp[i] = ap[i]; |
---|
| 2552 | |
---|
| 2553 | x.normalize(); |
---|
| 2554 | } |
---|
| 2555 | } |
---|
| 2556 | |
---|
| 2557 | void CyclicReduce(ZZ_pX& x, const ZZ_pX& a, long m) |
---|
| 2558 | |
---|
| 2559 | // computes x = a mod X^m-1 |
---|
| 2560 | |
---|
| 2561 | { |
---|
| 2562 | long n = deg(a); |
---|
| 2563 | long i, j; |
---|
| 2564 | ZZ_p accum; |
---|
| 2565 | |
---|
| 2566 | if (n < m) { |
---|
| 2567 | x = a; |
---|
| 2568 | return; |
---|
| 2569 | } |
---|
| 2570 | |
---|
| 2571 | if (&x != &a) |
---|
| 2572 | x.rep.SetLength(m); |
---|
| 2573 | |
---|
| 2574 | for (i = 0; i < m; i++) { |
---|
| 2575 | accum = a.rep[i]; |
---|
| 2576 | for (j = i + m; j <= n; j += m) |
---|
| 2577 | add(accum, accum, a.rep[j]); |
---|
| 2578 | x.rep[i] = accum; |
---|
| 2579 | } |
---|
| 2580 | |
---|
| 2581 | if (&x == &a) |
---|
| 2582 | x.rep.SetLength(m); |
---|
| 2583 | |
---|
| 2584 | x.normalize(); |
---|
| 2585 | } |
---|
| 2586 | |
---|
| 2587 | |
---|
| 2588 | |
---|
| 2589 | void InvTrunc(ZZ_pX& x, const ZZ_pX& a, long m) |
---|
| 2590 | { |
---|
| 2591 | if (m < 0) Error("InvTrunc: bad args"); |
---|
| 2592 | |
---|
| 2593 | if (m == 0) { |
---|
| 2594 | clear(x); |
---|
| 2595 | return; |
---|
| 2596 | } |
---|
| 2597 | |
---|
[09da99] | 2598 | if (NTL_OVERFLOW(m, 1, 0)) |
---|
[2cfffe] | 2599 | Error("overflow in InvTrunc"); |
---|
| 2600 | |
---|
| 2601 | if (&x == &a) { |
---|
| 2602 | ZZ_pX la; |
---|
| 2603 | la = a; |
---|
| 2604 | if (m > NTL_ZZ_pX_NEWTON_CROSSOVER && deg(a) > 0) |
---|
| 2605 | NewtonInvTrunc(x, la, m); |
---|
| 2606 | else |
---|
| 2607 | PlainInvTrunc(x, la, m); |
---|
| 2608 | } |
---|
| 2609 | else { |
---|
| 2610 | if (m > NTL_ZZ_pX_NEWTON_CROSSOVER && deg(a) > 0) |
---|
| 2611 | NewtonInvTrunc(x, a, m); |
---|
| 2612 | else |
---|
| 2613 | PlainInvTrunc(x, a, m); |
---|
| 2614 | } |
---|
| 2615 | } |
---|
| 2616 | |
---|
| 2617 | |
---|
| 2618 | |
---|
| 2619 | void build(ZZ_pXModulus& x, const ZZ_pX& f) |
---|
| 2620 | { |
---|
| 2621 | x.f = f; |
---|
| 2622 | x.n = deg(f); |
---|
| 2623 | |
---|
| 2624 | x.tracevec.SetLength(0); |
---|
| 2625 | |
---|
| 2626 | if (x.n <= 0) |
---|
| 2627 | Error("build: deg(f) must be at least 1"); |
---|
| 2628 | |
---|
| 2629 | if (x.n <= NTL_ZZ_pX_FFT_CROSSOVER + 1) { |
---|
| 2630 | x.UseFFT = 0; |
---|
| 2631 | return; |
---|
| 2632 | } |
---|
| 2633 | |
---|
| 2634 | x.UseFFT = 1; |
---|
| 2635 | |
---|
| 2636 | x.k = NextPowerOfTwo(x.n); |
---|
| 2637 | x.l = NextPowerOfTwo(2*x.n - 3); |
---|
| 2638 | ToFFTRep(x.FRep, f, x.k); |
---|
| 2639 | |
---|
| 2640 | ZZ_pX P1(INIT_SIZE, x.n+1), P2(INIT_SIZE, x.n); |
---|
| 2641 | |
---|
| 2642 | CopyReverse(P1, f, 0, x.n); |
---|
| 2643 | InvTrunc(P2, P1, x.n-1); |
---|
| 2644 | |
---|
| 2645 | CopyReverse(P1, P2, 0, x.n-2); |
---|
| 2646 | ToFFTRep(x.HRep, P1, x.l); |
---|
| 2647 | } |
---|
| 2648 | |
---|
| 2649 | ZZ_pXModulus::ZZ_pXModulus(const ZZ_pX& ff) |
---|
| 2650 | { |
---|
| 2651 | build(*this, ff); |
---|
| 2652 | } |
---|
| 2653 | |
---|
| 2654 | ZZ_pXMultiplier::ZZ_pXMultiplier(const ZZ_pX& b, const ZZ_pXModulus& F) |
---|
| 2655 | { |
---|
| 2656 | build(*this, b, F); |
---|
| 2657 | } |
---|
| 2658 | |
---|
| 2659 | void build(ZZ_pXMultiplier& x, const ZZ_pX& b, |
---|
| 2660 | const ZZ_pXModulus& F) |
---|
| 2661 | { |
---|
| 2662 | long db; |
---|
| 2663 | long n = F.n; |
---|
| 2664 | |
---|
| 2665 | if (n < 0) Error("build ZZ_pXMultiplier: uninitialized modulus"); |
---|
| 2666 | |
---|
| 2667 | x.b = b; |
---|
| 2668 | db = deg(b); |
---|
| 2669 | |
---|
| 2670 | if (db >= n) Error("build ZZ_pXMultiplier: deg(b) >= deg(f)"); |
---|
| 2671 | |
---|
| 2672 | if (!F.UseFFT || db <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2673 | x.UseFFT = 0; |
---|
| 2674 | return; |
---|
| 2675 | } |
---|
| 2676 | |
---|
| 2677 | x.UseFFT = 1; |
---|
| 2678 | |
---|
| 2679 | FFTRep R1(INIT_SIZE, F.l); |
---|
| 2680 | ZZ_pX P1(INIT_SIZE, n); |
---|
| 2681 | |
---|
| 2682 | |
---|
| 2683 | ToFFTRep(R1, b, F.l); |
---|
| 2684 | reduce(x.B2, R1, F.k); |
---|
| 2685 | mul(R1, R1, F.HRep); |
---|
| 2686 | FromFFTRep(P1, R1, n-1, 2*n-3); |
---|
| 2687 | ToFFTRep(x.B1, P1, F.l); |
---|
| 2688 | } |
---|
| 2689 | |
---|
| 2690 | |
---|
| 2691 | void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXMultiplier& B, |
---|
| 2692 | const ZZ_pXModulus& F) |
---|
| 2693 | { |
---|
| 2694 | |
---|
| 2695 | long n = F.n; |
---|
| 2696 | long da; |
---|
| 2697 | |
---|
| 2698 | da = deg(a); |
---|
| 2699 | |
---|
| 2700 | if (da >= n) |
---|
| 2701 | Error(" bad args to MulMod(ZZ_pX,ZZ_pX,ZZ_pXMultiplier,ZZ_pXModulus)"); |
---|
| 2702 | |
---|
| 2703 | if (da < 0) { |
---|
| 2704 | clear(x); |
---|
| 2705 | return; |
---|
| 2706 | } |
---|
| 2707 | |
---|
| 2708 | if (!B.UseFFT || !F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER) { |
---|
| 2709 | ZZ_pX P1; |
---|
| 2710 | mul(P1, a, B.b); |
---|
| 2711 | rem(x, P1, F); |
---|
| 2712 | return; |
---|
| 2713 | } |
---|
| 2714 | |
---|
| 2715 | ZZ_pX P1(INIT_SIZE, n), P2(INIT_SIZE, n); |
---|
| 2716 | FFTRep R1(INIT_SIZE, F.l), R2(INIT_SIZE, F.l); |
---|
| 2717 | |
---|
| 2718 | ToFFTRep(R1, a, F.l); |
---|
| 2719 | mul(R2, R1, B.B1); |
---|
| 2720 | FromFFTRep(P1, R2, n-1, 2*n-3); |
---|
| 2721 | |
---|
| 2722 | reduce(R1, R1, F.k); |
---|
| 2723 | mul(R1, R1, B.B2); |
---|
| 2724 | ToFFTRep(R2, P1, F.k); |
---|
| 2725 | mul(R2, R2, F.FRep); |
---|
| 2726 | sub(R1, R1, R2); |
---|
| 2727 | |
---|
| 2728 | FromFFTRep(x, R1, 0, n-1); |
---|
| 2729 | } |
---|
| 2730 | |
---|
| 2731 | |
---|
| 2732 | void PowerXMod(ZZ_pX& hh, const ZZ& e, const ZZ_pXModulus& F) |
---|
| 2733 | { |
---|
| 2734 | if (F.n < 0) Error("PowerXMod: uninitialized modulus"); |
---|
| 2735 | |
---|
| 2736 | if (IsZero(e)) { |
---|
| 2737 | set(hh); |
---|
| 2738 | return; |
---|
| 2739 | } |
---|
| 2740 | |
---|
| 2741 | long n = NumBits(e); |
---|
| 2742 | long i; |
---|
| 2743 | |
---|
| 2744 | ZZ_pX h; |
---|
| 2745 | |
---|
| 2746 | h.SetMaxLength(F.n); |
---|
| 2747 | set(h); |
---|
| 2748 | |
---|
| 2749 | for (i = n - 1; i >= 0; i--) { |
---|
| 2750 | SqrMod(h, h, F); |
---|
| 2751 | if (bit(e, i)) |
---|
| 2752 | MulByXMod(h, h, F); |
---|
| 2753 | } |
---|
| 2754 | |
---|
| 2755 | if (e < 0) InvMod(h, h, F); |
---|
| 2756 | |
---|
| 2757 | hh = h; |
---|
| 2758 | } |
---|
| 2759 | |
---|
| 2760 | |
---|
| 2761 | void PowerXPlusAMod(ZZ_pX& hh, const ZZ_p& a, const ZZ& e, const ZZ_pXModulus& F) |
---|
| 2762 | { |
---|
| 2763 | if (F.n < 0) Error("PowerXPlusAMod: uninitialized modulus"); |
---|
| 2764 | |
---|
| 2765 | if (IsZero(e)) { |
---|
| 2766 | set(hh); |
---|
| 2767 | return; |
---|
| 2768 | } |
---|
| 2769 | |
---|
| 2770 | ZZ_pX t1(INIT_SIZE, F.n), t2(INIT_SIZE, F.n); |
---|
| 2771 | long n = NumBits(e); |
---|
| 2772 | long i; |
---|
| 2773 | |
---|
| 2774 | ZZ_pX h; |
---|
| 2775 | |
---|
| 2776 | h.SetMaxLength(F.n); |
---|
| 2777 | set(h); |
---|
| 2778 | |
---|
| 2779 | for (i = n - 1; i >= 0; i--) { |
---|
| 2780 | SqrMod(h, h, F); |
---|
| 2781 | if (bit(e, i)) { |
---|
| 2782 | MulByXMod(t1, h, F); |
---|
| 2783 | mul(t2, h, a); |
---|
| 2784 | add(h, t1, t2); |
---|
| 2785 | } |
---|
| 2786 | } |
---|
| 2787 | |
---|
| 2788 | if (e < 0) InvMod(h, h, F); |
---|
| 2789 | |
---|
| 2790 | hh = h; |
---|
| 2791 | } |
---|
| 2792 | |
---|
| 2793 | |
---|
| 2794 | void PowerMod(ZZ_pX& h, const ZZ_pX& g, const ZZ& e, const ZZ_pXModulus& F) |
---|
| 2795 | { |
---|
| 2796 | if (deg(g) >= F.n) |
---|
| 2797 | Error("PowerMod: bad args"); |
---|
| 2798 | |
---|
| 2799 | if (IsZero(e)) { |
---|
| 2800 | set(h); |
---|
| 2801 | return; |
---|
| 2802 | } |
---|
| 2803 | |
---|
| 2804 | ZZ_pXMultiplier G; |
---|
| 2805 | |
---|
| 2806 | ZZ_pX res; |
---|
| 2807 | |
---|
| 2808 | long n = NumBits(e); |
---|
| 2809 | long i; |
---|
| 2810 | |
---|
| 2811 | build(G, g, F); |
---|
| 2812 | |
---|
| 2813 | res.SetMaxLength(F.n); |
---|
| 2814 | set(res); |
---|
| 2815 | |
---|
| 2816 | for (i = n - 1; i >= 0; i--) { |
---|
| 2817 | SqrMod(res, res, F); |
---|
| 2818 | if (bit(e, i)) |
---|
| 2819 | MulMod(res, res, G, F); |
---|
| 2820 | } |
---|
| 2821 | |
---|
| 2822 | if (e < 0) InvMod(res, res, F); |
---|
| 2823 | |
---|
| 2824 | h = res; |
---|
| 2825 | } |
---|
| 2826 | |
---|
| 2827 | |
---|
| 2828 | void NewtonInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m) |
---|
| 2829 | { |
---|
| 2830 | x.SetMaxLength(m); |
---|
| 2831 | long i, t, k; |
---|
| 2832 | |
---|
| 2833 | long log2_newton = NextPowerOfTwo(NTL_ZZ_pX_NEWTON_CROSSOVER)-1; |
---|
| 2834 | PlainInvTrunc(x, a, 1L << log2_newton); |
---|
| 2835 | |
---|
| 2836 | t = NextPowerOfTwo(m); |
---|
| 2837 | |
---|
| 2838 | FFTRep R1(INIT_SIZE, t), R2(INIT_SIZE, t); |
---|
| 2839 | ZZ_pX P1(INIT_SIZE, m/2); |
---|
| 2840 | |
---|
| 2841 | long a_len = min(m, a.rep.length()); |
---|
| 2842 | |
---|
| 2843 | ZZ_pXModRep a_rep; |
---|
| 2844 | ToZZ_pXModRep(a_rep, a, 0, a_len-1); |
---|
| 2845 | |
---|
| 2846 | k = 1L << log2_newton; |
---|
| 2847 | t = log2_newton; |
---|
| 2848 | |
---|
| 2849 | while (k < m) { |
---|
| 2850 | long l = min(2*k, m); |
---|
| 2851 | |
---|
| 2852 | ToFFTRep(R1, x, t+1); |
---|
| 2853 | ToFFTRep(R2, a_rep, t+1, 0, l-1); |
---|
| 2854 | mul(R2, R2, R1); |
---|
| 2855 | FromFFTRep(P1, R2, k, l-1); |
---|
| 2856 | |
---|
| 2857 | ToFFTRep(R2, P1, t+1); |
---|
| 2858 | mul(R2, R2, R1); |
---|
| 2859 | FromFFTRep(P1, R2, 0, l-k-1); |
---|
| 2860 | |
---|
| 2861 | x.rep.SetLength(l); |
---|
| 2862 | long y_len = P1.rep.length(); |
---|
| 2863 | for (i = k; i < l; i++) { |
---|
| 2864 | if (i-k >= y_len) |
---|
| 2865 | clear(x.rep[i]); |
---|
| 2866 | else |
---|
| 2867 | negate(x.rep[i], P1.rep[i-k]); |
---|
| 2868 | } |
---|
| 2869 | x.normalize(); |
---|
| 2870 | |
---|
| 2871 | t++; |
---|
| 2872 | k = l; |
---|
| 2873 | } |
---|
| 2874 | } |
---|
| 2875 | |
---|
| 2876 | |
---|
| 2877 | |
---|
| 2878 | void FFTDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 2879 | { |
---|
| 2880 | long n = deg(b); |
---|
| 2881 | long m = deg(a); |
---|
| 2882 | long k, l; |
---|
| 2883 | |
---|
| 2884 | if (m < n) { |
---|
| 2885 | clear(q); |
---|
| 2886 | r = a; |
---|
| 2887 | return; |
---|
| 2888 | } |
---|
| 2889 | |
---|
| 2890 | if (m >= 3*n) { |
---|
| 2891 | ZZ_pXModulus B; |
---|
| 2892 | build(B, b); |
---|
| 2893 | DivRem(q, r, a, B); |
---|
| 2894 | return; |
---|
| 2895 | } |
---|
| 2896 | |
---|
| 2897 | ZZ_pX P1, P2, P3; |
---|
| 2898 | |
---|
| 2899 | CopyReverse(P3, b, 0, n); |
---|
| 2900 | InvTrunc(P2, P3, m-n+1); |
---|
| 2901 | CopyReverse(P1, P2, 0, m-n); |
---|
| 2902 | |
---|
| 2903 | k = NextPowerOfTwo(2*(m-n)+1); |
---|
| 2904 | long k1 = NextPowerOfTwo(n); |
---|
| 2905 | long mx = max(k1, k); |
---|
| 2906 | |
---|
| 2907 | FFTRep R1(INIT_SIZE, mx), R2(INIT_SIZE, mx); |
---|
| 2908 | |
---|
| 2909 | ToFFTRep(R1, P1, k); |
---|
| 2910 | ToFFTRep(R2, a, k, n, m); |
---|
| 2911 | mul(R1, R1, R2); |
---|
| 2912 | FromFFTRep(P3, R1, m-n, 2*(m-n)); |
---|
| 2913 | |
---|
| 2914 | l = 1L << k1; |
---|
| 2915 | |
---|
| 2916 | |
---|
| 2917 | ToFFTRep(R1, b, k1); |
---|
| 2918 | ToFFTRep(R2, P3, k1); |
---|
| 2919 | mul(R1, R1, R2); |
---|
| 2920 | FromFFTRep(P1, R1, 0, n-1); |
---|
| 2921 | CyclicReduce(P2, a, l); |
---|
| 2922 | trunc(r, P2, n); |
---|
| 2923 | sub(r, r, P1); |
---|
| 2924 | q = P3; |
---|
| 2925 | } |
---|
| 2926 | |
---|
| 2927 | |
---|
| 2928 | |
---|
| 2929 | |
---|
| 2930 | void FFTDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 2931 | { |
---|
| 2932 | |
---|
| 2933 | long n = deg(b); |
---|
| 2934 | long m = deg(a); |
---|
| 2935 | long k; |
---|
| 2936 | |
---|
| 2937 | if (m < n) { |
---|
| 2938 | clear(q); |
---|
| 2939 | return; |
---|
| 2940 | } |
---|
| 2941 | |
---|
| 2942 | if (m >= 3*n) { |
---|
| 2943 | ZZ_pXModulus B; |
---|
| 2944 | build(B, b); |
---|
| 2945 | div(q, a, B); |
---|
| 2946 | return; |
---|
| 2947 | } |
---|
| 2948 | |
---|
| 2949 | ZZ_pX P1, P2, P3; |
---|
| 2950 | |
---|
| 2951 | CopyReverse(P3, b, 0, n); |
---|
| 2952 | InvTrunc(P2, P3, m-n+1); |
---|
| 2953 | CopyReverse(P1, P2, 0, m-n); |
---|
| 2954 | |
---|
| 2955 | k = NextPowerOfTwo(2*(m-n)+1); |
---|
| 2956 | |
---|
| 2957 | FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k); |
---|
| 2958 | |
---|
| 2959 | ToFFTRep(R1, P1, k); |
---|
| 2960 | ToFFTRep(R2, a, k, n, m); |
---|
| 2961 | mul(R1, R1, R2); |
---|
| 2962 | FromFFTRep(q, R1, m-n, 2*(m-n)); |
---|
| 2963 | } |
---|
| 2964 | |
---|
| 2965 | |
---|
| 2966 | |
---|
| 2967 | void FFTRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 2968 | { |
---|
| 2969 | long n = deg(b); |
---|
| 2970 | long m = deg(a); |
---|
| 2971 | long k, l; |
---|
| 2972 | |
---|
| 2973 | if (m < n) { |
---|
| 2974 | r = a; |
---|
| 2975 | return; |
---|
| 2976 | } |
---|
| 2977 | |
---|
| 2978 | if (m >= 3*n) { |
---|
| 2979 | ZZ_pXModulus B; |
---|
| 2980 | build(B, b); |
---|
| 2981 | rem(r, a, B); |
---|
| 2982 | return; |
---|
| 2983 | } |
---|
| 2984 | |
---|
| 2985 | ZZ_pX P1, P2, P3; |
---|
| 2986 | |
---|
| 2987 | CopyReverse(P3, b, 0, n); |
---|
| 2988 | InvTrunc(P2, P3, m-n+1); |
---|
| 2989 | CopyReverse(P1, P2, 0, m-n); |
---|
| 2990 | |
---|
| 2991 | k = NextPowerOfTwo(2*(m-n)+1); |
---|
| 2992 | long k1 = NextPowerOfTwo(n); |
---|
| 2993 | long mx = max(k, k1); |
---|
| 2994 | |
---|
| 2995 | FFTRep R1(INIT_SIZE, mx), R2(INIT_SIZE, mx); |
---|
| 2996 | |
---|
| 2997 | ToFFTRep(R1, P1, k); |
---|
| 2998 | ToFFTRep(R2, a, k, n, m); |
---|
| 2999 | mul(R1, R1, R2); |
---|
| 3000 | FromFFTRep(P3, R1, m-n, 2*(m-n)); |
---|
| 3001 | |
---|
| 3002 | l = 1L << k1; |
---|
| 3003 | |
---|
| 3004 | ToFFTRep(R1, b, k1); |
---|
| 3005 | ToFFTRep(R2, P3, k1); |
---|
| 3006 | mul(R1, R1, R2); |
---|
| 3007 | FromFFTRep(P3, R1, 0, n-1); |
---|
| 3008 | CyclicReduce(P2, a, l); |
---|
| 3009 | trunc(r, P2, n); |
---|
| 3010 | sub(r, r, P3); |
---|
| 3011 | } |
---|
| 3012 | |
---|
| 3013 | |
---|
| 3014 | void DivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 3015 | { |
---|
| 3016 | if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER) |
---|
| 3017 | FFTDivRem(q, r, a, b); |
---|
| 3018 | else |
---|
| 3019 | PlainDivRem(q, r, a, b); |
---|
| 3020 | } |
---|
| 3021 | |
---|
| 3022 | void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 3023 | { |
---|
| 3024 | if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER) |
---|
| 3025 | FFTDiv(q, a, b); |
---|
| 3026 | else |
---|
| 3027 | PlainDiv(q, a, b); |
---|
| 3028 | } |
---|
| 3029 | |
---|
| 3030 | void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_p& b) |
---|
| 3031 | { |
---|
| 3032 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 3033 | |
---|
| 3034 | inv(T, b); |
---|
| 3035 | mul(q, a, T); |
---|
| 3036 | } |
---|
| 3037 | |
---|
| 3038 | void div(ZZ_pX& q, const ZZ_pX& a, long b) |
---|
| 3039 | { |
---|
| 3040 | ZZ_pTemp TT; ZZ_p& T = TT.val(); |
---|
| 3041 | |
---|
| 3042 | T = b; |
---|
| 3043 | inv(T, T); |
---|
| 3044 | mul(q, a, T); |
---|
| 3045 | } |
---|
| 3046 | |
---|
| 3047 | |
---|
| 3048 | |
---|
| 3049 | void rem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b) |
---|
| 3050 | { |
---|
| 3051 | if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER) |
---|
| 3052 | FFTRem(r, a, b); |
---|
| 3053 | else |
---|
| 3054 | PlainRem(r, a, b); |
---|
| 3055 | } |
---|
| 3056 | |
---|
| 3057 | |
---|
| 3058 | long operator==(const ZZ_pX& a, long b) |
---|
| 3059 | { |
---|
| 3060 | if (b == 0) |
---|
| 3061 | return IsZero(a); |
---|
| 3062 | |
---|
| 3063 | if (b == 1) |
---|
| 3064 | return IsOne(a); |
---|
| 3065 | |
---|
| 3066 | long da = deg(a); |
---|
| 3067 | |
---|
| 3068 | if (da > 0) |
---|
| 3069 | return 0; |
---|
| 3070 | |
---|
| 3071 | ZZ_pTemp TT; ZZ_p& bb = TT.val(); |
---|
| 3072 | bb = b; |
---|
| 3073 | |
---|
| 3074 | if (da < 0) |
---|
| 3075 | return IsZero(bb); |
---|
| 3076 | |
---|
| 3077 | return a.rep[0] == bb; |
---|
| 3078 | } |
---|
| 3079 | |
---|
| 3080 | long operator==(const ZZ_pX& a, const ZZ_p& b) |
---|
| 3081 | { |
---|
| 3082 | if (IsZero(b)) |
---|
| 3083 | return IsZero(a); |
---|
| 3084 | |
---|
| 3085 | long da = deg(a); |
---|
| 3086 | |
---|
| 3087 | if (da != 0) |
---|
| 3088 | return 0; |
---|
| 3089 | |
---|
| 3090 | return a.rep[0] == b; |
---|
| 3091 | } |
---|
| 3092 | |
---|
| 3093 | void power(ZZ_pX& x, const ZZ_pX& a, long e) |
---|
| 3094 | { |
---|
| 3095 | if (e < 0) { |
---|
| 3096 | Error("power: negative exponent"); |
---|
| 3097 | } |
---|
| 3098 | |
---|
| 3099 | if (e == 0) { |
---|
| 3100 | x = 1; |
---|
| 3101 | return; |
---|
| 3102 | } |
---|
| 3103 | |
---|
| 3104 | if (a == 0 || a == 1) { |
---|
| 3105 | x = a; |
---|
| 3106 | return; |
---|
| 3107 | } |
---|
| 3108 | |
---|
| 3109 | long da = deg(a); |
---|
| 3110 | |
---|
| 3111 | if (da == 0) { |
---|
| 3112 | x = power(ConstTerm(a), e); |
---|
| 3113 | return; |
---|
| 3114 | } |
---|
| 3115 | |
---|
| 3116 | if (da > (NTL_MAX_LONG-1)/e) |
---|
| 3117 | Error("overflow in power"); |
---|
| 3118 | |
---|
| 3119 | ZZ_pX res; |
---|
| 3120 | res.SetMaxLength(da*e + 1); |
---|
| 3121 | res = 1; |
---|
| 3122 | |
---|
| 3123 | long k = NumBits(e); |
---|
| 3124 | long i; |
---|
| 3125 | |
---|
| 3126 | for (i = k - 1; i >= 0; i--) { |
---|
| 3127 | sqr(res, res); |
---|
| 3128 | if (bit(e, i)) |
---|
| 3129 | mul(res, res, a); |
---|
| 3130 | } |
---|
| 3131 | |
---|
| 3132 | x = res; |
---|
| 3133 | } |
---|
| 3134 | |
---|
| 3135 | void reverse(ZZ_pX& x, const ZZ_pX& a, long hi) |
---|
| 3136 | { |
---|
[09da99] | 3137 | if (hi < 0) { clear(x); return; } |
---|
| 3138 | if (NTL_OVERFLOW(hi, 1, 0)) |
---|
[2cfffe] | 3139 | Error("overflow in reverse"); |
---|
| 3140 | |
---|
| 3141 | if (&x == &a) { |
---|
| 3142 | ZZ_pX tmp; |
---|
| 3143 | CopyReverse(tmp, a, 0, hi); |
---|
| 3144 | x = tmp; |
---|
| 3145 | } |
---|
| 3146 | else |
---|
| 3147 | CopyReverse(x, a, 0, hi); |
---|
| 3148 | } |
---|
| 3149 | |
---|
| 3150 | NTL_END_IMPL |
---|