calc/lib/seedrandom.cal

/*
 * Copyright (c) 1996 Landon Curt Noll
 *
 * 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.
 *
 *	Landon Curt Noll
 *	http://reality.sgi.com/chongo/
 *
 *	chongo <was here> /\../\
 */

/*
 * seedrandom - seed the cryptographically strong Blum generator
 *
 *
 * The period of a  Blum generators with modulus 'n=p*q' (where p and
 * q are primes 3 mod 4) is:
 *
 *	lambda(n) = lcm(factors of p-1 & q-1)
 *
 * One can construct a generator with a maximal period when
 * 'p' and 'q' have the fewest possible factors in common.
 * The quickest way to select such primes is only use 'p'
 * and 'q' when '(p-1)/2' and '(q-1)/2' are both primes.
 * This function will seed the random() generator that uses
 * such primes.
 *
 * given:
 *	seed1 - a large random value (at least 10^20 and perhaps < 10^314)
 *	seed2 - a large random value (at least 10^20 and perhaps < 10^314)
 *	size - min Blum modulus as a power of 2 (at least 32, perhaps >= 512)
 *	trials - number of ptest() trials (default 25)
 *
 * returns:
 *	the previous random state
 *
 * NOTE: The [10^20, 10^314) range comes from the fact that the 13th internal
 *	 modulus is ~10^315.  We want the lower bound seed to be reasonably big.
 */
define seedrandom(seed1, seed2, size, trials)
{
	local p;		/* first Blum prime */
	local fp;		/* prime co-factor of p-1 */
	local sp;		/* min bit size of p */
	local q;		/* second Blum prime */
	local fq;		/* prime co-factor of q-1 */
	local sq;		/* min bit size of q */
	local n;		/* Blum modulus */
	local binsize;		/* smallest power of 2 > n=p*q */
	local r;		/* initial quadratic residue */
	local random_state;	/* the initial rand state */
	local random_junk;	/* rand state that is not needed */
	local old_state;	/* old random state to return */

	/*
	 * firewall
	 */
	if (!isint(seed1)) {
		quit "1st arg (seed1) is not an int";
	}
	if (!isint(seed2)) {
		quit "2nd arg (seed2) is not an int";
	}
	if (!isint(size)) {
		quit "3rd arg (size) is not an int";
	}
	if (!isint(trials)) {
		trials = 25;
	}
	if (digits(seed1) <= 20) {
		quit "1st arg (seed1) must be > 10^20 and perhaps < 10^314";
	}
	if (digits(seed2) <= 20) {
		quit "2nd arg (seed2) must be > 10^20 and perhaps < 10^314";
	}
	if (size < 32) {
		quit "3rd arg (size) needs to be >= 32 (perhaps >= 512)";
	}
	if (trials < 1) {
		quit "4th arg (trials) must be > 0";
	}

	/*
	 * determine the search parameters
	 */
	++size;		/* convert power of 2 to bit length */
	sp = int((size/2)-(size*0.03)+1);
	sq = size - sp;

	/*
	 * find the first Blum prime
	 */
	random_state = srandom(seed1, 13);
	do {
		do {
			fp = nextcand(2^sp+randombit(sp), 1, 1, 3, 4);
			p = 2*fp+1;
		} while (ptest(p,1,0) == 0);
	} while(ptest(p, trials) == 0 || ptest(fp, trials) == 0);
	if (config("lib_debug") & 3) {
		print "/* 1st Blum prime */ p=", p;
	}

	/*
	 * find the 2nd Blum prime
	 */
	random_junk = srandom(seed2, 13);
	do {
		do {
			fq = nextcand(2^sq+randombit(sq), 1, 1, 3, 4);
			q = 2*fq+1;
		} while (ptest(q,1,0) == 0);
	} while(ptest(q, trials) == 0 || ptest(fq, trials) == 0);
	if (config("lib_debug") & 3) {
		print "/* 2nd Blum prime */ q=", q;
	}

	/*
	 * seed the Blum generator
	 */
	n = p*q;				/* the Blum modulus */
	binsize = highbit(n)+1;			/* smallest power of 2 > p*q */
	r = pmod(rand(1<<ceil(binsize*4/5), 1<<(binsize-2)), 2, n);
	if (config("lib_debug") & 3) {
		print "/* seed quadratic residue */ r=", r;
		print "/* newn", binsize, "bit quadratic residue*/ newn=", n;
	}
	old_state = srandom(r, n);

	/*
	 * restore other states that we altered
	 */
	random_junk = srandom(random_state);

	/*
	 * return the previous random state
	 */
	return old_state;
}

if (config("lib_debug") & 3) {
    print "seedrandom(seed1, seed2, size [, trials]) defined";
}