zrand.c

Calc Software Package for Number Calc

/*
 * zrand - subtractive 100 shuffle generator
 *
 * Copyright (C) 1999,2004  Landon Curt Noll
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.9 $
 * @(#) $Id: zrand.c,v 29.9 2004/03/31 04:58:40 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/zrand.c,v $
 *
 * Under source code control:	1995/01/07 09:45:25
 * File existed as early as:	1994
 *
 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * AN OVERVIEW OF THE FUNCTIONS:
 *
 * This module contains an Subtractive 100 shuffle generator wrapped inside
 * of a shuffle generator.
 *
 *	We refer to this generator as the s100 generator.
 *
 *	rand	- s100 shuffle generator
 *	srand	- seed the s100 shuffle generator
 *
 *	This generator has two distinct parts, the s100 generator
 *	and the shuffle generator.
 *
 *	The subtractive 100 generator is described in Knuth's "The Art of
 *	Computer Programming - Seminumerical Algorithms", Vol 2, 3rd edition
 *	(1998), Section 3.6, page 186, formula (2).
 *
 *	The "use only the first 100 our of every 1009" is described in
 *	Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 *	Vol 2, 3rd edition (1998), Section 3.6, page 188".
 *
 *	The period and other properties of this generator make it very
 *	useful to 'seed' other generators.
 *
 *	The shuffle generator is described in Knuth's "The Art of Computer
 *	Programming - Seminumerical Algorithms", Vol 2, 3rd edition (1998),
 *	Section 3.2.2, page 34, Algorithm B.
 *
 *	The shuffle generator is fast and serves as a fairly good standard
 *	pseudo-random generator.   If you need a fast generator and do not
 *	need a cryptographically strong one, this generator is likely to do
 *	the job.
 *
 *	The shuffle generator is feed values by the subtractive 100 process.
 *
 ******************************************************************************
 *
 * 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 need
 *
 *	    By using the appropriate seed and other arguments, one may
 *	    increase the strength of the generator to suit the need of
 *	    the application.  One does not have just a few levels.
 *
 * Even though I have done my best to implement a good system, you still
 * must use these routines your own risk.
 *
 * Share and enjoy!  :-)
 */

/*
 * ON THE GENERATORS:
 *
 * The subtractive 100 generator has a good period, and is fast.  It is
 * reasonable as generators go, though there are better ones available.
 * The shuffle generator has a very good period, and is fast.  It is
 * fairly good as generators go, particularly when it is feed reasonably
 * random numbers.  Because of this, we use feed values from the subtractive
 * 100 process into the shuffle generator.
 *
 * The s100 generator uses 2 tables:
 *
 *	subtractive table - 100 entries of 64 bits used by the subtractive 100
 *			    part of the s100 generator
 *
 *	shuffle table - 256 entries of 64 bits used by the shuffle
 *			part of the s100 generator and feed by the
 *			subtractive table.
 *
 * Casual direct use of the shuffle generator may be acceptable.  If one
 * desires cryptographically strong random numbers, or if one is paranoid,
 * one should use the Blum generator instead (see zrandom.c).
 *
 * The s100 generator as the following calc interfaces:
 *
 *   rand(min,max)		(where min < max)
 *
 *	Print an s100 generator random value over interval [a,b).
 *
 *   rand()
 *
 *	Same as rand(0, 2^64).	Print 64 bits.
 *
 *   rand(lim)			(where 0 > lim)
 *
 *	Same as rand(0, lim).
 *
 *   randbit(x)			(where x > 0)
 *
 *	Same as rand(0, 2^x).  Print x bits.
 *
 *   randbit(skip)		(where skip < 0)
 *
 *	Skip random bits and return the bit skip count (-skip).
 */

/*
 * INITIALIZATION AND SEEDS:
 *
 * All generators come already seeded with precomputed initial constants.
 * Thus, it is not required to seed a generator before using it.
 *
 * The s100 generator may be initialized and seeded via srand().
 *
 * Using a seed of '0' will reload generators with their initial states.
 *
 *	srand(0)	restore subtractive 100 generator to the initial state
 *
 * The above single arg calls are fairly fast.
 *
 * Optimal seed range for the s100 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 * 100!) should be
 *	sufficient for most purposes.  An easy way to stay within this
 *	range to to use seeds that are between 21 and 178 digits, or
 *	64 to 588 bits long.
 *
 *	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.
 *
 ******************************************************************************
 *
 * srand(seed)
 *
 *    Seed the s100 generator.
 *
 *    seed != 0:
 *    ---------
 *	Any buffered random bits are flushed.  The subtractive table is loaded
 *	with the default subtractive table.  The low order 64 bits of seed is
 *	xor-ed against each table value.  The subtractive table is shuffled
 *	according to seed/2^64.
 *
 *	The following calc code produces the same effect:
 *
 *	    (* reload default subtractive table xored with low 64 seed bits *)
 *	    seed_xor = seed & ((1<<64)-1);
 *	    for (i=0; i < 100; ++i) {
 *		subtractive[i] = xor(default_subtractive[i], seed_xor);
 *	    }
 *
 *	    (* shuffle the subtractive table *)
 *	    seed >>= 64;
 *	    for (i=100; seed > 0 && i > 0; --i) {
 *		quomod(seed, i+1, seed, j);
 *		swap(subtractive[i], subtractive[j]);
 *	    }
 *
 *	Seed must be >= 0.  All seed values < 0 are reserved for future use.
 *
 *	The subtractive 100 pointers are reset to subtractive[36] and
 *	subtractive[99].  Last the shuffle table is loaded with successive
 *	values from the	subtractive 100 generator.
 *
 *    seed == 0:
 *    ---------
 *	Restore the initial state and modulus of the s100 generator.
 *	After this call, the s100 generator is restored to its initial
 *	state after calc started.
 *
 *	The subtractive 100 pointers are reset to subtractive[36] and
 *	subtractive[99].  Last the shuffle table is loaded with successive
 *	values from the	subtractive 100 generator.
 *
 ******************************************************************************
 *
 * srand(mat100)
 *
 *    Seed the s100 generator.
 *
 *    Any buffered random bits are flushed.  The subtractive table with the
 *    first 100 entries of the array mat100, mod 2^64.
 *
 *    The subtractive 100 pointers are reset to subtractive[36] and
 *    subtractive[99].  Last the shuffle table is loaded with successive
 *    values from the subtractive 100 generator.
 *
 ******************************************************************************
 *
 * srand()
 *
 *    Return current s100 generator state.  This call does not alter
 *    the generator state.
 *
 ******************************************************************************
 *
 * srand(state)
 *
 *    Restore the s100 state and return the previous state.  Note that
 *    the argument state is a rand state value (isrand(state) is true).
 *    Any internally buffered random bits are restored.
 *
 *    The states of the s100 generators can be saved by calling the seed
 *    function with no arguments, and later restored by calling the seed
 *    functions with that same return value.
 *
 *	rand_state = srand();
 *	... generate random bits ...
 *	prev_rand_state = srand(rand_state);
 *	... generate the same random bits ...
 *	srand() == prev_rand_state;		(* is true *)
 *
 *    Saving the state just after seeding a generator and restoring it later
 *    as a very fast way to reseed a generator.
 */

/*
 * TRUTH IN ADVERTISING:
 *
 * A "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.
 *
 * A final "truth in advertising" issue deals with how the magic numbers
 * found in this file were generated.  Detains can be found in the
 * various functions, while a overview can be found in the "SOURCE FOR
 * MAGIC NUMBERS" section below.
 */

/*
 * SOURCE OF MAGIC NUMBERS:
 *
 * Most of the magic constants used on this file ultimately are
 * based on LavaRnd.  LavaRnd produced them via a cryprographic
 * of the digitization of chaotic system that consisted of a noisy
 * digital camera and 6 Lava Lite(R) lamps.
 *
 * 	BTW: Lava Lite(R) is a trademark of Haggerty Enterprises, Inc.
 *
 * The first 100 groups of 64 bit bits were used to initialize init_s100.slot.
 *
 * The function, randreseed64(), uses 2 primes to scramble 64 bits
 * into 64 bits.  These primes were also extracted from the Rand
 * Book of Random Numbers.  See randreseed64() for details.
 *
 * The shuffle table size is longer than the 100 entries recommended
 * by Knuth.  We use a power of 2 shuffle table length so that the
 * shuffle process can select a table entry from a new subtractive 100
 * value by extracting its low order bits.  The value 256 is convenient
 * in that it is the size of a byte which allows for easy extraction.
 *
 * We use the upper byte of the subtractive 100 value to select the shuffle
 * table entry because it allows all of 64 bits to play a part in the
 * entry selection.  If we were to select a lower 8 bits in the 64 bit
 * value, carries that propagate above our 8 bits would not impact the
 * s100 generator output.
 *
 ******************************************************************************
 *
 * 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
 * if you that paranoid why would you use a non-cryprographically strong
 * pseudo-random number generator in the first place?  You would be
 * better off using the random() builtin function.
 *
 * The two constants that were picked from the Rand Book of Random Numbers
 * 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