zrandom.c

Calc Software Package for Number Calc

 *
 * since 'p' and 'q' are both odd, 'p-1' and 'q-1' have 2 as
 * a factor.  The calc resource file above ensures that '(p-1)/2' and
 * '(q-1)/2' are probable prime, thus maximizing the period
 * of the default generator to:
 *
 *	lambda(n) = lcm(2,2,fp,fq) = 2*fp*fq = ~2*(p/2)*(q/2) = ~n/2
 *
 * The process above resulted in a default generator Blum modulus n > 2^259
 * with period of at least 2^258 bits.	To be exact, the period of the
 * default Blum generator is:
 *
 *	    0x79560c818ab57cf1b9ebc309f68746881adc15e79c05e476f741e5f904b9beb1a
 *
 * which is approximately:
 *
 *	   ~8.781 * 10^77
 *
 * This period is more than long enough for computationally tractable tasks.
 *
 ***
 *
 * The 20 builtin generators, on the other hand, were selected
 * with more care.  While the lower order 20 generators have
 * factorable generators, the higher order are likely to be
 * be beyond the reach for a while.
 *
 * The lengths of the two Blum probable primes 'p' and 'q' used to make up
 * the 20 Blum modului 'n=p*q' differ slightly to avoid certain
 * 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 up to 6% of the product size to avoid such attacks.
 * Clearly one does not want the size of the two factors to differ
 * by a large percentage: p=3 and q large would result in a easy
 * to factor Blum modulus.  Thus we select sizes that differ by
 * up to 6% but not (significantly) greater than 6%.
 *
 * Using the '6% rule' above, a Blum modulus n=p*q > 2^1024 would have two
 * Blum factors p > 2^482 and q > 2^542.  This is because 482+542 = 1024.
 * The difference 542-482 is ~5.86% of 1024, and is the largest difference
 * that is < 6%.
 *
 ******************************************************************************
 *
 * 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.
 *
 ***
 *
 * One could take issue with the above resource file that produced a 260 bit
 * Blum modulus.  So if that bothers you, then seed your generator
 * with your own Blum modulus and initial quadratic residue.  And
 * if you are truly paranoid, why would you want to use the default seed,
 * which is well known?
 *
 ***
 *
 * If all the above fails to pacify the truly paranoid, then one may
 * select your own modulus and initial quadratic residue by calling:
 *
 *	srandom(r, n);
 *
 * Of course, you will need to select a correct Blum modulus 'n' as the
 * product of two Blum primes (both 3 mod 4) and with a long period (where
 * lcm(factors of one less than the two Blum primes) is large) and an
 * initial quadratic residue 'r' that is hard to guess selected at random.
 *
 * A simple way to seed the generator would be to:
 *
 *	srandom(ir, ip, iq, 25)
 *
 * where 'ip', 'iq' and 'ir' are large integers that are unlikely to be
 * 'guessed' and where numbers around the size of iq*ir are beyond
 * the current reach of the best factoring methods on the fastest
 * SGI/Cray supercomuters.
 *
 * Of course you can increase the '25' value if 1 of 4^25 odds of a
 * non-prime are too probable for you.
 *
 * The problem with the above call is that the period of the Blum generator
 * could be small if 'p' and 'q' have lots of small prime factors in common.
 *
 * A better way to do seed the Blum generator yourself is to use the
 * seedrandom(seed1, seed2, size [,tests]) function from "seedrandom.cal"
 * with the args:
 *
 *	seed1 - seed rand() to search for 'p', select from [2^64, 2^308)
 *	seed2 - seed rand() to search for 'q', select from [2^64, 2^308)
 *	size  - minimum bit size of the Blum modulus 'n=p*q'
 *	tests - optional arg for number of pseudo prime tests (default is 25)
 *
 * Last, one could use some external random source to select starting
 * search points for 'p', 'q' and the quadratic residue.  One way to
 * do this is:
 *
 *	fp = int((ip-1)/2);
 *	do {
 *	    fp = nextcand(fp+2, tests, 0, 3, 4);
 *	    p = 2*fp+1;
 *	} while (ptest(p, tests) == 0);
 *	fq = int((iq-1)/2);
 *	do {
 *	    fq = nextcand(fq+2, tests, 0, 3, 4);
 *	    q = 2*fq+1;
 *	} while (ptest(q, tests) == 0);
 *	srandom(pmod(ir,2,p*q), p*q);
 *
 * where 'tests' is the number of pseudo prime tests that a candidate must
 * pass before being considered a probable prime (must be >0, perhaps 25), and
 * where 'ip' is the initial search location for the Blum prime 'p', and
 * where 'iq' is the initial search location for the Blum prime 'q', and
 * where 'ir' is the initial Blum quadratic residue generator.	The 'ir'
 * value should be a random value in the range [2^(binsize*4/5), 2^(binsize-2))
 * where 2^(binsize-1) < n=p*q <= 2^binsize.
 *
 * Your external generator would need to generate 'ip', 'iq' and 'ir'.
 * While any value for 'ip' and 'iq will do (provided that their product
 * is large enough to meet your modulus needs), 'ir' should be selected
 * to avoid values near 0 or near 'n' (or ip*iq).
 *
 ***
 *
 * The Blum moduli used with the pre-defined generators (via the call
 * srandom(seed, 0<x<=16)) were generated using the above process.  The
 * 'ip', 'iq' and 'ir' values were produced by special purpose hardware
 * that produced cryptographically strong random numbers.  The Blum
 * primes used in these special pre-defined generators are unknown.
 *
 * Not being able to factor 'n=p*q' into 'p' and 'q' does not directly
 * improve the quality Blum generator.	On the other hand, it does
 * improve the security of it.
 *
 * I (Landon Curt Noll) did not keep the search values of these 20 special
 * pre-defined generators.  While some of the smaller Blum moduli is
 * within the range of some factoring methods, others are not.	As of
 * Feb 1997, the following is the estimate of what can factor the
 * pre-defined moduli:
 *
 *	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
 *
 * This is not to say that in the future things will not change.  One
 * can say that faster hardware will not help in the factoring of 1024+
 * bit Blum moduli found in 12 <= newn as well as in the default
 * Blum generator.  A combination of better algorithms, helped by faster
 * hardware will be needed.
 *
 * It is true that the the pre-defined moduli are 'magic'.  On the other
 * hand, there purpose was in part to produce users with pre-seeded
 * generators where the individual Blum primes are not well known.  If
 * this bothers you, don't use them and seed your own!
 *
 ***
 *
 * The output of the Blum generator has been proven to be cryptographically
 * strong.  I.e., there is absolutely no better way to predict the next
 * bit in the sequence than by tossing a coin (as with TRULY random numbers)
 * EVEN IF YOU HAVE A LARGE NUMBER OF PREVIOUSLY GENERATED BITS AND KNOW
 * A LARGE NUMBER OF BITS THAT FOLLOW THE NEXT BIT!  An adversary would
 * be far better advised to try to factor the modulus or exhaustively search
 * for the quadratic residue in use.
 *
 */


#include "zrandom.h"
#include "have_const.h"
#include "have_unused.h"


/*
 * current Blum generator state
 */
static RANDOM blum;


/*
 * Default Blum generator
 *
 * The init_blum generator is established via the srandom(0) call.
 *
 * The z_rdefault ZVALUE is the 'r' (quadratic residue) of init_blum.
 */
#if FULL_BITS == 64
static CONST HALF h_ndefvec[] = {
	(HALF)0xd6be18d9, (HALF)0xba55a689, (HALF)0x4a9027d9, (HALF)0x84f34b64,
	(HALF)0xd0e8d302, (HALF)0x3d78613e, (HALF)0x56af9e37, (HALF)0x2ac19031,
	(HALF)0xf
};
static CONST HALF h_rdefvec[] = {
	(HALF)0xd3e7b11e, (HALF)0x0f23fd7e, (HALF)0xa62546c9, (HALF)0x06c73ca5,
	(HALF)0x627d6265, (HALF)0x2f1e93a8, (HALF)0xff4b074e, (HALF)0x48b6d882,
	(HALF)0x7
};
#elif 2*FULL_BITS == 64
static CONST HALF h_ndefvec[] = {
	(HALF)0x18d9, (HALF)0xd6be, (HALF)0xa689, (HALF)0xba55,
	(HALF)0x27d9, (HALF)0x4a90, (HALF)0x4b64, (HALF)0x84f3,
	(HALF)0xd302, (HALF)0xd0e8, (HALF)0x613e, (HALF)0x3d78,
	(HALF)0x9e37, (HALF)0x56af, (HALF)0x9031, (HALF)0x2ac1,
	(HALF)0xf
};
static CONST HALF h_rdefvec[] = {
	(HALF)0xb11e, (HALF)0xd3e7, (HALF)0xfd7e, (HALF)0x0f23,
	(HALF)0x46c9, (HALF)0xa625, (HALF)0x3ca5, (HALF)0x06c7,
	(HALF)0x6265, (HALF)0x627d, (HALF)0x93a8, (HALF)0x2f1e,
	(HALF)0x074e, (HALF)0xff4b, (HALF)0xd882, (HALF)0x48b6,
	(HALF)0x7
};
#else
   /\../\	FULL_BITS must be 32 or 64	/\../\	 !!!
#endif
static CONST RANDOM init_blum = {1, 0, 8, (HALF)0, (HALF)0xff,
				 {(HALF *)h_ndefvec,
				  sizeof(h_ndefvec)/sizeof(HALF), 0},
				 {(HALF *)h_rdefvec,
				  sizeof(h_rdefvec)/sizeof(HALF), 0}};
#if FULL_BITS == 64
static CONST HALF h_rdefvec_2[] = {
	(HALF)0xd3e7b11e, (HALF)0x0f23fd7e, (HALF)0xa62546c9, (HALF)0x06c73ca5,
	(HALF)0x627d6265, (HALF)0x2f1e93a8, (HALF)0xff4b074e, (HALF)0x48b6d882,
	(HALF)0x7
};
#elif 2*FULL_BITS == 64
static CONST HALF h_rdefvec_2[] = {
	(HALF)0xb11e, (HALF)0xd3e7, (HALF)0xfd7e, (HALF)0x0f23,
	(HALF)0x46c9, (HALF)0xa625, (HALF)0x3ca5, (HALF)0x06c7,
	(HALF)0x6265, (HALF)0x627d, (HALF)0x93a8, (HALF)0x2f1e,
	(HALF)0x074e, (HALF)0xff4b, (HALF)0xd882, (HALF)0x48b6,
	(HALF)0x7
};
#else
   /\../\	FULL_BITS must be 32 or 64	/\../\	 !!!
#endif
static CONST ZVALUE z_rdefault = {
    (HALF *)h_rdefvec_2, sizeof(h_rdefvec_2)/sizeof(HALF), 0
};


/*
 * Pre-defined Blum generators
 *
 * These generators are seeded via the srandom(0, newn) where
 * 1 <= newn < BLUM_PREGEN.
 */
#if FULL_BITS == 64
static CONST HALF h_nvec01[] = {
	(HALF)0x83de9361, (HALF)0xf0db722d, (HALF)0x6fe328ca, (HALF)0x04944073,
	(HALF)0x5
};
static CONST HALF h_rvec01[] = {
	(HALF)0xa4cc42ec, (HALF)0x4e5dbb01, (HALF)0x11d952e7, (HALF)0xb226980f
};
static CONST HALF h_nvec02[] = {
	(HALF)0x353443f1, (HALF)0xeb286ea9, (HALF)0xdd374a18, (HALF)0x348a2555,
	(HALF)0x2c5
};
static CONST HALF h_rvec02[] = {
	(HALF)0x21e3a218, (HALF)0xe893616b, (HALF)0x6cd710e3, (HALF)0xf3d64344,
	(HALF)0x40
};
static CONST HALF h_nvec03[] = {
	(HALF)0x11d001f1, (HALF)0xf2ca661f, (HALF)0x3a81f1e0, (HALF)0x59d6ce4e,
	(HALF)0x0009cfd9
};
static CONST HALF h_rvec03[] = {
	(HALF)0xa0d7d76a, (HALF)0x3e142de2, (HALF)0xff5cea4f, (HALF)0xb44d9b64,
	(HALF)0xfae5
};
static CONST HALF h_nvec04[] = {
	(HALF)0xdfcc0751, (HALF)0x2decc680, (HALF)0x5df12a1a, (HALF)0x5c894ed7,
	(HALF)0x3070f924
};
static CONST HALF h_rvec04[] = {
	(HALF)0x4b984570, (HALF)0xa220ddba, (HALF)0xa2c0af8a, (HALF)0x131b2bdc,
	(HALF)0x0020c2d8
};
static CONST HALF h_nvec05[] = {
	(HALF)0x99166ef1, (HALF)0x8b99e5e7, (HALF)0x8769a010, (HALF)0x5d3fe111,
	(HALF)0x680bc2fa, (HALF)0x38f75aac, (HALF)0xdb81a85b, (HALF)0x109b1822,
	(HALF)0x2
};
static CONST HALF h_rvec05[] = {
	(HALF)0x59e2efa9, (HALF)0x0e6c77c8, (HALF)0x1e70aeed, (HALF)0x234f7b7d,
	(HALF)0x5f5df6db, (HALF)0xe821a960, (HALF)0xae33b792, (HALF)0x5e9b890e
};
static CONST HALF h_nvec06[] = {
	(HALF)0xe1ddf431, (HALF)0xd85557f1, (HALF)0x5ee732da, (HALF)0x3a38db77,
	(HALF)0x5c644026, (HALF)0xf2dbf218, (HALF)0x9ada2c79, (HALF)0x7bfd9d7d,
	(HALF)0xa
};
static CONST HALF h_rvec06[] = {
	(HALF)0xc9404daf, (HALF)0xc5dc2e80, (HALF)0x2c98eccf, (HALF)0xe1f3495d,
	(HALF)0xce1c925c, (HALF)0xe097aede, (HALF)0x88667154, (HALF)0x5e94a02f
};
static CONST HALF h_nvec07[] = {
	(HALF)0xcf9ec751, (HALF)0x602f9125, (HALF)0x52882e7f, (HALF)0x0dcf53ce,