![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to genprime.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [PGP] / pgp / src|
Return to genprime.c CVS log | Up to [PGP] / pgp / src |
1.1.1.7 ! root 1: /* genprime.c - C source code for generation of large primes ! 2: used by public-key key generation routines. ! 3: First version 17 Mar 87 ! 4: Last revised 2 Jun 91 by PRZ ! 5: 24 Apr 93 by CP ! 6: ! 7: (c) Copyright 1990-1994 by Philip Zimmermann. All rights reserved. ! 8: The author assumes no liability for damages resulting from the use ! 9: of this software, even if the damage results from defects in this ! 10: software. No warranty is expressed or implied. ! 11: ! 12: Note that while most PGP source modules bear Philip Zimmermann's ! 13: copyright notice, many of them have been revised or entirely written ! 14: by contributors who frequently failed to put their names in their ! 15: code. Code that has been incorporated into PGP from other authors ! 16: was either originally published in the public domain or is used with ! 17: permission from the various authors. ! 18: ! 19: PGP is available for free to the public under certain restrictions. ! 20: See the PGP User's Guide (included in the release package) for ! 21: important information about licensing, patent restrictions on ! 22: certain algorithms, trademarks, copyrights, and export controls. ! 23: ! 24: These functions are for the generation of large prime integers and ! 25: for other functions related to factoring and key generation for ! 26: many number-theoretic cryptographic algorithms, such as the NIST ! 27: Digital Signature Standard. ! 28: */ 1.1.1.6 root 29: 30: #define SHOWPROGRESS 31: 32: /* Define some error status returns for keygen... */ 33: #define NOPRIMEFOUND -14 /* slowtest probably failed */ 34: #define NOSUSPECTS -13 /* fastsieve probably failed */ 35: 36: 37: #ifdef MSDOS 38: #define poll_for_break() {while (kbhit()) getch();} 39: #endif 40: 41: #ifndef poll_for_break 1.1.1.7 ! root 42: #define poll_for_break() /* stub */ 1.1.1.6 root 43: #endif 44: 45: #ifdef SHOWPROGRESS 1.1.1.7 ! root 46: #include <stdio.h> /* needed for putchar() */ 1.1.1.6 root 47: #endif 48: 1.1.1.7 ! root 49: #ifdef EMBEDDED /* compiling for embedded target */ ! 50: #define _NOMALLOC /* defined if no malloc is available. */ ! 51: #endif /* EMBEDDED */ 1.1.1.6 root 52: 53: /* Decide whether malloc is available. Some embedded systems lack it. */ 1.1.1.7 ! root 54: #ifndef _NOMALLOC /* malloc library routine available */ ! 55: #include <stdlib.h> /* ANSI C library - for malloc() and free() */ ! 56: /* #include <alloc.h> *//* Borland Turbo C has malloc in <alloc.h> */ ! 57: #endif /* malloc available */ 1.1.1.6 root 58: 59: #include "mpilib.h" 60: #include "genprime.h" 61: #if defined(MSDOS) && !defined(__GO32__) 62: #include <conio.h> 63: #endif 64: 65: #include "random.h" 66: 67: 1.1.1.7 ! root 68: /* #define STRONGPRIMES *//* if defined, generate "strong" primes for key */ 1.1.1.6 root 69: /* 70: *"Strong" primes are no longer advantageous, due to the new 71: * elliptical curve method of factoring. Randomly selected primes 72: * are as good as any. See "Factoring", by Duncan A. Buell, Journal 73: * of Supercomputing 1 (1987), pages 191-216. 74: * This justifies disabling the lengthy search for strong primes. 75: * 76: * The advice about strong primes in the early RSA literature applies 77: * to 256-bit moduli where the attacks were the Pollard rho and P-1 78: * factoring algorithms. Later developments in factoring have entirely 79: * supplanted these methods. The later algorithms are always faster 80: * (so we need bigger primes), and don't care about STRONGPRIMES. 81: * 82: * The early literature was saying that you can get away with small 83: * moduli if you choose the primes carefully. The later developments 84: * say you can't get away with small moduli, period. And it doesn't 85: * matter how you choose the primes. 86: * 87: * It's just taking a heck of a long time for the advice on "strong primes" 88: * to disappear from the books. Authors keep going back to the original 89: * documents and repeating what they read there, even though it's out 90: * of date. 91: */ 92: 93: #define BLUM 94: /* If BLUM is defined, this looks for prines congruent to 3 modulo 4. 95: The product of two of these is a Blum integer. You can uniquely define 96: a square root Cmodulo a Blum integer, which leads to some extra 97: possibilities for encryption algorithms. This shrinks the key space by 98: 2 bits, which is not considered significant. 1.1.1.7 ! root 99: */ 1.1.1.6 root 100: 101: #ifdef STRONGPRIMES 102: 103: static boolean primetest(unitptr p); 104: /* Returns TRUE iff p is a prime. */ 105: 1.1.1.7 ! root 106: static int mp_sqrt(unitptr quotient, unitptr dividend); 1.1.1.6 root 107: /* Quotient is returned as the square root of dividend. */ 108: 109: #endif 110: 111: static int nextprime(unitptr p); 112: /* Find next higher prime starting at p, returning result in p. */ 113: 1.1.1.7 ! root 114: static void randombits(unitptr p, short nbits); 1.1.1.6 root 115: /* Make a random unit array p with nbits of precision. */ 116: 117: #ifdef DEBUG 118: #define DEBUGprintf1(x) printf(x) 119: #define DEBUGprintf2(x,y) printf(x,y) 120: #define DEBUGprintf3(x,y,z) printf(x,y,z) 121: #else 122: #define DEBUGprintf1(x) 123: #define DEBUGprintf2(x,y) 124: #define DEBUGprintf3(x,y,z) 125: #endif 126: 127: 1.1.1.7 ! root 128: /* primetable is a table of 16-bit prime numbers used for sieving ! 129: and for other aspects of public-key cryptographic key generation */ 1.1.1.6 root 130: 1.1.1.7 ! root 131: static word16 primetable[] = ! 132: { ! 133: 2, 3, 5, 7, 11, 13, 17, 19, ! 134: 23, 29, 31, 37, 41, 43, 47, 53, ! 135: 59, 61, 67, 71, 73, 79, 83, 89, ! 136: 97, 101, 103, 107, 109, 113, 127, 131, ! 137: 137, 139, 149, 151, 157, 163, 167, 173, ! 138: 179, 181, 191, 193, 197, 199, 211, 223, ! 139: 227, 229, 233, 239, 241, 251, 257, 263, ! 140: 269, 271, 277, 281, 283, 293, 307, 311, ! 141: #ifndef EMBEDDED /* not embedded, use larger table */ ! 142: 313, 317, 331, 337, 347, 349, 353, 359, ! 143: 367, 373, 379, 383, 389, 397, 401, 409, ! 144: 419, 421, 431, 433, 439, 443, 449, 457, ! 145: 461, 463, 467, 479, 487, 491, 499, 503, ! 146: 509, 521, 523, 541, 547, 557, 563, 569, ! 147: 571, 577, 587, 593, 599, 601, 607, 613, ! 148: 617, 619, 631, 641, 643, 647, 653, 659, ! 149: 661, 673, 677, 683, 691, 701, 709, 719, ! 150: 727, 733, 739, 743, 751, 757, 761, 769, ! 151: 773, 787, 797, 809, 811, 821, 823, 827, ! 152: 829, 839, 853, 857, 859, 863, 877, 881, ! 153: 883, 887, 907, 911, 919, 929, 937, 941, ! 154: 947, 953, 967, 971, 977, 983, 991, 997, ! 155: 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, ! 156: 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, ! 157: 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, ! 158: 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, ! 159: 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, ! 160: 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, ! 161: 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, ! 162: 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, ! 163: 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, ! 164: 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, ! 165: 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, ! 166: 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, ! 167: 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, ! 168: 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, ! 169: 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, ! 170: 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, ! 171: 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, ! 172: 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, ! 173: 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, ! 174: 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, ! 175: 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, ! 176: 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, ! 177: 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, ! 178: 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, ! 179: 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, ! 180: 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, ! 181: 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, ! 182: 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, ! 183: 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, ! 184: 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, ! 185: 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, ! 186: 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, ! 187: 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, ! 188: 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, ! 189: 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, ! 190: 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, ! 191: 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, ! 192: 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, ! 193: 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, ! 194: 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, ! 195: 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, ! 196: 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, ! 197: 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, ! 198: 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, ! 199: 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, ! 200: 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, ! 201: 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, ! 202: 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, ! 203: 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, ! 204: 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, ! 205: 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, ! 206: 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, ! 207: 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, ! 208: 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, ! 209: 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, ! 210: 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, ! 211: 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, ! 212: 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, ! 213: 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, ! 214: 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, ! 215: 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, ! 216: 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, ! 217: 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, ! 218: 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, ! 219: 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, ! 220: 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, ! 221: 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, ! 222: 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, ! 223: 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, ! 224: 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, ! 225: 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, ! 226: 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, ! 227: 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, ! 228: 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, ! 229: 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, ! 230: 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, ! 231: 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, ! 232: 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, ! 233: 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, ! 234: 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, ! 235: 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, ! 236: 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, ! 237: 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, ! 238: 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, ! 239: 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, ! 240: 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, ! 241: 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, ! 242: 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, ! 243: 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, ! 244: 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, ! 245: 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, ! 246: 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, ! 247: 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, ! 248: 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, ! 249: 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, ! 250: 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, ! 251: 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, ! 252: 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, ! 253: 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, ! 254: 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, ! 255: 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, ! 256: 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, ! 257: 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, ! 258: 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, ! 259: 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, ! 260: 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, ! 261: 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, ! 262: 8167, 8171, 8179, 8191, ! 263: #endif /* not EMBEDDED, use larger table */ ! 264: 0}; /* null-terminated list, with only one null at end */ 1.1.1.6 root 265: 266: 267: 268: #ifdef UNIT8 269: static word16 bottom16(unitptr r) 270: /* Called from nextprime and primetest. Returns low 16 bits of r. */ 271: { 1.1.1.7 ! root 272: make_lsbptr(r, (global_precision - ((2 / BYTES_PER_UNIT) - 1))); ! 273: return *(word16 *) r; ! 274: } /* bottom16 */ ! 275: #else /* UNIT16 or UNIT32 */ 1.1.1.6 root 276: #define bottom16(r) ((word16) lsunit(r)) 1.1.1.7 ! root 277: /* or UNIT32 could mask off lower 16 bits, instead of typecasting it. */ ! 278: #endif /* UNIT16 or UNIT32 */ 1.1.1.6 root 279: 280: 281: /* 282: * This routine tests p for primality by applying Fermat's theorem: 283: * For any x, if ((x**(p-1)) mod p) != 1, then p is not prime. 284: * By trying a few values for x, we can determine if p is "probably" prime. 285: * 286: * Because this test is so slow, it is recommended that p be sieved first 287: * to weed out numbers that are obviously not prime. 288: * 289: * Contrary to what you may have read in the literature, empirical evidence 290: * shows this test weeds out a LOT more than 50% of the composite candidates 291: * for each trial x. Each test catches nearly all the composites. 292: */ 1.1.1.7 ! root 293: static boolean slowtest(unitptr p) 1.1.1.6 root 294: { 1.1.1.7 ! root 295: unit x[MAX_UNIT_PRECISION], is_one[MAX_UNIT_PRECISION]; ! 296: unit pminus1[MAX_UNIT_PRECISION]; ! 297: short i; ! 298: ! 299: mp_move(pminus1, p); ! 300: mp_dec(pminus1); ! 301: ! 302: for (i = 0; i < 4; i++) { /* Just do a few tests. */ ! 303: poll_for_break(); /* polls keyboard, allows ctrl-C to abort program */ ! 304: mp_init(x, primetable[i]); /* Use any old random trial x */ ! 305: /* if ((x**(p-1)) mod p) != 1, then p is not prime */ ! 306: if (mp_modexp(is_one, x, pminus1, p) < 0) /* modexp error? */ ! 307: return FALSE; /* error means return not prime status */ ! 308: if (testne(is_one, 1)) /* then p is not prime */ ! 309: return FALSE; /* return not prime status */ 1.1.1.6 root 310: #ifdef SHOWPROGRESS 1.1.1.7 ! root 311: putchar('+'); /* let user see how we are progressing */ ! 312: fflush(stdout); ! 313: #endif /* SHOWPROGRESS */ ! 314: } ! 315: ! 316: /* If it gets to this point, it's very likely that p is prime */ ! 317: mp_burn(x); /* burn the evidence on the stack... */ ! 318: mp_burn(is_one); ! 319: mp_burn(pminus1); ! 320: return TRUE; ! 321: } /* slowtest -- fermattest */ 1.1.1.6 root 322: 323: 324: #ifdef STRONGPRIMES 325: 326: static boolean primetest(unitptr p) 327: /* 328: * Returns TRUE iff p is a prime. 329: * If p doesn't pass through the sieve, then p is definitely NOT a prime. 1.1.1.7 ! root 330: * If p is small enough for the sieve to prove primality or not, 1.1.1.6 root 331: * and p passes through the sieve, then p is definitely a prime. 332: * If p is large and p passes through the sieve and may be a prime, 333: * then p is further tested for primality with a slower test. 334: */ 335: { 1.1.1.7 ! root 336: short i; ! 337: static word16 lastprime = 0; /* last prime in primetable */ ! 338: word16 sqrt_p; /* to limit sieving past sqrt(p), for small p's */ ! 339: ! 340: if (!lastprime) { /* lastprime still undefined. So define it. */ ! 341: /* executes this code only once, then skips it next time */ ! 342: for (i = 0; primetable[i]; i++); /* seek end of primetable */ ! 343: lastprime = primetable[i - 1]; /* get last prime in table */ ! 344: } ! 345: if (significance(p) <= (2 / BYTES_PER_UNIT)) /* if p <= 16 bits */ ! 346: /* p may be in primetable. Search it. */ ! 347: if (bottom16(p) <= lastprime) ! 348: for (i = 0; primetable[i]; i++) { ! 349: /* scan until null-terminator */ ! 350: if (primetable[i] == bottom16(p)) ! 351: return TRUE; /* yep, definitely a prime. */ ! 352: if (primetable[i] > bottom16(p)) /* we missed. */ ! 353: return FALSE; /* definitely NOT a prime. */ ! 354: } /* We got past the whole primetable without a hit. */ ! 355: /* p is bigger than any prime in primetable, so let's sieve. */ ! 356: if (!(lsunit(p) & 1)) /* if least significant bit is 0... */ ! 357: return FALSE; /* divisible by 2, not prime */ ! 358: ! 359: if (mp_tstminus(p)) /* error if p<0 */ ! 360: return FALSE; /* not prime if p<0 */ ! 361: ! 362: /* ! 363: * Optimization for small (32 bits or less) p's. ! 364: * If p is small, compute sqrt_p = sqrt(p), or else ! 365: * if p is >32 bits then just set sqrt_p to something ! 366: * at least as big as the largest primetable entry. ! 367: */ ! 368: if (significance(p) <= (4 / BYTES_PER_UNIT)) { /* if p <= 32 bits */ ! 369: unit sqrtp[MAX_UNIT_PRECISION]; ! 370: /* Just sieve up to sqrt(p) */ ! 371: if (mp_sqrt(sqrtp, p) == 0) /* 0 means p is a perfect square */ ! 372: return FALSE; /* perfect square is not a prime */ ! 373: /* we know that sqrtp <= 16 bits because p <= 32 bits */ ! 374: sqrt_p = bottom16(sqrtp); ! 375: } else { ! 376: /* p > 32 bits, so obviate sqrt(p) test. */ ! 377: sqrt_p = lastprime; /* ensures that we do ENTIRE sieve. */ ! 378: } ! 379: ! 380: /* p is assumed odd, so begin sieve at 3 */ ! 381: for (i = 1; primetable[i]; i++) { ! 382: /* Compute p mod (primetable[i]). If it divides evenly... */ ! 383: if (mp_shortmod(p, primetable[i]) == 0) ! 384: return FALSE; /* then p is definitely NOT prime */ ! 385: if (primetable[i] > sqrt_p) /* fully sieved p? */ ! 386: return TRUE; /* yep, fully passed sieve, definitely a prime. */ ! 387: } ! 388: /* It passed the sieve, so p is a suspected prime. */ ! 389: ! 390: /* Now try slow complex primality test on suspected prime. */ ! 391: return slowtest(p); /* returns TRUE or FALSE */ ! 392: } /* primetest */ 1.1.1.6 root 393: 394: #endif 395: 396: /* 397: * Used in conjunction with fastsieve. Builds a table of remainders 398: * relative to the random starting point p, so that fastsieve can 399: * sequentially sieve for suspected primes quickly. Call buildsieve 400: * once, then call fastsieve for consecutive prime candidates. 401: * Note that p must be odd, because the sieve begins at 3. 402: */ 1.1.1.7 ! root 403: static void buildsieve(unitptr p, word16 remainders[]) 1.1.1.6 root 404: { 1.1.1.7 ! root 405: short i; ! 406: for (i = 1; primetable[i]; i++) { ! 407: remainders[i] = mp_shortmod(p, primetable[i]); ! 408: } ! 409: } /* buildsieve */ 1.1.1.6 root 410: 411: /* 1.1.1.7 ! root 412: Fast prime sieving algorithm by Philip Zimmermann, March 1987. ! 413: This is the fastest algorithm I know of for quickly sieving for ! 414: large (hundreds of bits in length) "random" probable primes, because ! 415: it uses only single-precision (16-bit) arithmetic. Because rigorous ! 416: prime testing algorithms are very slow, it is recommended that ! 417: potential prime candidates be quickly passed through this fast ! 418: sieving algorithm first to weed out numbers that are obviously not ! 419: prime. ! 420: ! 421: This algorithm is optimized to search sequentially for a large prime ! 422: from a random starting point. For generalized nonsequential prime ! 423: testing, the slower conventional sieve should be used, as given ! 424: in primetest(p). ! 425: ! 426: This algorithm requires a fixed table (called primetable) of the ! 427: first hundred or so small prime numbers. ! 428: First we select a random odd starting point (p) for our prime ! 429: search. Then we build a table of 16-bit remainders calculated ! 430: from that initial p. This table of 16-bit remainders is exactly ! 431: the same length as the table of small 16-bit primes. Each ! 432: remainders table entry contains the remainder of p divided by the ! 433: corresponding primetable entry. Then we begin sequentially testing ! 434: all odd integers, starting from the initial odd random p. The ! 435: distance we have searched from the huge random starting point p is ! 436: a small 16-bit number, pdelta. If pdelta plus a remainders table ! 437: entry is evenly divisible by the corresponding primetable entry, ! 438: then p+pdelta is factorable by that primetable entry, which means ! 439: p+pdelta is not prime. ! 440: */ 1.1.1.6 root 441: 1.1.1.7 ! root 442: /* Fastsieve is used for searching sequentially from a random starting ! 443: point for a suspected prime. Requires that buildsieve be called ! 444: first, to build a table of remainders relative to the random starting ! 445: point p. ! 446: Returns true iff pdelta passes through the sieve and thus p+pdelta ! 447: may be a prime. Note that p must be odd, because the sieve begins ! 448: at 3. ! 449: */ 1.1.1.6 root 450: static boolean fastsieve(word16 pdelta, word16 remainders[]) 1.1.1.7 ! root 451: { ! 452: short i; ! 453: for (i = 1; primetable[i]; i++) { ! 454: /* ! 455: * If pdelta plus a remainders table entry is evenly ! 456: * divisible by the corresponding primetable entry, ! 457: * then p+pdelta is factorable by that primetable entry, ! 458: * which means p+pdelta is not prime. ! 459: */ ! 460: if ((pdelta + remainders[i]) % primetable[i] == 0) ! 461: return FALSE; /* then p+pdelta is not prime */ ! 462: } ! 463: /* It passed the sieve. It is now a suspected prime. */ ! 464: return TRUE; ! 465: } /* fastsieve */ 1.1.1.6 root 466: 467: 1.1.1.7 ! root 468: #define numberof(x) (sizeof(x)/sizeof(x[0])) /* number of table entries */ 1.1.1.6 root 469: 470: 471: static int nextprime(unitptr p) 472: /* 473: * Find next higher prime starting at p, returning result in p. 474: * Uses fast prime sieving algorithm to search sequentially. 475: * Returns 0 for normal completion status, < 0 for failure status. 476: */ 477: { 1.1.1.7 ! root 478: word16 pdelta, range; ! 479: short oldprecision; ! 480: short i, suspects; ! 481: ! 482: /* start search at candidate p */ ! 483: mp_inc(p); /* remember, it's the NEXT prime from p, noninclusive. */ ! 484: if (significance(p) <= 1) { ! 485: /* ! 486: * p might be smaller than the largest prime in primetable. ! 487: * We can't sieve for primes that are already in primetable. ! 488: * We will have to directly search the table. ! 489: */ ! 490: /* scan until null-terminator */ ! 491: for (i = 0; primetable[i]; i++) { ! 492: if (primetable[i] >= lsunit(p)) { ! 493: mp_init(p, primetable[i]); ! 494: return 0; /* return next higher prime from primetable */ ! 495: } ! 496: } /* We got past the whole primetable without a hit. */ ! 497: } /* p is bigger than any prime in primetable, so let's sieve. */ ! 498: if (mp_tstminus(p)) { /* error if p<0 */ ! 499: mp_init(p, 2); /* next prime >0 is 2 */ ! 500: return 0; /* normal completion status */ ! 501: } 1.1.1.6 root 502: #ifndef BLUM 1.1.1.7 ! root 503: lsunit(p) |= 1; /* set candidate's lsb - make it odd */ 1.1.1.6 root 504: #else 1.1.1.7 ! root 505: lsunit(p) |= 3; /* Make candidate ==3 mod 4 */ 1.1.1.6 root 506: #endif 507: 1.1.1.7 ! root 508: /* Adjust the global_precision downward to the optimum size for p... */ ! 509: oldprecision = global_precision; /* save global_precision */ ! 510: /* We will need 2-3 extra bits of precision for the falsekeytest. */ ! 511: set_precision(bits2units(countbits(p) + 4 + SLOP_BITS)); ! 512: /* Rescale p to global_precision we just defined */ ! 513: rescale(p, oldprecision, global_precision); 1.1.1.6 root 514: 1.1.1.7 ! root 515: { 1.1.1.6 root 516: #ifdef _NOMALLOC /* No malloc and free functions available. Use stack. */ 1.1.1.7 ! root 517: word16 remainders[numberof(primetable)]; ! 518: #else /* malloc available, we can conserve stack space. */ ! 519: word16 *remainders; ! 520: /* Allocate some memory for the table of remainders: */ ! 521: remainders = (word16 *) malloc(sizeof(primetable)); ! 522: #endif /* malloc available */ ! 523: ! 524: /* Build remainders table relative to initial p: */ ! 525: buildsieve(p, remainders); ! 526: pdelta = 0; /* offset from initial random p */ ! 527: /* Sieve preparation complete. Now for some fast fast sieving... */ ! 528: /* slowtest will not be called unless fastsieve is true */ 1.1.1.6 root 529: 1.1.1.7 ! root 530: /* range is how far to search before giving up. */ 1.1.1.6 root 531: #ifndef BLUM 1.1.1.7 ! root 532: range = 4 * units2bits(global_precision); 1.1.1.6 root 533: #else 1.1.1.7 ! root 534: /* Twice as many because step size is twice as large, */ ! 535: range = 8 * units2bits(global_precision); 1.1.1.6 root 536: #endif 1.1.1.7 ! root 537: suspects = 0; /* number of suspected primes and slowtest trials */ ! 538: for (;;) { ! 539: /* equivalent to: if (primetest(p)) break; */ ! 540: if (fastsieve(pdelta, remainders)) { /* found suspected prime */ ! 541: suspects++; /* tally for statistical purposes */ 1.1.1.6 root 542: #ifdef SHOWPROGRESS 1.1.1.7 ! root 543: putchar('.'); /* let user see how we are progressing */ ! 544: fflush(stdout); ! 545: #endif /* SHOWPROGRESS */ ! 546: if (slowtest(p)) ! 547: break; /* found a prime */ ! 548: } 1.1.1.6 root 549: #ifndef BLUM 1.1.1.7 ! root 550: pdelta += 2; /* try next odd number */ 1.1.1.6 root 551: #else 1.1.1.7 ! root 552: pdelta += 4; ! 553: mp_inc(p); ! 554: mp_inc(p); 1.1.1.6 root 555: #endif 1.1.1.7 ! root 556: mp_inc(p); ! 557: mp_inc(p); 1.1.1.6 root 558: 1.1.1.7 ! root 559: if (pdelta > range) /* searched too many candidates? */ ! 560: break; /* something must be wrong--bail out of search */ 1.1.1.6 root 561: 1.1.1.7 ! root 562: } /* while (TRUE) */ 1.1.1.6 root 563: 564: #ifdef SHOWPROGRESS 1.1.1.7 ! root 565: putchar(' '); /* let user see how we are progressing */ ! 566: #endif /* SHOWPROGRESS */ 1.1.1.6 root 567: 1.1.1.7 ! root 568: for (i = 0; primetable[i]; i++) /* scan until null-terminator */ ! 569: remainders[i] = 0; /* don't leave remainders exposed in RAM */ 1.1.1.6 root 570: #ifndef _NOMALLOC 1.1.1.7 ! root 571: free(remainders); /* free allocated memory */ ! 572: #endif /* not _NOMALLOC */ ! 573: } ! 574: ! 575: set_precision(oldprecision); /* restore precision */ ! 576: ! 577: if (pdelta > range) { /* searched too many candidates? */ ! 578: if (suspects < 1) /* unreasonable to have found no suspects */ ! 579: return NOSUSPECTS; /* fastsieve failed, probably */ ! 580: return NOPRIMEFOUND; /* return error status */ ! 581: } ! 582: return 0; /* return normal completion status */ 1.1.1.6 root 583: 1.1.1.7 ! root 584: } /* nextprime */ 1.1.1.6 root 585: 586: 587: /* We will need a series of truly random bits for key generation. 588: In most implementations, our random number supply is derived from 589: random keyboard delays rather than a hardware random number 590: chip. So we will have to ensure we have a large enough pool of 591: accumulated random numbers from the keyboard. trueRandAccum() 592: performs this operation. 1.1.1.7 ! root 593: */ 1.1.1.6 root 594: 595: /* Fills 1 unit with random bytes, and returns unit. */ 1.1.1.7 ! root 596: static unit randomunit(void) 1.1.1.6 root 597: { 1.1.1.7 ! root 598: unit u = 0; ! 599: byte i; ! 600: i = BYTES_PER_UNIT; ! 601: do ! 602: u = (u << 8) + trueRandByte(); ! 603: while (--i != 0); ! 604: return u; ! 605: } /* randomunit */ 1.1.1.6 root 606: 607: /* 608: * Make a random unit array p with nbits of precision. Used mainly to 609: * generate large random numbers to search for primes. 610: */ 1.1.1.7 ! root 611: static void randombits(unitptr p, short nbits) 1.1.1.6 root 612: { 1.1.1.7 ! root 613: mp_init(p, 0); ! 614: make_lsbptr(p, global_precision); 1.1.1.6 root 615: 1.1.1.7 ! root 616: /* Add whole units of randomness */ ! 617: while (nbits >= UNITSIZE) { ! 618: *post_higherunit(p) = randomunit(); ! 619: nbits -= UNITSIZE; ! 620: } ! 621: ! 622: /* Add most-significant partial unit (if any) */ ! 623: if (nbits) ! 624: *p = randomunit() & (power_of_2(nbits) - 1); ! 625: } /* randombits */ 1.1.1.6 root 626: 627: /* 628: * Makes a "random" prime p with nbits significant bits of precision. 629: * Since these primes are used to compute a modulus of a guaranteed 630: * length, the top 2 bits of the prime are set to 1, so that the 631: * product of 2 primes (the modulus) is of a deterministic length. 632: * Returns 0 for normal completion status, < 0 for failure status. 633: */ 1.1.1.7 ! root 634: int randomprime(unitptr p, short nbits) 1.1.1.6 root 635: { 1.1.1.7 ! root 636: DEBUGprintf2("\nGenerating a %d-bit random prime. ", nbits); ! 637: /* Get an initial random candidate p to start search. */ ! 638: randombits(p, nbits - 2); /* 2 less random bits for nonrandom top bits */ ! 639: /* To guarantee exactly nbits of significance, set the top 2 bits to 1 */ ! 640: mp_setbit(p, nbits - 1); /* highest bit is nonrandom */ ! 641: mp_setbit(p, nbits - 2); /* next highest bit is also nonrandom */ ! 642: return nextprime(p); /* search for next higher prime ! 643: from starting point p */ ! 644: } /* randomprime */ ! 645: 1.1.1.6 root 646: 1.1.1.7 ! root 647: #ifdef STRONGPRIMES /* generate "strong" primes for keys */ 1.1.1.6 root 648: 1.1.1.7 ! root 649: #define log_1stprime 6 /* log base 2 of firstprime */ 1.1.1.6 root 650: 1.1.1.7 ! root 651: /* 1st primetable entry used by tryprime */ ! 652: #define firstprime (1<<log_1stprime) 1.1.1.6 root 653: 654: /* This routine attempts to generate a prime p such that p-1 has prime p1 655: as its largest factor. Prime p will have no more than maxbits bits of 656: significance. Prime p1 must be less than maxbits-log_1stprime in length. 657: This routine is called only from goodprime. 1.1.1.7 ! root 658: */ ! 659: static boolean tryprime(unitptr p, unitptr p1, short maxbits) 1.1.1.6 root 660: { 1.1.1.7 ! root 661: int i; ! 662: unit i2[MAX_UNIT_PRECISION]; ! 663: /* Generate p such that p = (i*2*p1)+1, for i=1,2,3,5,7,11,13,17... ! 664: and test p for primality for each small prime i. ! 665: It's better to start i at firstprime rather than at 1, ! 666: because then p grows slower in significance. ! 667: Start looking for small primes that are > firstprime... ! 668: */ ! 669: if ((countbits(p1) + log_1stprime) >= maxbits) { ! 670: DEBUGprintf1("\007[Error: overconstrained prime]"); 1.1.1.6 root 671: return FALSE; /* failed to make a good prime */ 1.1.1.7 ! root 672: } ! 673: for (i = 0; primetable[i]; i++) { ! 674: if (primetable[i] < firstprime) ! 675: continue; ! 676: /* note that mp_init doesn't extend sign bit for >32767 */ ! 677: mp_init(i2, primetable[i] << 1); ! 678: mp_mult(p, p1, i2); ! 679: mp_inc(p); ! 680: if (countbits(p) > maxbits) ! 681: break; ! 682: DEBUGprintf1("."); ! 683: if (primetest(p)) ! 684: return TRUE; ! 685: } ! 686: return FALSE; /* failed to make a good prime */ ! 687: } /* tryprime */ 1.1.1.6 root 688: 689: 690: /* 691: * Make a "strong" prime p with at most maxbits and at least minbits of 692: * significant bits of precision. This algorithm is called to generate 693: * a high-quality prime p for key generation purposes. It must have 694: * special characteristics for making a modulus n that is hard to factor. 695: * Returns 0 for normal completion status, < 0 for failure status. 696: */ 1.1.1.7 ! root 697: int goodprime(unitptr p, short maxbits, short minbits) 1.1.1.6 root 698: { 1.1.1.7 ! root 699: unit p1[MAX_UNIT_PRECISION]; ! 700: short oldprecision, midbits; ! 701: int status; ! 702: ! 703: mp_init(p, 0); ! 704: /* Adjust the global_precision downward to the optimum size for p... */ ! 705: oldprecision = global_precision; /* save global_precision */ ! 706: /* We will need 2-3 extra bits of precision for the falsekeytest. */ ! 707: set_precision(bits2units(maxbits + 4 + SLOP_BITS)); ! 708: /* rescale p to global_precision we just defined */ ! 709: rescale(p, oldprecision, global_precision); ! 710: ! 711: minbits -= 2 * log_1stprime; /* length of p" */ ! 712: midbits = (maxbits + minbits) / 2; /* length of p' */ ! 713: DEBUGprintf3("\nGenerating a %d-%d bit refined prime. ", ! 714: minbits + 2 * log_1stprime, maxbits); ! 715: do { 1.1.1.6 root 716: do { 1.1.1.7 ! root 717: status = randomprime(p, minbits - 1); ! 718: if (status < 0) ! 719: return status; /* failed to find a random prime */ ! 720: DEBUGprintf2("\n(p\042=%d bits)", countbits(p)); ! 721: } while (!tryprime(p1, p, midbits)); ! 722: DEBUGprintf2("(p'=%d bits)", countbits(p1)); ! 723: } while (!tryprime(p, p1, maxbits)); ! 724: DEBUGprintf2("\n\007(p=%d bits) ", countbits(p)); ! 725: mp_burn(p1); /* burn the evidence on the stack */ ! 726: set_precision(oldprecision); /* restore precision */ ! 727: return 0; /* normal completion status */ ! 728: } /* goodprime */ 1.1.1.6 root 729: 1.1.1.7 ! root 730: #endif /* STRONGPRIMES */ 1.1.1.6 root 731: 732: 733: #define iplus1 ( i==2 ? 0 : i+1 ) /* used by Euclid algorithms */ 734: #define iminus1 ( i==0 ? 2 : i-1 ) /* used by Euclid algorithms */ 735: 736: /* Computes greatest common divisor via Euclid's algorithm. */ 1.1.1.7 ! root 737: void mp_gcd(unitptr result, unitptr a, unitptr n) 1.1.1.6 root 738: { 1.1.1.7 ! root 739: short i; ! 740: unit gcopies[3][MAX_UNIT_PRECISION]; 1.1.1.6 root 741: #define g(i) ( &(gcopies[i][0]) ) 1.1.1.7 ! root 742: mp_move(g(0), n); ! 743: mp_move(g(1), a); 1.1.1.6 root 744: 1.1.1.7 ! root 745: i = 1; ! 746: while (testne(g(i), 0)) { ! 747: mp_mod(g(iplus1), g(iminus1), g(i)); ! 748: i = iplus1; ! 749: } ! 750: mp_move(result, g(iminus1)); ! 751: mp_burn(g(iminus1)); /* burn the evidence on the stack... */ ! 752: mp_burn(g(iplus1)); ! 753: #undef g ! 754: } /* mp_gcd */ 1.1.1.6 root 755: 756: /* 757: * Euclid's algorithm extended to compute multiplicative inverse. 758: * Computes x such that a*x mod n = 1, where 0<a<n 759: * 760: * The variable u is unnecessary for the algorithm, but is 761: * included in comments for mathematical clarity. 762: */ 1.1.1.7 ! root 763: void mp_inv(unitptr x, unitptr a, unitptr n) 1.1.1.6 root 764: { 1.1.1.7 ! root 765: short i; ! 766: unit y[MAX_UNIT_PRECISION], temp[MAX_UNIT_PRECISION]; ! 767: unit gcopies[3][MAX_UNIT_PRECISION], vcopies[3][MAX_UNIT_PRECISION]; 1.1.1.6 root 768: #define g(i) ( &(gcopies[i][0]) ) 769: #define v(i) ( &(vcopies[i][0]) ) 1.1.1.7 ! root 770: /* unit ucopies[3][MAX_UNIT_PRECISION]; */ 1.1.1.6 root 771: /* #define u(i) ( &(ucopies[i][0]) ) */ 1.1.1.7 ! root 772: mp_move(g(0), n); ! 773: mp_move(g(1), a); ! 774: /* mp_init(u(0),1); mp_init(u(1),0); */ ! 775: mp_init(v(0), 0); ! 776: mp_init(v(1), 1); ! 777: i = 1; ! 778: while (testne(g(i), 0)) { ! 779: /* we know that at this point, g(i) = u(i)*n + v(i)*a */ ! 780: mp_udiv(g(iplus1), y, g(iminus1), g(i)); ! 781: mp_mult(temp, y, v(i)); ! 782: mp_move(v(iplus1), v(iminus1)); ! 783: mp_sub(v(iplus1), temp); ! 784: /* mp_mult(temp,y,u(i)); mp_move(u(iplus1),u(iminus1)); ! 785: mp_sub(u(iplus1),temp); */ ! 786: i = iplus1; ! 787: } ! 788: mp_move(x, v(iminus1)); ! 789: if (mp_tstminus(x)) ! 790: mp_add(x, n); ! 791: mp_burn(g(iminus1)); /* burn the evidence on the stack... */ ! 792: mp_burn(g(iplus1)); ! 793: mp_burn(v(0)); ! 794: mp_burn(v(1)); ! 795: mp_burn(v(2)); ! 796: mp_burn(y); ! 797: mp_burn(temp); 1.1.1.6 root 798: #undef g 799: #undef v 1.1.1.7 ! root 800: } /* mp_inv */ 1.1.1.6 root 801: 802: #ifdef STRONGPRIMES 803: 1.1.1.7 ! root 804: /* mp_sqrt - returns square root of a number. ! 805: returns -1 for error, 0 for perfect square, 1 for not perfect square. ! 806: Not used by any RSA-related functions. Used by factoring algorithms. ! 807: This version needs optimization. ! 808: by Charles W. Merritt July 15, 1989, refined by PRZ. ! 809: ! 810: These are notes on computing the square root the manual old-fashioned ! 811: way. This is the basis of the fast sqrt algorithm mp_sqrt below: ! 812: ! 813: 1) Separate the number into groups (periods) of two digits each, ! 814: beginning with units or at the decimal point. ! 815: 2) Find the greatest perfect square in the left hand period & write ! 816: its square root as the first figure of the required root. Subtract ! 817: the square of this number from the left hand period. Annex to the ! 818: remainder the next group so as to form a dividend. ! 819: 3) Double the root already found and write it as a partial divisor at ! 820: the left of the new dividend. Annex one zero digit, making a trial ! 821: divisor, and divide the new dividend by the trial divisor. ! 822: 4) Write the quotient in the root as the trial term and also substitute ! 823: this quotient for the annexed zero digit in the partial divisor, ! 824: making the latter complete. ! 825: 5) Multiply the complete divisor by the figure just obtained and, ! 826: if possible, subtract the product from the last remainder. ! 827: If this product is too large, the trial term of the quotient ! 828: must be replaced by the next smaller number and the operations ! 829: preformed as before. ! 830: (IN BINARY, OUR TRIAL TERM IS ALWAYS 1 AND WE USE IT OR DO NOTHING.) ! 831: 6) Proceed in this manner until all periods are used. ! 832: If there is still a remainder, it's not a perfect square. ! 833: */ ! 834: 1.1.1.6 root 835: /* Quotient is returned as the square root of dividend. */ 1.1.1.7 ! root 836: static int mp_sqrt(unitptr quotient, unitptr dividend) 1.1.1.6 root 837: { 1.1.1.7 ! root 838: register short next2bits; /* "period", or group of 2 bits of dividend */ ! 839: register unit dvdbitmask, qbitmask; ! 840: unit remainder[MAX_UNIT_PRECISION]; ! 841: unit rjq[MAX_UNIT_PRECISION], divisor[MAX_UNIT_PRECISION]; ! 842: unsigned int qbits, qprec, dvdbits, dprec, oldprecision; ! 843: int notperfect; ! 844: ! 845: mp_init(quotient, 0); ! 846: if (mp_tstminus(dividend)) { /* if dividend<0, return error */ ! 847: mp_dec(quotient); /* quotient = -1 */ ! 848: return -1; ! 849: } ! 850: /* normalize and compute number of bits in dividend first */ ! 851: init_bitsniffer(dividend, dvdbitmask, dprec, dvdbits); ! 852: /* init_bitsniffer returns a 0 if dvdbits is 0 */ ! 853: if (dvdbits == 1) { ! 854: mp_init(quotient, 1); /* square root of 1 is 1 */ ! 855: return 0; ! 856: } ! 857: /* rescale quotient to half the precision of dividend */ ! 858: qbits = (dvdbits + 1) >> 1; ! 859: qprec = bits2units(qbits); ! 860: rescale(quotient, global_precision, qprec); ! 861: make_msbptr(quotient, qprec); ! 862: qbitmask = power_of_2((qbits - 1) & (UNITSIZE - 1)); ! 863: ! 864: /* ! 865: * Set smallest optimum precision for this square root. ! 866: * The low-level primitives are affected by the call to set_precision. ! 867: * Even though the dividend precision is bigger than the precision ! 868: * we will use, no low-level primitives will be used on the dividend. ! 869: * They will be used on the quotient, remainder, and rjq, which are ! 870: * smaller precision. ! 871: */ ! 872: oldprecision = global_precision; /* save global_precision */ ! 873: set_precision(bits2units(qbits + 3)); /* 3 bits of precision slop */ ! 874: ! 875: /* special case: sqrt of 1st 2 (binary) digits of dividend ! 876: is 1st (binary) digit of quotient. This is always 1. */ ! 877: stuff_bit(quotient, qbitmask); ! 878: bump_bitsniffer(quotient, qbitmask); ! 879: mp_init(rjq, 1); /* rjq is Right Justified Quotient */ ! 880: ! 881: if (!(dvdbits & 1)) { ! 882: /* even number of bits in dividend */ ! 883: next2bits = 2; ! 884: bump_bitsniffer(dividend, dvdbitmask); ! 885: dvdbits--; ! 886: if (sniff_bit(dividend, dvdbitmask)) ! 887: next2bits++; ! 888: bump_bitsniffer(dividend, dvdbitmask); ! 889: dvdbits--; ! 890: } else { ! 891: /* odd number of bits in dividend */ ! 892: next2bits = 1; ! 893: bump_bitsniffer(dividend, dvdbitmask); ! 894: dvdbits--; ! 895: } ! 896: ! 897: mp_init(remainder, next2bits - 1); ! 898: ! 899: /* dvdbits is guaranteed to be even at this point */ ! 900: ! 901: while (dvdbits) { ! 902: next2bits = 0; ! 903: if (sniff_bit(dividend, dvdbitmask)) ! 904: next2bits = 2; ! 905: bump_bitsniffer(dividend, dvdbitmask); ! 906: dvdbits--; ! 907: if (sniff_bit(dividend, dvdbitmask)) ! 908: next2bits++; ! 909: bump_bitsniffer(dividend, dvdbitmask); ! 910: dvdbits--; ! 911: mp_rotate_left(remainder, (boolean) ((next2bits & 2) != 0)); ! 912: mp_rotate_left(remainder, (boolean) ((next2bits & 1) != 0)); 1.1.1.6 root 913: 914: /* 1.1.1.7 ! root 915: * "divisor" is trial divisor, complete divisor is 4*rjq ! 916: * or 4*rjq+1. ! 917: * Subtract divisor times its last digit from remainder. ! 918: * If divisor ends in 1, remainder -= divisor*1, ! 919: * or if divisor ends in 0, remainder -= divisor*0 (do nothing). ! 920: * Last digit of divisor inflates divisor as large as possible ! 921: * yet still subtractable from remainder. 1.1.1.6 root 922: */ 1.1.1.7 ! root 923: mp_move(divisor, rjq); /* divisor = 4*rjq+1 */ ! 924: mp_rotate_left(divisor, 0); ! 925: mp_rotate_left(divisor, 1); ! 926: if (mp_compare(remainder, divisor) >= 0) { ! 927: mp_sub(remainder, divisor); ! 928: stuff_bit(quotient, qbitmask); ! 929: mp_rotate_left(rjq, 1); 1.1.1.6 root 930: } else { 1.1.1.7 ! root 931: mp_rotate_left(rjq, 0); 1.1.1.6 root 932: } 1.1.1.7 ! root 933: bump_bitsniffer(quotient, qbitmask); ! 934: } ! 935: notperfect = testne(remainder, 0); /* not a perfect square? */ ! 936: set_precision(oldprecision); /* restore original precision */ ! 937: return notperfect; /* normal return */ ! 938: } /* mp_sqrt */ 1.1.1.6 root 939: #endif 940: 1.1.1.7 ! root 941: /*------------------- End of genprime.c -----------------------------*/
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.