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 /* 17 * The purpose of this example is to implement Pollard's rho 18 * algorithm to find non-trivial factors of a composite natural 19 * number. 20 * The prime numbers decomposition of the natural number is 21 * recovered through repeated Pollard's rho. Primality checking 22 * is performed using a Miller-Rabin test. 23 * 24 * WARNING: the code in this example is only here to illustrate 25 * how to use the NN layer API. This code has not been designed 26 * for production purposes (e.g. no effort has been made to make 27 * it constant time). 28 * 29 * 30 */ 31 32 /* We include the NN layer API header */ 33 #include <libecc/libarith.h> 34 35 /* Declare our Miller-Rabin test implemented 36 * in another module. 37 */ 38 ATTRIBUTE_WARN_UNUSED_RET int miller_rabin(nn_src_t n, const unsigned int t, int *check); 39 40 ATTRIBUTE_WARN_UNUSED_RET int pollard_rho(nn_t d, nn_src_t n, const word_t c); 41 /* Pollard's rho main function, as described in 42 * "Handbook of Applied Cryptography". 43 * 44 * Pollard's rho: 45 * ============== 46 * See "Handbook of Applied Cryptography", alorithm 3.9: 47 * 48 * Algorithm Pollard’s rho algorithm for factoring integers 49 * INPUT: a composite integer n that is not a prime power. 50 * OUTPUT: a non-trivial factor d of n. 51 * 1. Set a←2, b←2. 52 * 2. For i = 1, 2, ... do the following: 53 * 2.1 Compute a←a^2 + 1 mod n, b←b^2 + 1 mod n, b←b^2 + 1 mod n. 54 * 2.2 Compute d = gcd(a − b, n). 55 * 2.3 If 1 < d < n then return(d) and terminate with success. 56 * 2.4 If d = n then terminate the algorithm with failure (see Note 3.12). 57 */ 58 int pollard_rho(nn_t d, nn_src_t n, const word_t c) 59 { 60 int ret, cmp, cmp1, cmp2; 61 /* Temporary a and b variables */ 62 nn a, b, tmp, one, c_bignum; 63 a.magic = b.magic = tmp.magic = one.magic = c_bignum.magic = WORD(0); 64 65 /* Initialize variables */ 66 ret = nn_init(&a, 0); EG(ret, err); 67 ret = nn_init(&b, 0); EG(ret, err); 68 ret = nn_init(&tmp, 0); EG(ret, err); 69 ret = nn_init(&one, 0); EG(ret, err); 70 ret = nn_init(&c_bignum, 0); EG(ret, err); 71 ret = nn_init(d, 0); EG(ret, err); 72 73 MUST_HAVE((c > 0), ret, err); 74 75 /* Zeroize the output */ 76 ret = nn_zero(d); EG(ret, err); 77 ret = nn_one(&one); EG(ret, err); 78 /* 1. Set a←2, b←2. */ 79 ret = nn_set_word_value(&a, WORD(2)); EG(ret, err); 80 ret = nn_set_word_value(&b, WORD(2)); EG(ret, err); 81 ret = nn_set_word_value(&c_bignum, c); EG(ret, err); 82 83 /* For i = 1, 2, . . . do the following: */ 84 while (1) { 85 int sign; 86 /* 2.1 Compute a←a^2 + c mod n */ 87 ret = nn_sqr(&a, &a); EG(ret, err); 88 ret = nn_add(&a, &a, &c_bignum); EG(ret, err); 89 ret = nn_mod(&a, &a, n); EG(ret, err); 90 /* 2.1 Compute b←b^2 + c mod n twice in a row */ 91 ret = nn_sqr(&b, &b); EG(ret, err); 92 ret = nn_add(&b, &b, &c_bignum); EG(ret, err); 93 ret = nn_mod(&b, &b, n); EG(ret, err); 94 ret = nn_sqr(&b, &b); EG(ret, err); 95 ret = nn_add(&b, &b, &c_bignum); EG(ret, err); 96 ret = nn_mod(&b, &b, n); EG(ret, err); 97 /* 2.2 Compute d = gcd(a − b, n) */ 98 ret = nn_cmp(&a, &b, &cmp); EG(ret, err); 99 if (cmp >= 0) { 100 ret = nn_sub(&tmp, &a, &b); EG(ret, err); 101 } else { 102 ret = nn_sub(&tmp, &b, &a); EG(ret, err); 103 } 104 ret = nn_gcd(d, &tmp, n, &sign); EG(ret, err); 105 ret = nn_cmp(d, n, &cmp1); EG(ret, err); 106 ret = nn_cmp(d, &one, &cmp2); EG(ret, err); 107 if ((cmp1 < 0) && (cmp2 > 0)) { 108 ret = 0; 109 goto err; 110 } 111 ret = nn_cmp(d, n, &cmp); EG(ret, err); 112 if (cmp == 0) { 113 ret = -1; 114 goto err; 115 } 116 } 117 err: 118 /* Uninitialize local variables */ 119 nn_uninit(&a); 120 nn_uninit(&b); 121 nn_uninit(&tmp); 122 nn_uninit(&one); 123 nn_uninit(&c_bignum); 124 125 return ret; 126 } 127 128 ATTRIBUTE_WARN_UNUSED_RET int find_divisors(nn_src_t in); 129 /* Maximum number of divisors we support */ 130 #define MAX_DIVISORS 10 131 /* Function to find prime divisors of the NN input */ 132 int find_divisors(nn_src_t in) 133 { 134 int n_divisors_found, i, found, ret, check, cmp; 135 nn n; 136 nn divisors[MAX_DIVISORS]; 137 word_t c; 138 139 n.magic = WORD(0); 140 for(i = 0; i < MAX_DIVISORS; i++){ 141 divisors[i].magic = WORD(0); 142 } 143 144 ret = nn_check_initialized(in); EG(ret, err); 145 146 ext_printf("=================\n"); 147 nn_print("Finding factors of:", in); 148 149 /* First, check primality of the input */ 150 ret = miller_rabin(in, 10, &check); EG(ret, err); 151 if (check) { 152 ext_printf("The number is probably prime, leaving ...\n"); 153 ret = -1; 154 goto err; 155 } 156 ext_printf("The number is composite, performing Pollard's rho\n"); 157 158 ret = nn_init(&n, 0); EG(ret, err); 159 ret = nn_copy(&n, in); EG(ret, err); 160 for (i = 0; i < MAX_DIVISORS; i++) { 161 ret = nn_init(&(divisors[i]), 0); EG(ret, err); 162 } 163 164 n_divisors_found = 0; 165 c = 0; 166 while (1) { 167 c++; 168 ret = pollard_rho(&(divisors[n_divisors_found]), &n, c); 169 if (ret) { 170 continue; 171 } 172 found = 0; 173 for (i = 0; i < n_divisors_found; i++) { 174 ret = nn_cmp(&(divisors[n_divisors_found]), &(divisors[i]), &cmp); EG(ret, err); 175 if (cmp == 0) { 176 found = 1; 177 } 178 } 179 if (found == 0) { 180 nn q, r; 181 ret = nn_init(&q, 0); EG(ret, err); 182 ret = nn_init(&r, 0); EG(ret, err); 183 ext_printf("Pollard's rho succeded\n"); 184 nn_print("d:", &(divisors[n_divisors_found])); 185 /* 186 * Now we can launch the algorithm again on n / d 187 * to find new divisors. If n / d is prime, we are done! 188 */ 189 ret = nn_divrem(&q, &r, &n, &(divisors[n_divisors_found])); EG(ret, err); 190 /* 191 * Check n / d primality with Miller-Rabin (security 192 * parameter of 10) 193 */ 194 ret = miller_rabin(&q, 10, &check); EG(ret, err); 195 if (check == 1) { 196 nn_print("Last divisor is prime:", &q); 197 nn_uninit(&q); 198 nn_uninit(&r); 199 break; 200 } 201 nn_print("Now performing Pollard's rho on:", &q); 202 ret = nn_copy(&n, &q); EG(ret, err); 203 nn_uninit(&q); 204 nn_uninit(&r); 205 c = 0; 206 n_divisors_found++; 207 if (n_divisors_found == MAX_DIVISORS) { 208 ext_printf 209 ("Max divisors reached, leaving ...\n"); 210 break; 211 } 212 } 213 } 214 ret = 0; 215 216 err: 217 ext_printf("=================\n"); 218 nn_uninit(&n); 219 for (i = 0; i < MAX_DIVISORS; i++) { 220 nn_uninit(&(divisors[i])); 221 } 222 return ret; 223 } 224 225 #ifdef NN_EXAMPLE 226 int main(int argc, char *argv[]) 227 { 228 int ret; 229 230 /* Fermat F5 = 2^32 + 1 = 641 x 6700417 */ 231 const unsigned char fermat_F5[] = { 0x01, 0x00, 0x00, 0x00, 0x01 }; 232 /* Fermat F6 = 2^64 + 1 = 274177 x 67280421310721 */ 233 const unsigned char fermat_F6[] = 234 { 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x01 }; 235 nn n; 236 n.magic = WORD(0); 237 238 FORCE_USED_VAR(argc); 239 FORCE_USED_VAR(argv); 240 241 ret = nn_init(&n, 0); EG(ret, err); 242 /* Execute factorization on F5 */ 243 ret = nn_init_from_buf(&n, fermat_F5, sizeof(fermat_F5)); EG(ret, err); 244 ret = find_divisors(&n); EG(ret, err); 245 /* Execute factorization on F6 */ 246 ret = nn_init_from_buf(&n, fermat_F6, sizeof(fermat_F6)); EG(ret, err); 247 ret = find_divisors(&n); EG(ret, err); 248 /* Execute factorization on a random 80 bits number */ 249 ret = nn_one(&n); EG(ret, err); 250 /* Compute 2**80 = 0x1 << 80 */ 251 ret = nn_lshift(&n, &n, 80); EG(ret, err); 252 ret = nn_get_random_mod(&n, &n); EG(ret, err); 253 ret = find_divisors(&n); EG(ret, err); 254 255 return 0; 256 err: 257 return -1; 258 } 259 #endif 260