+/*
+ * 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<b, into random(b,a) */
+ range = b-a+1;
+ offset = a;
+ }
+ rangebits = highbit(range-1)+1;
+ /*
+ * At this point, we seek a pseudo-random number [0,range) to which
+ * we will add offset to produce a number [offset,range+offset).
+ */
+
+ /*
+ * loop until we get a value that is in range
+ *
+ * We will obtain pseudo-random values over the range [0,2^rangebits)
+ * where 2^rangebits >= 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";
+}
projects / unix-history / commitdiff