1 /*
2 * Copyright (C) 2017 - This file is part of libecc project
3 *
4 * Authors:
5 * Ryad BENADJILA <ryadbenadjila@gmail.com>
6 * Arnaud EBALARD <arnaud.ebalard@ssi.gouv.fr>
7 * Jean-Pierre FLORI <jean-pierre.flori@ssi.gouv.fr>
8 *
9 * Contributors:
10 * Nicolas VIVET <nicolas.vivet@ssi.gouv.fr>
11 * Karim KHALFALLAH <karim.khalfallah@ssi.gouv.fr>
12 *
13 * This software is licensed under a dual BSD and GPL v2 license.
14 * See LICENSE file at the root folder of the project.
15 */
16 /* We include the NN layer API header */
17 #include <libecc/libarith.h>
18
19 ATTRIBUTE_WARN_UNUSED_RET int miller_rabin(nn_src_t n, const unsigned int t, int *res);
20
21 /* Miller-Rabin primality test.
22 * See "Handbook of Applied Cryptography", alorithm 4.24:
23 *
24 * Algorithm: Miller-Rabin probabilistic primality test
25 * MILLER-RABIN(n,t)
26 * INPUT: an odd integer n ≥ 3 and security parameter t ≥ 1.
27 * OUTPUT: an answer “prime” or “composite” to the question: “Is n prime?”
28 * 1. Write n − 1 = 2**s x r such that r is odd.
29 * 2. For i from 1 to t do the following:
30 * 2.1 Choose a random integer a, 2 ≤ a ≤ n − 2.
31 * 2.2 Compute y = a**r mod n using Algorithm 2.143.
32 * 2.3 If y != 1 and y != n − 1 then do the following:
33 * j←1.
34 * While j ≤ s − 1 and y != n − 1 do the following:
35 * Compute y←y2 mod n.
36 * If y = 1 then return(“composite”).
37 * j←j + 1.
38 * If y != n − 1 then return (“composite”).
39 * 3. Return(“maybe prime”).
40 *
41 * The Miller-Rabin test can give false positives when
42 * answering "maybe prime", but is always right when answering
43 * "composite".
44 */
miller_rabin(nn_src_t n,const unsigned int t,int * res)45 int miller_rabin(nn_src_t n, const unsigned int t, int *res)
46 {
47 int ret, iszero, cmp, isodd, cmp1, cmp2;
48 unsigned int i;
49 bitcnt_t k;
50 /* Temporary NN variables */
51 nn s, q, r, d, a, y, j, one, two, tmp;
52 s.magic = q.magic = r.magic = d.magic = a.magic = y.magic = j.magic = WORD(0);
53 one.magic = two.magic = tmp.magic = WORD(0);
54
55 ret = nn_check_initialized(n); EG(ret, err);
56 MUST_HAVE((res != NULL), ret, err);
57 (*res) = 0;
58
59 /* Initialize our local NN variables */
60 ret = nn_init(&s, 0); EG(ret, err);
61 ret = nn_init(&q, 0); EG(ret, err);
62 ret = nn_init(&r, 0); EG(ret, err);
63 ret = nn_init(&d, 0); EG(ret, err);
64 ret = nn_init(&a, 0); EG(ret, err);
65 ret = nn_init(&y, 0); EG(ret, err);
66 ret = nn_init(&j, 0); EG(ret, err);
67 ret = nn_init(&one, 0); EG(ret, err);
68 ret = nn_init(&two, 0); EG(ret, err);
69 ret = nn_init(&tmp, 0); EG(ret, err);
70
71 /* Security parameter t must be >= 1 */
72 MUST_HAVE((t >= 1), ret, err);
73
74 /* one = 1 */
75 ret = nn_one(&one); EG(ret, err);
76 /* two = 2 */
77 ret = nn_set_word_value(&two, WORD(2)); EG(ret, err);
78
79 /* If n = 0, this is not a prime */
80 ret = nn_iszero(n, &iszero); EG(ret, err);
81 if (iszero) {
82 ret = 0;
83 (*res) = 0;
84 goto err;
85 }
86 /* If n = 1, this is not a prime */
87 ret = nn_cmp(n, &one, &cmp); EG(ret, err);
88 if (cmp == 0) {
89 ret = 0;
90 (*res) = 0;
91 goto err;
92 }
93 /* If n = 2, this is a prime number */
94 ret = nn_cmp(n, &two, &cmp); EG(ret, err);
95 if (cmp == 0) {
96 ret = 0;
97 (*res) = 1;
98 goto err;
99 }
100 /* If n = 3, this is a prime number */
101 ret = nn_copy(&tmp, n); EG(ret, err);
102 ret = nn_dec(&tmp, &tmp); EG(ret, err);
103 ret = nn_cmp(&tmp, &two, &cmp); EG(ret, err);
104 if (cmp == 0) {
105 ret = 0;
106 (*res) = 1;
107 goto err;
108 }
109
110 /* If n >= 4 is even, this is not a prime */
111 ret = nn_isodd(n, &isodd); EG(ret, err);
112 if (!isodd) {
113 ret = 0;
114 (*res) = 0;
115 goto err;
116 }
117
118 /* n − 1 = 2^s x r, repeatedly try to divide n-1 by 2 */
119 /* s = 0 and r = n-1 */
120 ret = nn_zero(&s); EG(ret, err);
121 ret = nn_copy(&r, n); EG(ret, err);
122 ret = nn_dec(&r, &r); EG(ret, err);
123 while (1) {
124 ret = nn_divrem(&q, &d, &r, &two); EG(ret, err);
125 ret = nn_inc(&s, &s); EG(ret, err);
126 ret = nn_copy(&r, &q); EG(ret, err);
127 /* If r is odd, we have finished our division */
128 ret = nn_isodd(&r, &isodd); EG(ret, err);
129 if (isodd) {
130 break;
131 }
132 }
133 /* 2. For i from 1 to t do the following: */
134 for (i = 1; i <= t; i++) {
135 bitcnt_t blen;
136 /* 2.1 Choose a random integer a, 2 ≤ a ≤ n − 2 */
137 ret = nn_copy(&tmp, n); EG(ret, err);
138 ret = nn_dec(&tmp, &tmp); EG(ret, err);
139 ret = nn_zero(&a); EG(ret, err);
140 ret = nn_cmp(&a, &two, &cmp); EG(ret, err);
141 while (cmp < 0) {
142 ret = nn_get_random_mod(&a, &tmp); EG(ret, err);
143 ret = nn_cmp(&a, &two, &cmp); EG(ret, err);
144 }
145 /* A very loose (and NOT robust) implementation of
146 * modular exponentiation with square and multiply
147 * to compute y = a**r (n)
148 * WARNING: NOT to be used in production code!
149 */
150 ret = nn_one(&y); EG(ret, err);
151 ret = nn_bitlen(&r, &blen); EG(ret, err);
152 for (k = 0; k < blen; k++) {
153 u8 bit;
154 ret = nn_getbit(&r, k, &bit); EG(ret, err);
155 if (bit) {
156 /* Warning: the multiplication is not modular, we
157 * have to take care of our size here!
158 */
159 MUST_HAVE((NN_MAX_BIT_LEN >=
160 (WORD_BITS * (y.wlen + a.wlen))), ret, err);
161 ret = nn_mul(&y, &y, &a); EG(ret, err);
162 ret = nn_mod(&y, &y, n); EG(ret, err);
163 }
164 MUST_HAVE((NN_MAX_BIT_LEN >= (2 * WORD_BITS * a.wlen)), ret, err);
165 ret = nn_sqr(&a, &a); EG(ret, err);
166 ret = nn_mod(&a, &a, n); EG(ret, err);
167 }
168 /* 2.3 If y != 1 and y != n − 1 then do the following
169 * Note: tmp still contains n - 1 here.
170 */
171 ret = nn_cmp(&y, &one, &cmp1); EG(ret, err);
172 ret = nn_cmp(&y, &tmp, &cmp2); EG(ret, err);
173 if ((cmp1 != 0) && (cmp2 != 0)) {
174 /* j←1. */
175 ret = nn_one(&j); EG(ret, err);
176 /* While j ≤ s − 1 and y != n − 1 do the following: */
177 ret = nn_cmp(&j, &s, &cmp1); EG(ret, err);
178 ret = nn_cmp(&y, &tmp, &cmp2); EG(ret, err);
179 while ((cmp1 < 0) && (cmp2 != 0)) {
180 /* Compute y←y2 mod n. */
181 MUST_HAVE((NN_MAX_BIT_LEN >=
182 (2 * WORD_BITS * y.wlen)), ret, err);
183 ret = nn_sqr(&y, &y); EG(ret, err);
184 ret = nn_mod(&y, &y, n); EG(ret, err);
185 /* If y = 1 then return(“composite”). */
186 ret = nn_cmp(&y, &one, &cmp); EG(ret, err);
187 if (cmp == 0) {
188 ret = 0;
189 (*res) = 0;
190 goto err;
191 }
192 /* j←j + 1. */
193 ret = nn_inc(&j, &j); EG(ret, err);
194 ret = nn_cmp(&j, &s, &cmp1); EG(ret, err);
195 ret = nn_cmp(&y, &tmp, &cmp2); EG(ret, err);
196 }
197 /* If y != n − 1 then return (“composite”). */
198 ret = nn_cmp(&y, &tmp, &cmp); EG(ret, err);
199 if (cmp != 0) {
200 ret = 0;
201 (*res) = 0;
202 goto err;
203 }
204 }
205 /* 3. Return(“maybe prime”). */
206 ret = 0;
207 (*res) = 1;
208 }
209
210 err:
211 nn_uninit(&s);
212 nn_uninit(&q);
213 nn_uninit(&r);
214 nn_uninit(&d);
215 nn_uninit(&a);
216 nn_uninit(&y);
217 nn_uninit(&j);
218 nn_uninit(&one);
219 nn_uninit(&two);
220 nn_uninit(&tmp);
221
222 return ret;
223 }
224