/*
 * seedrandom - seed the cryptographically strong Blum generator
 *
 * Copyright (C) 1999  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.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * @(#) $Revision: 30.1 $
 * @(#) $Id: seedrandom.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/seedrandom.cal,v $
 *
 * Under source code control:	1996/01/01 08:21:00
 * File existed as early as:	1996
 *
 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * The period of 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("resource_debug") & 8) {
		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("resource_debug") & 8) {
		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("resource_debug") & 8) {
		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("resource_debug") & 3) {
    print "seedrandom(seed1, seed2, size [, trials]) defined";
}