Return to genprime.c CVS log

Up to [PGP] / pgp / src

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

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

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