zrandom.c

Calc Software Package for Number Calc

 *	      to having their Blum moduli factored, depending in their size,
 *	      by small PCs in a reasonable to large supercomputers/highly
 *	      parallel processors over a long time.  Their value lies in their
 *	      speed relative the the default Blum generator.  As of Jan 1997,
 *	      the Blum moduli associated with 13 <= newn < 20 appear to
 *	      be well beyond the scope of hardware and algorithms,
 *	      and 9 <= newn < 12 might be factorable with extreme difficulty.
 *
 *	The following table may be useful as a guide for how easy it
 *	is to factor the modulus:
 *
 *	1 <= newn <= 4	PC using ECM in a short amount of time
 *	5 <= newn <= 8	Workstation using MPQS in a short amount of time
 *	8 <= newn <= 12 High end supercomputer or high parallel processor
 *			using state of the art factoring over a long time
 *     12 <= newn <= 16 Beyond Feb 1997 systems and factoring methods
 *     17 <= newn <= 20 Well beyond Feb 1997 systems and factoring methods
 *
 *	      See the section titled 'FOR THE PARANOID' for more details.
 *
 *    seed < 0, 0 < newn <= 20:
 *    -------------------------
 *	Reserved for future use.
 *
 ******************************************************************************
 *
 * srandom(seed, ip, iq, trials)
 *
 *    We attempt to set the Blum modulus and quadratic residue.
 *    Any internally buffered random bits are flushed.
 *
 *    Use the ip and iq args as starting points for Blum primes.
 *    The trials arg determines how many ptest cycles are performed
 *    on each candidate.
 *
 *    The new quadratic residue is set according to the seed
 *    arg as defined below.
 *
 *    seed >= 2^32, ip >= 2^16, iq >= 2^16:
 *    -------------------------------------
 *	Set the Blum modulus by searching from the ip and iq search points.
 *
 *	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:
 *
 *	  p = nextcand(ip-2, trials, 0, 3, 4);	  (* find the 1st Blum prime *)
 *	  q = nextcand(iq-2, trials, 0, 3, 4);	  (* find the 2nd Blum prime *)
 *	  n = p * q;			(* 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 *)
 *	  srandom(r, n);
 *
 *    any seed, ip <= 2^16 or iq <= 2^16:
 *    -----------------------------------
 *	Reserved for future use.
 *
 *    0 < seed < 2^32, any ip, any iq:
 *    --------------------------------
 *	Reserved for future use.
 *
 *    seed == 0, ip > 2^16, iq > 2^16:
 *    --------------------------------
 *	Set the Blum modulus by searching from the ip and iq search points.
 *	If trials is omitted, 1 is assumed.
 *
 *	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, ip, iq, trials)
 *
 *	or in other words:
 *
 *	  (* trials, if omitted, is assumed to be 1 *)
 *	  p = nextcand(ip-2, trials, 0, 3, 4);	  (* find the 1st Blum prime *)
 *	  q = nextcand(iq-2, trials, 0, 3, 4);	  (* find the 2nd Blum prime *)
 *	  n = p * q;			(* n is the new Blum modulus *)
 *	  r = default_residue;		(* as used by the initial state *)
 *	  do {
 *	      last_r = r;
 *	      r = pmod(last_r, 2, n);
 *	  } while (r > last_r);		(* r is the new quadratic residue *)
 *
 *    seed < 0, any ip, any iq:
 *    -------------------------
 *	Reserved for future use.
 *
 ******************************************************************************
 *
 * srandom()
 *
 *    Return current Blum generator state.  This call does not alter
 *    the generator state.
 *
 ******************************************************************************
 *
 * srandom(state)
 *
 *	Restore the Blum state and return the previous state.  Note that
 *	the argument state is a random state value (israndom(state) is true).
 *	Any internally buffered random bits are restored.
 *
 *	The states of the Blum 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.
 *
 *	    random_state = srandom();
 *	    ... generate random bits ...
 *	    prev_random_state = srandom(random_state);
 *	    ... generate the same random bits ...
 *	    srandom() == prev_random_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:
 *
 * When the call:
 *
 *	srandom(seed, nextcand(ip,25,0,3,4)*nextcand(iq,25,0,3,4))
 *
 * probable primes from nextcand are used.  We use the word
 * 'probable' because of an extremely extremely small chance that a
 * composite (a non-prime) may be returned.
 *
 * We use the builtin function nextcand in its 5 arg form:
 *
 *	nextcand(value, 25, 0, 3, 4)
 *
 * The odds that a number returned by the above call is not prime is
 * less than 1 in 4^25.	 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 nextcand() 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.
 *
 * Even though Blum generator is 'pseudo-random', there is no statistical
 * test, which runs in polynomial time, that can distinguish the Blum
 * generator from a truly random source.  See the comment under
 * the "Blum-Blum-Shub generator" section above.
 *
 * 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:
 *
 * When seeding the Blum generator, we disallow seeds < 2^32 in an
 * effort to avoid trivial seed values such as 0, 1 or other small values.
 * The 2^32 lower bound limit was also selected because it provides a
 * large reserved value space for special seeds.  Currently the
 * [1,2^32) seed range is reserved for future use.
 *
 ***
 *
 * When using the 2 arg srandom(seed, newn), we require newn > 2^32
 * to avoid trivial Blum moduli.  The 2^32 lower bound limit was also
 * selected because it provides a large reserved value space for special
 * moduli.  Currently the [21,2^32) newn range is reserved for future use.
 *
 * When using the 3 or 4 arg srandom(seed, ip, iq [, trials]) form,
 * we require ip>2^16 and ip>2^16.  This ensures that the resulting
 * Blum modulus is > 2^32.
 *
 ***
 *
 * Taking some care to select a good initial residue helps eliminate cheap
 * search attacks.  It is true that a subsequent residue could be one of the
 * residues that we would first avoid.	However such an occurrence will
 * happen after the generator is well underway and any such seed information
 * has been lost.
 *
 * In an effort to avoid trivial seed values, we force the seed arg
 * to srandom() to be > 2^32.  We then square this value mod the
 * Blum modulus until it is less than the previous value.  This ensures
 * that the previous seed value is large enough that its square is > Blum
 * modulus, and this the square mod Blum modulus is non-trivial.
 *
 ***
 *
 * The size of default Blum modulus 'n=p*q' was taken to be > 2^259, or
 * 260 bits (79 digits) long.  A modulus > 2^256 will generate 8 bits
 * per crank of the generator.	The period of this generator is long
 * enough to be reasonable, and the modulus is small enough to be fast.
 *
 * The default Blum modulus is not a secure modulus because it can
 * be factored with ease.  As if Feb 1997,  the upper reach of the
 * state of the art for factoring general numbers was around 2^512.
 * Clearly factoring a 260 bit number if well within the reach of even
 * a low life Pentium.
 *
 * The fact that the default modulus can be factored with ease is
 * not a drawback.  Afterall, if we are going to keep to the goal
 * of disclosing the source magic numbers, we need to disclose how
 * the Blum Modulus was produced ... including its factors.  Knowing
 * the facotrs of the Blum modulus does not reduce its quality,
 * only the ability for someone to determine where you are in the
 * sequence.  But heck, the default seed is well known anyway so
 * there is no real loss if the factors are also known.
 *
 ***
 *
 * The default Blum modulus is the product of two Blum probable primes
 * that were selected by the Rand Book of Random Numbers.  Using the '6% rule',
 * a default Blum modulus n=p*q > 2^256 would be satisfied if p were
 * 38 decimal digits and q were 42 decimal digits in length.  We restate
 * the sizes in decimal digits because the Rand Book of Random Numbers is a
 * book of decimal digits.  Using the first 38 rand digits as a
 * starting search point for 'p', and the next 42 digits for a starting
 * search point for 'q'.
 *
 *	(*
 *	 * setup the search points (lines split for readability)
 *	 *)
 *	ip = 10097325337652013586346735487680959091;
 *	iq = 173929274945375420480564894742962480524037;
 *
 *	(*
 *	 * find the first Blum prime
 *	 *)
 *	fp = int((ip-1)/2);
 *	do {
 *		fp = nextcand(fp+2, 25, 0, 3, 4);
 *		p = 2*fp+1;
 *	} while (ptest(p, 25) == 0);
 *
 *	(*
 *	 * find the 2nd Blum prime
 *	 *)
 *	fq = int((iq-1)/2);
 *	do {
 *		fq = nextcand(fq+2, 25, 0, 3, 4);
 *		q = 2*fq+1;
 *	} while (ptest(q, 25) == 0);
 *
 * The above resource file produces the Blum probable primes and initial
 * quadratic residue (line wrapped for readability):
 *
 *	    p= 0x798ac934c7a3318ad446190f3474e57
 *
 *	    q= 0x1ff21d7e1dd7d5965e224d485d84c3ef44f
 *
 * These Blum primes were found after 1.81s of CPU time on a 195 Mhz IP28
 * R10000 version 2.5 processor.  The first Blum prime 'p' was 31716 higher
 * than the initial search value 'ip'.	The second Blum prime 'q' was 18762
 * higher than the initial starting 'iq'.
 *
 * The product of the two Blum primes results in a 260 bit Blum modulus of:
 *
 *	n = 0xf2ac1903156af9e373d78613ed0e8d30284f34b644a9027d9ba55a689d6be18d9
 *
 * The selection if the initial quadratic residue comes from the next
 * unused digits of the Rand Book of Random Numbers.  Now the two initial
 * search values 'ip' and 'iq' used above needed the first 38 digits and
 * the next 42 digits.	Thus we will skip the first 38+42=80 digits
 * and begin to build in initial search value for a quadratic residue (most
 * significant digit first) from the Rand Book of Numbers digits until we
 * have a value whose square mod n > 4th power mod n.  In other words, we
 * need to build ir up from new Rand Book of Random Numbers digits until
 * we find a value in which srandom(ir), for the Blum Modulus 'n' produces
 * an initial quadratic residue on the first loop.
 *
 * Clearly we need to find an ir that is > sqrt(n).  The first ir:
 *
 *	ir = 2063610402008229166508422689531964509303
 *
 * fails the single loop criteria.  So we add the next digit:
 *
 *	ir = 20636104020082291665084226895319645093032
 *
 * Here we find that:
 *
 *	pmod(ir,2,n) > pmod(pmod(ir,2,n),2,n)
 *
 * Thus, for thw Blum modulus 'n', the method outlined for srandom(ir) yields
 * the initial quadratic residue of:
 *
 *    r = 0x748b6d882ff4b074e2f1e93a8627d626506c73ca5a62546c90f23fd7ed3e7b11e
 *
 ***
 *
 * In the above process of finding the Blum primes used in the default
 * Blum modulus, we selected primes of the form:
 *
 *	2*x + 1		x is also prime
 *
 * because Blum generators with modulus 'n=p*q' have a period:
 *
 *	lambda(n) = lcm(factors of p-1 & q-1)