1/*
2 * Licensed Materials - Property of IBM
3 *
4 * trousers - An open source TCG Software Stack
5 *
6 * (C) Copyright International Business Machines Corp. 2004, 2005
7 *
8 */
9
10#include <stdlib.h>
11#include <stdio.h>
12#include <string.h>
13
14#include "bi.h"
15#include "list.h"
16#include "tsplog.h"
17
18static unsigned long *primes;
19static int primes_length;
20
21/* Generates a random number of bit_length bit length. The first two bits and the last bit of this
22 * number are always set, therefore the number is odd and >= (2^(bit_length-1)+2^(bit_length-2)+1)
23 *
24 * bit_length: The length of the number to be generated, in bits
25 * return: a random number of bitLength bit length with first and last bits set
26 */
27void
28random_odd_bi(bi_ptr bi, int bit_length)
29{
30	if (bit_length > 0) {
31		bi_urandom(bi, bit_length);
32		//bi_generate_prime(bi, bit_length);
33		bi_setbit(bi, 0);
34		bi_setbit(bi, bit_length - 1);
35		bi_setbit(bi, bit_length - 2);
36	}
37}
38
39/* This method generates small prime numbers up to a specified bounds using the Sieve of
40 * Eratosthenes algorithm.
41 *
42 * prime_bound: the upper bound for the primes to be generated
43 * starting_prime: the first prime in the list of primes that is returned
44 * return: list of primes up to the specified bound. Each prime is of type bi_ptr
45 */
46void
47generate_small_primes(int prime_bound, int starting_prime)
48{
49	list_ptr res;
50	int length;
51	int *is_primes;
52	int i;
53	int k;
54	int prime;
55	node_t *current;
56
57	primes_length = 0;
58	res = list_new();
59	if (allocs == NULL) {
60		LogError("malloc of list failed");
61		return;
62	}
63	if ((prime_bound <= 1) || (starting_prime > prime_bound))
64		return;
65
66	if (starting_prime <= 2) {
67		starting_prime = 2;
68		list_add(res, bi_2);
69	}
70	length = (prime_bound - 1) >> 1; // length = (prime_bound -1) / 2;
71	is_primes = (int *)malloc(sizeof(int)*length);
72	if (is_primes == NULL) {
73		LogError("malloc of %zd bytes failed", sizeof(int) * length);
74		return;
75	}
76
77	for (i = 0; i < length; i++)
78		is_primes[i] = 1;
79
80	for (i = 0; i < length; i++) {
81		if (is_primes[i] == 1) {
82			prime = 2 * i + 3;
83			for (k = i + prime; k < length; k+= prime)
84				is_primes[k] = 0;
85
86			if (prime >= starting_prime) {
87				list_add(res, (void *)prime);
88				primes_length++;
89			}
90		}
91	}
92	// converti the list to a table
93	current = res->head; // go to first node
94	primes  = (unsigned long *)malloc(sizeof(unsigned long) * primes_length);
95	if (primes == NULL) {
96		LogError("malloc of %d bytes failed",
97			sizeof(unsigned long)*primes_length);
98		return;
99	}
100
101	i = 0;
102	while (current != NULL) {
103		primes[i++] = (unsigned long)current->obj;
104		current = current->next; // traverse through the list
105	}
106
107	free(is_primes);
108	list_freeall(res);
109}
110
111void
112prime_init()
113{
114	generate_small_primes(16384, 3);
115}
116
117/* Test whether the provided pDash or p = 2*pDash + 1 are divisible by any of the small primes
118 * saved in the listOfSmallPrimes. A limit for the largest prime to be tested against can be
119 * specified, but it will be ignored if it exeeds the number of precalculated primes.
120 *
121 * p_dash: the number to be tested (p_dash)
122 * prime_bound: the limit for the small primes to be tested against.
123 */
124static int
125test_small_prime_factors(const bi_ptr p_dash, const unsigned long prime_bound)
126{
127	int sievePassed = 1;
128	unsigned long r;
129	unsigned long small_prime;
130	bi_t temp; bi_new(temp);
131
132	small_prime = 1;
133	int i = 0;
134	while (i < primes_length && small_prime < prime_bound ) {
135		small_prime = primes[i++];
136		// r = p_dash % small_prime
137		bi_mod_si(temp, p_dash, small_prime);
138		r = bi_get_si(temp);
139		// test if pDash = 0 (mod smallPrime)
140		if (r == 0) {
141			sievePassed = 0;
142			break;
143		}
144		// test if p = 0 (mod smallPrime) (or r == smallPrime - r - 1)
145		if (r == (small_prime - r - 1)) {
146			sievePassed = 0;
147			break;
148		}
149	}
150	bi_free(temp);
151
152	return sievePassed;
153}
154
155/* Tests if a is a Miller-Rabin witness for n
156 *
157 * a: number which is supposed to be the witness
158 * n: number to be tested against
159 * return: true if a is Miller-Rabin witness for n, false otherwise
160 */
161int
162is_miller_rabin_witness(const bi_ptr a, const bi_ptr n)
163{
164	bi_t n_1;
165	bi_t temp;
166	bi_t _2_power_t;
167	bi_t u;
168	bi_t x0;
169	bi_t x1;
170	int t = -1;
171	int i;
172
173	bi_new(n_1);
174	bi_new(temp);
175	bi_new(_2_power_t);
176	bi_new(u);
177
178	// n1 = n - 1
179	bi_sub_si(n_1, n, 1);
180
181	// test if n-1 = 2^t*u with t >= 1 && u even
182	do {
183		t++;
184		// _2_power_t = bi_1 << t  ( ==  2 ^ t)
185		bi_shift_left(_2_power_t, bi_1, t);
186		// u = n_1 / (2 ^ t)
187		bi_div(u, n_1, _2_power_t);
188	} while (bi_equals_si(bi_mod(temp, u, bi_2), 0));
189
190	bi_new(x0);
191	bi_new(x1);
192
193	// x1 = (a ^ u ) % n
194	bi_mod_exp(x1, a, u, n);
195
196	// finished to use u, _2_power_t and temp
197	bi_free(u);
198	bi_free(_2_power_t);
199	bi_free(temp);
200	for (i = 0; i < t; i++) {
201		bi_set(x0, x1);
202
203		// x1 = (x0 ^ 2) % n
204		bi_mod_exp(x1, x0, bi_2, n);
205		if (bi_equals_si(x1, 1) && !bi_equals_si(x0, 1) && !bi_equals(x0, n_1) != 0) {
206			bi_free(x0);
207			bi_free(x1);
208			bi_free(n_1);
209			return 1;
210		}
211	}
212
213	bi_free(x0);
214	bi_free(x1);
215	bi_free(n_1);
216
217	if (!bi_equals(x1, bi_1))
218		return 1;
219
220	return 0;
221}
222
223bi_ptr
224compute_trivial_safe_prime(bi_ptr result, int bit_length)
225{
226	LogDebugFn("Enter");
227	do {
228		bi_generate_prime(result, bit_length-1);
229		bi_shift_left(result, result, 1);	// result := result << 1
230		bi_add_si(result, result, 1);	// result := result -1
231		if (getenv("TSS_DEBUG_OFF") == NULL) {
232			printf(".");
233			fflush(stdout);
234		}
235	} while (bi_is_probable_prime(result)==0);
236
237	return result;
238}
239
240/* The main method to compute a random safe prime of the specified bit length.
241 * IMPORTANT: The computer prime will have two first bits and the last bit set to 1 !!
242 * i.e. > (2^(bitLength-1)+2^(bitLength-2)+1). This is done to be sure that if two primes of
243 * bitLength n are multiplied, the result will have the bitLenght of 2*n exactly This
244 * implementation uses the algorithm proposed by Ronald Cramer and Victor Shoup in "Signature
245 * Schemes Based on the strong RSA Assumption" May 9, 2000.
246 *
247 * bitLength: the bit length of the safe prime to be computed.
248 * return: a number which is considered to be safe prime
249 */
250bi_ptr
251compute_safe_prime(bi_ptr p, int bit_length)
252{
253	bi_ptr p_dash;
254	bi_ptr temp_p;
255	bi_ptr p_minus_1;
256	int stop;
257	unsigned long prime_bound;
258
259	LogDebug("compute Safe Prime: length: %d bits\n", bit_length);
260
261	p_dash = bi_new_ptr();
262	temp_p = bi_new_ptr();
263	p_minus_1 = bi_new_ptr();
264
265	/* some heuristic checks to limit the number of small primes to check against and the
266	 * number of Miller-Rabin primality tests at the end */
267	if (bit_length <= 256) {
268		prime_bound = 768;
269	} else if (bit_length <= 512) {
270		prime_bound = 3072;
271	} else if (bit_length <= 768) {
272		prime_bound = 6144;
273	} else if (bit_length <= 1024) {
274		prime_bound = 1024;
275	} else {
276		prime_bound = 16384;
277	}
278
279	do {
280		stop = 0;
281		/* p_dash = generated random with basic bit settings (odd) */
282		random_odd_bi(p_dash, bit_length - 1);
283
284		if (test_small_prime_factors(p_dash, prime_bound) == 0)  {
285			LogDebugFn("1");
286			continue;
287		}
288		/* test if p_dash or p are divisible by some small primes */
289		if (is_miller_rabin_witness(bi_2, p_dash)) {
290			LogDebugFn("2");
291			continue;
292		}
293		/* test if 2^(pDash) = +1/-1 (mod p)
294		 * bi can not handle negative operation, we compare to (p-1) instead of -1
295		 * calculate p = 2*pDash+1 -> (pDash << 1) + 1
296		 */
297		bi_shift_left(p, p_dash, 1);
298		bi_add(p, p, bi_1);
299
300		// p_minus_1:= p - 1
301		bi_sub(p_minus_1, p, bi_1);
302
303		//  temp_p := ( 2 ^ p_dash ) mod p
304		bi_mod_exp(temp_p, bi_2, p_dash, p);
305		if (!bi_equals_si(temp_p, 1)  && !bi_equals(temp_p, p_minus_1) ) {
306			LogDebugFn("3");
307			continue;
308		}
309
310		// test if pDash or p are divisible by some small primes
311		if (is_miller_rabin_witness(bi_2, p_dash)) {
312			LogDebugFn("4");
313			continue;
314		}
315		// test the library dependent probable_prime
316		if (bi_is_probable_prime(p_dash))
317			stop = 1;
318	} while (stop == 0);
319
320	bi_free(p_minus_1);
321	bi_free(temp_p);
322	bi_free(p_dash);
323
324	LogDebug("found Safe Prime: %s bits", bi_2_hex_char(p));
325
326	return p;
327}
328