zrandom.c

Calc Software Package for Number Calc

/*
 * zrandom - Blum-Blum-Shub pseudo-random 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.8 $
 * @(#) $Id: zrandom.c,v 29.8 2004/03/31 04:58:40 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/zrandom.c,v $
 *
 * Under source code control:	1997/02/15 04:01:56
 * File existed as early as:	1997
 *
 * 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 a Blum-Blum-Shub generator:
 *
 *	We refer to this generator as the Blum generator.
 *
 *	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.
 *
 *	    Sciences 28, pp. 270-299.
 *
 *	    Bruce Schneier, "Applied Cryptography", John Wiley & Sons,
 *	    1st edition (1994), pp 365-366.
 *
 *	This generator is considered 'strong' in that it passes all
 *	polynomial-time statistical tests.  The sequences produced
 *	are random in an absolutely precise way.  There is absolutely
 *	no better way to predict the sequence than by tossing a coin
 *	(as with TRULY random numbers) EVEN IF YOU KNOW THE MODULUS!
 *	Furthermore, having a large chunk of output from the sequence
 *	does not help.	The BITS THAT FOLLOW OR PRECEDE A SEQUENCE
 *	ARE UNPREDICTABLE!
 *
 *	Of course the Blum modulus should have a long period.  The default
 *	Blum modulus as well as the compiled in Blum moduli have very long
 *	periods.  When using your own Blum modulus, a little care is needed
 *	to avoid generators with very short periods.  (see below)
 *
 *	To compromise the generator, an adversary must either factor the
 *	modulus or perform an exhaustive search just to determine the next
 *	(or previous) bit.  If we make the modulus hard to factor
 *	(such as the product of two large well chosen primes) breaking
 *	the sequence could be intractable for todays computers and methods.
 *
 ******************************************************************************
 *
 * 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 GENERATOR:
 *
 * The Blum 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 Blum generator is not too slow, though
 * seeding the generator via srandom(seed,plen,qlen) can be slow.
 * Shortcuts and pre-defined generators have been provided for this reason.
 * Use of Blum should be more than acceptable for many applications.
 *
 * The Blum generator as the following calc interfaces:
 *
 *   random(min, max)		(where min < max)
 *
 *	Print a Blum generator random value over interval [min,max).
 *
 *   random()
 *
 *	Same as random(0, 2^64).  Print 64 bits.
 *
 *   random(lim)		(where 0 > lim)
 *
 *	Same as random(0, lim).
 *
 *   randombit(x)		(where x > 0)
 *
 *	Same as random(0, 2^x).	 Print x bits.
 *
 *   randombit(skip)		(where skip < 0)
 *
 *	Skip 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 Blum generator may be initialized and seeded via srandom().
 *
 * Using a seed of '0' will reload generators with their initial states.
 *
 *	srandom(0)	restore Blum generator to the initial state
 *
 * The above single arg calls are fairly fast.
 *
 * The call:
 *
 *	srandom(seed, newn)
 *
 * is fast when the config value "srandom" is 0, 1 or 2.
 *
 * Optimal seed range for the Blum generator:
 *
 *	There is no limit on the size of a seed.  On the other hand,
 *	in most cases the seed is taken modulo the Blum modulus.
 *	Using a seed that is too small (except for 0) results in
 *	an internal generator be used to increase its size.
 *
 *	It is faster to use seeds that are in the half open internal
 *	[sqrt(n), n) where n is the Blum modulus.
 *
 *	The default Blum modulus is 260 bits long, so when using a the
 *	single arg call, a seed of between 128 and 256 bits is reasonable.
 *
 ******************************************************************************
 *
 * srandom(seed)
 *
 *    We attempt to set the quadratic residue and possibly the Blum modulus.
 *
 *    Any internally buffered random bits are flushed.
 *
 *    The Blum modulus is only set if seed == 0.
 *
 *    The initial quadratic residue is set according to the seed
 *    arg as defined below.
 *
 *    seed >= 2^32:
 *    -------------
 *	Use seed to compute a new quadratic residue for use with
 *	the current Blum modulus.  We will successively square mod Blum
 *	modulus until we get a smaller value (modulus wrap).
 *
 *	The follow calc resource file produces an equivalent effect:
 *
 *	  n = default_modulus;		(* n is the new Blum modulus *)
 *	  r = seed;
 *	  do {
 *	      last_r = r;
 *	      r = pmod(last_r, 2, n);
 *	  } while (r > last_r);		(* r is the new quadratic residue *)
 *
 *	NOTE: The Blum modulus is not set by this call.
 *
 *    0 < seed < 2^32:
 *    ----------------
 *	Reserved for future use.
 *
 *    seed == 0:
 *    ----------
 *	Restore the initial state and modulus of the Blum generator.
 *	After this call, the Blum  generator is restored to its initial
 *	state after calc started.
 *
 *	The Blum prime factors of the modulus have been disclosed (see
 *	"SOURCE OF MAGIC NUMBERS" below).  If you want to use moduli that
 *	have not been disclosed, use srandom(seed, newn) with the
 *	appropriate args as noted below.
 *
 *	The follow calc resource file produces an equivalent effect:
 *
 *	  n = default_modulus;	(* as used by the initial state *)
 *	  r = default_residue;	(* as used by the initial state *)
 *
 *	NOTE: The Blum modulus is changed by this call.
 *
 *    seed < 0:
 *    ---------
 *	Reserved for future use.
 *
 ******************************************************************************
 *
 * srandom(seed, newn)
 *
 *    We attempt to set the Blum modulus and quadratic residue.
 *    Any internally buffered random bits are flushed.
 *
 *    If newn == 1 mod 4, then we will assume that it is the
 *    product of two Blum primes (primes == 3 mod 4) and use it
 *    as the Blum modulus.
 *
 *    The new quadratic residue is set according to the seed
 *    arg as defined below.
 *
 *    seed >= 2^32, newn >= 2^32:
 *    ---------------------------
 *	Assuming that 'newn' == 3 mod 4, then we will use it as
 *	the Blum modulus.
 *
 *	We will use the seed arg to compute a new quadratic residue.
 *	We will successively square it mod Blum modulus until we get
 *	a smaller value (modulus wrap).
 *
 *	The follow calc resource file produces an equivalent effect:
 *
 *	  if (newn % 4 == 1) {
 *	      n = newn;			(* n is the new Blum modulus *)
 *	      r = seed;
 *	      do {
 *		  last_r = r;
 *		  r = pmod(last_r, 2, n);
 *	      } while (r > last_r);	(* r is the new quadratic residue *)
 *	  } else {
 *	      quit "newn (2nd arg) must be 3 mod 4";
 *	  }
 *
 *    0 < seed < 2^32, newn >= 2^32:
 *    ------------------------------
 *	Reserved for future use.
 *
 *    seed == 0, newn >= 2^32:
 *    ------------------------
 *	Assuming that 'newn' == 3 mod 4, then we will use it as
 *	the Blum modulus.
 *
 *	The initial quadratic residue will be as if the default initial
 *	quadratic residue arg was given.
 *
 *	The follow calc resource file produces an equivalent effect:
 *
 *	  srandom(default_residue, newn)
 *
 *	or in other words:
 *
 *	  if (newn % 4 == 1) {
 *	      n = newn;			(* n is the new Blum modulus *)
 *	      r = default_residue;	(* as used by the initial state *)
 *	      do {
 *		  last_r = r;