1 | |
---|
2 | |
---|
3 | #include <NTL/ZZ.h> |
---|
4 | #include <NTL/vec_ZZ.h> |
---|
5 | |
---|
6 | |
---|
7 | #include <NTL/new.h> |
---|
8 | |
---|
9 | |
---|
10 | |
---|
11 | NTL_START_IMPL |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | |
---|
16 | const ZZ& ZZ::zero() |
---|
17 | { |
---|
18 | static ZZ z; |
---|
19 | return z; |
---|
20 | } |
---|
21 | |
---|
22 | |
---|
23 | const ZZ& ZZ_expo(long e) |
---|
24 | { |
---|
25 | static ZZ expo_helper; |
---|
26 | conv(expo_helper, e); |
---|
27 | return expo_helper; |
---|
28 | } |
---|
29 | |
---|
30 | |
---|
31 | |
---|
32 | |
---|
33 | void AddMod(ZZ& x, const ZZ& a, long b, const ZZ& n) |
---|
34 | { |
---|
35 | static ZZ B; |
---|
36 | conv(B, b); |
---|
37 | AddMod(x, a, B, n); |
---|
38 | } |
---|
39 | |
---|
40 | |
---|
41 | void SubMod(ZZ& x, const ZZ& a, long b, const ZZ& n) |
---|
42 | { |
---|
43 | static ZZ B; |
---|
44 | conv(B, b); |
---|
45 | SubMod(x, a, B, n); |
---|
46 | } |
---|
47 | |
---|
48 | void SubMod(ZZ& x, long a, const ZZ& b, const ZZ& n) |
---|
49 | { |
---|
50 | static ZZ A; |
---|
51 | conv(A, a); |
---|
52 | SubMod(x, A, b, n); |
---|
53 | } |
---|
54 | |
---|
55 | |
---|
56 | |
---|
57 | // ****** input and output |
---|
58 | |
---|
59 | static long iodigits = 0; |
---|
60 | static long ioradix = 0; |
---|
61 | |
---|
62 | // iodigits is the greatest integer such that 10^{iodigits} < NTL_WSP_BOUND |
---|
63 | // ioradix = 10^{iodigits} |
---|
64 | |
---|
65 | static void InitZZIO() |
---|
66 | { |
---|
67 | long x; |
---|
68 | |
---|
69 | x = (NTL_WSP_BOUND-1)/10; |
---|
70 | iodigits = 0; |
---|
71 | ioradix = 1; |
---|
72 | |
---|
73 | while (x) { |
---|
74 | x = x / 10; |
---|
75 | iodigits++; |
---|
76 | ioradix = ioradix * 10; |
---|
77 | } |
---|
78 | |
---|
79 | if (iodigits <= 0) Error("problem with I/O"); |
---|
80 | } |
---|
81 | |
---|
82 | // The class _ZZ_local_stack should be defined in an empty namespace, |
---|
83 | // but since I don't want to rely on namespaces, we just give it a funny |
---|
84 | // name to avoid accidental name clashes. |
---|
85 | |
---|
86 | struct _ZZ_local_stack { |
---|
87 | long top; |
---|
88 | long alloc; |
---|
89 | long *elts; |
---|
90 | |
---|
91 | _ZZ_local_stack() { top = -1; alloc = 0; elts = 0; } |
---|
92 | ~_ZZ_local_stack() { } |
---|
93 | |
---|
94 | long pop() { return elts[top--]; } |
---|
95 | long empty() { return (top == -1); } |
---|
96 | void push(long x); |
---|
97 | }; |
---|
98 | |
---|
99 | void _ZZ_local_stack::push(long x) |
---|
100 | { |
---|
101 | if (alloc == 0) { |
---|
102 | alloc = 100; |
---|
103 | elts = (long *) NTL_MALLOC(alloc, sizeof(long), 0); |
---|
104 | } |
---|
105 | |
---|
106 | top++; |
---|
107 | |
---|
108 | if (top + 1 > alloc) { |
---|
109 | alloc = 2*alloc; |
---|
110 | elts = (long *) NTL_REALLOC(elts, alloc, sizeof(long), 0); |
---|
111 | } |
---|
112 | |
---|
113 | if (!elts) { |
---|
114 | Error("out of space in ZZ output"); |
---|
115 | } |
---|
116 | |
---|
117 | elts[top] = x; |
---|
118 | } |
---|
119 | |
---|
120 | |
---|
121 | long GCD(long a, long b) |
---|
122 | { |
---|
123 | long u, v, t, x; |
---|
124 | |
---|
125 | if (a < 0) { |
---|
126 | if (a < -NTL_MAX_LONG) Error("GCD: integer overflow"); |
---|
127 | a = -a; |
---|
128 | } |
---|
129 | |
---|
130 | if (b < 0) { |
---|
131 | if (b < -NTL_MAX_LONG) Error("GCD: integer overflow"); |
---|
132 | b = -b; |
---|
133 | } |
---|
134 | |
---|
135 | |
---|
136 | if (b==0) |
---|
137 | x = a; |
---|
138 | else { |
---|
139 | u = a; |
---|
140 | v = b; |
---|
141 | do { |
---|
142 | t = u % v; |
---|
143 | u = v; |
---|
144 | v = t; |
---|
145 | } while (v != 0); |
---|
146 | |
---|
147 | x = u; |
---|
148 | } |
---|
149 | |
---|
150 | return x; |
---|
151 | } |
---|
152 | |
---|
153 | |
---|
154 | |
---|
155 | void XGCD(long& d, long& s, long& t, long a, long b) |
---|
156 | { |
---|
157 | long u, v, u0, v0, u1, v1, u2, v2, q, r; |
---|
158 | |
---|
159 | long aneg = 0, bneg = 0; |
---|
160 | |
---|
161 | if (a < 0) { |
---|
162 | if (a < -NTL_MAX_LONG) Error("XGCD: integer overflow"); |
---|
163 | a = -a; |
---|
164 | aneg = 1; |
---|
165 | } |
---|
166 | |
---|
167 | if (b < 0) { |
---|
168 | if (b < -NTL_MAX_LONG) Error("XGCD: integer overflow"); |
---|
169 | b = -b; |
---|
170 | bneg = 1; |
---|
171 | } |
---|
172 | |
---|
173 | u1=1; v1=0; |
---|
174 | u2=0; v2=1; |
---|
175 | u = a; v = b; |
---|
176 | |
---|
177 | while (v != 0) { |
---|
178 | q = u / v; |
---|
179 | r = u % v; |
---|
180 | u = v; |
---|
181 | v = r; |
---|
182 | u0 = u2; |
---|
183 | v0 = v2; |
---|
184 | u2 = u1 - q*u2; |
---|
185 | v2 = v1- q*v2; |
---|
186 | u1 = u0; |
---|
187 | v1 = v0; |
---|
188 | } |
---|
189 | |
---|
190 | if (aneg) |
---|
191 | u1 = -u1; |
---|
192 | |
---|
193 | if (bneg) |
---|
194 | v1 = -v1; |
---|
195 | |
---|
196 | d = u; |
---|
197 | s = u1; |
---|
198 | t = v1; |
---|
199 | } |
---|
200 | |
---|
201 | |
---|
202 | long InvMod(long a, long n) |
---|
203 | { |
---|
204 | long d, s, t; |
---|
205 | |
---|
206 | XGCD(d, s, t, a, n); |
---|
207 | if (d != 1) Error("InvMod: inverse undefined"); |
---|
208 | if (s < 0) |
---|
209 | return s + n; |
---|
210 | else |
---|
211 | return s; |
---|
212 | } |
---|
213 | |
---|
214 | |
---|
215 | long PowerMod(long a, long ee, long n) |
---|
216 | { |
---|
217 | long x, y; |
---|
218 | |
---|
219 | unsigned long e; |
---|
220 | |
---|
221 | if (ee < 0) |
---|
222 | e = - ((unsigned long) ee); |
---|
223 | else |
---|
224 | e = ee; |
---|
225 | |
---|
226 | x = 1; |
---|
227 | y = a; |
---|
228 | while (e) { |
---|
229 | if (e & 1) x = MulMod(x, y, n); |
---|
230 | y = MulMod(y, y, n); |
---|
231 | e = e >> 1; |
---|
232 | } |
---|
233 | |
---|
234 | if (ee < 0) x = InvMod(x, n); |
---|
235 | |
---|
236 | return x; |
---|
237 | } |
---|
238 | |
---|
239 | long ProbPrime(long n, long NumTests) |
---|
240 | { |
---|
241 | long m, x, y, z; |
---|
242 | long i, j, k; |
---|
243 | |
---|
244 | if (n <= 1) return 0; |
---|
245 | |
---|
246 | |
---|
247 | if (n == 2) return 1; |
---|
248 | if (n % 2 == 0) return 0; |
---|
249 | |
---|
250 | if (n == 3) return 1; |
---|
251 | if (n % 3 == 0) return 0; |
---|
252 | |
---|
253 | if (n == 5) return 1; |
---|
254 | if (n % 5 == 0) return 0; |
---|
255 | |
---|
256 | if (n == 7) return 1; |
---|
257 | if (n % 7 == 0) return 0; |
---|
258 | |
---|
259 | if (n >= NTL_SP_BOUND) { |
---|
260 | return ProbPrime(to_ZZ(n), NumTests); |
---|
261 | } |
---|
262 | |
---|
263 | m = n - 1; |
---|
264 | k = 0; |
---|
265 | while((m & 1) == 0) { |
---|
266 | m = m >> 1; |
---|
267 | k++; |
---|
268 | } |
---|
269 | |
---|
270 | // n - 1 == 2^k * m, m odd |
---|
271 | |
---|
272 | for (i = 0; i < NumTests; i++) { |
---|
273 | do { |
---|
274 | x = RandomBnd(n); |
---|
275 | } while (x == 0); |
---|
276 | // x == 0 is not a useful candidtae for a witness! |
---|
277 | |
---|
278 | |
---|
279 | if (x == 0) continue; |
---|
280 | z = PowerMod(x, m, n); |
---|
281 | if (z == 1) continue; |
---|
282 | |
---|
283 | j = 0; |
---|
284 | do { |
---|
285 | y = z; |
---|
286 | z = MulMod(y, y, n); |
---|
287 | j++; |
---|
288 | } while (j != k && z != 1); |
---|
289 | |
---|
290 | if (z != 1 || y != n-1) return 0; |
---|
291 | } |
---|
292 | |
---|
293 | return 1; |
---|
294 | } |
---|
295 | |
---|
296 | |
---|
297 | long MillerWitness(const ZZ& n, const ZZ& x) |
---|
298 | { |
---|
299 | ZZ m, y, z; |
---|
300 | long j, k; |
---|
301 | |
---|
302 | if (x == 0) return 0; |
---|
303 | |
---|
304 | add(m, n, -1); |
---|
305 | k = MakeOdd(m); |
---|
306 | // n - 1 == 2^k * m, m odd |
---|
307 | |
---|
308 | PowerMod(z, x, m, n); |
---|
309 | if (z == 1) return 0; |
---|
310 | |
---|
311 | j = 0; |
---|
312 | do { |
---|
313 | y = z; |
---|
314 | SqrMod(z, y, n); |
---|
315 | j++; |
---|
316 | } while (j != k && z != 1); |
---|
317 | |
---|
318 | if (z != 1) return 1; |
---|
319 | add(y, y, 1); |
---|
320 | if (y != n) return 1; |
---|
321 | return 0; |
---|
322 | } |
---|
323 | |
---|
324 | |
---|
325 | // ComputePrimeBound computes a reasonable bound for trial |
---|
326 | // division in the Miller-Rabin test. |
---|
327 | // It is computed a bit on the "low" side, since being a bit |
---|
328 | // low doesn't hurt much, but being too high can hurt a lot. |
---|
329 | |
---|
330 | static |
---|
331 | long ComputePrimeBound(long bn) |
---|
332 | { |
---|
333 | long wn = (bn+NTL_ZZ_NBITS-1)/NTL_ZZ_NBITS; |
---|
334 | |
---|
335 | long fn; |
---|
336 | |
---|
337 | if (wn <= 36) |
---|
338 | fn = wn/4 + 1; |
---|
339 | else |
---|
340 | fn = long(1.67*sqrt(double(wn))); |
---|
341 | |
---|
342 | long prime_bnd; |
---|
343 | |
---|
344 | if (NumBits(bn) + NumBits(fn) > NTL_SP_NBITS) |
---|
345 | prime_bnd = NTL_SP_BOUND; |
---|
346 | else |
---|
347 | prime_bnd = bn*fn; |
---|
348 | |
---|
349 | return prime_bnd; |
---|
350 | } |
---|
351 | |
---|
352 | |
---|
353 | long ProbPrime(const ZZ& n, long NumTrials) |
---|
354 | { |
---|
355 | if (n <= 1) return 0; |
---|
356 | |
---|
357 | if (n.SinglePrecision()) { |
---|
358 | return ProbPrime(to_long(n), NumTrials); |
---|
359 | } |
---|
360 | |
---|
361 | |
---|
362 | long prime_bnd = ComputePrimeBound(NumBits(n)); |
---|
363 | |
---|
364 | |
---|
365 | PrimeSeq s; |
---|
366 | long p; |
---|
367 | |
---|
368 | p = s.next(); |
---|
369 | while (p && p < prime_bnd) { |
---|
370 | if (rem(n, p) == 0) |
---|
371 | return 0; |
---|
372 | |
---|
373 | p = s.next(); |
---|
374 | } |
---|
375 | |
---|
376 | ZZ W; |
---|
377 | W = 2; |
---|
378 | |
---|
379 | // first try W == 2....the exponentiation |
---|
380 | // algorithm runs slightly faster in this case |
---|
381 | |
---|
382 | if (MillerWitness(n, W)) |
---|
383 | return 0; |
---|
384 | |
---|
385 | |
---|
386 | long i; |
---|
387 | |
---|
388 | for (i = 0; i < NumTrials; i++) { |
---|
389 | do { |
---|
390 | RandomBnd(W, n); |
---|
391 | } while (W == 0); |
---|
392 | // W == 0 is not a useful candidate for a witness! |
---|
393 | |
---|
394 | if (MillerWitness(n, W)) |
---|
395 | return 0; |
---|
396 | } |
---|
397 | |
---|
398 | return 1; |
---|
399 | } |
---|
400 | |
---|
401 | |
---|
402 | void RandomPrime(ZZ& n, long l, long NumTrials) |
---|
403 | { |
---|
404 | if (l <= 1) |
---|
405 | Error("RandomPrime: l out of range"); |
---|
406 | |
---|
407 | if (l == 2) { |
---|
408 | if (RandomBnd(2)) |
---|
409 | n = 3; |
---|
410 | else |
---|
411 | n = 2; |
---|
412 | |
---|
413 | return; |
---|
414 | } |
---|
415 | |
---|
416 | do { |
---|
417 | RandomLen(n, l); |
---|
418 | if (!IsOdd(n)) add(n, n, 1); |
---|
419 | } while (!ProbPrime(n, NumTrials)); |
---|
420 | } |
---|
421 | |
---|
422 | void NextPrime(ZZ& n, const ZZ& m, long NumTrials) |
---|
423 | { |
---|
424 | ZZ x; |
---|
425 | |
---|
426 | if (m <= 2) { |
---|
427 | n = 2; |
---|
428 | return; |
---|
429 | } |
---|
430 | |
---|
431 | x = m; |
---|
432 | |
---|
433 | while (!ProbPrime(x, NumTrials)) |
---|
434 | add(x, x, 1); |
---|
435 | |
---|
436 | n = x; |
---|
437 | } |
---|
438 | |
---|
439 | long NextPrime(long m, long NumTrials) |
---|
440 | { |
---|
441 | long x; |
---|
442 | |
---|
443 | if (m <= 2) |
---|
444 | return 2; |
---|
445 | |
---|
446 | x = m; |
---|
447 | |
---|
448 | while (x < NTL_SP_BOUND && !ProbPrime(x, NumTrials)) |
---|
449 | x++; |
---|
450 | |
---|
451 | if (x >= NTL_SP_BOUND) |
---|
452 | Error("NextPrime: no more primes"); |
---|
453 | |
---|
454 | return x; |
---|
455 | } |
---|
456 | |
---|
457 | |
---|
458 | |
---|
459 | long NextPowerOfTwo(long m) |
---|
460 | { |
---|
461 | long k; |
---|
462 | unsigned long n, um; |
---|
463 | |
---|
464 | if (m < 0) return 0; |
---|
465 | |
---|
466 | um = m; |
---|
467 | n = 1; |
---|
468 | k = 0; |
---|
469 | |
---|
470 | while (n < um) { |
---|
471 | n = n << 1; |
---|
472 | k++; |
---|
473 | } |
---|
474 | |
---|
475 | if (k >= NTL_BITS_PER_LONG-1) |
---|
476 | Error("NextPowerOfTwo: overflow"); |
---|
477 | |
---|
478 | return k; |
---|
479 | } |
---|
480 | |
---|
481 | |
---|
482 | |
---|
483 | long NumBits(long a) |
---|
484 | { |
---|
485 | unsigned long aa; |
---|
486 | if (a < 0) |
---|
487 | aa = - ((unsigned long) a); |
---|
488 | else |
---|
489 | aa = a; |
---|
490 | |
---|
491 | long k = 0; |
---|
492 | while (aa) { |
---|
493 | k++; |
---|
494 | aa = aa >> 1; |
---|
495 | } |
---|
496 | |
---|
497 | return k; |
---|
498 | } |
---|
499 | |
---|
500 | |
---|
501 | long bit(long a, long k) |
---|
502 | { |
---|
503 | unsigned long aa; |
---|
504 | if (a < 0) |
---|
505 | aa = - ((unsigned long) a); |
---|
506 | else |
---|
507 | aa = a; |
---|
508 | |
---|
509 | if (k < 0 || k >= NTL_BITS_PER_LONG) |
---|
510 | return 0; |
---|
511 | else |
---|
512 | return long((aa >> k) & 1); |
---|
513 | } |
---|
514 | |
---|
515 | |
---|
516 | |
---|
517 | long divide(ZZ& q, const ZZ& a, const ZZ& b) |
---|
518 | { |
---|
519 | static ZZ qq, r; |
---|
520 | |
---|
521 | if (IsZero(b)) { |
---|
522 | if (IsZero(a)) { |
---|
523 | clear(q); |
---|
524 | return 1; |
---|
525 | } |
---|
526 | else |
---|
527 | return 0; |
---|
528 | } |
---|
529 | |
---|
530 | |
---|
531 | if (IsOne(b)) { |
---|
532 | q = a; |
---|
533 | return 1; |
---|
534 | } |
---|
535 | |
---|
536 | DivRem(qq, r, a, b); |
---|
537 | if (!IsZero(r)) return 0; |
---|
538 | q = qq; |
---|
539 | return 1; |
---|
540 | } |
---|
541 | |
---|
542 | long divide(const ZZ& a, const ZZ& b) |
---|
543 | { |
---|
544 | static ZZ r; |
---|
545 | |
---|
546 | if (IsZero(b)) return IsZero(a); |
---|
547 | if (IsOne(b)) return 1; |
---|
548 | |
---|
549 | rem(r, a, b); |
---|
550 | return IsZero(r); |
---|
551 | } |
---|
552 | |
---|
553 | long divide(ZZ& q, const ZZ& a, long b) |
---|
554 | { |
---|
555 | static ZZ qq; |
---|
556 | |
---|
557 | if (!b) { |
---|
558 | if (IsZero(a)) { |
---|
559 | clear(q); |
---|
560 | return 1; |
---|
561 | } |
---|
562 | else |
---|
563 | return 0; |
---|
564 | } |
---|
565 | |
---|
566 | if (b == 1) { |
---|
567 | q = a; |
---|
568 | return 1; |
---|
569 | } |
---|
570 | |
---|
571 | long r = DivRem(qq, a, b); |
---|
572 | if (r) return 0; |
---|
573 | q = qq; |
---|
574 | return 1; |
---|
575 | } |
---|
576 | |
---|
577 | long divide(const ZZ& a, long b) |
---|
578 | { |
---|
579 | if (!b) return IsZero(a); |
---|
580 | if (b == 1) { |
---|
581 | return 1; |
---|
582 | } |
---|
583 | |
---|
584 | long r = rem(a, b); |
---|
585 | return (r == 0); |
---|
586 | } |
---|
587 | |
---|
588 | |
---|
589 | |
---|
590 | long RandomPrime_long(long l, long NumTrials) |
---|
591 | { |
---|
592 | if (l <= 1 || l >= NTL_BITS_PER_LONG) |
---|
593 | Error("RandomPrime: length out of range"); |
---|
594 | |
---|
595 | long n; |
---|
596 | do { |
---|
597 | n = RandomLen_long(l); |
---|
598 | } while (!ProbPrime(n, NumTrials)); |
---|
599 | |
---|
600 | return n; |
---|
601 | } |
---|
602 | |
---|
603 | |
---|
604 | PrimeSeq::PrimeSeq() |
---|
605 | { |
---|
606 | movesieve = 0; |
---|
607 | movesieve_mem = 0; |
---|
608 | pshift = -1; |
---|
609 | pindex = -1; |
---|
610 | exhausted = 0; |
---|
611 | } |
---|
612 | |
---|
613 | PrimeSeq::~PrimeSeq() |
---|
614 | { |
---|
615 | if (movesieve_mem) |
---|
616 | free(movesieve_mem); |
---|
617 | } |
---|
618 | |
---|
619 | long PrimeSeq::next() |
---|
620 | { |
---|
621 | if (exhausted) { |
---|
622 | return 0; |
---|
623 | } |
---|
624 | |
---|
625 | if (pshift < 0) { |
---|
626 | shift(0); |
---|
627 | return 2; |
---|
628 | } |
---|
629 | |
---|
630 | for (;;) { |
---|
631 | char *p = movesieve; |
---|
632 | long i = pindex; |
---|
633 | |
---|
634 | while ((++i) < NTL_PRIME_BND) { |
---|
635 | if (p[i]) { |
---|
636 | pindex = i; |
---|
637 | return pshift + 2 * i + 3; |
---|
638 | } |
---|
639 | } |
---|
640 | |
---|
641 | long newshift = pshift + 2*NTL_PRIME_BND; |
---|
642 | |
---|
643 | if (newshift > 2 * NTL_PRIME_BND * (2 * NTL_PRIME_BND + 1)) { |
---|
644 | /* end of the road */ |
---|
645 | exhausted = 1; |
---|
646 | return 0; |
---|
647 | } |
---|
648 | |
---|
649 | shift(newshift); |
---|
650 | } |
---|
651 | } |
---|
652 | |
---|
653 | static char *lowsieve = 0; |
---|
654 | |
---|
655 | void PrimeSeq::shift(long newshift) |
---|
656 | { |
---|
657 | long i; |
---|
658 | long j; |
---|
659 | long jstep; |
---|
660 | long jstart; |
---|
661 | long ibound; |
---|
662 | char *p; |
---|
663 | |
---|
664 | if (!lowsieve) |
---|
665 | start(); |
---|
666 | |
---|
667 | pindex = -1; |
---|
668 | exhausted = 0; |
---|
669 | |
---|
670 | if (newshift < 0) { |
---|
671 | pshift = -1; |
---|
672 | return; |
---|
673 | } |
---|
674 | |
---|
675 | if (newshift == pshift) return; |
---|
676 | |
---|
677 | pshift = newshift; |
---|
678 | |
---|
679 | if (pshift == 0) { |
---|
680 | movesieve = lowsieve; |
---|
681 | } |
---|
682 | else { |
---|
683 | if (!movesieve_mem) { |
---|
684 | movesieve_mem = (char *) NTL_MALLOC(NTL_PRIME_BND, 1, 0); |
---|
685 | if (!movesieve_mem) |
---|
686 | Error("out of memory in PrimeSeq"); |
---|
687 | } |
---|
688 | |
---|
689 | p = movesieve = movesieve_mem; |
---|
690 | for (i = 0; i < NTL_PRIME_BND; i++) |
---|
691 | p[i] = 1; |
---|
692 | |
---|
693 | jstep = 3; |
---|
694 | ibound = pshift + 2 * NTL_PRIME_BND + 1; |
---|
695 | for (i = 0; jstep * jstep <= ibound; i++) { |
---|
696 | if (lowsieve[i]) { |
---|
697 | if (!((jstart = (pshift + 2) / jstep + 1) & 1)) |
---|
698 | jstart++; |
---|
699 | if (jstart <= jstep) |
---|
700 | jstart = jstep; |
---|
701 | jstart = (jstart * jstep - pshift - 3) / 2; |
---|
702 | for (j = jstart; j < NTL_PRIME_BND; j += jstep) |
---|
703 | p[j] = 0; |
---|
704 | } |
---|
705 | jstep += 2; |
---|
706 | } |
---|
707 | } |
---|
708 | } |
---|
709 | |
---|
710 | |
---|
711 | void PrimeSeq::start() |
---|
712 | { |
---|
713 | long i; |
---|
714 | long j; |
---|
715 | long jstep; |
---|
716 | long jstart; |
---|
717 | long ibnd; |
---|
718 | char *p; |
---|
719 | |
---|
720 | p = lowsieve = (char *) NTL_MALLOC(NTL_PRIME_BND, 1, 0); |
---|
721 | if (!p) |
---|
722 | Error("out of memory in PrimeSeq"); |
---|
723 | |
---|
724 | for (i = 0; i < NTL_PRIME_BND; i++) |
---|
725 | p[i] = 1; |
---|
726 | |
---|
727 | jstep = 1; |
---|
728 | jstart = -1; |
---|
729 | ibnd = (SqrRoot(2 * NTL_PRIME_BND + 1) - 3) / 2; |
---|
730 | for (i = 0; i <= ibnd; i++) { |
---|
731 | jstart += 2 * ((jstep += 2) - 1); |
---|
732 | if (p[i]) |
---|
733 | for (j = jstart; j < NTL_PRIME_BND; j += jstep) |
---|
734 | p[j] = 0; |
---|
735 | } |
---|
736 | } |
---|
737 | |
---|
738 | void PrimeSeq::reset(long b) |
---|
739 | { |
---|
740 | if (b > (2*NTL_PRIME_BND+1)*(2*NTL_PRIME_BND+1)) { |
---|
741 | exhausted = 1; |
---|
742 | return; |
---|
743 | } |
---|
744 | |
---|
745 | if (b <= 2) { |
---|
746 | shift(-1); |
---|
747 | return; |
---|
748 | } |
---|
749 | |
---|
750 | if ((b & 1) == 0) b++; |
---|
751 | |
---|
752 | shift(((b-3) / (2*NTL_PRIME_BND))* (2*NTL_PRIME_BND)); |
---|
753 | pindex = (b - pshift - 3)/2 - 1; |
---|
754 | } |
---|
755 | |
---|
756 | long Jacobi(const ZZ& aa, const ZZ& nn) |
---|
757 | { |
---|
758 | ZZ a, n; |
---|
759 | long t, k; |
---|
760 | long d; |
---|
761 | |
---|
762 | a = aa; |
---|
763 | n = nn; |
---|
764 | t = 1; |
---|
765 | |
---|
766 | while (a != 0) { |
---|
767 | k = MakeOdd(a); |
---|
768 | d = trunc_long(n, 3); |
---|
769 | if ((k & 1) && (d == 3 || d == 5)) t = -t; |
---|
770 | |
---|
771 | if (trunc_long(a, 2) == 3 && (d & 3) == 3) t = -t; |
---|
772 | swap(a, n); |
---|
773 | rem(a, a, n); |
---|
774 | } |
---|
775 | |
---|
776 | if (n == 1) |
---|
777 | return t; |
---|
778 | else |
---|
779 | return 0; |
---|
780 | } |
---|
781 | |
---|
782 | |
---|
783 | void SqrRootMod(ZZ& x, const ZZ& aa, const ZZ& nn) |
---|
784 | { |
---|
785 | if (aa == 0 || aa == 1) { |
---|
786 | x = aa; |
---|
787 | return; |
---|
788 | } |
---|
789 | |
---|
790 | // at this point, we must have nn >= 5 |
---|
791 | |
---|
792 | if (trunc_long(nn, 2) == 3) { // special case, n = 3 (mod 4) |
---|
793 | ZZ n, a, e, z; |
---|
794 | |
---|
795 | n = nn; |
---|
796 | a = aa; |
---|
797 | |
---|
798 | add(e, n, 1); |
---|
799 | RightShift(e, e, 2); |
---|
800 | |
---|
801 | PowerMod(z, a, e, n); |
---|
802 | x = z; |
---|
803 | |
---|
804 | return; |
---|
805 | } |
---|
806 | |
---|
807 | ZZ n, m; |
---|
808 | int h, nlen; |
---|
809 | |
---|
810 | n = nn; |
---|
811 | nlen = NumBits(n); |
---|
812 | |
---|
813 | sub(m, n, 1); |
---|
814 | h = MakeOdd(m); // h >= 2 |
---|
815 | |
---|
816 | |
---|
817 | if (nlen > 50 && h < SqrRoot(nlen)) { |
---|
818 | long i, j; |
---|
819 | ZZ a, b, a_inv, c, r, m1, d; |
---|
820 | |
---|
821 | a = aa; |
---|
822 | InvMod(a_inv, a, n); |
---|
823 | |
---|
824 | if (h == 2) |
---|
825 | b = 2; |
---|
826 | else { |
---|
827 | do { |
---|
828 | RandomBnd(b, n); |
---|
829 | } while (Jacobi(b, n) != -1); |
---|
830 | } |
---|
831 | |
---|
832 | |
---|
833 | PowerMod(c, b, m, n); |
---|
834 | |
---|
835 | add(m1, m, 1); |
---|
836 | RightShift(m1, m1, 1); |
---|
837 | PowerMod(r, a, m1, n); |
---|
838 | |
---|
839 | for (i = h-2; i >= 0; i--) { |
---|
840 | SqrMod(d, r, n); |
---|
841 | MulMod(d, d, a_inv, n); |
---|
842 | for (j = 0; j < i; j++) |
---|
843 | SqrMod(d, d, n); |
---|
844 | if (!IsOne(d)) |
---|
845 | MulMod(r, r, c, n); |
---|
846 | SqrMod(c, c, n); |
---|
847 | } |
---|
848 | |
---|
849 | x = r; |
---|
850 | return; |
---|
851 | } |
---|
852 | |
---|
853 | |
---|
854 | |
---|
855 | |
---|
856 | |
---|
857 | long i, k; |
---|
858 | ZZ ma, t, u, v, e; |
---|
859 | ZZ t1, t2, t3, t4; |
---|
860 | |
---|
861 | n = nn; |
---|
862 | NegateMod(ma, aa, n); |
---|
863 | |
---|
864 | // find t such that t^2 - 4*a is not a square |
---|
865 | |
---|
866 | MulMod(t1, ma, 4, n); |
---|
867 | do { |
---|
868 | RandomBnd(t, n); |
---|
869 | SqrMod(t2, t, n); |
---|
870 | AddMod(t2, t2, t1, n); |
---|
871 | } while (Jacobi(t2, n) != -1); |
---|
872 | |
---|
873 | // compute u*X + v = X^{(n+1)/2} mod f, where f = X^2 - t*X + a |
---|
874 | |
---|
875 | add(e, n, 1); |
---|
876 | RightShift(e, e, 1); |
---|
877 | |
---|
878 | u = 0; |
---|
879 | v = 1; |
---|
880 | |
---|
881 | k = NumBits(e); |
---|
882 | |
---|
883 | for (i = k - 1; i >= 0; i--) { |
---|
884 | add(t2, u, v); |
---|
885 | sqr(t3, t2); // t3 = (u+v)^2 |
---|
886 | sqr(t1, u); |
---|
887 | sqr(t2, v); |
---|
888 | sub(t3, t3, t1); |
---|
889 | sub(t3, t3, t2); // t1 = u^2, t2 = v^2, t3 = 2*u*v |
---|
890 | rem(t1, t1, n); |
---|
891 | mul(t4, t1, t); |
---|
892 | add(t4, t4, t3); |
---|
893 | rem(u, t4, n); |
---|
894 | |
---|
895 | mul(t4, t1, ma); |
---|
896 | add(t4, t4, t2); |
---|
897 | rem(v, t4, n); |
---|
898 | |
---|
899 | if (bit(e, i)) { |
---|
900 | MulMod(t1, u, t, n); |
---|
901 | AddMod(t1, t1, v, n); |
---|
902 | MulMod(v, u, ma, n); |
---|
903 | u = t1; |
---|
904 | } |
---|
905 | |
---|
906 | } |
---|
907 | |
---|
908 | x = v; |
---|
909 | } |
---|
910 | |
---|
911 | |
---|
912 | |
---|
913 | // Chinese Remaindering. |
---|
914 | // |
---|
915 | // This version in new to v3.7, and is significantly |
---|
916 | // simpler and faster than the previous version. |
---|
917 | // |
---|
918 | // This function takes as input g, a, G, p, |
---|
919 | // such that a > 0, 0 <= G < p, and gcd(a, p) = 1. |
---|
920 | // It computes a' = a*p and g' such that |
---|
921 | // * g' = g (mod a); |
---|
922 | // * g' = G (mod p); |
---|
923 | // * -a'/2 < g' <= a'/2. |
---|
924 | // It then sets g := g' and a := a', and returns 1 iff g has changed. |
---|
925 | // |
---|
926 | // Under normal use, the input value g satisfies -a/2 < g <= a/2; |
---|
927 | // however, this was not documented or enforced in earlier versions, |
---|
928 | // so to maintain backward compatability, no restrictions are placed |
---|
929 | // on g. This routine runs faster, though, if -a/2 < g <= a/2, |
---|
930 | // and the first thing the routine does is to make this condition |
---|
931 | // hold. |
---|
932 | // |
---|
933 | // Also, under normal use, both a and p are odd; however, the routine |
---|
934 | // will still work even if this is not so. |
---|
935 | // |
---|
936 | // The routine is based on the following simple fact. |
---|
937 | // |
---|
938 | // Let -a/2 < g <= a/2, and let h satisfy |
---|
939 | // * g + a h = G (mod p); |
---|
940 | // * -p/2 < h <= p/2. |
---|
941 | // Further, if p = 2*h and g > 0, set |
---|
942 | // g' := g - a h; |
---|
943 | // otherwise, set |
---|
944 | // g' := g + a h. |
---|
945 | // Then g' so defined satisfies the above requirements. |
---|
946 | // |
---|
947 | // It is trivial to see that g's satisfies the congruence conditions. |
---|
948 | // The only thing is to check that the "balancing" condition |
---|
949 | // -a'/2 < g' <= a'/2 also holds. |
---|
950 | |
---|
951 | |
---|
952 | long CRT(ZZ& gg, ZZ& a, long G, long p) |
---|
953 | { |
---|
954 | if (p >= NTL_SP_BOUND) { |
---|
955 | ZZ GG, pp; |
---|
956 | conv(GG, G); |
---|
957 | conv(pp, p); |
---|
958 | return CRT(gg, a, GG, pp); |
---|
959 | } |
---|
960 | |
---|
961 | long modified = 0; |
---|
962 | |
---|
963 | static ZZ g; |
---|
964 | |
---|
965 | if (!CRTInRange(gg, a)) { |
---|
966 | modified = 1; |
---|
967 | ZZ a1; |
---|
968 | rem(g, gg, a); |
---|
969 | RightShift(a1, a, 1); |
---|
970 | if (g > a1) sub(g, g, a); |
---|
971 | } |
---|
972 | else |
---|
973 | g = gg; |
---|
974 | |
---|
975 | |
---|
976 | long p1; |
---|
977 | p1 = p >> 1; |
---|
978 | |
---|
979 | long a_inv; |
---|
980 | a_inv = rem(a, p); |
---|
981 | a_inv = InvMod(a_inv, p); |
---|
982 | |
---|
983 | long h; |
---|
984 | h = rem(g, p); |
---|
985 | h = SubMod(G, h, p); |
---|
986 | h = MulMod(h, a_inv, p); |
---|
987 | if (h > p1) |
---|
988 | h = h - p; |
---|
989 | |
---|
990 | if (h != 0) { |
---|
991 | modified = 1; |
---|
992 | |
---|
993 | if (!(p & 1) && g > 0 && (h == p1)) |
---|
994 | MulSubFrom(g, a, h); |
---|
995 | else |
---|
996 | MulAddTo(g, a, h); |
---|
997 | } |
---|
998 | |
---|
999 | mul(a, a, p); |
---|
1000 | gg = g; |
---|
1001 | |
---|
1002 | return modified; |
---|
1003 | } |
---|
1004 | |
---|
1005 | long CRT(ZZ& gg, ZZ& a, const ZZ& G, const ZZ& p) |
---|
1006 | { |
---|
1007 | long modified = 0; |
---|
1008 | |
---|
1009 | ZZ g; |
---|
1010 | |
---|
1011 | if (!CRTInRange(gg, a)) { |
---|
1012 | modified = 1; |
---|
1013 | ZZ a1; |
---|
1014 | rem(g, gg, a); |
---|
1015 | RightShift(a1, a, 1); |
---|
1016 | if (g > a1) sub(g, g, a); |
---|
1017 | } |
---|
1018 | else |
---|
1019 | g = gg; |
---|
1020 | |
---|
1021 | |
---|
1022 | ZZ p1; |
---|
1023 | RightShift(p1, p, 1); |
---|
1024 | |
---|
1025 | ZZ a_inv; |
---|
1026 | rem(a_inv, a, p); |
---|
1027 | InvMod(a_inv, a_inv, p); |
---|
1028 | |
---|
1029 | ZZ h; |
---|
1030 | rem(h, g, p); |
---|
1031 | SubMod(h, G, h, p); |
---|
1032 | MulMod(h, h, a_inv, p); |
---|
1033 | if (h > p1) |
---|
1034 | sub(h, h, p); |
---|
1035 | |
---|
1036 | if (h != 0) { |
---|
1037 | modified = 1; |
---|
1038 | ZZ ah; |
---|
1039 | mul(ah, a, h); |
---|
1040 | |
---|
1041 | if (!IsOdd(p) && g > 0 && (h == p1)) |
---|
1042 | sub(g, g, ah); |
---|
1043 | else |
---|
1044 | add(g, g, ah); |
---|
1045 | } |
---|
1046 | |
---|
1047 | mul(a, a, p); |
---|
1048 | gg = g; |
---|
1049 | |
---|
1050 | return modified; |
---|
1051 | } |
---|
1052 | |
---|
1053 | |
---|
1054 | |
---|
1055 | void sub(ZZ& x, const ZZ& a, long b) |
---|
1056 | { |
---|
1057 | static ZZ B; |
---|
1058 | conv(B, b); |
---|
1059 | sub(x, a, B); |
---|
1060 | } |
---|
1061 | |
---|
1062 | void sub(ZZ& x, long a, const ZZ& b) |
---|
1063 | { |
---|
1064 | static ZZ A; |
---|
1065 | conv(A, a); |
---|
1066 | sub(x, A, b); |
---|
1067 | } |
---|
1068 | |
---|
1069 | |
---|
1070 | void power2(ZZ& x, long e) |
---|
1071 | { |
---|
1072 | if (e < 0) Error("power2: negative exponent"); |
---|
1073 | set(x); |
---|
1074 | LeftShift(x, x, e); |
---|
1075 | } |
---|
1076 | |
---|
1077 | |
---|
1078 | void conv(ZZ& x, const char *s) |
---|
1079 | { |
---|
1080 | long c; |
---|
1081 | long cval; |
---|
1082 | long sign; |
---|
1083 | long ndigits; |
---|
1084 | long acc; |
---|
1085 | long i = 0; |
---|
1086 | |
---|
1087 | static ZZ a; |
---|
1088 | |
---|
1089 | if (!s) Error("bad ZZ input"); |
---|
1090 | |
---|
1091 | if (!iodigits) InitZZIO(); |
---|
1092 | |
---|
1093 | a = 0; |
---|
1094 | |
---|
1095 | c = s[i]; |
---|
1096 | while (IsWhiteSpace(c)) { |
---|
1097 | i++; |
---|
1098 | c = s[i]; |
---|
1099 | } |
---|
1100 | |
---|
1101 | if (c == '-') { |
---|
1102 | sign = -1; |
---|
1103 | i++; |
---|
1104 | c = s[i]; |
---|
1105 | } |
---|
1106 | else |
---|
1107 | sign = 1; |
---|
1108 | |
---|
1109 | cval = CharToIntVal(c); |
---|
1110 | if (cval < 0 || cval > 9) Error("bad ZZ input"); |
---|
1111 | |
---|
1112 | ndigits = 0; |
---|
1113 | acc = 0; |
---|
1114 | while (cval >= 0 && cval <= 9) { |
---|
1115 | acc = acc*10 + cval; |
---|
1116 | ndigits++; |
---|
1117 | |
---|
1118 | if (ndigits == iodigits) { |
---|
1119 | mul(a, a, ioradix); |
---|
1120 | add(a, a, acc); |
---|
1121 | ndigits = 0; |
---|
1122 | acc = 0; |
---|
1123 | } |
---|
1124 | |
---|
1125 | i++; |
---|
1126 | c = s[i]; |
---|
1127 | cval = CharToIntVal(c); |
---|
1128 | } |
---|
1129 | |
---|
1130 | if (ndigits != 0) { |
---|
1131 | long mpy = 1; |
---|
1132 | while (ndigits > 0) { |
---|
1133 | mpy = mpy * 10; |
---|
1134 | ndigits--; |
---|
1135 | } |
---|
1136 | |
---|
1137 | mul(a, a, mpy); |
---|
1138 | add(a, a, acc); |
---|
1139 | } |
---|
1140 | |
---|
1141 | if (sign == -1) |
---|
1142 | negate(a, a); |
---|
1143 | |
---|
1144 | x = a; |
---|
1145 | } |
---|
1146 | |
---|
1147 | |
---|
1148 | |
---|
1149 | void bit_and(ZZ& x, const ZZ& a, long b) |
---|
1150 | { |
---|
1151 | static ZZ B; |
---|
1152 | conv(B, b); |
---|
1153 | bit_and(x, a, B); |
---|
1154 | } |
---|
1155 | |
---|
1156 | void bit_or(ZZ& x, const ZZ& a, long b) |
---|
1157 | { |
---|
1158 | static ZZ B; |
---|
1159 | conv(B, b); |
---|
1160 | bit_or(x, a, B); |
---|
1161 | } |
---|
1162 | |
---|
1163 | void bit_xor(ZZ& x, const ZZ& a, long b) |
---|
1164 | { |
---|
1165 | static ZZ B; |
---|
1166 | conv(B, b); |
---|
1167 | bit_xor(x, a, B); |
---|
1168 | } |
---|
1169 | |
---|
1170 | |
---|
1171 | long power_long(long a, long e) |
---|
1172 | { |
---|
1173 | if (e < 0) Error("power_long: negative exponent"); |
---|
1174 | |
---|
1175 | if (e == 0) return 1; |
---|
1176 | |
---|
1177 | if (a == 1) return 1; |
---|
1178 | if (a == -1) { |
---|
1179 | if (e & 1) |
---|
1180 | return -1; |
---|
1181 | else |
---|
1182 | return 1; |
---|
1183 | } |
---|
1184 | |
---|
1185 | // no overflow check --- result is computed correctly |
---|
1186 | // modulo word size |
---|
1187 | |
---|
1188 | unsigned long res = 1; |
---|
1189 | unsigned long aa = a; |
---|
1190 | long i; |
---|
1191 | |
---|
1192 | for (i = 0; i < e; i++) |
---|
1193 | res *= aa; |
---|
1194 | |
---|
1195 | return to_long(res); |
---|
1196 | } |
---|
1197 | |
---|
1198 | // RANDOM NUMBER GENERATION |
---|
1199 | |
---|
1200 | // Idea for this PRNG. Iteratively hash seed using md5 |
---|
1201 | // to get 256 bytes to initialize arc4. |
---|
1202 | // Then use arc4 to get a pseudo-random byte stream. |
---|
1203 | |
---|
1204 | // I've taken care that the pseudo-random numbers generated by |
---|
1205 | // the routines RandomBnd, RandomBits, and RandomLen |
---|
1206 | // are completely platform independent. |
---|
1207 | |
---|
1208 | // I make use of the md5 compression function, |
---|
1209 | // which I've modified to work on 64-bit machines |
---|
1210 | |
---|
1211 | |
---|
1212 | /* |
---|
1213 | * BEGIN RSA's md5 stuff |
---|
1214 | * |
---|
1215 | */ |
---|
1216 | |
---|
1217 | /* |
---|
1218 | ********************************************************************** |
---|
1219 | ** md5.c ** |
---|
1220 | ** RSA Data Security, Inc. MD5 Message Digest Algorithm ** |
---|
1221 | ** Created: 2/17/90 RLR ** |
---|
1222 | ** Revised: 1/91 SRD,AJ,BSK,JT Reference C Version ** |
---|
1223 | ********************************************************************** |
---|
1224 | */ |
---|
1225 | |
---|
1226 | /* |
---|
1227 | ********************************************************************** |
---|
1228 | ** Copyright (C) 1990, RSA Data Security, Inc. All rights reserved. ** |
---|
1229 | ** ** |
---|
1230 | ** License to copy and use this software is granted provided that ** |
---|
1231 | ** it is identified as the "RSA Data Security, Inc. MD5 Message ** |
---|
1232 | ** Digest Algorithm" in all material mentioning or referencing this ** |
---|
1233 | ** software or this function. ** |
---|
1234 | ** ** |
---|
1235 | ** License is also granted to make and use derivative works ** |
---|
1236 | ** provided that such works are identified as "derived from the RSA ** |
---|
1237 | ** Data Security, Inc. MD5 Message Digest Algorithm" in all ** |
---|
1238 | ** material mentioning or referencing the derived work. ** |
---|
1239 | ** ** |
---|
1240 | ** RSA Data Security, Inc. makes no representations concerning ** |
---|
1241 | ** either the merchantability of this software or the suitability ** |
---|
1242 | ** of this software for any particular purpose. It is provided "as ** |
---|
1243 | ** is" without express or implied warranty of any kind. ** |
---|
1244 | ** ** |
---|
1245 | ** These notices must be retained in any copies of any part of this ** |
---|
1246 | ** documentation and/or software. ** |
---|
1247 | ********************************************************************** |
---|
1248 | */ |
---|
1249 | |
---|
1250 | |
---|
1251 | #if (NTL_BITS_PER_LONG <= 32) |
---|
1252 | #define TRUNC32(x) (x) |
---|
1253 | #else |
---|
1254 | #define TRUNC32(x) ((x) & ((1UL << 32)-1UL)) |
---|
1255 | #endif |
---|
1256 | |
---|
1257 | /* F, G and H are basic MD5 functions: selection, majority, parity */ |
---|
1258 | #define F(x, y, z) (((x) & (y)) | ((~x) & (z))) |
---|
1259 | #define G(x, y, z) (((x) & (z)) | ((y) & (~z))) |
---|
1260 | #define H(x, y, z) ((x) ^ (y) ^ (z)) |
---|
1261 | #define I(x, y, z) (TRUNC32((y) ^ ((x) | (~z)))) |
---|
1262 | |
---|
1263 | /* ROTATE_LEFT rotates x left n bits */ |
---|
1264 | #define ROTATE_LEFT(x, n) (TRUNC32(((x) << (n)) | ((x) >> (32-(n))))) |
---|
1265 | |
---|
1266 | /* FF, GG, HH, and II transformations for rounds 1, 2, 3, and 4 */ |
---|
1267 | /* Rotation is separate from addition to prevent recomputation */ |
---|
1268 | #define FF(a, b, c, d, x, s, ac) \ |
---|
1269 | {(a) = TRUNC32((a) + F((b), (c), (d)) + (x) + (ac)); \ |
---|
1270 | (a) = ROTATE_LEFT((a), (s)); \ |
---|
1271 | (a) = TRUNC32((a) + (b)); \ |
---|
1272 | } |
---|
1273 | #define GG(a, b, c, d, x, s, ac) \ |
---|
1274 | {(a) = TRUNC32((a) + G((b), (c), (d)) + (x) + (ac)); \ |
---|
1275 | (a) = ROTATE_LEFT((a), (s)); \ |
---|
1276 | (a) = TRUNC32((a) + (b)); \ |
---|
1277 | } |
---|
1278 | #define HH(a, b, c, d, x, s, ac) \ |
---|
1279 | {(a) = TRUNC32((a) + H((b), (c), (d)) + (x) + (ac)); \ |
---|
1280 | (a) = ROTATE_LEFT((a), (s)); \ |
---|
1281 | (a) = TRUNC32((a) + (b)); \ |
---|
1282 | } |
---|
1283 | #define II(a, b, c, d, x, s, ac) \ |
---|
1284 | {(a) = TRUNC32((a) + I((b), (c), (d)) + (x) + (ac)); \ |
---|
1285 | (a) = ROTATE_LEFT((a), (s)); \ |
---|
1286 | (a) = TRUNC32((a) + (b)); \ |
---|
1287 | } |
---|
1288 | |
---|
1289 | |
---|
1290 | |
---|
1291 | static |
---|
1292 | void MD5_default_IV(unsigned long *buf) |
---|
1293 | { |
---|
1294 | buf[0] = 0x67452301UL; |
---|
1295 | buf[1] = 0xefcdab89UL; |
---|
1296 | buf[2] = 0x98badcfeUL; |
---|
1297 | buf[3] = 0x10325476UL; |
---|
1298 | } |
---|
1299 | |
---|
1300 | |
---|
1301 | |
---|
1302 | /* Basic MD5 step. Transform buf based on in. |
---|
1303 | */ |
---|
1304 | |
---|
1305 | static |
---|
1306 | void MD5_compress(unsigned long *buf, unsigned long *in) |
---|
1307 | { |
---|
1308 | unsigned long a = buf[0], b = buf[1], c = buf[2], d = buf[3]; |
---|
1309 | |
---|
1310 | /* Round 1 */ |
---|
1311 | #define S11 7 |
---|
1312 | #define S12 12 |
---|
1313 | #define S13 17 |
---|
1314 | #define S14 22 |
---|
1315 | FF ( a, b, c, d, in[ 0], S11, 3614090360UL); /* 1 */ |
---|
1316 | FF ( d, a, b, c, in[ 1], S12, 3905402710UL); /* 2 */ |
---|
1317 | FF ( c, d, a, b, in[ 2], S13, 606105819UL); /* 3 */ |
---|
1318 | FF ( b, c, d, a, in[ 3], S14, 3250441966UL); /* 4 */ |
---|
1319 | FF ( a, b, c, d, in[ 4], S11, 4118548399UL); /* 5 */ |
---|
1320 | FF ( d, a, b, c, in[ 5], S12, 1200080426UL); /* 6 */ |
---|
1321 | FF ( c, d, a, b, in[ 6], S13, 2821735955UL); /* 7 */ |
---|
1322 | FF ( b, c, d, a, in[ 7], S14, 4249261313UL); /* 8 */ |
---|
1323 | FF ( a, b, c, d, in[ 8], S11, 1770035416UL); /* 9 */ |
---|
1324 | FF ( d, a, b, c, in[ 9], S12, 2336552879UL); /* 10 */ |
---|
1325 | FF ( c, d, a, b, in[10], S13, 4294925233UL); /* 11 */ |
---|
1326 | FF ( b, c, d, a, in[11], S14, 2304563134UL); /* 12 */ |
---|
1327 | FF ( a, b, c, d, in[12], S11, 1804603682UL); /* 13 */ |
---|
1328 | FF ( d, a, b, c, in[13], S12, 4254626195UL); /* 14 */ |
---|
1329 | FF ( c, d, a, b, in[14], S13, 2792965006UL); /* 15 */ |
---|
1330 | FF ( b, c, d, a, in[15], S14, 1236535329UL); /* 16 */ |
---|
1331 | |
---|
1332 | /* Round 2 */ |
---|
1333 | #define S21 5 |
---|
1334 | #define S22 9 |
---|
1335 | #define S23 14 |
---|
1336 | #define S24 20 |
---|
1337 | GG ( a, b, c, d, in[ 1], S21, 4129170786UL); /* 17 */ |
---|
1338 | GG ( d, a, b, c, in[ 6], S22, 3225465664UL); /* 18 */ |
---|
1339 | GG ( c, d, a, b, in[11], S23, 643717713UL); /* 19 */ |
---|
1340 | GG ( b, c, d, a, in[ 0], S24, 3921069994UL); /* 20 */ |
---|
1341 | GG ( a, b, c, d, in[ 5], S21, 3593408605UL); /* 21 */ |
---|
1342 | GG ( d, a, b, c, in[10], S22, 38016083UL); /* 22 */ |
---|
1343 | GG ( c, d, a, b, in[15], S23, 3634488961UL); /* 23 */ |
---|
1344 | GG ( b, c, d, a, in[ 4], S24, 3889429448UL); /* 24 */ |
---|
1345 | GG ( a, b, c, d, in[ 9], S21, 568446438UL); /* 25 */ |
---|
1346 | GG ( d, a, b, c, in[14], S22, 3275163606UL); /* 26 */ |
---|
1347 | GG ( c, d, a, b, in[ 3], S23, 4107603335UL); /* 27 */ |
---|
1348 | GG ( b, c, d, a, in[ 8], S24, 1163531501UL); /* 28 */ |
---|
1349 | GG ( a, b, c, d, in[13], S21, 2850285829UL); /* 29 */ |
---|
1350 | GG ( d, a, b, c, in[ 2], S22, 4243563512UL); /* 30 */ |
---|
1351 | GG ( c, d, a, b, in[ 7], S23, 1735328473UL); /* 31 */ |
---|
1352 | GG ( b, c, d, a, in[12], S24, 2368359562UL); /* 32 */ |
---|
1353 | |
---|
1354 | /* Round 3 */ |
---|
1355 | #define S31 4 |
---|
1356 | #define S32 11 |
---|
1357 | #define S33 16 |
---|
1358 | #define S34 23 |
---|
1359 | HH ( a, b, c, d, in[ 5], S31, 4294588738UL); /* 33 */ |
---|
1360 | HH ( d, a, b, c, in[ 8], S32, 2272392833UL); /* 34 */ |
---|
1361 | HH ( c, d, a, b, in[11], S33, 1839030562UL); /* 35 */ |
---|
1362 | HH ( b, c, d, a, in[14], S34, 4259657740UL); /* 36 */ |
---|
1363 | HH ( a, b, c, d, in[ 1], S31, 2763975236UL); /* 37 */ |
---|
1364 | HH ( d, a, b, c, in[ 4], S32, 1272893353UL); /* 38 */ |
---|
1365 | HH ( c, d, a, b, in[ 7], S33, 4139469664UL); /* 39 */ |
---|
1366 | HH ( b, c, d, a, in[10], S34, 3200236656UL); /* 40 */ |
---|
1367 | HH ( a, b, c, d, in[13], S31, 681279174UL); /* 41 */ |
---|
1368 | HH ( d, a, b, c, in[ 0], S32, 3936430074UL); /* 42 */ |
---|
1369 | HH ( c, d, a, b, in[ 3], S33, 3572445317UL); /* 43 */ |
---|
1370 | HH ( b, c, d, a, in[ 6], S34, 76029189UL); /* 44 */ |
---|
1371 | HH ( a, b, c, d, in[ 9], S31, 3654602809UL); /* 45 */ |
---|
1372 | HH ( d, a, b, c, in[12], S32, 3873151461UL); /* 46 */ |
---|
1373 | HH ( c, d, a, b, in[15], S33, 530742520UL); /* 47 */ |
---|
1374 | HH ( b, c, d, a, in[ 2], S34, 3299628645UL); /* 48 */ |
---|
1375 | |
---|
1376 | /* Round 4 */ |
---|
1377 | #define S41 6 |
---|
1378 | #define S42 10 |
---|
1379 | #define S43 15 |
---|
1380 | #define S44 21 |
---|
1381 | II ( a, b, c, d, in[ 0], S41, 4096336452UL); /* 49 */ |
---|
1382 | II ( d, a, b, c, in[ 7], S42, 1126891415UL); /* 50 */ |
---|
1383 | II ( c, d, a, b, in[14], S43, 2878612391UL); /* 51 */ |
---|
1384 | II ( b, c, d, a, in[ 5], S44, 4237533241UL); /* 52 */ |
---|
1385 | II ( a, b, c, d, in[12], S41, 1700485571UL); /* 53 */ |
---|
1386 | II ( d, a, b, c, in[ 3], S42, 2399980690UL); /* 54 */ |
---|
1387 | II ( c, d, a, b, in[10], S43, 4293915773UL); /* 55 */ |
---|
1388 | II ( b, c, d, a, in[ 1], S44, 2240044497UL); /* 56 */ |
---|
1389 | II ( a, b, c, d, in[ 8], S41, 1873313359UL); /* 57 */ |
---|
1390 | II ( d, a, b, c, in[15], S42, 4264355552UL); /* 58 */ |
---|
1391 | II ( c, d, a, b, in[ 6], S43, 2734768916UL); /* 59 */ |
---|
1392 | II ( b, c, d, a, in[13], S44, 1309151649UL); /* 60 */ |
---|
1393 | II ( a, b, c, d, in[ 4], S41, 4149444226UL); /* 61 */ |
---|
1394 | II ( d, a, b, c, in[11], S42, 3174756917UL); /* 62 */ |
---|
1395 | II ( c, d, a, b, in[ 2], S43, 718787259UL); /* 63 */ |
---|
1396 | II ( b, c, d, a, in[ 9], S44, 3951481745UL); /* 64 */ |
---|
1397 | |
---|
1398 | buf[0] = TRUNC32(buf[0] + a); |
---|
1399 | buf[1] = TRUNC32(buf[1] + b); |
---|
1400 | buf[2] = TRUNC32(buf[2] + c); |
---|
1401 | buf[3] = TRUNC32(buf[3] + d); |
---|
1402 | } |
---|
1403 | |
---|
1404 | |
---|
1405 | /* |
---|
1406 | * END RSA's md5 stuff |
---|
1407 | * |
---|
1408 | */ |
---|
1409 | |
---|
1410 | |
---|
1411 | static |
---|
1412 | void words_from_bytes(unsigned long *txtl, unsigned char *txtc, long n) |
---|
1413 | { |
---|
1414 | long i; |
---|
1415 | unsigned long v; |
---|
1416 | |
---|
1417 | for (i = 0; i < n; i++) { |
---|
1418 | v = txtc[4*i]; |
---|
1419 | v += ((unsigned long) (txtc[4*i+1])) << 8; |
---|
1420 | v += ((unsigned long) (txtc[4*i+2])) << 16; |
---|
1421 | v += ((unsigned long) (txtc[4*i+3])) << 24; |
---|
1422 | txtl[i] = v; |
---|
1423 | } |
---|
1424 | } |
---|
1425 | |
---|
1426 | static |
---|
1427 | void bytes_from_words(unsigned char *txtc, unsigned long *txtl, long n) |
---|
1428 | { |
---|
1429 | long i; |
---|
1430 | unsigned long v; |
---|
1431 | |
---|
1432 | for (i = 0; i < n; i++) { |
---|
1433 | v = txtl[i]; |
---|
1434 | txtc[4*i] = v & 255; |
---|
1435 | v = v >> 8; |
---|
1436 | txtc[4*i+1] = v & 255; |
---|
1437 | v = v >> 8; |
---|
1438 | txtc[4*i+2] = v & 255; |
---|
1439 | v = v >> 8; |
---|
1440 | txtc[4*i+3] = v & 255; |
---|
1441 | } |
---|
1442 | } |
---|
1443 | |
---|
1444 | |
---|
1445 | static |
---|
1446 | void MD5_compress1(unsigned long *buf, unsigned char *in, long n) |
---|
1447 | { |
---|
1448 | unsigned long txtl[16]; |
---|
1449 | unsigned char txtc[64]; |
---|
1450 | long i, j, k; |
---|
1451 | |
---|
1452 | if (n < 0) n = 0; |
---|
1453 | |
---|
1454 | i = 0; |
---|
1455 | while (i < n) { |
---|
1456 | k = n-i; |
---|
1457 | if (k > 64) k = 64; |
---|
1458 | for (j = 0; j < k; j++) |
---|
1459 | txtc[j] = in[i+j]; |
---|
1460 | for (; j < 64; j++) |
---|
1461 | txtc[j] = 0; |
---|
1462 | words_from_bytes(txtl, txtc, 16); |
---|
1463 | MD5_compress(buf, txtl); |
---|
1464 | i += k; |
---|
1465 | } |
---|
1466 | } |
---|
1467 | |
---|
1468 | |
---|
1469 | // the "cipherpunk" version of arc4 |
---|
1470 | |
---|
1471 | struct _ZZ_arc4_key |
---|
1472 | { |
---|
1473 | unsigned char state[256]; |
---|
1474 | unsigned char x; |
---|
1475 | unsigned char y; |
---|
1476 | }; |
---|
1477 | |
---|
1478 | |
---|
1479 | static inline |
---|
1480 | void swap_byte(unsigned char *a, unsigned char *b) |
---|
1481 | { |
---|
1482 | unsigned char swapByte; |
---|
1483 | |
---|
1484 | swapByte = *a; |
---|
1485 | *a = *b; |
---|
1486 | *b = swapByte; |
---|
1487 | } |
---|
1488 | |
---|
1489 | static |
---|
1490 | void prepare_key(unsigned char *key_data_ptr, |
---|
1491 | long key_data_len, _ZZ_arc4_key *key) |
---|
1492 | { |
---|
1493 | unsigned char index1; |
---|
1494 | unsigned char index2; |
---|
1495 | unsigned char* state; |
---|
1496 | long counter; |
---|
1497 | |
---|
1498 | state = &key->state[0]; |
---|
1499 | for(counter = 0; counter < 256; counter++) |
---|
1500 | state[counter] = counter; |
---|
1501 | key->x = 0; |
---|
1502 | key->y = 0; |
---|
1503 | index1 = 0; |
---|
1504 | index2 = 0; |
---|
1505 | for(counter = 0; counter < 256; counter++) |
---|
1506 | { |
---|
1507 | index2 = (key_data_ptr[index1] + state[counter] + index2) & 255; |
---|
1508 | swap_byte(&state[counter], &state[index2]); |
---|
1509 | |
---|
1510 | index1 = (index1 + 1) % key_data_len; |
---|
1511 | } |
---|
1512 | } |
---|
1513 | |
---|
1514 | |
---|
1515 | |
---|
1516 | static |
---|
1517 | void arc4(unsigned char *buffer_ptr, long buffer_len, _ZZ_arc4_key *key) |
---|
1518 | { |
---|
1519 | unsigned char x; |
---|
1520 | unsigned char y; |
---|
1521 | unsigned char* state; |
---|
1522 | unsigned char xorIndex; |
---|
1523 | long counter; |
---|
1524 | |
---|
1525 | x = key->x; |
---|
1526 | y = key->y; |
---|
1527 | |
---|
1528 | state = &key->state[0]; |
---|
1529 | for(counter = 0; counter < buffer_len; counter ++) |
---|
1530 | { |
---|
1531 | x = (x + 1) & 255; |
---|
1532 | y = (state[x] + y) & 255; |
---|
1533 | swap_byte(&state[x], &state[y]); |
---|
1534 | |
---|
1535 | xorIndex = (state[x] + state[y]) & 255; |
---|
1536 | |
---|
1537 | buffer_ptr[counter] = state[xorIndex]; |
---|
1538 | } |
---|
1539 | key->x = x; |
---|
1540 | key->y = y; |
---|
1541 | } |
---|
1542 | |
---|
1543 | // global state information for PRNG |
---|
1544 | |
---|
1545 | static long ran_initialized = 0; |
---|
1546 | static _ZZ_arc4_key ran_key; |
---|
1547 | |
---|
1548 | static unsigned long default_md5_tab[16] = { |
---|
1549 | 744663023UL, 1011602954UL, 3163087192UL, 3383838527UL, |
---|
1550 | 3305324122UL, 3197458079UL, 2266495600UL, 2760303563UL, |
---|
1551 | 346234297UL, 1919920720UL, 1896169861UL, 2192176675UL, |
---|
1552 | 2027150322UL, 2090160759UL, 2134858730UL, 1131796244UL |
---|
1553 | }; |
---|
1554 | |
---|
1555 | |
---|
1556 | |
---|
1557 | static |
---|
1558 | void build_arc4_tab(unsigned char *seed_bytes, const ZZ& s) |
---|
1559 | { |
---|
1560 | long nb = NumBytes(s); |
---|
1561 | |
---|
1562 | unsigned char *txt; |
---|
1563 | |
---|
1564 | typedef unsigned char u_char; |
---|
1565 | txt = NTL_NEW_OP u_char[nb + 68]; |
---|
1566 | if (!txt) Error("out of memory"); |
---|
1567 | |
---|
1568 | BytesFromZZ(txt + 4, s, nb); |
---|
1569 | |
---|
1570 | bytes_from_words(txt + nb + 4, default_md5_tab, 16); |
---|
1571 | |
---|
1572 | unsigned long buf[4]; |
---|
1573 | |
---|
1574 | unsigned long i; |
---|
1575 | for (i = 0; i < 16; i++) { |
---|
1576 | MD5_default_IV(buf); |
---|
1577 | bytes_from_words(txt, &i, 1); |
---|
1578 | |
---|
1579 | MD5_compress1(buf, txt, nb + 68); |
---|
1580 | |
---|
1581 | bytes_from_words(seed_bytes + 16*i, buf, 4); |
---|
1582 | } |
---|
1583 | |
---|
1584 | delete [] txt; |
---|
1585 | } |
---|
1586 | |
---|
1587 | |
---|
1588 | void SetSeed(const ZZ& s) |
---|
1589 | { |
---|
1590 | unsigned char seed_bytes[256]; |
---|
1591 | |
---|
1592 | build_arc4_tab(seed_bytes, s); |
---|
1593 | prepare_key(seed_bytes, 256, &ran_key); |
---|
1594 | |
---|
1595 | // we discard the first 1024 bytes of the arc4 stream, as this is |
---|
1596 | // recommended practice. |
---|
1597 | |
---|
1598 | arc4(seed_bytes, 256, &ran_key); |
---|
1599 | arc4(seed_bytes, 256, &ran_key); |
---|
1600 | arc4(seed_bytes, 256, &ran_key); |
---|
1601 | arc4(seed_bytes, 256, &ran_key); |
---|
1602 | |
---|
1603 | ran_initialized = 1; |
---|
1604 | } |
---|
1605 | |
---|
1606 | static |
---|
1607 | void ran_bytes(unsigned char *bytes, long n) |
---|
1608 | { |
---|
1609 | if (!ran_initialized) SetSeed(ZZ::zero()); |
---|
1610 | arc4(bytes, n, &ran_key); |
---|
1611 | } |
---|
1612 | |
---|
1613 | |
---|
1614 | unsigned long RandomWord() |
---|
1615 | { |
---|
1616 | unsigned char buf[NTL_BITS_PER_LONG/8]; |
---|
1617 | long i; |
---|
1618 | unsigned long res; |
---|
1619 | |
---|
1620 | ran_bytes(buf, NTL_BITS_PER_LONG/8); |
---|
1621 | |
---|
1622 | res = 0; |
---|
1623 | for (i = NTL_BITS_PER_LONG/8 - 1; i >= 0; i--) { |
---|
1624 | res = res << 8; |
---|
1625 | res = res | buf[i]; |
---|
1626 | } |
---|
1627 | |
---|
1628 | return res; |
---|
1629 | } |
---|
1630 | |
---|
1631 | long RandomBits_long(long l) |
---|
1632 | { |
---|
1633 | if (l <= 0) return 0; |
---|
1634 | if (l >= NTL_BITS_PER_LONG) |
---|
1635 | Error("RandomBits: length too big"); |
---|
1636 | |
---|
1637 | unsigned char buf[NTL_BITS_PER_LONG/8]; |
---|
1638 | unsigned long res; |
---|
1639 | long i; |
---|
1640 | |
---|
1641 | long nb = (l+7)/8; |
---|
1642 | ran_bytes(buf, nb); |
---|
1643 | |
---|
1644 | res = 0; |
---|
1645 | for (i = nb - 1; i >= 0; i--) { |
---|
1646 | res = res << 8; |
---|
1647 | res = res | buf[i]; |
---|
1648 | } |
---|
1649 | |
---|
1650 | return long(res & ((1UL << l)-1UL)); |
---|
1651 | } |
---|
1652 | |
---|
1653 | unsigned long RandomBits_ulong(long l) |
---|
1654 | { |
---|
1655 | if (l <= 0) return 0; |
---|
1656 | if (l > NTL_BITS_PER_LONG) |
---|
1657 | Error("RandomBits: length too big"); |
---|
1658 | |
---|
1659 | unsigned char buf[NTL_BITS_PER_LONG/8]; |
---|
1660 | unsigned long res; |
---|
1661 | long i; |
---|
1662 | |
---|
1663 | long nb = (l+7)/8; |
---|
1664 | ran_bytes(buf, nb); |
---|
1665 | |
---|
1666 | res = 0; |
---|
1667 | for (i = nb - 1; i >= 0; i--) { |
---|
1668 | res = res << 8; |
---|
1669 | res = res | buf[i]; |
---|
1670 | } |
---|
1671 | |
---|
1672 | if (l < NTL_BITS_PER_LONG) |
---|
1673 | res = res & ((1UL << l)-1UL); |
---|
1674 | |
---|
1675 | return res; |
---|
1676 | } |
---|
1677 | |
---|
1678 | long RandomLen_long(long l) |
---|
1679 | { |
---|
1680 | if (l <= 0) return 0; |
---|
1681 | if (l == 1) return 1; |
---|
1682 | if (l >= NTL_BITS_PER_LONG) |
---|
1683 | Error("RandomLen: length too big"); |
---|
1684 | |
---|
1685 | return RandomBits_long(l-1) + (1L << (l-1)); |
---|
1686 | } |
---|
1687 | |
---|
1688 | |
---|
1689 | void RandomBits(ZZ& x, long l) |
---|
1690 | { |
---|
1691 | if (l <= 0) { |
---|
1692 | x = 0; |
---|
1693 | return; |
---|
1694 | } |
---|
1695 | |
---|
1696 | if (NTL_OVERFLOW(l, 1, 0)) |
---|
1697 | Error("RandomBits: length too big"); |
---|
1698 | |
---|
1699 | long nb = (l+7)/8; |
---|
1700 | |
---|
1701 | static unsigned char *buf = 0; |
---|
1702 | static long buf_len = 0; |
---|
1703 | |
---|
1704 | if (nb > buf_len) { |
---|
1705 | if (buf) delete [] buf; |
---|
1706 | buf_len = ((nb + 1023)/1024)*1024; // allocate in 1024-byte lots |
---|
1707 | typedef unsigned char u_char; |
---|
1708 | buf = NTL_NEW_OP u_char[buf_len]; |
---|
1709 | if (!buf) Error("out of memory"); |
---|
1710 | } |
---|
1711 | |
---|
1712 | ran_bytes(buf, nb); |
---|
1713 | |
---|
1714 | static ZZ res; |
---|
1715 | |
---|
1716 | ZZFromBytes(res, buf, nb); |
---|
1717 | trunc(res, res, l); |
---|
1718 | |
---|
1719 | x = res; |
---|
1720 | } |
---|
1721 | |
---|
1722 | void RandomLen(ZZ& x, long l) |
---|
1723 | { |
---|
1724 | if (l <= 0) { |
---|
1725 | x = 0; |
---|
1726 | return; |
---|
1727 | } |
---|
1728 | |
---|
1729 | if (l == 1) { |
---|
1730 | x = 1; |
---|
1731 | return; |
---|
1732 | } |
---|
1733 | |
---|
1734 | if (NTL_OVERFLOW(l, 1, 0)) |
---|
1735 | Error("RandomLen: length too big"); |
---|
1736 | |
---|
1737 | // pre-allocate space to avoid two allocations |
---|
1738 | long nw = (l + NTL_ZZ_NBITS - 1)/NTL_ZZ_NBITS; |
---|
1739 | x.SetSize(nw); |
---|
1740 | |
---|
1741 | RandomBits(x, l-1); |
---|
1742 | SetBit(x, l-1); |
---|
1743 | } |
---|
1744 | |
---|
1745 | |
---|
1746 | const long RandomBndExcess = 8; |
---|
1747 | |
---|
1748 | |
---|
1749 | void RandomBnd(ZZ& x, const ZZ& bnd) |
---|
1750 | { |
---|
1751 | if (bnd <= 1) { |
---|
1752 | x = 0; |
---|
1753 | return; |
---|
1754 | } |
---|
1755 | |
---|
1756 | long k = NumBits(bnd); |
---|
1757 | |
---|
1758 | if (weight(bnd) == 1) { |
---|
1759 | RandomBits(x, k-1); |
---|
1760 | return; |
---|
1761 | } |
---|
1762 | |
---|
1763 | long l = k + RandomBndExcess; |
---|
1764 | |
---|
1765 | static ZZ t, r, t1; |
---|
1766 | |
---|
1767 | do { |
---|
1768 | RandomBits(t, l); |
---|
1769 | rem(r, t, bnd); |
---|
1770 | sub(t1, bnd, r); |
---|
1771 | add(t, t, t1); |
---|
1772 | } while (NumBits(t) > l); |
---|
1773 | |
---|
1774 | x = r; |
---|
1775 | } |
---|
1776 | |
---|
1777 | long RandomBnd(long bnd) |
---|
1778 | { |
---|
1779 | if (bnd <= 1) return 0; |
---|
1780 | |
---|
1781 | long k = NumBits(bnd); |
---|
1782 | |
---|
1783 | if (((bnd - 1) & bnd) == 0) |
---|
1784 | return RandomBits_long(k-1); |
---|
1785 | |
---|
1786 | long l = k + RandomBndExcess; |
---|
1787 | |
---|
1788 | if (l > NTL_BITS_PER_LONG-2) { |
---|
1789 | static ZZ Bnd, res; |
---|
1790 | |
---|
1791 | Bnd = bnd; |
---|
1792 | RandomBnd(res, Bnd); |
---|
1793 | return to_long(res); |
---|
1794 | } |
---|
1795 | |
---|
1796 | long t, r; |
---|
1797 | |
---|
1798 | do { |
---|
1799 | t = RandomBits_long(l); |
---|
1800 | r = t % bnd; |
---|
1801 | } while (t + bnd - r > (1L << l)); |
---|
1802 | |
---|
1803 | return r; |
---|
1804 | } |
---|
1805 | |
---|
1806 | |
---|
1807 | |
---|
1808 | |
---|
1809 | // More prime generation stuff... |
---|
1810 | |
---|
1811 | static |
---|
1812 | double Log2(double x) |
---|
1813 | { |
---|
1814 | static double log2 = log(2.0); |
---|
1815 | return log(x)/log2; |
---|
1816 | } |
---|
1817 | |
---|
1818 | // Define p(k,t) to be the conditional probability that a random, odd, k-bit |
---|
1819 | // number is composite, given that it passes t iterations of the |
---|
1820 | // Miller-Rabin test. |
---|
1821 | // This routine returns 0 or 1, and if it returns 1 then |
---|
1822 | // p(k,t) <= 2^{-n}. |
---|
1823 | // This basically encodes the estimates of Damgard, Landrock, and Pomerance; |
---|
1824 | // it uses floating point arithmetic, but is coded in such a way |
---|
1825 | // that its results should be correct, assuming that the log function |
---|
1826 | // is computed with reasonable precision. |
---|
1827 | // |
---|
1828 | // It is assumed that k >= 3 and t >= 1; if this does not hold, |
---|
1829 | // then 0 is returned. |
---|
1830 | |
---|
1831 | static |
---|
1832 | long ErrBoundTest(long kk, long tt, long nn) |
---|
1833 | |
---|
1834 | { |
---|
1835 | const double fudge = (1.0 + 1024.0/NTL_FDOUBLE_PRECISION); |
---|
1836 | const double log2_3 = Log2(3.0); |
---|
1837 | const double log2_7 = Log2(7.0); |
---|
1838 | const double log2_20 = Log2(20.0); |
---|
1839 | |
---|
1840 | double k = kk; |
---|
1841 | double t = tt; |
---|
1842 | double n = nn; |
---|
1843 | |
---|
1844 | if (k < 3 || t < 1) return 0; |
---|
1845 | if (n < 1) return 1; |
---|
1846 | |
---|
1847 | // the following test is largely academic |
---|
1848 | if (9*t > NTL_FDOUBLE_PRECISION) Error("ErrBoundTest: t too big"); |
---|
1849 | |
---|
1850 | double log2_k = Log2(k); |
---|
1851 | |
---|
1852 | if ((n + log2_k)*fudge <= 2*t) |
---|
1853 | return 1; |
---|
1854 | |
---|
1855 | if ((2*log2_k + 4.0 + n)*fudge <= 2*sqrt(k)) |
---|
1856 | return 2; |
---|
1857 | |
---|
1858 | if ((t == 2 && k >= 88) || (3 <= t && 9*t <= k && k >= 21)) { |
---|
1859 | if ((1.5*log2_k + t + 4.0 + n)*fudge <= 0.5*Log2(t) + 2*(sqrt(t*k))) |
---|
1860 | return 3; |
---|
1861 | } |
---|
1862 | |
---|
1863 | if (k <= 9*t && 4*t <= k && k >= 21) { |
---|
1864 | if ( ((log2_3 + log2_7 + log2_k + n)*fudge <= log2_20 + 5*t) && |
---|
1865 | ((log2_3 + (15.0/4.0)*log2_k + n)*fudge <= log2_7 + k/2 + 2*t) && |
---|
1866 | ((2*log2_3 + 2 + log2_k + n)*fudge <= k/4 + 3*t) ) |
---|
1867 | return 4; |
---|
1868 | } |
---|
1869 | |
---|
1870 | if (4*t >= k && k >= 21) { |
---|
1871 | if (((15.0/4.0)*log2_k + n)*fudge <= log2_7 + k/2 + 2*t) |
---|
1872 | return 5; |
---|
1873 | } |
---|
1874 | |
---|
1875 | return 0; |
---|
1876 | } |
---|
1877 | |
---|
1878 | |
---|
1879 | void GenPrime(ZZ& n, long k, long err) |
---|
1880 | { |
---|
1881 | if (k <= 1) Error("GenPrime: bad length"); |
---|
1882 | |
---|
1883 | if (k > (1L << 20)) Error("GenPrime: length too large"); |
---|
1884 | |
---|
1885 | if (err < 1) err = 1; |
---|
1886 | if (err > 512) err = 512; |
---|
1887 | |
---|
1888 | if (k == 2) { |
---|
1889 | if (RandomBnd(2)) |
---|
1890 | n = 3; |
---|
1891 | else |
---|
1892 | n = 2; |
---|
1893 | |
---|
1894 | return; |
---|
1895 | } |
---|
1896 | |
---|
1897 | |
---|
1898 | long t; |
---|
1899 | |
---|
1900 | t = 1; |
---|
1901 | while (!ErrBoundTest(k, t, err)) |
---|
1902 | t++; |
---|
1903 | |
---|
1904 | RandomPrime(n, k, t); |
---|
1905 | } |
---|
1906 | |
---|
1907 | |
---|
1908 | long GenPrime_long(long k, long err) |
---|
1909 | { |
---|
1910 | if (k <= 1) Error("GenPrime: bad length"); |
---|
1911 | |
---|
1912 | if (k >= NTL_BITS_PER_LONG) Error("GenPrime: length too large"); |
---|
1913 | |
---|
1914 | if (err < 1) err = 1; |
---|
1915 | if (err > 512) err = 512; |
---|
1916 | |
---|
1917 | if (k == 2) { |
---|
1918 | if (RandomBnd(2)) |
---|
1919 | return 3; |
---|
1920 | else |
---|
1921 | return 2; |
---|
1922 | } |
---|
1923 | |
---|
1924 | long t; |
---|
1925 | |
---|
1926 | t = 1; |
---|
1927 | while (!ErrBoundTest(k, t, err)) |
---|
1928 | t++; |
---|
1929 | |
---|
1930 | return RandomPrime_long(k, t); |
---|
1931 | } |
---|
1932 | |
---|
1933 | |
---|
1934 | void GenGermainPrime(ZZ& n, long k, long err) |
---|
1935 | { |
---|
1936 | if (k <= 1) Error("GenGermainPrime: bad length"); |
---|
1937 | |
---|
1938 | if (k > (1L << 20)) Error("GenGermainPrime: length too large"); |
---|
1939 | |
---|
1940 | if (err < 1) err = 1; |
---|
1941 | if (err > 512) err = 512; |
---|
1942 | |
---|
1943 | if (k == 2) { |
---|
1944 | if (RandomBnd(2)) |
---|
1945 | n = 3; |
---|
1946 | else |
---|
1947 | n = 2; |
---|
1948 | |
---|
1949 | return; |
---|
1950 | } |
---|
1951 | |
---|
1952 | |
---|
1953 | long prime_bnd = ComputePrimeBound(k); |
---|
1954 | |
---|
1955 | if (NumBits(prime_bnd) >= k/2) |
---|
1956 | prime_bnd = (1L << (k/2-1)); |
---|
1957 | |
---|
1958 | |
---|
1959 | ZZ two; |
---|
1960 | two = 2; |
---|
1961 | |
---|
1962 | ZZ n1; |
---|
1963 | |
---|
1964 | |
---|
1965 | PrimeSeq s; |
---|
1966 | |
---|
1967 | ZZ iter; |
---|
1968 | iter = 0; |
---|
1969 | |
---|
1970 | |
---|
1971 | for (;;) { |
---|
1972 | iter++; |
---|
1973 | |
---|
1974 | RandomLen(n, k); |
---|
1975 | if (!IsOdd(n)) add(n, n, 1); |
---|
1976 | |
---|
1977 | s.reset(3); |
---|
1978 | long p; |
---|
1979 | |
---|
1980 | long sieve_passed = 1; |
---|
1981 | |
---|
1982 | p = s.next(); |
---|
1983 | while (p && p < prime_bnd) { |
---|
1984 | long r = rem(n, p); |
---|
1985 | |
---|
1986 | if (r == 0) { |
---|
1987 | sieve_passed = 0; |
---|
1988 | break; |
---|
1989 | } |
---|
1990 | |
---|
1991 | // test if 2*r + 1 = 0 (mod p) |
---|
1992 | if (r == p-r-1) { |
---|
1993 | sieve_passed = 0; |
---|
1994 | break; |
---|
1995 | } |
---|
1996 | |
---|
1997 | p = s.next(); |
---|
1998 | } |
---|
1999 | |
---|
2000 | if (!sieve_passed) continue; |
---|
2001 | |
---|
2002 | |
---|
2003 | if (MillerWitness(n, two)) continue; |
---|
2004 | |
---|
2005 | // n1 = 2*n+1 |
---|
2006 | mul(n1, n, 2); |
---|
2007 | add(n1, n1, 1); |
---|
2008 | |
---|
2009 | |
---|
2010 | if (MillerWitness(n1, two)) continue; |
---|
2011 | |
---|
2012 | // now do t M-R iterations...just to make sure |
---|
2013 | |
---|
2014 | // First compute the appropriate number of M-R iterations, t |
---|
2015 | // The following computes t such that |
---|
2016 | // p(k,t)*8/k <= 2^{-err}/(5*iter^{1.25}) |
---|
2017 | // which suffices to get an overall error probability of 2^{-err}. |
---|
2018 | // Note that this method has the advantage of not requiring |
---|
2019 | // any assumptions on the density of Germain primes. |
---|
2020 | |
---|
2021 | long err1 = max(1, err + 7 + (5*NumBits(iter) + 3)/4 - NumBits(k)); |
---|
2022 | long t; |
---|
2023 | t = 1; |
---|
2024 | while (!ErrBoundTest(k, t, err1)) |
---|
2025 | t++; |
---|
2026 | |
---|
2027 | ZZ W; |
---|
2028 | long MR_passed = 1; |
---|
2029 | |
---|
2030 | long i; |
---|
2031 | for (i = 1; i <= t; i++) { |
---|
2032 | do { |
---|
2033 | RandomBnd(W, n); |
---|
2034 | } while (W == 0); |
---|
2035 | // W == 0 is not a useful candidate witness! |
---|
2036 | |
---|
2037 | if (MillerWitness(n, W)) { |
---|
2038 | MR_passed = 0; |
---|
2039 | break; |
---|
2040 | } |
---|
2041 | } |
---|
2042 | |
---|
2043 | if (MR_passed) break; |
---|
2044 | } |
---|
2045 | } |
---|
2046 | |
---|
2047 | long GenGermainPrime_long(long k, long err) |
---|
2048 | { |
---|
2049 | if (k >= NTL_BITS_PER_LONG-1) |
---|
2050 | Error("GenGermainPrime_long: length too long"); |
---|
2051 | |
---|
2052 | ZZ n; |
---|
2053 | GenGermainPrime(n, k, err); |
---|
2054 | return to_long(n); |
---|
2055 | } |
---|
2056 | |
---|
2057 | |
---|
2058 | NTL_END_IMPL |
---|