/* * cryrand - cryptographically strong pseudo-random number generator library */ /* * Copyright (c) 1994 by Landon Curt Noll. All Rights Reserved. * * Permission to use, copy, modify, and distribute this software and * its documentation for any purpose and without fee is hereby granted, * provided that the above copyright, this permission notice and text * this comment, and the disclaimer below appear in all of the following: * * supporting documentation * source copies * source works derived from this source * binaries derived from this source or from derived source * * LANDON CURT NOLL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO * EVENT SHALL LANDON CURT NOLL BE LIABLE FOR ANY SPECIAL, INDIRECT OR * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF * USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR * PERFORMANCE OF THIS SOFTWARE. * * chongo was here /\../\ */ /* * These routines are written in the calc language. At the time of this * notice, calc was at version 2.9.2 (We refer to calc, as in the C * program, not the Emacs subsystem). * * Calc is available by anonymous ftp from ftp.uu.net under the * directory /pub/calc. * * If you can't get calc any other way, EMail a request to my EMail * address as shown below. * * Comments, suggestions, bug fixes and questions about these routines * are welcome. Send EMail to the address given below. * * Happy bit twiddling, * * Landon Curt Noll * * chongo@toad.com * ...!{pyramid,sun,uunet}!hoptoad!chongo */ /* * AN OVERVIEW OF THE FUNCTIONS: * * This calc library contains several pseudo-random number generators: * * additive 55: * * a55rand - external interface to the additive 55 generator * sa55rand - seed the additive 55 generator * * This is a generator based on the additive 55 generator as * described in Knuth's "The Art of Computer Programming - * Seminumerical Algorithms", vol 2, 2nd edition (1981), * Section 3.2.2, page 27, Algorithm A. * * The period and other properties of this generator make it very * useful to 'seed' other generators. * * This generator is used by other other generators to produce * various internal values. In fact, the lower 64 bits of seed * given to other other generators is passed to sa55rand(). * * If you need a fast generator and do not need a crypto strong one, * you should consider using the shuffle generator instead. * * shuffle: * * shufrand - generate a pseudo-random value via a shuffle process * sshufrand - seed the shufrand generator * rand - generate a pseudo-random value over a given range * srand - seed the random stream generator * * This is a generator based on the shuffle generator as described in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B. * * The shuffle generator is fast and serves as a fairly good standard * pseudo-random generator. * * The shuffle generator is feed random values by the additive 55 * generator via a55rand(). Calling a55rand() or sa55rand() will * affect the shuffle generator output. * * The rand function is really another interface to the shuffle * generator. Unlike shufrand(), rand() can return a value of any * given size. The value returned by rand() may be confined to * either a half open [0,a) (0 <= value < a) or closed interval * [a,b] (a <= value <= b). By default, a 64 bit value is returned. * * Calling srand() simply calls sshufrand() with the same seed. * * crypto: * * cryrand - produce a cryptographically strong pseudo-random number * scryrand - seed the crypto generator * random - produce a cryptographically strong pseudo-random number * over a given range * srandom - seed random * * This generator is described in the papers: * * Blum, Blum, and Shub, "Comparison of Two Pseudorandom Number * Generators", in Chaum, D. et. al., "Advances in Cryptology: * Proceedings Crypto 82", pp. 61-79, Plenum Press, 1983. * * Blum, Blum, and Shub, "A Simple Unpredictable Pseudo-Random * Number Generator", SIAM Journal of Computing, v. 15, n. 2, * 1986, pp. 364-383. * * U. V. Vazirani and V. V. Vazirani, "Trapdoor Pseudo-Random * Number Generators with Applications to Protocol Design", * Proceedings of the 24th IEEE Symposium on the Foundations * of Computer Science, 1983, pp. 23-30. * * U. V. Vazirani and V. V. Vazirani, "Efficient and Secure * Pseudo-Random Number Generation", Proceedings of the 24th * IEEE Symposium on the Foundations of Computer Science, * 1984, pp. 458-463. * * U. V. Vazirani and V. V. Vazirani, "Efficient and Secure * Pseudo-Random Number Generation", Advances in Cryptology - * Proceedings of CRYPTO '84, Berlin: Springer-Verlag, 1985, * pp. 193-202. * * "Probabilistic Encryption", Journal of Computer & System * Sciences 28, pp. 270-299. * * We also refer to this generator as the 'Blum' generator. * * This generator is considered 'strong' in that it passes all * polynomial-time statistical tests. * * The crypto generator is not as fast as most generators, though * it is not painfully slow either. * * One may fully seed this generator via scryrand(). Calling * scryrand() with 1 or 3 arguments will result in the additive * 55 generator being seeded with the same seed. Calling * scryrand() with 4 arguments, where the first argument * is >= 0 will also result in the additive 55 generator * being seeded with the same seed. * * The random() generator is really another interface to the * crypto generator. Unlike cryrand(), random() can return a * value confined to either a half open (0 <= value < a) or closed * interval (a <= value <= b). By default, a 64 bit value is * returned. * * Calling srandom() simply calls scryrand(seed). The additive * 55 generator will be seeded with the same seed. * * As a bonus, the function 'nxtprime' is provided to produce a probable * prime number. * * All generators come already seeded with precomputed initial constants. * Thus, it is not required to seed a generator before using it. * * Using a seed of '0' will reload generators with their initial states. * In the case of scryrand(), lengths of -1 must also be supplied. * * sa55rand(0) initializes only additive 55 * sshufrand(0) initializes additive 55 and shuffle * srand(0) initializes additive 55 and shuffle * scryrand(0,-1,-1) initializes all generators * scryrand(0) initializes all generators * srandom(0) initializes all generators * randstate(0) initializes all generators * * All of the above single arg calls are fairly fast. In fact, passing * seeding with a non-zero seed, in the above cases, where seed is * not excessively large (many bits long), is also reasonably fast. * * The call: * * scryrand(-1, ip, iq, ir) * * is fast because no checking is performed on the 'ip', 'iq', or 'ir' * when seed is -1. NOTE: One must ensure that 'ip', 'iq', are valid * Blum primes and that 'ir' is a quadratic residue of their product! * * A call of scryrand(seed,len1,len2), with len1,len2 > 4, (3 args) * is a somewhat slow as the length args increase. This type of * call selects 2 primes that are used internally by the crypto * generator. A call of scryrand(seed,ip,iq,ir) where seed >= 0 * is as slow as the 3 arg case. See scryrand() for more information. * * A call of scryrand() (no args) may be used to quickly change the * internal state of the crypto and random generators. Only one major * internal crypto generator value (a quadratic residue) is randomly * selected via the additive 55 generator. The other 2 major internal * values (the 2 Blum primes) are preserved. In this form, the additive * 55 is not seeded. * * Calling scryrand(seed,[len1,len2]) (1 or 3 args), or calling * srandom(seed) will leave the additive 55 and shuffle generator in a * seeded state as if srand(seed) has been called. Calling * scryrand(seed,ip,iq,ir) (4 args), with seed>0 will also leave * the additive 55 generator in the same scryrand(seed) state. * * Calling scryrand() (no args) will not seed the additive * 55 or shuffle generator before or afterwards. The additive 55 * and shuffle generators will be changed as a side effect of that call. * * Calling srandom(seed) produces the same results as calling scryrand(seed). * * The states of all generators (additive 55, shuffle and crypto) can be * saved and restored via the randstate() function. Saving the state just * after seeding a generator and restoring it later as a very fast way * to reseed a generator. * * TRUTH IN ADVERTISING: * * The word 'probable', in reference to the nxtprime() function, is used * because of an extremely extremely small chance that a composite (a * non-prime) may be returned. In no cases will a prime be skipped. * For our purposes, this is sufficient as the chance of returning a * composite is much smaller than the chance that a hardware glitch * will cause nxtprime() to return a bogus result. * * Another "truth in advertising" issue is the use of the term * 'pseudo-random'. All deterministic generators are pseudo random. * This is opposed to true random generators based on some special * physical device. * * The crypto generator is 'pseudo-random'. There is no statistical * test, which runs in polynomial time, that can distinguish the crypto * generator from a truly random source. * * A final "truth in advertising" issue deals with how the magic numbers * found in this library were generated. Detains can be found in the * various functions, while a overview can be found in the SOURCE FOR * MAGIC NUMBERS section below. * **** * * ON THE GENERATORS: * * The additive 55 generator has a good period, and is fast. It is * reasonable as generators go, though there are better ones available. * We use it in seeding the crypto generator as its period and * other statistical properties are good enough for our purposes. * * This shuffle generator has a very good period, and is fast. It is * fairly good as generators go, and is better than the additive 55 * generator. Casual direct use of the shuffle generator may be * acceptable. Because of this, the interface to the shuffle generator, * but not the additive 55 generator, is advertised when this file is * loaded. * * The shuffle generator functions, shufrand() and rand() use the same * seed and tables. The shuffle generator shuffles the values produced * by the additive 55 generator. Calling or seeding the additive 55 * generator will affect the output of the shuffle generator. * * The crypto generator is the best generator in this package. It * produces a cryptographically strong pseudo-random bit sequence. * Internally, a fixed number of bits are generated after each * generator iteration. Any unused bits are saved for the next call * to the generator. The crypto generator is not too slow, though * seeding the generator from scratch is slow. Shortcuts and * pre-computer seeds have been provided for this reason. Use of * crypto should be more than acceptable for many applications. * * The crypto seed is in 4 parts, the additive 55 seed (lower 64 * bits of seed), the shuffle seed (all but the lower 64 bits of * seed), and two lengths. The two lengths specifies the minimum * bit size of two primes used internal to the crypto generator. * Not specifying the lengths, or using -1 will cause crypto to * use the default minimum lengths of 248 and 264 bits, respectively. * * The random() function just another interface to the crypto * generator. Like rand(), random() provides an interval capability * (closed or open) as well as a 64 bit default return value. * The random() function as good as crypto, and produces numbers * that are equally cryptographically strong. One may use the * seed functions srandom() or scryrand() for either the random() * or cryrand() functions. * * The seed for all of the generators may be of any size. Only the * lower 64 bits of seed affect the additive 55 generator. Bits * beyond the lower 64 bits affect the shuffle generators. Excessively * large values of seed will result in increased memory usage as * well as a larger seed time for the shuffle and crypto generators. * See REGARDING SEEDS below, for more information. * * One may save and restore the state of all generators via the * randstate() function. * **** * * REGARDING SEEDS: * * Because the generators are interrelated, seeding one generator * will directly or indirectly affect the other generators. Seeding * the shuffle generator seeds the additive 55 generator. Seeding * the crypto generator seeds the shuffle generator. * * The seed of '0' implies that a generator should be seeded back * to its initial default state. * * For the moment, seeds < -1 are reserved for future use. The * value of -1 is used as a special indicator to the fourth form * of scryrand(), and it not a real seed. * * A seed may be of any size. The additive 55 generator uses only * the lower 64 bits, while the shuffle generator uses bytes beyond * the lower 64 bits. * * To help make the generator produced by seed S, significantly * different from S+1, seeds are scrambled prior to use. The * function randreseed64() maps [0,2^64) into [0,2^64) in a 1-to-1 * and onto fashion. * * The purpose of the randreseed64() is not to add security. It * simply helps remove the human perception of the relationship * between the seed and the production of the generator. * * The randreseed64() process does not reduce the security of the * generators. Every seed is converted into a different unique seed. * No seed is ignored or favored. * * There is no limit on the size of a seed. On the other hand, * extremely large seeds require large tables and long seed times. * Using a seed in the range of [2^64, 2^64 * 128!) should be * sufficient for most purposes. An easy way to stay within this * range to to use seeds that are between 21 and 215 digits, or 64 to * 780 bits long. * **** * * SOURCE OF MAGIC NUMBERS: * * Most of the magic constants used on this library ultimately are * based on the Rand book of random numbers. The Rand book contains * 10^6 decimal digits, generated by a physical process. This book, * produced by the Rand corporation in the 1950's is considered * a standard against which other generators may be measured. * * The Rand book of numbers was groups into groups of 20 digits. * The first 55 groups < 2^64 were used to initialize add55_init_tbl. * The size of 20 digits was used because 2^64 is 20 digits long. * The restriction of < 2^64 was used to prevent modulus biasing. * (see the note on modulus biasing in rand()). * * The additive 55 generator during seeding is used 128 times to help * remove the initial seed state from the initial values produced. * The loop count of 128 was a power of 2 that permits each of the * 55 table entries to be processed at least twice. * * The function, randreseed64(), uses 4 primes to scramble 64 bits * into 64 bits. These primes were also extracted from the Rand * book of numbers. See sshufrand() for details. * * The default shuffle table size of 128 entries is the power of 2 * that is longer than the 100 entries recommended by Knuth for * the shuffle algorithm (see the text cited in shufrand()). * We use a power of 2 shuffle table length so that the shuffle * process can select a table entry from a new additive 55 value * by extracting its top most bits. * * The quadratic residue search performed by cryres() starts at * a value that is in the interval [2^sqrpow,n-100], where '2^sqrpow' * is the smallest power of 2 >= 'n^(3/4)' where 'n=p*q'. We also * reject any initial residue whose square (mod n) does not fit * this same restriction. Finally, we reject any residue that * is within 100 of its square (mod n). * * The use of 'n^(3/4)' insures that the quadratic residue is * large, but not too large. We want to avoid residues that are * near 0 or that are near 'n'. Such residues are trivial or * semi-trivial. Applying the same restriction to the square * of the initial residue avoid initial residues near 'sqrt(n)'. * Such residues are trivial or semi-trivial as well. * * Avoiding residues whose squares (mod n) are not within 100 of * itself helps avoid selecting residue sequences (repeated * squaring mod n) that initally do not change very much. * Such residues might be somewhat trivial, so we play it safe. * * Taking some care to select a good initial residue helps * eliminate cheep search attacks. It is true that a subsequent * residue could be one of the residues that we would initially * avoid. However such an occurance will happen after the * generator is well underway and any such information * has been lost. * * The use of '100' in the above initial residue selection is * somewhat arbitrary. It could be argued that a value as * small as 10 are sufficient. The value '100' was selected * because it is the first 3 digits of the Rand Book of Numbers. * We used 3 digits instead of 2 or 1 because '10' was too close * for comfort and '1' was clearly too small. * * Because of the initial 'n-100' upper bound part of the initial * residue selection range, the smallest Blum prime that may be * used is 19. The first 3 Blum primes 3, 7, and 11 cannot be * used. The largest value of 'n' that is a product of those * Blum primes is 121. The 'n-100' value (21) is already smaller * than the smallest power of 2 >= 'n^(3/4)'. The next Blum prime, * 19, produces the smallest value of 'n' (19*19=361) for which * one can find an initial residue that can satisfy the above. * By not considering Blum primes that are less than 5 bits long, * we avoid the smaller problem values. * * The final magic numbers: '248' and '264' are the exponents the * largest powers of 2 that are < the two default Blum primes 'p' * and 'q' used by the crypto generator. The values of '248' and * '264' implies that the product n=p*q > 2^512. Each iteration * of the crypto generator produces log2(log2(n=p*q)) random bits. * The crypto generator is the most efficient when n is slightly > * 2^(2^b). The product n > 2^(2^9)) produces 9 bits as efficiently * as possible under the crypto generator process. * * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly * improve the crypto generator. On the other hand, it can't hurt. * The two len values differ slightly to avoid factorization attacks * that work on numbers that are a perfect square, or where the two * primes are nearly the same. I elected to have the sizes differ * by 3% of the product size. The difference between '248' and * '264', is '16', which is ~3.125% of '512'. Now 3% of '512' is * '~15.36', and the next largest whole number is '16'. * * The product n=p*q > 2^512 implies a product if at least 155 digits. * A product of two primes that is at least 155 digits is slightly * beyond Number Theory and computer power of Nov 1992, though this * will likely change in the future. * * Again, the ability (or lack thereof) to factor 'n=p*q' does not * directly relate to the strength of the crypto generator. We * selected n=p*q > 2^512 mainly because '512 was a power of 2 and * only slightly because it is up in the range where it is difficult * to factor. * **** * * FOR THE PARANOID: * * The truly paranoid might suggest that my claims in the MAGIC NUMBERS * section are a lie intended to entrap people. Well they are not, but * you need not take my word for it. * * The random numbers from the Rand book of random numbers can be * verified by anyone who obtains the book. As these numbers were * created before I (Landon Curt Noll) was born (you can look up * my birth record if you want), I claim to have no possible influence * on their generation. * * There is a very slight chance that the electronic copy of the * Rand book that I was given access to differs from the printed text. * I am willing to provide access to this electronic copy should * anyone wants to compare it to the printed text. * * One might object to the complexity of the seed scramble/mapping * via the randreseed64() function. The randreseed64() function maps: * * 1 ==> 4967126403401436567 * * calling srand(1) with the randreseed64() process would be the same * as calling srand(4967126403401436567) without it. No extra security * is gained or reduced by using the randreseed64() process. The meaning * of seeds are exchanged, but not lost or favored (used by more than * one input seed). * * One could take issue with my selection of the default sizes '248' * and '264'. As far as I know, 155 digits (512 bits) is beyond the * state of the art of Number Theory and Computation as of 01 Sep 92. * It will likely be true that 155 digit products of two primes could * come within reach in the next few years, but so what? If you are * truly paranoid, why would you want to use the default seed, which * is well known? * * The paranoid today might consider using the lengths of at least '345' * and '325' will produce a product of two primes that is 202 digits. * (the 2nd and 3rd args of scryrand > 345 & >325 respectively) Factoring * 200+ digit product of two primes is well beyond the current hopes of * Number Theory and Computer power, though even this limit may be passed * someday. * * One might ask if value of '100' is too small with respect to the * initial residue selection. Showing that '100' is too small would * be difficult. Even if one could make that case, the chance that * a 'problem' initial reside would be used would be very very small * for non-trivial values of 'p' and 'q'. * * If all the above fails to pacify the truly paranoid, then one may * select by some independent means, 2 Blum primes (primes mod 4 == 3, * p < q), and a quadratic residue if p*q. Then by calling: * * scryrand(-1, p, q, r) * * and then calling cryrand() or random(), one may bypass all magic * numbers and use the pure generator. * * Note that randstate() may also be used by the truly paranoid. * Even though it holds state for the other generators, their states * are independent. * **** * * GOALS: * * The goals of this package are: * * all magic numbers are explained * * I distrust systems with constants (magic numbers) and tables * that have no justification (e.g., DES). I believe that I have * done my best to justify all of the magic numbers used. * * full documentation * * You have this source file, plus background publications, * what more could you ask? * * large selection of seeds * * Seeds are not limited to a small number of bits. A seed * may be of any size. * * the strength of the generators may be tuned to meet the application need * * By using the appropriate seed arguments, one may increase * the strength of the generator to suit the need of the * application. One does not have just a few levels. * * This calc lib file is intended for demonstration purposes. Writing * a C program (with possible assembly or libmp assist) would produce * a faster generator. * * Even though I have done my best to implement a good system, you still * must use these routines your own risk. * * Share and enjoy! :-) */ /* * These constants are used by all of the generators in various direct and * indirect forms. */ static cry_seed = 0; /* master seed */ static cry_mask = 0xffffffffffffffff; /* 64 bit work mask */ /* * Initial magic values for the additive 55 generator - used by sa55rand() * * These values were generated from the Rand book of random numbers. * In groups of 20 digits, we took the first 55 groups < 2^64. Leading * 0 digits were dropped off to avoid confusion with octal values. */ static mat add55_init_tbl[55] = { 10097325337652013586, 8422689531964509303, 9376707153831131165, 12807999708015736147, 12171768336606574717, 2051656926866574818, 3529647783580834282, 13746700781847540610, 17468509505804776974, 14385537637435099817, 14225685144642756788, 11100020401286074697, 7207317906119690446, 15474452669527079953, 16868487670307112059, 4493524947524633824, 13021248927856520106, 15956600001874392423, 1758753794041921585, 1540354560501451176, 5335129695612719255, 9973334408846123356, 2295368703230757546, 15020099946907494138, 10518216150184876938, 9188200973282539527, 4220863048338987374, 682273982071453295, 7706178130835869910, 4618975533122308420, 397583911260717646, 5686731560708285046, 10123916228549657560, 1304775865627110086, 15501295782182641134, 3061180729620744156, 6958929830512809719, 10850627469959910507, 13499063195307571839, 6410193623982098952, 4111084083850807341, 17719042079595449953, 5462692006544395659, 18288274374963224041, 8337656769629990836, 7477446061798548911, 9815931464890815877, 6913451974267278601, 11883095286301198901, 14974403441045516019, 14210337129134237821, 12883973436502761184, 4285013921797415077, 16435915296724552670, 3742838738308012451 }; /* * additive 55 tables - used by a55rand() and sa55rand() */ static add55_j = 23; /* the first walking table pointer */ static add55_k = 54; /* the second walking table pointer */ static add55_seed64 = 0; /* lower 64 bits of the reseeded seed */ static mat add55_tbl[55]; /* additive 55 state table */ /* * cryobj - cryptographic pseudo-random state object */ obj cryobj { \ p, /* first Blum prime (prime 3 mod 4) */ \ q, /* second Blum prime (prime 3 mod 4) */ \ r, /* quadratic residue of n=p*q */ \ exp, /* used in computing crypto good bits */ \ left, /* bits unused from the last cryrand() call */ \ bitcnt, /* left contains bitcnt crypto good bits */ \ a55j, /* 1st additive 55 table pointer */ \ a55k, /* 2nd additive 55 table pointer */ \ seed, /* last seed set by sa55rand() or 0 */ \ shufy, /* Y (previous a55rand output for shuffle) */ \ shufsz, /* size of the shuffle table */ \ shuftbl, /* a matrix of shufsz entries */ \ a55tbl /* additive 55 generator state table */ \ } /* * initial cryptographic pseudo-random values - used by scryrand() * * These values are what the crypto generator is initialized with * with this library first read. These values may be reproduced the * hard way by calling scryrand(0,248,264) or scryrand(0,-1,-1). * * We will build up these values a piece at a time to avoid long lines * that are difficult to send via EMail. */ /* p, first Blum prime */ static cryrand_init_p = 0x1aa9e726a735044; cryrand_init_p <<= 200; cryrand_init_p |= 0x73b7457c5297122310880fcbfa8d4e38380a00396d533300bb; /* q, second Blum prime */ static cryrand_init_q = 0xa62ee0481aa12059c3; cryrand_init_q <<= 200; cryrand_init_q |= 0x79ef44e72ff58ea70cacabbe9d264a0b16db72117d96f77e17; /* quadratic residue of n=p*q */ static cryrand_init_r = 0x3d214853f9a5bb4b12f467c9052129a9; cryrand_init_r <<= 200; cryrand_init_r |= 0xd215cc6b3c2eae8c7090591b16d8044a886b3c6a58759b1a07; cryrand_init_r <<= 200; cryrand_init_r |= 0x2b50a914b42e1b6f9703be86742837c4bc637804c2dc668c5b; /* * cryptographic pseudo-random values - used by cryrand() and scryrand() */ /* p, first Blum prime */ static cryrand_p = cryrand_init_p; /* q, second Blum prime */ static cryrand_q = cryrand_init_q; /* n = p*q */ static cryrand_n = cryrand_p*cryrand_q; /* quad residue of n */ static cryrand_r = cryrand_init_r; /* this cryrand() running exp used in computing crypto good bits */ static cryrand_exp = cryrand_r; /* good crypto bits unused from the last cryrand() call */ static cryrand_left = 0; /* the value cryrand_left contains cryrand_bitcnt crypto good bits */ static cryrand_bitcnt = 0; /* * shufrand - shuffle generator constants */ static shuf_size = 128; /* entries in shuffle table */ static shuf_shift = (64-highbit(shuf_size)); /* shift a55 value to get tbl indx */ static shuf_y; /* Y value (previous output) */ static mat shuf_tbl[shuf_size]; /* shuffle state table */ /* * products of consecutive primes - used by nxtprime() * * We compute these products now to avoid recalculating them on each call. * These values help weed out numbers that have a prime factor < 1000. */ static nxtprime_pfact10 = pfact(10); static nxtprime_pfact100 = pfact(100)/nxtprime_pfact10; static nxtprime_pfact1000 = pfact(1000)/nxtprime_pfact100; /* * a55rand - additive 55 pseudo-random number generator * * returns: * A number in the half open interval [0,2^64) * * This function implements the additive 55 pseudo-random number generator. * * This is a generator based on the additive 55 generator as described in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A. * * This function is called from the shuffle generator function shufrand(). * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you need a fast generator and do not need a crypto strong one, * you should consider using the shuffle generator (see shufrand() * and rand()). Direct use of this function is not recommended. */ define a55rand() { local ret; /* the pseudo-random number to return */ /* generate the next value using the additive 55 generator method */ ret = add55_tbl[add55_k] + add55_tbl[add55_j]; ret &= cry_mask; add55_tbl[add55_k] = ret; /* post-process out data with the seed */ ret = xor(ret, add55_seed64); /* step the pointers */ --add55_j; if (add55_j < 0) { add55_j = 54; } --add55_k; if (add55_k < 0) { add55_k = 54; } /* return the new pseudo-random number */ return(ret); } /* * sa55rand - initialize the additive 55 pseudo-random generator * * given: * seed * * returns: * old_seed * * This function seeds the additive 55 pseudo-random generator. * * This is a generator based on the additive 55 generator as described in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.2.2, page 27, Algorithm A. * * Unlike Knuth's description, we will let a seed post process our data. * * We do not apply the seed processing to the additive 55 table * data as this would disturb its pseudo-random generator properties. * Instead, we xor the output with the low 64 bits of seed. * * If seed == 0: * * This function produces values in exactly the same way * described by Knuth. * * else seed > 0: * * Each value produced is xor-ed by a 64 bit value. This value * is the result of randreseed64(rand), which produces a 64 * bit value. * * else seed == -1: * * This is a reserved seed for sshufrand(0) and srand(0). One should * not directly call srand(-1). * * else: * * Reserved for future use. * * Anyone comfortable with seed == 0 should also be comfortable with * non-zero seeds. A non-zero seeded sequence will produce a values * that have the exact same pseudo-random properties as the algorithm * described by Knuth. I.e., the sequence, while different, is as good * (or bad) as the sequence produced by a seed of 0. * * This function updates both the cry_seed and a55_seed64 global values. * * This function is called from the shuffle generator seed function sshufrand(). * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you need a fast generator and do not need a crypto strong one, * you should consider using the shuffle generator (see sshufrand() * and srand()). Direct use of this function is not recommended. */ define sa55rand(seed) { local oldseed; /* previous seed */ local junk; /* discards the first few numbers */ local j; /* firewall */ if (!isint(seed)) { quit "bad arg: arg must be an integer"; } if (seed < -1) { quit "bad arg: seed < 0 is reserved for future use"; } /* misc table setup */ oldseed = cry_seed; /* remember the previous seed */ cry_seed = seed; /* save the new seed */ if (cry_seed == -1) { /* since -1 was a special case, pretend it really was zero */ cry_seed = 0; } add55_tbl = add55_init_tbl; /* reload the table */ add55_j = 23; /* reset first walking table pointer */ add55_k = 54; /* reset second walking table pointer */ /* obtain our 64 bit xor seed */ add55_seed64 = randreseed64(cry_seed); /* spin the pseudo-random number generator a while */ if (seed == 0) { /* we will act as if sshufrand(0) or srand(0) had been called */ for (j=0; j < 513; ++j) { junk = a55rand(); } } else { for (j=0; j < 128; ++j) { junk = a55rand(); } } /* return the old seed */ return(oldseed); } /* * shufrand - implement the shuffle pseudo-random number generator * * returns: * A number in the half open interval [0,2^64) * * This function implements the shuffle number generator. * * This is a generator based on the shuffle generator as described in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B. * * The function rand() calls this function. * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you do not need a crypto strong pseudo-random generator, * this function may very well serve your needs. */ define shufrand() { local j; /* table index to replace */ /* * obtain a new random value * determine the table entry to shuffle * shuffle out the value we will return */ shuf_y = shuf_tbl[j = (shuf_y >> shuf_shift)]; /* shuffle in the new random value */ shuf_tbl[j] = a55rand(); /* return the shuffled out value */ return (shuf_y); } /* * sshufrand - seed the shuffle pseudo-random generator * * given: * a seed * * returns: * the previous seed * * This function implements the shuffle number generator. * * This is a generator based on the shuffle generator as described in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B. * * The low 64 bits of seed are used by the additive 55 generator. * See the sa55rand() function for details. The remaining bits of seed * are used to perform an initial shuffle on the shuffle state table. * The size of the seed also determines the size of the shuffle table. * * The shuffle table size is always a power of 2, and is at least 128 * entries long. Let the table size be: * * shuf_size = 2^shuf_pow * * The number of ways one could shuffle that table is: * * (2^shuf_pow)! * * We select the smallest 'shuf_pow' (and thus the size of the shuffle table) * such that the following are true: * * (2^shuf_pow)! >= (seed / 2^64) and 2^shuf_pow >= 128 * * Given that we now have the table size of 'shuf_size', we must go about * loading the table and shuffling it. * * Loading is easy, we will generate random values via the additive 55 * generator (a55rand()) and load them into successive entries. * * We enumerate the (2^shuf_pow)! values via: * * shuf_seed = 2*(U[2] + 3*(U[3] + 4*(U[4] + ... * + (U[shuf_pow-1]*(shuf_pow-1)) ... ))) * 0 <= U[j] < j * * We swap the swap table entries shuf_tbl[U[j]] & shuf_tbl[j-1] for all * 1 < j < shuf_pow. * * We will convert 'seed / 2^64' into shuf_seed, by applying the 64 bit * scramble function on 64 bit chunks of 'seed / 2^64'. * * The function srand() calls this function. * * There is no limit on the size of a seed. On the other hand, * extremely large seeds require large tables and long seed times. * Using a seed in the range of [2^64, 2^64 * 128!) should be * sufficient for most purposes. An easy way to stay within this * range to to use seeds that are between 21 and 215 digits long, or * 64 to 780 bits long. * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you do not need a crypto strong pseudo-random generator, * this function may very well serve your needs. */ define sshufrand(seed) { local shuf_pow; /* power of two - log2(shuf_size) */ local shuf_seed; /* seed bits beyond low 64 bits */ local oldseed; /* previous seed */ local shift; /* shift factor to access 64 bit chunks */ local swap_indx; /* what to swap shuf_tbl[0] with */ local rval; /* random value form additive 55 generator */ local j; /* firewall */ if (!isint(seed)) { quit "bad arg: seed must be an integer"; } if (seed < 0) { quit "bad arg: seed < 0 is reserved for future use"; } /* * seed the additive 55 generator */ if (seed == 0) { /* allow sshufrand(0) and srand(0) to arrive at the same state */ oldseed = sa55rand(-1); } else { oldseed = sa55rand(seed); } /* * form the shuffle table size and constants */ shuf_seed = seed >> 64; for (shuf_pow = 7; shuf_seed > (j=fact(1<<(shuf_pow))) && shuf_pow < 64; \ ++shuf_pow) { } shuf_size = (1 << shuf_pow); shuf_shift = 64 - shuf_pow; /* reallocate the shuffle table */ mat shuf_tbl[shuf_size]; /* * scramble the seed above the low 64 bits */ if (shuf_seed > 0) { j = 0; for (shift=0; shift < highbit(shuf_seed)+1; shift += 64) { j |= (randreseed64(shuf_seed >> shift) << shift); } shuf_seed = j; } /* * load the shuffle table */ for (j=0; j < shuf_size; ++j) { shuf_tbl[j] = a55rand(); } shuf_y = a55rand(); /* get the next Y value */ /* * We will shuffle based the process outlined in: * * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.4.2, page 139, Algorithm P. * * Here, we will let j run over the range [0,shuf_size) instead of * [shuf_size,0) as outlined in algorithm P. We will also generate * U values from shuf_seed. */ j = 0; while (shuf_seed > 0 && ++j < shuf_size) { /* determine what index we will swap with the '0' index */ quomod(shuf_seed, (j+1), shuf_seed, swap_indx); /* swap table entries, if needed */ if (swap_indx != j) { swap(shuf_tbl[j], shuf_tbl[swap_indx]); } } /* * run the generator for twice the table size */ for (j=0; j < shuf_size*2; ++j) { rval = shufrand(); } /* return the old seed */ return (oldseed); } /* * rand - generate a pseudo-random value over a given range via additive 55 * * usage: * rand() - generate a pseudo-random integer >=0 and < 2^64 * rand(a) - generate a pseudo-random integer >=0 and < a * rand(a,b) - generate a pseudo-random integer >=a and <= b * * returns: * a large pseudo-random integer over a give range (see usage) * * When no arguments are given, a pseudo-random number in the half open * interval [0,2^64) is produced. This form is identical to calling * shufrand(). * * When 1 argument is given, a pseudo-random number in the half open interval * [0,a) is produced. * * When 2 arguments are given, a pseudo-random number in the closed interval * [a,b] is produced. * * The input values a and b, if given, must be integers. * * This generator is simply a different interface to the shuffle generator. * calling shufrand(), or seeding via sshufrand() will affect the output * of this function. * * NOTE: Unlike cryrand(), this function does not preserve unused bits for * use by the next function call. * * NOTE: The Un*x rand() function returns only 16 bit or 31 bits, while we * return a number of any given size (default is 64 bits). * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you do not need a crypto strong pseudo-random generator * this function may very well serve your needs. */ define rand(a,b) { local range; /* we must generate [0,range) first */ local offset; /* what to add to get a adjusted range */ local ret; /* pseudo-random value */ local fullwords; /* number of full 64 bit chunks in ret */ local finalmask; /* mask of bits in final chunk of range */ local j; /* * setup and special cases */ /* deal with the rand() case */ if (isnull(a) && isnull(b)) { /* no args means same range as additive 55 generator */ return(a55rand()); } /* firewall - args, if given must be in integers */ if (!isint(a) || (!isnull(b) && !isint(b))) { quit "bad args: args, if given, must be integers"; } /* convert rand(x) into rand(0,x-1) */ if (isnull(b)) { /* convert call into a closed interval */ b = a-1; a = 0; /* firewall - rand(0) should act like rand(0,0) */ if (b == -1) { return(0); } } /* determine the range and offset */ if (a >= b) { /* deal with the case of rand(a,a) */ if (a == b) { /* not very random, but it is true! */ return(a); } range = a-b+1; offset = b; } else { /* convert rand(a,b), where a < b into rand(b,a) */ range = b-a+1; offset = a; } /* * At this point, we seek a pseudo-random number [0,range) to which * we will add offset to produce a number [offset,range+offset). */ fullwords = highbit(range-1)//64; finalmask = (1 << (1+(highbit(range-1)%64)))-1; /* * loop until we get a value that is in range * * A note in modulus biasing: * * We will not fall into the trap of thinking that we can simply take * a value mod 'range'. Consider the case where 'range' is '80' * and we are given pseudo-random numbers [0,100). If we took them * mod 80, then the numbers [0,20) would be produced more frequently * because the numbers [81,100) mod 80 wrap back into [0,20). */ do { /* load up all lower full 64 bit chunks with pseudo-random bits */ ret = 0; for (j=0; j < fullwords; ++j) { ret = ((ret << 64) | shufrand()); } /* load the highest chunk */ ret <<= (highbit(finalmask)+1); ret |= (shufrand() >> (64-highbit(finalmask)-1)); } while (ret >= range); /* * return the adjusted range value */ return(ret+offset); } /* * srand - seed the pseudo-random additive 55 generator * * input: * seed * * returns: * old_seed * * This function actually seeds the shuffle generator (and indirectly * the additive 55 generator used by rand() and a55rand(). * * See sshufrand() and sa55rand() for information about a seed. * * There is no limit on the size of a seed. On the other hand, * extremely large seeds require large tables and long seed times. * Using a seed in the range of [2^64, 2^64 * 128!) should be * sufficient for most purposes. An easy way to stay within this * range to to use seeds that are between 21 and 215 digits long, or * 64 to 780 bits long. * * NOTE: This is NOT Blum's method. This is used by Blum's method to * generate some internal values. * * NOTE: If you do not need a crypto strong pseudo-random generator * this function may very well serve your needs. */ define srand(seed) { if (!isint(seed)) { quit "bad arg: seed must be an integer"; } if (seed < 0) { quit "bad arg: seed < 0 is reserved for future use"; } return(sshufrand(seed)); } /* * cryrand - cryptographically strong pseudo-random number generator * * usage: * cryrand(len) * * given: * len number of pseudo-random bits to generate * * returns: * a cryptographically strong pseudo-random number of len bits * * Internally, bits are produced log2(log2(n=p*q)) at a time. If a * call to this function does not exhaust all of the collected bits, * the unused bits will be saved away and used at a later call. * Setting the seed via scryrand() or srandom() will clear out all * unused bits. Thus: * * scryrand(0); <-- restore generator to initial state * cryrand(16); <-- 16 bits * * will produce the same value as: * * scryrand(0); <-- restore generator to initial state * cryrand(4)<<12 | cryrand(12); <-- 4+12 = 16 bits * * and will produce the same value as: * * scryrand(0); <-- restore generator to initial state * cryrand(3)<<13 | cryrand(7)<<6 | cryrand(6); <-- 3+7+6 = 16 bits * * The crypto generator is not as fast as most generators, though it is not * painfully slow either. * * NOTE: This function is the Blum cryptographically strong * pseudo-random number generator. */ define cryrand(len) { local goodbits; /* the number of good bits generated each pass */ local goodmask; /* mask for the low order good bits */ local randval; /* pseudo-random value being generated */ /* * firewall */ if (!isint(len) || len < 1) { quit "bad arg: len must be an integer > 0"; } /* * Determine how many bits may be generated each pass. * * The result by Alexi et. al., says that the log2(log2(n=p*q)) * least significant bits are secure, where log2(x) is log base 2. */ goodbits = highbit(highbit(cryrand_n)); goodmask = (1 << goodbits)-1; /* * If we have bits left over from the previous call, collect * them now. */ if (cryrand_bitcnt > 0) { /* case where the left over bits are enough for this call */ if (len <= cryrand_bitcnt) { /* we need only len bits */ randval = (cryrand_left >> (cryrand_bitcnt-len)); /* save the unused bits for later use */ cryrand_left &= ((1 << (cryrand_bitcnt-len))-1); /* save away the number of bits that we will not use */ cryrand_bitcnt -= len; /* return our complete result */ return(randval); /* case where we need more than just the left over bits */ } else { /* clear out the number of left over bits */ len -= cryrand_bitcnt; cryrand_bitcnt = 0; /* collect all of the left over bits for now */ randval = cryrand_left; } /* case where we have no previously left over bits */ } else { randval = 0; } /* * Pump out len cryptographically strong pseudo-random bits, * 'goodbits' at a time using Blum's process. */ while (len >= goodbits) { /* generate the bits */ cryrand_exp = (cryrand_exp^2) % cryrand_n; randval <<= goodbits; randval |= (cryrand_exp & goodmask); /* reduce the need count */ len -= goodbits; } /* if needed, save the unused bits for later use */ if (len > 0) { /* generate the bits */ cryrand_exp = (cryrand_exp^2) % cryrand_n; randval <<= len; randval |= ((cryrand_exp&goodmask) >> (goodbits-len)); /* save away the number of bits that we will not use */ cryrand_left = cryrand_exp & ((1 << (goodbits-len))-1); cryrand_bitcnt = goodbits-len; } /* * return our pseudo-random bits */ return(randval); } /* * scryrand - seed the cryptographically strong pseudo-random number generator * * usage: * scryrand(seed) * scryrand() * scryrand(seed,len1,len2) * scryrand(seed,ip,iq,ir) * * input: * [seed pseudo-random seed * [len1 len2] minimum bit length of the Blum primes 'p' and 'q' * -1 => default lengths * [ip iq ir] Initial search values for Blum primes 'p', 'q' and * a quadratic residue 'r' * * returns: * the previous seed * * * This function will seed and setup the generator needed to produce * cryptographically strong pseudo-random numbers. See the function * a55rand() and sshufrand() for information about how 'seed' works. * * The first form of this function are fairly fast if the seed is not * excessively large. The second form is also fairly fast if the internal * primes are not too large. The third form, can take a long time to call. * (see below) The fourth form, if the 'seed' arg is not -1, can take * as long as the third form to call. If the fourth form is called with * a 'seed' arg of -1, then it is fairly fast. * * Calling scryrand() with 1 or 3 args (first and third forms), or * calling srandom(), or calling scryrand() with 4 args with the first * arg >0, will leave the shuffle generator in a seeded state as if * sshufrand(seed) has been called. * * Calling scryrand() with no args will not seed the shuffle generator, * before or afterwards, however the shuffle generator will have been * changed as a side effect of that call. * * Calling scryrand() with 4 args where the first arg is 0 or '-1' * will not change the other generators. * * * First form of call: scryrand(seed) * * The first form of this function will seed the shuffle generator * (via srand). The default precomputed constants will be used. * * * Second form of call: scryrand() * * Only a new quadratic residue of n=p*q is recomputed. The previous prime * values are kept. * * Unlike the first and second forms of this function, the shuffle * generator function is not seeded before or after the call. The * current state is used to generate a new quadratic residue of n=p*q. * * * Third form of call: scryrand(seed,len1,len2) * * In the third form, 'len1' and 'len2' guide this function in selecting * internally used prime numbers. The larger the lengths, the longer * the time this function will take. The impact on execution time of * cryrand() and random() may also be noticed, but not as much. * * If a length is '-1', then the default lengths (248 for len1, and 264 * for len2) are used. The call scryrand(0,-1,-1) recreates the initial * crypto state the slow and hard way. (use scryrand(0) or srandom(0)) * * This function can take a long time to call given reasonable values * of len1 and len2. On an R3000, the time to seed was: * * Approx value digits seed time * of len1+len2 in n=p*q in sec * ------------ -------- ------ * 8 3 0.53 * 16 5 0.54 * 32 10 0.79 * 64 20 1.17 * 128 39 2.89 * 200 61 4.68 * 256 78 7.49 * 322 100 12.47 * 464 140 35.56 * 512 155 53.57 * 664 200 83.97 * 830 250 122.93 * 996 300 242.49 * 1024 309 295.66 * 1328 400 663.44 * 1586 478 2002.10 * 1660 500 1643.45 (Faster mult/square methods kick in * 1992 600 2885.81 in certain cases. Type help config * 2048 617 1578.06 in calc for more details.) * * NOTE: The small lengths above are given for comparison * purposes and are NOT recommended for actual use. * * NOTE: Generating crypto pseudo-random numbers is MUCH * faster than seeding a crypto generator. * * NOTE: This calc lib file is intended for demonstration * purposes. Writing a C program (with possible assembly * or libmp assist) would produce a faster generator. * * * Fourth form of call: scryrand(seed,ip,iq,ir) * * In the fourth form, 'ip', 'iq' and 'ir' serve as initial search * values for the two Blum primes 'p' and 'q' and an associated * quadratic residue 'r' respectively. Unlike the 3rd form, where * lengths are given, the fourth form allows one to specify minimum * search values. * * The 'seed' value is interpreted as follows: * * If seed > 0: * * Seed and use the shuffle generator to generate 3 jump values * that are in the range '[0,ip)', '[0,iq)' and '[0,ir)' respectively. * Start searching for legal 'p', 'q' and 'r' values by adding * the jump values to their respective argument values. * * If seed == 0: * * Start searching for legal 'p', 'q' and 'r' values from * 'ip', 'iq' and 'ir' respectively. * * This form does not change/seed the other generators. * * If seed == -1: * * Let 'p' == 'ip', 'q' == 'iq' and 'r' == 'ir'. Do not check * if the value given are legal Blum primes or an associated * quadratic residue respectively. * * This form does not change/seed the other generators. * * WARNING: No checks are performed on the args passed. * Passing improper values will likely produce * poor results, or worse! * * * It should be noted that calling scryrand() while using the default * primes took only 0.04 seconds. Calling scryrand(0,-1,-1) took * 47.19 seconds. * * The paranoid, when giving explicit lengths, should keep in mind that * len1 and len2 are the largest powers of 2 that are less than the two * probable primes ('p' and 'q'). These two primes will be used * internally to cryrand(). For simplicity, we refer to len1 and len2 * as bit lengths, even though they are actually 1 less then the * minimum possible prime length. * * The actual lengths may exceed the lengths by slightly more than 3%. * Furthermore, part of the strength of this generator rests on the * difficultly to factor 'p*q'. Thus one should select 'len1' and 'len2' * (from which 'p' and 'q' are selected) such that factoring a 'len1+len2' * bit number is difficult. * * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly * improve the crypto generator. On the other hand, it can't hurt. * * There is no limit on the size of a seed. On the other hand, * extremely large seeds require large tables and long seed times. * Using a seed in the range of [2^64, 2^64 * 128!) should be * sufficient for most purposes. An easy way to stay within this * range to to use seeds that are between 21 and 215 digits long, or * 64 to 780 bits long. * * NOTE: This function will clear any internally buffer bits. See * cryrand() for details. * * NOTE: This function seeds the Blum cryptographically strong * pseudo-random number generator. */ define scryrand(seed,len1,len2,arg4) { local rval; /* a temporary pseudo-random value */ local oldseed; /* the previous seed */ local newres; /* the new quad res */ local ip; /* initial Blum prime 'p' search value */ local iq; /* initial Blum prime 'q' search value */ local ir; /* initial quadratic residue search value */ local sqir; /* square of ir mod n */ local minres; /* minimum residue allowed */ local maxres; /* maximum residue allowed */ /* * firewall - avoid bogus args and very trivial lengths */ /* catch the case of no args - compute a different quadratic residue */ if (isnull(seed) && isnull(len1) && isnull(len2)) { /* generate the next quadratic residue */ do { newres = cryres(cryrand_n); } while (newres == cryrand_r); cryrand_r = newres; cryrand_exp = cryrand_r; /* clear the internal bits */ cryrand_left = 0; cryrand_bitcnt = 0; /* return the current seed early */ return (cry_seed); } if (!isint(seed)) { quit "bad arg: seed arg (1st) must be an integer"; } if (param(0) == 4) { if (seed < -1) { quit "bad arg: with 4 args: a seed < -1 is reserved for future use"; } } else if (param(0) > 0 && seed < 0) { quit "bad arg: a seed arg (1st) < 0 is reserved for future use"; } /* * 4 arg case: select or search for 'p', 'q' and 'r' from given values */ if (param(0) == 4) { /* set initial values */ ip = len1; iq = len2; ir = arg4; /* * Unless prohibited by a seed of -1, force minimum values on * 'ip', 'iq' and 'ir'. */ if (seed >= 0) { /* * Due to the initial quadratic residue selection process, * the smallest Blum prime that is usable is 19. This * in turn implies that the smallest 'n' is 19*19 = 361. * This in turn imples that the smallest initial residue * that is possible is 128 (for the value 'n = 23*23 = 529). */ if (!isint(ip) || ip < 19) { ip = 19; } if (!isint(iq) || iq < 19) { iq = 19; } if (!isint(ir) || ir < 128) { ir = 128; } } /* * Determine our prime search points * * Unless we have a seed <= 0, we will add a random 64 bit * value to the initial search points. */ if (seed > 0) { /* add in a random value */ oldseed = srand(seed); ip += rand(ip); iq += rand(iq); } /* * force ip <= iq */ if (ip > iq) { swap(ip, iq); } /* * find the first Blum prime 'p' */ if (seed >= 0) { cryrand_p = nxtprime(ip,3,4); } else { /* seed == -1: assume 'ip' is a Blum prime */ cryrand_p = ip; } /* * find the 2nd Blum prime 'q' > 'p', if needed */ if (seed >= 0) { if (iq <= cryrand_p) { iq = cryrand_p + 2; } cryrand_q = nxtprime(iq,3,4); } else { /* seed == -1: assume 'iq' is a Blum prime */ cryrand_q = iq; } /* remember our p*q Blum prime product as well */ cryrand_n = cryrand_p*cryrand_q; /* * search for a quadratic residue */ if (seed >= 0) { /* * add in a random value to 'ir' if seeded * * Unless we have a seed <= 0, we will add a random 64 bit * value to the initial search point. */ if (seed > 0) { ir += rand(ir); } /* * increment until we find a quadratic value * * p is prime and J(r,p) == 1 ==> r is a quadratic residue of p * q is prime and J(q,p) == 1 ==> r is a quadratic residue of q * * J(r,p)*J(r,q) == J(r,p*q) * J(r,p) and J(q,p) == 1 ==> J(r,p*q) == 1 * J(r,p*q) == 1 ==> r is a quadratic residue of p*q * * We could simply compute the square of a value mod n like * we do in cryres(), but here we want to climb a little higher * than the ir value given. We will start sequentially searching * values starting at the initial search point 'ir', while at * same time confining our search to the interval [2^sqrpow,n-100], * where 2^sqrpow is the smallest power of 2 >= n^(3/4). This * range helps us avoid selecting trivial residues. * * We will also reject any quadratic residue whose square mod n * is outside of the [2^sqrpow,n-100] range, or whose square mod n * is within 100 of itself. */ minres = 2^(highbit(floor(power(cryrand_n,0.75)))+1); maxres = cryrand_n - 100; --ir; do { /* consdier the next residue that is in the allowed range */ ++ir; if (ir < minres || ir > maxres) { ir = minres; } sqir = pmod(ir, 2, cryrand_n); } while (jacobi(ir,cryrand_p) != 1 || \ jacobi(ir,cryrand_q) != 1 || \ sqir < minres || sqir > maxres || \ abs(sqir-ir) <= 100); } cryrand_r = ir; /* * clear the previously unused cryrand bits & other misc setup */ cryrand_left = 0; cryrand_bitcnt = 0; cryrand_exp = cryrand_r; /* * reseed the generator, if needed * * The crypto generator no longer needs the additive 55 and shuffle * generators. We however, restore the additive 55 and shuffle * generators back to its seeded state in order to be sure that it * will be left in the same state. * * This will make more reproducible, calls to the additive 55 and * shuffle generator; or more reproducible, calls to this function * without args. */ if (seed > 0) { ir = srand(seed); /* ignore this return value */ return(oldseed); } else { /* no seed change, return the current seed */ return (cry_seed); } } /* * If not the 4 arg case: * * convert explicit -1 args into default values * convert missing args into -1 values (use precomputed tables) */ if ((isint(len1) && len1 != -1 && len1 < 5) || (isint(len2) && len2 != -1 && len2 < 5) || (!isint(len1) && isint(len2)) || (isint(len1) && !isint(len2))) { quit "bad args: 2 & 3: if 2nd is given, must be -1 or ints > 4"; } if (isint(len1) && len1 == -1) { len1 = 248; /* default len1 value */ } if (isint(len2) && len2 == -1) { len2 = 264; /* default len2 value */ } if (!isint(len1) && !isint(len2)) { /* from here down, -1 means use precomputed values */ len1 = -1; len2 = -1; } /* * force len1 <= len2 */ if (len1 > len2) { swap(len1, len2); } /* * seed the generator */ oldseed = srand(seed); /* * generate p and q Blum primes * * The Blum process requires the primes to be of the form 3 mod 4. * We also generate n=p*q for future reference. * * We make sure that the lengths are the minimum lengths possible. * We want some range to select a random prime from, so we * go at least 3 bits higher, and as much as 3% plus 3 bits * higher. Since the section is a random, how high really * does not matter that much, but we want to avoid going to * an extreme to keep the execution time from getting too long. * * Finally, we generate a quadratic residue of n=p*q. */ if (len1 > 0) { /* generate a pseudo-random prime ~len1 bits long */ rval = rand(2^(len1-1), 2^((int(len1*1.03))+3)); cryrand_p = nxtprime(rval,3,4); } else { /* use precomputed 'p' value */ cryrand_p = cryrand_init_p; } if (len2 > 0) { /* generate a pseudo-random prime ~len1 bits long */ rval = rand(2^(len2-1), 2^((int(len2*1.03))+3)); cryrand_q = nxtprime(rval,3,4); } else { /* use precomputed 'q' value */ cryrand_q = cryrand_init_q; } /* * find the quadratic residue */ cryrand_n = cryrand_p*cryrand_q; if (len1 == 248 && len2 == 264 && seed == 0) { cryrand_r = cryrand_init_r; } else { cryrand_r = cryres(cryrand_n); } /* * clear the previously unused cryrand bits & other misc setup */ cryrand_left = 0; cryrand_bitcnt = 0; cryrand_exp = cryrand_r; /* * reseed the generator * * The crypto generator no longer needs the additive 55 and shuffle * generators. We however, restore the additive 55 and shuffle * generators back to its seeded state in order to be sure that it * will be left in the same state. * * This will make more reproducible, calls to the additive 55 and * shuffle generator; or more reproducible, calls to this function * without args. */ /* we do not care about this old seed */ rval = srand(seed); /* * return the old seed */ return(oldseed); } /* * random - a cryptographically strong pseudo-random number generator * * usage: * random() - generate a pseudo-random integer >=0 and < 2^64 * random(a) - generate a pseudo-random integer >=0 and < a * random(a,b) - generate a pseudo-random integer >=a and <= b * * returns: * a large cryptographically strong pseudo-random number (see usage) * * This function is just another interface to the crypto generator. * (see the cryrand() function). * * When no arguments are given, a pseudo-random number in the half open * interval [0,2^64) is produced. This form is identical to calling * cryrand(64). * * When 1 argument is given, a pseudo-random number in the half open interval * [0,a) is produced. * * When 2 arguments are given, a pseudo-random number in the closed interval * [a,b] is produced. * * This generator uses the crypto to return a large pseudo-random number. * * The input values a and b, if given, must be integers. * * Internally, bits are produced log2(log2(n=p*q)) at a time. If a * call to this function does not exhaust all of the collected bits, * the unused bits will be saved away and used at a later call. * Setting the seed via scryrand(), srandom() or cryrand(len,1) * will clear out all unused bits. * * NOTE: The BSD random() function returns only 31 bits, while we return 64. * * NOTE: This function is the Blum cryptographically strong * pseudo-random number generator. */ define random(a,b) { local range; /* we must generate [0,range) first */ local offset; /* what to add to get a adjusted range */ local rangebits; /* the number of bits in range */ local ret; /* pseudo-random bit value */ /* * setup and special cases */ /* deal with the rand() case */ if (isnull(a) && isnull(b)) { /* no args means return 64 bits */ return(cryrand(64)); } /* firewall - args, if given must be in integers */ if (!isint(a) || (!isnull(b) && !isint(b))) { quit "bad args: args, if given, must be integers"; } /* convert rand(x) into rand(0,x-1) */ if (isnull(b)) { /* convert call into a closed interval */ b = a-1; a = 0; /* firewall - rand(0) should act like rand(0,0) */ if (b == -1) { return(0); } } /* determine the range and offset */ if (a >= b) { /* deal with the case of rand(a,a) */ if (a == b) { /* not very random, but it is true! */ return(a); } range = a-b+1; offset = b; } else { /* convert random(a,b), where a= range and 2^(rangebits-1) < range. We * will ignore any results that are > the range that we want. * * A note in modulus biasing: * * We will not fall into the trap of thinking that we can simply take * a value mod 'range'. Consider the case where 'range' is '80' * and we are given pseudo-random numbers [0,100). If we took them * mod 80, then the numbers [0,20) would be produced more often * because the numbers [81,100) mod 80 wrap back into [0,20). */ do { /* obtain a pseudo-random value */ ret = cryrand(rangebits); } while (ret >= range); /* * return the adjusted range value */ return(ret+offset); } /* * srandom - seed the cryptographically strong pseudo-random number generator * * given: * seed a random number seed * * returns: * the previous seed * * This function is just another interface to the crypto generator. * (see the scryrand() function). * * This function makes indirect use of the additive 55 and shuffle * generator. * * There is no limit on the size of a seed. On the other hand, * extremely large seeds require large tables and long seed times. * Using a seed in the range of [2^64, 2^64 * 128!) should be * sufficient for most purposes. An easy way to stay within this * range to to use seeds that are between 21 and 215 digits long, or * 64 to 780 bits long. * * NOTE: Calling this function will clear any internally buffer bits. * See cryrand() for details. * * NOTE: This function seeds the Blum cryptographically strong * pseudo-random number generator. */ define srandom(seed) { if (!isint(seed)) { quit "bad arg: seed must be an integer"; } if (seed < 0) { quit "bad arg: seed < 0 is reserved for future use"; } return(scryrand(seed)); } /* * randstate - set/get the state of all of the generators * * usage: * randstate() return the current state * randstate(0) return the previous state, set the default state * randstate(cobj) return the previous state, set a new state * * In the first form: randstate() * * This function returns an cryobj object containing information * about the current state of all of the generators. * * In the second form: randstate(0) * * This function sets all of the generators to the default initial * state (i.e., the state when this library was loaded). * * This function returns an cryobj object containing information * about the previous state of all of the generators. * * In the third form: randstate(cobj) * * This function sets all of the generators to the state as found * in the cryobj object. * * This function returns an cryobj object containing information * about the previous state of all of the generators. * * This function may be used to save and restore cryrand() & random() * generator states. For example: * * state = randstate() <-- save the current state * random() <-- print the next 64 bits * randstate(state) <-- restore previous state * random() <-- print the same 64 bits * * One may quickly reseed a generator. For example: * * srandom(1,330,350) <-- seed the generator * seed1state = randstate() <-- remember this 1st seeded state * random() <-- print 1st 64 bits seed 1 generator * srandom(2,331,351) <-- seed the generator again * seed2state = randstate() <-- remember this 2nd seeded state * random() <-- print 1st 64 bits seed 2 generator * * randstate(seed1state) <-- reseed to the 1st seeded state * random() <-- reprint 1st 64 bits seed 1 generator * randstate(seed2state) <-- reseed to the 2nd seeded state * random() <-- reprint 1st 64 bits seed 2 generator * * oldstate = randstate(0) <-- seed to the default generator * random() <-- print 1st 64 bits from default * a55rand() <-- print 1st 64 bits a55 generator * prevstate = randstate(oldstate) <-- restore seed 2 generator * random() <-- print 2nd 64 bits seed 2 generator * randstate(prevstate) <-- restore def generator in progress * random() <-- print 2nd 64 bits default generator * a55rand() <-- print 2nd 64 bits a55 generator * * given: * cobj if a cryobj object, use that object to set the current state * if 0, set to the default state * * return: * return the state of the crypto pseudo-random number generator in * the form of an cryobj object, as it was prior to this call * * NOTE: No checking is performed on the data the 3rd form (cryobj object * arg) is used. The user must ensure that the arg represents a valid * generator state. * * NOTE: When using the second form (passing an integer arg), only 0 is * defined. All other integer values are reserved for future use. */ define randstate(arg) { /* declare our objects */ local obj cryobj x; /* firewall comparator */ local obj cryobj prev; /* previous states of the generators */ local junk; /* dummy holder of random values */ /* firewall */ if (!isint(arg) && !istype(arg,x) && !isnull(arg)) { quit "bad arg: argument must be integer, an cryobj object or missing"; } if (isint(arg) && arg != 0) { quit "bad arg: non-zero integer arguments are reserved for future use"; } /* * save the current state */ prev.p = cryrand_p; prev.q = cryrand_q; prev.r = cryrand_r; prev.exp = cryrand_exp; prev.left = cryrand_left; prev.bitcnt = cryrand_bitcnt; prev.a55j = add55_j; prev.a55k = add55_k; prev.seed = cry_seed; prev.shufy = shuf_y; prev.shufsz = shuf_size; prev.shuftbl = shuf_tbl; prev.a55tbl = add55_tbl; if (isnull(x)) { /* if no args, just return current state */ return (prev); } /* * deal with the cryobj arg - set the state */ if (istype(arg, x)) { /* set the state from this object */ cryrand_p = arg.p; cryrand_q = arg.q; cryrand_n = cryrand_p*cryrand_q; cryrand_r = arg.r; cryrand_exp = arg.exp; cryrand_left = arg.left; cryrand_bitcnt = arg.bitcnt; add55_j = arg.a55j; add55_k = arg.a55k; cry_seed = arg.seed; add55_seed64 = randreseed64(cry_seed); shuf_y = arg.shufy; shuf_size = arg.shufsz; shuf_shift = (64-highbit(shuf_size)); shuf_tbl = arg.shuftbl; add55_tbl = arg.a55tbl; /* * deal with the 0 integer arg - set the default initial state */ } else if (isint(arg) && arg == 0) { cryrand_p = cryrand_init_p; cryrand_q = cryrand_init_q; cryrand_n = cryrand_p * cryrand_q; cryrand_r = cryrand_init_r; cryrand_exp = cryrand_r; cryrand_left = 0; cryrand_bitcnt = 0; cry_seed = 0; cry_seed = sshufrand(0); } /* * return the previous state */ return (prev); } /* * nxtprime - find a probable prime >= n_arg * * usage: * nxtprime(n_arg) * nxtprime(n_arg, modval, modulus) * * given: * n_arg lower bound of the search * [modval modulus] if given, look for numbers mod modulus == modval * * returns: * A number is that is very likely prime. * * In the first case 'nxtprime(n_arg)', this function returns a probable * prime >= n_arg. In the second case 'nxtprime(n_arg, v, u)', this * function returns a probable prime >= n_arg that is also == v mod u. * * This function will not skip over a prime, through there is a * extremely unlikely chance that it will return a composite. * The odds that a number returned by this function is not prime * are 1 in 4^50. The failure rate of this function is many orders * or magnitude lower than the failure rate due to a hardware error. * * NOTE: This function can take a long time, given a large value of n_arg. */ define nxtprime(n_arg, modval, modulus) { local modgcd; /* gcd(modulus,modval) */ local n; /* value >= n_arg that is being tested */ local j; /* * firewall */ if (!isint(n_arg)) { quit "bad args: 1st arg must be an integer"; } if (!isnull(modulus) && !isint(modval)) { quit "bad args: 3rd arg, if 2nd arg is given, must be an integer"; } if (!isnull(modulus) && (!isint(modulus) || modulus <= 0)) { quit "bad args: 3nd arg, if given, must be an integer > 0"; } /* * get values < 3 out of the way */ n = n_arg; if (n < 3) { /* get the even prime out of the way, if possible */ if (isnull(modulus) || modulus == 1 || (n%modulus == modval%modulus)) { /* * 2 is the greatest odd prime, because * 2 is the least even prime :-) */ return(2); } /* we have eliminated everything < 3 */ n = 3; } /* * convert nxtprime(n) to nxtprime(n,1,2) * convert nxtprime(n,x,1) to nxtprime(n,1,2) * convert nxtprime(n,a,b) to nxtprime(n,a mod b,b) */ if (isnull(modulus) || modulus < 2) { modulus = 2; modval = 1; } modval %= modulus; /* * catch cases where no more primes == 'modval' mod 'modulus' exist */ modgcd = gcd(modval,modulus); if (modgcd > 1) { /* if beyond the modgcd, then no primes can exist */ if (n > modgcd) { print "n_arg:",n_arg," modval:",modval," modulus:",modulus; quit "no such prime of that form exists"; } /* else n <= modgcd, then our only chance is if modgcd is prime */ /* reject if modgcd has an obvious prime factor */ if (modgcd > 10 && gcd(modgcd,nxtprime_pfact10) != 1) { print "n_arg:",n_arg," modval:",modval," modulus:",modulus; quit "no such prime of that form exists"; } if (modgcd > 100 && gcd(modgcd,nxtprime_pfact100) != 1) { print "n_arg:",n_arg," modval:",modval," modulus:",modulus; quit "no such prime of that form exists"; } if (modgcd > 1000 && gcd(modgcd,nxtprime_pfact1000) != 1) { print "n_arg:",n_arg," modval:",modval," modulus:",modulus; quit "no such prime of that form exists"; } /* do 50 probable prime tests, for a 1 in 4^50 false prime rate */ if (!ptest(modgcd,50)) { print "n_arg:",n_arg," modval:",modval," modulus:",modulus; quit "no such prime of that form exists"; } /* modgcd is the only prime >= n */ return(modgcd); } /* * bump n up to the next possible case * * n will be an odd number == 'modval' mod 'modulus' */ if (n%modulus != modval) { j = n - (n%modulus) + modval; if (j < n) { n = j+modulus; } else { n = j; } } if (n%2 == 0) { n += modulus; } /* look for a prime */ n = n-modulus; do { /* try the next integer */ n = n+modulus; /* reject if it has an obvious prime factor */ if (n > 10 && gcd(n,nxtprime_pfact10) != 1) { continue; } if (n > 100 && gcd(n,nxtprime_pfact100) != 1) { continue; } if (n > 1000 && gcd(n,nxtprime_pfact1000) != 1) { continue; } /* do 50 probable prime tests */ if (!ptest(n,50)) { continue; } /* n is very likely a prime number */ return(n); } while (1); } /* * cryobj - how to initialize a cryobj object * * given: * p first Blum prime (prime 3 mod 4) * q second Blum prime (prime 3 mod 4) * r quadratic residue of n=p*q * exp used in computing crypto good bits * left bits unused from the last cryrand() call * bitcnt left contains bitcnt crypto good bits * a55j 1st additive 55 table pointer * a55k 2nd additive 55 table pointer * seed last seed set by sa55rand() or 0 * shufy Y (previous a55rand() output for shuffle) * shufsz size of the shuffle table * shuftbl a matrix of shufsz entries * a55tbl additive 55 generator state table * * return: * an cryobj object * * NOTE: This function, by convention, returns an cryobj object. */ define cryobj(p,q,r,exp,left,bitcnt,a55j,a55k,seed,shufy,shufsz,shuftbl,a55tbl) { /* declare our objects */ local obj cryobj x; /* firewall */ if (!isint(p) || !isint(q) || !isint(r) || !isint(exp) || \ !isint(left) || !isint(bitcnt) || !isint(a55j) || \ !isint(a55k) || !isint(seed) || !isint(shufy) || !isint(shufsz)) { quit "bad args: first 11 args must be integers"; } if (!ismat(shuftbl) || matdim(shuftbl) != 1 || \ matmin(shuftbl,1) != 0 || matmax(shuftbl,1) != shuf_size-1) { quit "bad arg: 12th is not a mat[0:shuf_size-1]"; } if (!ismat(a55tbl) || matdim(a55tbl) != 1 || \ matmin(a55tbl,1) != 0 || matmax(a55tbl,1) != 54) { quit "bad arg: 13th is not a mat[0:54]"; } /* initialize object with default startup values */ x.p = p; x.q = q; x.r = r; x.exp = exp; x.left = left; x.bitcnt = bitcnt; x.a55j = a55j; x.a55k = a55k; x.seed = seed; x.shufy = shuf_y; x.shufsz = shuf_size; x.shuftbl = shuf_tbl; x.a55tbl = a55tbl; /* return the initialized object */ return (x); } /* * cmpobj - compare two cryrand objects * * usage: * a an cryobj object * b an cryobj object * * NOTE: This function is intended for debug purposes. */ define cmpobj(a,b) { local obj cryobj x; /* firewall comparator */ local shufsiz; local i; /* firewall */ if (!istype(a, x)) { quit "bad arg: 1st arg is not an cryobj object"; } if (!istype(b, x)) { quit "bad arg: 2nd arg is not an cryobj object"; } if (!ismat(a.shuftbl) || matdim(a.shuftbl) != 1 || \ matmin(a.shuftbl,1) != 0 || matmax(a.shuftbl,1) != a.shufsz-1) { quit "bad arg: 1st arg is not a mat[0:shuf_size-1]"; } if (!ismat(b.shuftbl) || matdim(b.shuftbl) != 1 || \ matmin(b.shuftbl,1) != 0 || matmax(b.shuftbl,1) != b.shufsz-1) { quit "bad arg: 2nd arg is not a mat[0:shuf_size-1]"; } if (!ismat(a.a55tbl) || matdim(a.a55tbl) != 1 || \ matmin(a.a55tbl,1) != 0 || matmax(a.a55tbl,1) != 54) { quit "bad arg: 1st arg is not a mat[0:54]"; } if (!ismat(b.a55tbl) || matdim(b.a55tbl) != 1 || \ matmin(b.a55tbl,1) != 0 || matmax(b.a55tbl,1) != 54) { quit "bad arg: 2nd arg is not a mat[0:54]"; } /* compare values */ if (a.p != b.p) { print "a.p - b.p:", a.p - b.p; } if (a.q != b.q) { print "a.q - b.q:", a.q - b.q; } if (a.r != b.r) { print "a.r - b.r:", a.r - b.r; } if (a.exp != b.exp) { print "a.exp - b.exp:", a.exp - b.exp; } if (a.left != b.left) { print "a.left - b.left:", a.left - b.left; } if (a.bitcnt != b.bitcnt) { print "a.bitcnt - b.bitcnt:", a.bitcnt - b.bitcnt; } if (a.a55j != b.a55j) { print "a.a55j - b.a55j:", a.a55j - b.a55j; } if (a.a55k != b.a55k) { print "a.a55k - b.a55j:", a.a55k - b.a55k; } if (a.seed != b.seed) { print "a.seed - b.seed:", a.seed - b.seed; } if (a.shufy != b.shufy) { print "a.shufy - b.shufy:", a.shufy - b.shufy; } if (a.shufsz != b.shufsz) { print "a.shufsz - b.shufsz:", a.shufsz - b.shufsz; shufsiz = min(a.shufsz, b.shufsz); } else { shufsiz = a.shufsz; } for (i=0; i < shufsiz; ++i) { if (a.shuftbl[i] != b.shuftbl[i]) { print "a.shuftbl[" : i : "] - b.shuftbl[" : i : "]:", \ a.shuftbl[i] - b.shuftbl[i]; } } if (a.shufsz > shufsiz) { print " skipping a.shuftbl[" : shufsiz : ".." : a.shufsz-1 : "]"; } else if (b.shufsz > shufsiz) { print " skipping b.shuftbl[" : shufsiz : ".." : b.shufsz-1 : "]"; } for (i=0; i < 55; ++i) { if (a.a55tbl[i] != b.a55tbl[i]) { print "a.a55tbl[" : i : "] - b.a55tbl[" : i : "]:", \ a.a55tbl[i] - b.a55tbl[i]; } } } /* * cryobj_print - print the value of a cryobj object * * usage: * a an cryobj object * * NOTE: This function is called automatically when an cryobj object * is displayed. */ define cryobj_print(a) { /* declare our objects */ local obj cryobj x; /* firewall comparator */ /* firewall */ if (!istype(a, x)) { quit "bad arg: arg is not an cryobj object"; } if (!ismat(a.shuftbl) || matdim(a.shuftbl) != 1 || \ matmin(a.shuftbl,1) != 0 || matmax(a.shuftbl,1) != a.shufsz-1) { quit "bad arg: arg is not a mat[0:shuf_size-1]"; } if (!ismat(a.a55tbl) || matdim(a.a55tbl) != 1 || \ matmin(a.a55tbl,1) != 0 || matmax(a.a55tbl,1) != 54) { quit "bad arg: arg is not a mat[0:54]"; } /* print the value */ print "cryobj(" : a.p : "," : a.q : "," : a.r : "," : a.exp : "," : \ a.left : "," : a.bitcnt : "," : a.a55j : "," : a.a55k : "," : \ a.seed : "," : a.shufy : "," : a.shufsz : \ ",[" : a.shuftbl[0] : "," : a.shuftbl[1] : "," : \ a.shuftbl[2] : ",...," : a.shuftbl[52] : "," : \ a.shuftbl[53] : "," : a.shuftbl[54] : "]" : \ ",[" : a.a55tbl[0] : "," : a.a55tbl[1] : "," : \ a.a55tbl[2] : ",...," : a.a55tbl[52] : "," : \ a.a55tbl[53] : "," : a.a55tbl[54] : "])" : ; } /* * cryres - find a pseudo-random quadratic residue for scryrand() and cryrand() * * given: * n product of two Blum primes * * returns: * a number that is a quadratic residue of n=p*q * * This function is returns the pseudo-random quadratic residue of * the product of two primes. Normally this function is called * only by the crypto generator. * * NOTE: No check is made to ensure that the values passed are prime. */ define cryres(n) { local quadres; /* quadratic residue of n */ local sqquadres; /* square of quadres mod n */ local minres; /* minimum residue allowed */ local maxres; /* maximum residue allowed */ local j; /* * firewall */ if (!isint(n)) { quit "bad arg: must an integer"; } /* * find a pseudo-random quadratic residue of n = p*q * * We will start sequentially searching for quadratic residue * values starting at the initial search point 'ir', while at * same time confining our search to the interval [2^sqrpow,n-100], * where 2^sqrpow is the smallest power of 2 >= n^(3/4). This * range helps us avoid selecting trivial residues. * * We will also reject any quadratic residue whose square mod n * is outside of the [2^sqrpow,n-100] range, or whose square mod n * is within 100 of itself. */ minres = 2^(highbit(floor(power(n,0.75)))+1); maxres = n - 100; do { /* form a quadratic residue */ quadres = pmod(rand(minres,maxres+1), 2, n); sqquadres = pmod(quadres, 2, n); } while (quadres < minres || quadres > maxres || \ sqquadres < minres || sqquadres > maxres || \ abs(sqquadres-quadres) <= 100); /* * return the quadratic residue of n */ return (quadres); } /* * randreseed64 - scramble a 64 bit seed * * given: * a 64 bit seed * * returns: * a 64 scrambled seed, or 0 if seed was 0 * * It is 'nice' when a seed of "n" produces a 'significantly different' * sequence than a seed of "n+1". Generators, by convention, assign * special significance to the seed of '0'. It is an unfortunate that * people often pick small seed values, particularly when large seed * are of significance to the generators found in this file. * * If 'seed' is 0, then 0 is returned. If 'seed' is non-zero, we will * produce a different and unique new scrambled 'seed'. This scrambling * will effectively eliminate the human factors and perceptions that * are noted above. * * It should be noted that the purpose of this process to scramble a seed * ONLY. We do not care if these generators produce good random numbers. * We only want to help eliminate the human factors and perceptions * noted above. * * This function scrambles the low 64 bits of a seed, by mapping [0,2^64) * into [0,2^64). This map is one-to-one and onto. Mapping is performed * using a linear congruence generator of the form: * * X1 <-- (a*X0 + c) mod m * * The generator are based on the linear congruential generators found in * Knuth's "The Art of Computer Programming - Seminumerical Algorithms", * vol 2, 2nd edition (1981), Section 3.6, pages 170-171. * * Because we process 64 bits we will take: * * m = 2^64 (based on note ii) * * We will scan the Rand book of numbers to look for an 'a' such that: * * a mod 8 == 5 (based on note iii) * 0.01*m < a < 0.99*m (based on note iv) * 0.01*2^64 < a < 0.99*2^64 * * To help keep the generators independent, we want: * * a is prime * * The choice of an adder 'c' is considered immaterial according (based * in note v). Knuth suggests 'c==1' or 'c==a'. We elect to select 'c' * using the same process as we used to select 'a'. The choice is * 'immaterial' after all, and as long as: * * gcd(c, m) == 1 (based on note v) * gcd(c, 2^64) == 1 * * the concerns are met. It can be shown that if we have: * * gcd(a, c) == 1 * * then the adders and multipliers will be more independent. * * We will obtain the values 'a' and 'c for our generator from the * Rand book of numbers. Because m=2^64 is 20 decimal digits long, we * will search the Rand book of numbers 20 at a time. We will skip any * of the 55 values that were used to initialize the additive 55 generators. * The values obtained from the Rand book are: * * a = 6316878969928993981 * c = 1363042948800878693 * * As we stated before, we must map 0 ==> 0. The linear congruence * generator would normally map as follows: * * 0 ==> 1363042948800878693 (0 ==> c) * * To determine which value maps back into 0, we compute: * * (-c*minv(a,m)) % m * * and thus we find that the congruence generator would also normally map: * * 10239951819489363767 ==> 0 * * To overcome this, and preserve the 1-to-1 and onto map, we force: * * 0 ==> 0 * 10239951819489363767 ==> 1363042948800878693 * * To repeat, this function converts a values into a seed value. With the * except of 'seed == 0', every value is mapped into a unique seed value. * This mapping need not be complex, random or secure. All we attempt * to do here is to allow humans who pick small or successive seed values * to obtain reasonably different sequences from the generators below. * * NOTE: This is NOT a pseudo random number generator. This function is * intended to be used internally by sa55rand() and sshufrand(). */ define randreseed64(seed) { local ret; /* return value */ local a; /* generator 0 multiplier */ local c; /* generator 0 adder */ /* firewall */ if (!isint(seed)) { quit "bad args: seed must be an integer"; } if (seed < 0) { quit "bad arg: seed < 0 is reserved for future use"; } /* if seed is 0, we will return 0 */ if (seed == 0) { return (0); } /* * process the low 64 bits of the seed */ a = 6316878969928993981; c = 1363042948800878693; ret = (((seed & cry_mask)*a + c) & cry_mask); /* * Normally, the above equation would map: * * f(0) == 1363042948800878693 * f(10239951819489363767) == 0 * * However, we have already forced f(0) == 0. To preserve the * 1-to-1 and onto map property, we force: * * f(10239951819489363767) ==> 1363042948800878693 */ if (ret == 0) { ret = c; } /* return the scrambled value */ return (ret); } /* * Initial read execution code */ cry_seed = srandom(0); /* pre-initialize the tables */ cryrand_init_r = cryrand_r; global cryrand_ver = "10.7 23:47:35 13-Mar-1994"; global lib_debug; if (lib_debug >= 0) { print "cryrand_ver:", cryrand_ver; print "shufrand() defined"; print "sshufrand(seed) defined"; print "rand([a, [b]]) defined"; print "srand(seed) defined"; print "cryrand([a, [b]]) defined"; print "scryrand([seed, [len1, len2]]) defined"; print "scryrand(seed, ip, iq, ir) defined"; print "random([a, [b]]) defined"; print "srandom(seed) defined"; print "obj cryobj defined"; print "randstate([cryobj | 0]) defined"; print "nxtprime(n, [val, modulus]) defined"; }