4.4BSD/usr/src/contrib/calc-1.26.4/lib/cryrand.cal
/*
* cryrand - cryptographically strong pseudo-random number generator library
*/
/*
* Copyright (c) 1993 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 1.24.7 (We refer to calc, as in the C
* program, not the Emacs subsystem).
*
* Calc is available by anonymous ftp from ftp.uu.net (137.39.1.9)
* 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().
*
* 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.
*
* "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. This means that there
* is no statistical test, which runs in polynomial time, that
* can distinguish the crypto generator from a truly random source.
*
* 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.
*
* 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)
* sshufrand(0)
* srand(0)
* scryrand(0,-1,-1)
* srandom(0)
* scryrand(-1,ip,iq,ir)
*
* All of the above 'seed 0' 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 4 arg scryrand() call is fairly fast because no checking is
* performed on the 'ip', 'iq', or 'ir' values.
*
* 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() with 1 or 3 args, or calling srandom() will leave
* the additive 55 generator in a seeded state as if sa55rand(seed) has
* been called. Calling scryrand() with 4 args, where the first arg is
* >0 will also leave the additive 55 generator in the same scryrand(seed)
* state. Calling scryrand() with no args will not seed the additive
* 55 generator before or afterwards, however the additive 55 and shuffle
* generators will be changed as a side effect of that call.
*
* 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.
*
* As a bonus, the function 'nxtprime' is provided to produce a probable
* prime number.
*
* 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-3), where 2^sqrpow
* is the first power of 2 > n^(3/4) and n = p*q. We look at 1009
* random numbers in this interval for a quadratic residue. If
* none are found, sqrpow is decremented by 1 if sqrpow > 2.
* In practice, decrementing sqrpow never happens.
*
* 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 the square root of 'n', and are nearly 'n' itself (hence
* the n-3 upper limit) as such residues are trivial or semi-trivial.
*
* The '1009' trials is simply an attempt to look for a while before
* picking a new range. The number '1009' comes from the first 4
* digits of the rand book of numbers.
*
* Keeping sqrpow > 2 means that we do not look for quadratic
* residues that are below 4. This is in keeping with the
* n-3 in the [2^sqrpow,n-3) interval.
*
* 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.
*
* 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.
*/
cry_seed = 0; /* master seed */
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.
*/
mat add55_init_tbl[55];
add55_init_tbl[0] = 10097325337652013586;
add55_init_tbl[1] = 8422689531964509303;
add55_init_tbl[2] = 9376707153831131165;
add55_init_tbl[3] = 12807999708015736147;
add55_init_tbl[4] = 12171768336606574717;
add55_init_tbl[5] = 2051656926866574818;
add55_init_tbl[6] = 3529647783580834282;
add55_init_tbl[7] = 13746700781847540610;
add55_init_tbl[8] = 17468509505804776974;
add55_init_tbl[9] = 14385537637435099817;
add55_init_tbl[10] = 14225685144642756788;
add55_init_tbl[11] = 11100020401286074697;
add55_init_tbl[12] = 7207317906119690446;
add55_init_tbl[13] = 15474452669527079953;
add55_init_tbl[14] = 16868487670307112059;
add55_init_tbl[15] = 4493524947524633824;
add55_init_tbl[16] = 13021248927856520106;
add55_init_tbl[17] = 15956600001874392423;
add55_init_tbl[18] = 1758753794041921585;
add55_init_tbl[19] = 1540354560501451176;
add55_init_tbl[20] = 5335129695612719255;
add55_init_tbl[21] = 9973334408846123356;
add55_init_tbl[22] = 2295368703230757546;
add55_init_tbl[23] = 15020099946907494138;
add55_init_tbl[24] = 10518216150184876938;
add55_init_tbl[25] = 9188200973282539527;
add55_init_tbl[26] = 4220863048338987374;
add55_init_tbl[27] = 682273982071453295;
add55_init_tbl[28] = 7706178130835869910;
add55_init_tbl[29] = 4618975533122308420;
add55_init_tbl[30] = 397583911260717646;
add55_init_tbl[31] = 5686731560708285046;
add55_init_tbl[32] = 10123916228549657560;
add55_init_tbl[33] = 1304775865627110086;
add55_init_tbl[34] = 15501295782182641134;
add55_init_tbl[35] = 3061180729620744156;
add55_init_tbl[36] = 6958929830512809719;
add55_init_tbl[37] = 10850627469959910507;
add55_init_tbl[38] = 13499063195307571839;
add55_init_tbl[39] = 6410193623982098952;
add55_init_tbl[40] = 4111084083850807341;
add55_init_tbl[41] = 17719042079595449953;
add55_init_tbl[42] = 5462692006544395659;
add55_init_tbl[43] = 18288274374963224041;
add55_init_tbl[44] = 8337656769629990836;
add55_init_tbl[45] = 7477446061798548911;
add55_init_tbl[46] = 9815931464890815877;
add55_init_tbl[47] = 6913451974267278601;
add55_init_tbl[48] = 11883095286301198901;
add55_init_tbl[49] = 14974403441045516019;
add55_init_tbl[50] = 14210337129134237821;
add55_init_tbl[51] = 12883973436502761184;
add55_init_tbl[52] = 4285013921797415077;
add55_init_tbl[53] = 16435915296724552670;
add55_init_tbl[54] = 3742838738308012451;
/*
* additive 55 tables - used by a55rand() and sa55rand()
*/
add55_j = 23; /* the first walking table pointer */
add55_k = 54; /* the second walking table pointer */
add55_seed64 = 0; /* lower 64 bits of the reseeded seed */
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 */
cryrand_init_p = 0x1aa9e726a735044;
cryrand_init_p <<= 200;
cryrand_init_p |= 0x73b7457c5297122310880fcbfa8d4e38380a00396d533300bb;
/* q, second Blum prime */
cryrand_init_q = 0xa62ee0481aa12059c3;
cryrand_init_q <<= 200;
cryrand_init_q |= 0x79ef44e72ff58ea70cacabbe9d264a0b16db72117d96f77e17;
/* quadratic residue of n=p*q */
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 */
cryrand_p = cryrand_init_p;
/* q, second Blum prime */
cryrand_q = cryrand_init_q;
/* n = p*q */
cryrand_n = cryrand_p*cryrand_q;
/* quad residue of n */
cryrand_r = cryrand_init_r;
/* this cryrand() running exp used in computing crypto good bits */
cryrand_exp = cryrand_r;
/* good crypto bits unused from the last cryrand() call */
cryrand_left = 0;
/* the value cryrand_left contains cryrand_bitcnt crypto good bits */
cryrand_bitcnt = 0;
/*
* shufrand - shuffle generator constants
*/
shuf_size = 128; /* entries in shuffle table */
shuf_shift = (64-highbit(shuf_size)); /* shift a55 value to get tbl indx */
global shuf_y; /* Y value (previous output) */
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.
*/
nxtprime_pfact10 = pfact(10);
nxtprime_pfact100 = pfact(100)/nxtprime_pfact10;
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 are looking for a fast generator and do not need a crypto
* strong generator, you should not use this generator. Use of
* the shuffle generator functions shufrand() and rand() should
* be considered instead.
*/
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:
*
* 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.
*
* 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 are looking for a fast generator and do not need a crypto
* strong generator, you should not use this generator. Use of the
* shuffle generator seed functions sshufrand() and srand() should
* be considered instead.
*/
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 < 0) {
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 */
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(seed);
/* spin the pseudo-random number generator a while */
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
*/
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 a Pyramid 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 rlim; /* high quadratic res search value */
/*
* 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_p, cryrand_q);
} 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 (4 && 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 if -1, force minimum values on
* 'ip', 'iq' and 'ir'.
*/
if (seed >= 0) {
if (!isint(ip) || ip < 3) {
/* smallest Blum prime */
ip = 3;
}
if (!isint(iq) || iq < 3) {
/* smallest Blum prime */
iq = 3;
}
if (!isint(ir) || ir < 4) {
/* cryres() never selects a value < 4 */
ir = 4;
}
}
/*
* Determine our prime 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 */
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
*
* Try to avoid trivial quadratic residues by forcing the search
* over the interval [4, n-3).
*/
rlim = cryrand_n-3;
/* if ir is too big, drop down to the bottom of the range */
if (ir >= rlim) {
ir = 4;
}
while (jacobi(ir,cryrand_p) != 1 || jacobi(ir,cryrand_q) != 1) {
/* if ir is too big, drop down to the bottom of the range */
if (++ir >= rlim) {
ir = 4;
}
}
}
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 < 4) ||
(isint(len2) && len2 != -1 && len2 < 4) ||
(!isint(len1) && isint(len2)) ||
(isint(len1) && !isint(len2))) {
quit "bad args: 2 & 3: if 2nd is given, must be -1 or ints > 3";
}
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;
cryrand_r = cryres(cryrand_p, cryrand_q);
/*
* 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)
{
local x; /* firewall comparator */
local prev; /* previous states of the generators */
local junk; /* dummy holder of random values */
/* declare our objects */
obj cryobj x; /* firewall comparator */
obj cryobj prev; /* previous states of the generators */
/* 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));
mat shuf_tbl[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)
{
local x; /* the object being created */
/* declare our objects */
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);
}
/*
* 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)
{
local x; /* firewall comparator */
/* declare our objects */
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: 12th 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: 13th 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:
* p first prime
* q second prime
*
* 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 primes.
*
* NOTE: This is NOT Blum's method. This is used by Blum's method to
* generate some internal values.
*/
define
cryres(p,q)
{
local rval; /* a temporary pseudo-random value */
local sqrpow; /* a power of 2 starting > n^(3/4) */
local quadres; /* 0 or a quadratic residue of n = p*q */
local n; /* n=p*q */
local j;
/*
* firewall
*/
if (!isint(p) || !isint(q) || p < 3 || q < 3) {
quit "bad args: p and q must be integers > 2";
}
/*
* find a pseudo-random quadratic residue of n = p*q
*
* To find a pseudo-random quadratic residue, we will generate
* pseudo-random numbers to try in the interval [2^sqrpow,n-3),
* where 2^sqrpow is the first power of 2 > n^(3/4). This range
* helps us avoid selecting trivial residues abs(quadres mod n=p*q)
* is tiny. We will never select a residue < 4.
*
* When we fail to find a quadratic residue after 1009 tries, we will
* drop sqrpow by 1 and start at another pseudo-random location.
*
* It is very unlikely that we will need to search more than a
* few numbers to find a quadratic residue of n = p*q.
*
* We will reject any quadratic residue that is a square of
* itself, mod n=p*q.
*/
n = p*q;
quadres = 0;
sqrpow = highbit(n)//4*3;
do {
/* generate a pseudo-random starting point */
rval = rand(2^sqrpow, n-4);
/* look for 1009 times or until >= cryrand_n */
for (j=0; j < 1009; ++j, ++rval) {
/*
* check if we have a quadratic residue of n = p*q
*
* 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
*/
if (jacobi(rval,p) == 1 && jacobi(rval,q) == 1) {
quadres = rval;
break;
}
}
/* setup for a new pseudo-random starting point if we found nothing */
if (quadres == 0) {
/* reduce sqrpow if reasonable */
if (sqrpow > 2) {
--sqrpow;
}
}
} while (quadres == 0 || ((quadres^2)%n) == quadres);
/*
* return the quadratic residue of n = p*q
*/
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 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 = srand(0); /* pre-initialize the tables */
global lib_debug;
if (!isnum(lib_debug) || lib_debug>0) print "ver: 10.5 06:19:34 5/9/93";
if (!isnum(lib_debug) || lib_debug>0) print "shufrand() defined";
if (!isnum(lib_debug) || lib_debug>0) print "sshufrand(seed) defined";
if (!isnum(lib_debug) || lib_debug>0) print "rand([a, [b]]) defined";
if (!isnum(lib_debug) || lib_debug>0) print "srand(seed) defined";
if (!isnum(lib_debug) || lib_debug>0) print "cryrand([a, [b]]) defined";
if (!isnum(lib_debug) || lib_debug>0) print "scryrand([seed, [len1, len2]]) defined";
if (!isnum(lib_debug) || lib_debug>0) print "scryrand(seed, ip, iq, ir) defined";
if (!isnum(lib_debug) || lib_debug>0) print "random([a, [b]]) defined"
if (!isnum(lib_debug) || lib_debug>0) print "srandom(seed) defined"
if (!isnum(lib_debug) || lib_debug>0) print "obj cryobj defined";
if (!isnum(lib_debug) || lib_debug>0) print "randstate([cryobj | 0]) defined"
if (!isnum(lib_debug) || lib_debug>0) print "nxtprime(n, [val, modulus]) defined";