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 */
pollard_rho(nn_t d,nn_src_t n,const word_t c)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 */
find_divisors(nn_src_t in)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
main(int argc,char * argv[])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