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 */ 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