1/*
2 * seedrandom - seed the cryptographically strong Blum generator
3 *
4 * Copyright (C) 1999  Landon Curt Noll
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1996/01/01 08:21:00
21 * File existed as early as:	1996
22 *
23 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
24 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
25 */
26
27/*
28 * The period of Blum generators with modulus 'n=p*q' (where p and
29 * q are primes 3 mod 4) is:
30 *
31 *	lambda(n) = lcm(factors of p-1 & q-1)
32 *
33 * One can construct a generator with a maximal period when
34 * 'p' and 'q' have the fewest possible factors in common.
35 * The quickest way to select such primes is only use 'p'
36 * and 'q' when '(p-1)/2' and '(q-1)/2' are both primes.
37 * This function will seed the random() generator that uses
38 * such primes.
39 *
40 * given:
41 *	seed1 - a large random value (at least 10^20 and perhaps < 10^314)
42 *	seed2 - a large random value (at least 10^20 and perhaps < 10^314)
43 *	size - min Blum modulus as a power of 2 (at least 32, perhaps >= 512)
44 *	trials - number of ptest() trials (default 25)
45 *
46 * returns:
47 *	the previous random state
48 *
49 * NOTE: The [10^20, 10^314) range comes from the fact that the 13th internal
50 *	 modulus is ~10^315.  We want the lower bound seed to be reasonably big.
51 */
52
53
54define seedrandom(seed1, seed2, size, trials)
55{
56	local p;		/* first Blum prime */
57	local fp;		/* prime co-factor of p-1 */
58	local sp;		/* min bit size of p */
59	local q;		/* second Blum prime */
60	local fq;		/* prime co-factor of q-1 */
61	local sq;		/* min bit size of q */
62	local n;		/* Blum modulus */
63	local binsize;		/* smallest power of 2 > n=p*q */
64	local r;		/* initial quadratic residue */
65	local random_state;	/* the initial rand state */
66	local random_junk;	/* rand state that is not needed */
67	local old_state;	/* old random state to return */
68
69	/*
70	 * firewall
71	 */
72	if (!isint(seed1)) {
73		quit "1st arg (seed1) is not an int";
74	}
75	if (!isint(seed2)) {
76		quit "2nd arg (seed2) is not an int";
77	}
78	if (!isint(size)) {
79		quit "3rd arg (size) is not an int";
80	}
81	if (!isint(trials)) {
82		trials = 25;
83	}
84	if (digits(seed1) <= 20) {
85		quit "1st arg (seed1) must be > 10^20 and perhaps < 10^314";
86	}
87	if (digits(seed2) <= 20) {
88		quit "2nd arg (seed2) must be > 10^20 and perhaps < 10^314";
89	}
90	if (size < 32) {
91		quit "3rd arg (size) needs to be >= 32 (perhaps >= 512)";
92	}
93	if (trials < 1) {
94		quit "4th arg (trials) must be > 0";
95	}
96
97	/*
98	 * determine the search parameters
99	 */
100	++size;		/* convert power of 2 to bit length */
101	sp = int((size/2)-(size*0.03)+1);
102	sq = size - sp;
103
104	/*
105	 * find the first Blum prime
106	 */
107	random_state = srandom(seed1, 13);
108	do {
109		do {
110			fp = nextcand(2^sp+randombit(sp), 1, 1, 3, 4);
111			p = 2*fp+1;
112		} while (ptest(p,1,0) == 0);
113	} while(ptest(p, trials) == 0 || ptest(fp, trials) == 0);
114	if (config("resource_debug") & 8) {
115		print "/* 1st Blum prime */ p=", p;
116	}
117
118	/*
119	 * find the 2nd Blum prime
120	 */
121	random_junk = srandom(seed2, 13);
122	do {
123		do {
124			fq = nextcand(2^sq+randombit(sq), 1, 1, 3, 4);
125			q = 2*fq+1;
126		} while (ptest(q,1,0) == 0);
127	} while(ptest(q, trials) == 0 || ptest(fq, trials) == 0);
128	if (config("resource_debug") & 8) {
129		print "/* 2nd Blum prime */ q=", q;
130	}
131
132	/*
133	 * seed the Blum generator
134	 */
135	n = p*q;				/* the Blum modulus */
136	binsize = highbit(n)+1;			/* smallest power of 2 > p*q */
137	r = pmod(rand(1<<ceil(binsize*4/5), 1<<(binsize-2)), 2, n);
138	if (config("resource_debug") & 8) {
139		print "/* seed quadratic residue */ r=", r;
140		print "/* newn", binsize, "bit quadratic residue*/ newn=", n;
141	}
142	old_state = srandom(r, n);
143
144	/*
145	 * restore other states that we altered
146	 */
147	random_junk = srandom(random_state);
148
149	/*
150	 * return the previous random state
151	 */
152	return old_state;
153}
154
155if (config("resource_debug") & 3) {
156    print "seedrandom(seed1, seed2, size [, trials]) defined";
157}
158