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 #include <libecc/fp/fp_sqrt.h>
17 #include <libecc/nn/nn_add.h>
18 #include <libecc/nn/nn_logical.h>
19
20 /*
21 * Compute the legendre symbol of an element of Fp:
22 *
23 * Legendre(a) = a^((p-1)/2) (p) = { -1, 0, 1 }
24 *
25 */
legendre(fp_src_t a)26 ATTRIBUTE_WARN_UNUSED_RET static int legendre(fp_src_t a)
27 {
28 int ret, iszero, cmp;
29 fp pow; /* The result if the exponentiation is in Fp */
30 fp one; /* The element 1 in the field */
31 nn exp; /* The power exponent is in NN */
32 pow.magic = one.magic = WORD(0);
33 exp.magic = WORD(0);
34
35 /* Initialize elements */
36 ret = fp_check_initialized(a); EG(ret, err);
37 ret = fp_init(&pow, a->ctx); EG(ret, err);
38 ret = fp_init(&one, a->ctx); EG(ret, err);
39 ret = nn_init(&exp, 0); EG(ret, err);
40
41 /* Initialize our variables from the Fp context of the
42 * input a.
43 */
44 ret = fp_init(&pow, a->ctx); EG(ret, err);
45 ret = fp_init(&one, a->ctx); EG(ret, err);
46 ret = nn_init(&exp, 0); EG(ret, err);
47
48 /* one = 1 in Fp */
49 ret = fp_one(&one); EG(ret, err);
50
51 /* Compute the exponent (p-1)/2
52 * The computation is done in NN, and the division by 2
53 * is performed using a right shift by one
54 */
55 ret = nn_dec(&exp, &(a->ctx->p)); EG(ret, err);
56 ret = nn_rshift(&exp, &exp, 1); EG(ret, err);
57
58 /* Compute a^((p-1)/2) in Fp using our fp_pow
59 * API.
60 */
61 ret = fp_pow(&pow, a, &exp); EG(ret, err);
62
63 ret = fp_iszero(&pow, &iszero); EG(ret, err);
64 ret = fp_cmp(&pow, &one, &cmp); EG(ret, err);
65 if (iszero) {
66 ret = 0;
67 } else if (cmp == 0) {
68 ret = 1;
69 } else {
70 ret = -1;
71 }
72
73 err:
74 /* Cleaning */
75 fp_uninit(&pow);
76 fp_uninit(&one);
77 nn_uninit(&exp);
78
79 return ret;
80 }
81
82 /*
83 * We implement the Tonelli-Shanks algorithm for finding
84 * square roots (quadratic residues) modulo a prime number,
85 * i.e. solving the equation:
86 * x^2 = n (p)
87 * where p is a prime number. This can be seen as an equation
88 * over the finite field Fp where a and x are elements of
89 * this finite field.
90 * Source: https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
91 * All ≡ are taken to mean (mod p) unless stated otherwise.
92 * Input : p an odd prime, and an integer n .
93 * Step 0. Check that n is indeed a square : (n | p) must be ≡ 1
94 * Step 1. [Factors out powers of 2 from p-1] Define q -odd- and s such as p-1 = q * 2^s
95 * - if s = 1 , i.e p ≡ 3 (mod 4) , output the two solutions r ≡ +/- n^((p+1)/4) .
96 * Step 2. Select a non-square z such as (z | p) = -1 , and set c ≡ z^q .
97 * Step 3. Set r ≡ n ^((q+1)/2) , t ≡ n^q, m = s .
98 * Step 4. Loop.
99 * - if t ≡ 1 output r, p-r .
100 * - Otherwise find, by repeated squaring, the lowest i , 0 < i < m , such as t^(2^i) ≡ 1
101 * - Let b ≡ c^(2^(m-i-1)), and set r ≡ r*b, t ≡ t*b^2 , c ≡ b^2 and m = i.
102 *
103 * NOTE: the algorithm is NOT constant time.
104 *
105 * The outputs, sqrt1 and sqrt2 ARE initialized by the function.
106 * The function returns 0 on success, -1 on error (in which case values of sqrt1 and sqrt2
107 * must not be considered).
108 *
109 * Aliasing is supported.
110 *
111 */
fp_sqrt(fp_t sqrt1,fp_t sqrt2,fp_src_t n)112 int fp_sqrt(fp_t sqrt1, fp_t sqrt2, fp_src_t n)
113 {
114 int ret, iszero, cmp, isodd;
115 nn q, s, one_nn, two_nn, m, i, tmp_nn;
116 fp z, t, b, r, c, one_fp, tmp_fp, __n;
117 fp_t _n = &__n;
118 q.magic = s.magic = one_nn.magic = two_nn.magic = m.magic = WORD(0);
119 i.magic = tmp_nn.magic = z.magic = t.magic = b.magic = WORD(0);
120 r.magic = c.magic = one_fp.magic = tmp_fp.magic = __n.magic = WORD(0);
121
122 ret = nn_init(&q, 0); EG(ret, err);
123 ret = nn_init(&s, 0); EG(ret, err);
124 ret = nn_init(&tmp_nn, 0); EG(ret, err);
125 ret = nn_init(&one_nn, 0); EG(ret, err);
126 ret = nn_init(&two_nn, 0); EG(ret, err);
127 ret = nn_init(&m, 0); EG(ret, err);
128 ret = nn_init(&i, 0); EG(ret, err);
129 ret = fp_init(&z, n->ctx); EG(ret, err);
130 ret = fp_init(&t, n->ctx); EG(ret, err);
131 ret = fp_init(&b, n->ctx); EG(ret, err);
132 ret = fp_init(&r, n->ctx); EG(ret, err);
133 ret = fp_init(&c, n->ctx); EG(ret, err);
134 ret = fp_init(&one_fp, n->ctx); EG(ret, err);
135 ret = fp_init(&tmp_fp, n->ctx); EG(ret, err);
136
137 /* Handle input aliasing */
138 ret = fp_copy(_n, n); EG(ret, err);
139
140 /* Initialize outputs */
141 ret = fp_init(sqrt1, _n->ctx); EG(ret, err);
142 ret = fp_init(sqrt2, _n->ctx); EG(ret, err);
143
144 /* one_nn = 1 in NN */
145 ret = nn_one(&one_nn); EG(ret, err);
146 /* two_nn = 2 in NN */
147 ret = nn_set_word_value(&two_nn, WORD(2)); EG(ret, err);
148
149 /* If our p prime of Fp is 2, then return the input as square roots */
150 ret = nn_cmp(&(_n->ctx->p), &two_nn, &cmp); EG(ret, err);
151 if (cmp == 0) {
152 ret = fp_copy(sqrt1, _n); EG(ret, err);
153 ret = fp_copy(sqrt2, _n); EG(ret, err);
154 ret = 0;
155 goto err;
156 }
157
158 /* Square root of 0 is 0 */
159 ret = fp_iszero(_n, &iszero); EG(ret, err);
160 if (iszero) {
161 ret = fp_zero(sqrt1); EG(ret, err);
162 ret = fp_zero(sqrt2); EG(ret, err);
163 ret = 0;
164 goto err;
165 }
166 /* Step 0. Check that n is indeed a square : (n | p) must be ≡ 1 */
167 if (legendre(_n) != 1) {
168 /* a is not a square */
169 ret = -1;
170 goto err;
171 }
172 /* Step 1. [Factors out powers of 2 from p-1] Define q -odd- and s such as p-1 = q * 2^s */
173 /* s = 0 */
174 ret = nn_zero(&s); EG(ret, err);
175 /* q = p - 1 */
176 ret = nn_copy(&q, &(_n->ctx->p)); EG(ret, err);
177 ret = nn_dec(&q, &q); EG(ret, err);
178 while (1) {
179 /* i is used as a temporary unused variable here */
180 ret = nn_divrem(&tmp_nn, &i, &q, &two_nn); EG(ret, err);
181 ret = nn_inc(&s, &s); EG(ret, err);
182 ret = nn_copy(&q, &tmp_nn); EG(ret, err);
183 /* If r is odd, we have finished our division */
184 ret = nn_isodd(&q, &isodd); EG(ret, err);
185 if (isodd) {
186 break;
187 }
188 }
189 /* - if s = 1 , i.e p ≡ 3 (mod 4) , output the two solutions r ≡ +/- n^((p+1)/4) . */
190 ret = nn_cmp(&s, &one_nn, &cmp); EG(ret, err);
191 if (cmp == 0) {
192 ret = nn_inc(&tmp_nn, &(_n->ctx->p)); EG(ret, err);
193 ret = nn_rshift(&tmp_nn, &tmp_nn, 2); EG(ret, err);
194 ret = fp_pow(sqrt1, _n, &tmp_nn); EG(ret, err);
195 ret = fp_neg(sqrt2, sqrt1); EG(ret, err);
196 ret = 0;
197 goto err;
198 }
199 /* Step 2. Select a non-square z such as (z | p) = -1 , and set c ≡ z^q . */
200 ret = fp_zero(&z); EG(ret, err);
201 while (legendre(&z) != -1) {
202 ret = fp_inc(&z, &z); EG(ret, err);
203 }
204 ret = fp_pow(&c, &z, &q); EG(ret, err);
205 /* Step 3. Set r ≡ n ^((q+1)/2) , t ≡ n^q, m = s . */
206 ret = nn_inc(&tmp_nn, &q); EG(ret, err);
207 ret = nn_rshift(&tmp_nn, &tmp_nn, 1); EG(ret, err);
208 ret = fp_pow(&r, _n, &tmp_nn); EG(ret, err);
209 ret = fp_pow(&t, _n, &q); EG(ret, err);
210 ret = nn_copy(&m, &s); EG(ret, err);
211 ret = fp_one(&one_fp); EG(ret, err);
212
213 /* Step 4. Loop. */
214 while (1) {
215 /* - if t ≡ 1 output r, p-r . */
216 ret = fp_cmp(&t, &one_fp, &cmp); EG(ret, err);
217 if (cmp == 0) {
218 ret = fp_copy(sqrt1, &r); EG(ret, err);
219 ret = fp_neg(sqrt2, sqrt1); EG(ret, err);
220 ret = 0;
221 goto err;
222 }
223 /* - Otherwise find, by repeated squaring, the lowest i , 0 < i < m , such as t^(2^i) ≡ 1 */
224 ret = nn_one(&i); EG(ret, err);
225 ret = fp_copy(&tmp_fp, &t); EG(ret, err);
226 while (1) {
227 ret = fp_sqr(&tmp_fp, &tmp_fp); EG(ret, err);
228 ret = fp_cmp(&tmp_fp, &one_fp, &cmp); EG(ret, err);
229 if (cmp == 0) {
230 break;
231 }
232 ret = nn_inc(&i, &i); EG(ret, err);
233 ret = nn_cmp(&i, &m, &cmp); EG(ret, err);
234 if (cmp == 0) {
235 /* i has reached m, that should not happen ... */
236 ret = -2;
237 goto err;
238 }
239 }
240 /* - Let b ≡ c^(2^(m-i-1)), and set r ≡ r*b, t ≡ t*b^2 , c ≡ b^2 and m = i. */
241 ret = nn_sub(&tmp_nn, &m, &i); EG(ret, err);
242 ret = nn_dec(&tmp_nn, &tmp_nn); EG(ret, err);
243 ret = fp_copy(&b, &c); EG(ret, err);
244 ret = nn_iszero(&tmp_nn, &iszero); EG(ret, err);
245 while (!iszero) {
246 ret = fp_sqr(&b, &b); EG(ret, err);
247 ret = nn_dec(&tmp_nn, &tmp_nn); EG(ret, err);
248 ret = nn_iszero(&tmp_nn, &iszero); EG(ret, err);
249 }
250 /* r ≡ r*b */
251 ret = fp_mul(&r, &r, &b); EG(ret, err);
252 /* c ≡ b^2 */
253 ret = fp_sqr(&c, &b); EG(ret, err);
254 /* t ≡ t*b^2 */
255 ret = fp_mul(&t, &t, &c); EG(ret, err);
256 /* m = i */
257 ret = nn_copy(&m, &i); EG(ret, err);
258 }
259
260 err:
261 /* Uninitialize local variables */
262 nn_uninit(&q);
263 nn_uninit(&s);
264 nn_uninit(&tmp_nn);
265 nn_uninit(&one_nn);
266 nn_uninit(&two_nn);
267 nn_uninit(&m);
268 nn_uninit(&i);
269 fp_uninit(&z);
270 fp_uninit(&t);
271 fp_uninit(&b);
272 fp_uninit(&r);
273 fp_uninit(&c);
274 fp_uninit(&one_fp);
275 fp_uninit(&tmp_fp);
276 fp_uninit(&__n);
277
278 PTR_NULLIFY(_n);
279
280 return ret;
281 }
282