xref: /plan9/sys/src/libsec/port/genstrongprime.c (revision 80ee5cbfe36716af62da8896207e9763b8e3d760)
1*80ee5cbfSDavid du Colombier #include "os.h"
2*80ee5cbfSDavid du Colombier #include <mp.h>
3*80ee5cbfSDavid du Colombier #include <libsec.h>
4*80ee5cbfSDavid du Colombier 
5*80ee5cbfSDavid du Colombier // Gordon's algorithm for generating a strong prime
6*80ee5cbfSDavid du Colombier //	Menezes et al () Handbook, p.150
7*80ee5cbfSDavid du Colombier void
genstrongprime(mpint * p,int n,int accuracy)8*80ee5cbfSDavid du Colombier genstrongprime(mpint *p, int n, int accuracy)
9*80ee5cbfSDavid du Colombier {
10*80ee5cbfSDavid du Colombier 	mpint *s, *t, *r, *i;
11*80ee5cbfSDavid du Colombier 
12*80ee5cbfSDavid du Colombier 	if(n < 64)
13*80ee5cbfSDavid du Colombier 		n = 64;
14*80ee5cbfSDavid du Colombier 
15*80ee5cbfSDavid du Colombier 	s = mpnew(n/2);
16*80ee5cbfSDavid du Colombier 	genprime(s, (n/2)-16, accuracy);
17*80ee5cbfSDavid du Colombier 	t = mpnew(n/2);
18*80ee5cbfSDavid du Colombier 	genprime(t, n-mpsignif(s)-32, accuracy);
19*80ee5cbfSDavid du Colombier 
20*80ee5cbfSDavid du Colombier 	// first r = 2it + 1 that's prime
21*80ee5cbfSDavid du Colombier 	i = mpnew(16);
22*80ee5cbfSDavid du Colombier 	r = mpnew(0);
23*80ee5cbfSDavid du Colombier 	itomp(0x8000, i);
24*80ee5cbfSDavid du Colombier 	mpleft(t, 1, t);	// 2t
25*80ee5cbfSDavid du Colombier 	mpmul(i, t, r);		// 2it
26*80ee5cbfSDavid du Colombier 	mpadd(r, mpone, r);	// 2it + 1
27*80ee5cbfSDavid du Colombier 	for(;;){
28*80ee5cbfSDavid du Colombier 		if(probably_prime(r, 18))
29*80ee5cbfSDavid du Colombier 			break;
30*80ee5cbfSDavid du Colombier 		mpadd(r, t, r);	// r += 2t
31*80ee5cbfSDavid du Colombier 	}
32*80ee5cbfSDavid du Colombier 
33*80ee5cbfSDavid du Colombier 	// p0 = 2(s**(r-2) mod r)s - 1
34*80ee5cbfSDavid du Colombier 	itomp(2, p);
35*80ee5cbfSDavid du Colombier 	mpsub(r, p, p);
36*80ee5cbfSDavid du Colombier 	mpexp(s, p, r, p);
37*80ee5cbfSDavid du Colombier 	mpmul(s, p, p);
38*80ee5cbfSDavid du Colombier 	mpleft(p, 1, p);
39*80ee5cbfSDavid du Colombier 	mpsub(p, mpone, p);
40*80ee5cbfSDavid du Colombier 
41*80ee5cbfSDavid du Colombier 	// first p = p0 + 2irs that's prime
42*80ee5cbfSDavid du Colombier 	itomp(0x8000, i);
43*80ee5cbfSDavid du Colombier 	mpleft(r, 1, r);	// 2r
44*80ee5cbfSDavid du Colombier 	mpmul(r, s, r);		// 2rs
45*80ee5cbfSDavid du Colombier 	mpmul(r, i, i);		// 2irs
46*80ee5cbfSDavid du Colombier 	mpadd(p, i, p);		// p0 + 2irs
47*80ee5cbfSDavid du Colombier 	for(;;){
48*80ee5cbfSDavid du Colombier 		if(probably_prime(p, accuracy))
49*80ee5cbfSDavid du Colombier 			break;
50*80ee5cbfSDavid du Colombier 		mpadd(p, r, p); // p += 2rs
51*80ee5cbfSDavid du Colombier 	}
52*80ee5cbfSDavid du Colombier 
53*80ee5cbfSDavid du Colombier 	mpfree(i);
54*80ee5cbfSDavid du Colombier 	mpfree(s);
55*80ee5cbfSDavid du Colombier 	mpfree(r);
56*80ee5cbfSDavid du Colombier 	mpfree(t);
57*80ee5cbfSDavid du Colombier }
58