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