Return to genprime.c CVS log

Up to [PGP] / pgp / src

Annotation of pgp/src/genprime.c, revision 1.1.1.9

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.