i62_modpow2.c
1 /* 2 * Copyright (c) 2017 Thomas Pornin <pornin@bolet.org> 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining 5 * a copy of this software and associated documentation files (the 6 * "Software"), to deal in the Software without restriction, including 7 * without limitation the rights to use, copy, modify, merge, publish, 8 * distribute, sublicense, and/or sell copies of the Software, and to 9 * permit persons to whom the Software is furnished to do so, subject to 10 * the following conditions: 11 * 12 * The above copyright notice and this permission notice shall be 13 * included in all copies or substantial portions of the Software. 14 * 15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 16 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 17 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 18 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS 19 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN 20 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 21 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 22 * SOFTWARE. 23 */ 24 25 #include "inner.h" 26 27 #if BR_INT128 || BR_UMUL128 28 29 #if BR_INT128 30 31 /* 32 * Compute x*y+v1+v2. Operands are 64-bit, and result is 128-bit, with 33 * high word in "hi" and low word in "lo". 34 */ 35 #define FMA1(hi, lo, x, y, v1, v2) do { \ 36 unsigned __int128 fmaz; \ 37 fmaz = (unsigned __int128)(x) * (unsigned __int128)(y) \ 38 + (unsigned __int128)(v1) + (unsigned __int128)(v2); \ 39 (hi) = (uint64_t)(fmaz >> 64); \ 40 (lo) = (uint64_t)fmaz; \ 41 } while (0) 42 43 /* 44 * Compute x1*y1+x2*y2+v1+v2. Operands are 64-bit, and result is 128-bit, 45 * with high word in "hi" and low word in "lo". 46 * 47 * Callers should ensure that the two inner products, and the v1 and v2 48 * operands, are multiple of 4 (this is not used by this specific definition 49 * but may help other implementations). 50 */ 51 #define FMA2(hi, lo, x1, y1, x2, y2, v1, v2) do { \ 52 unsigned __int128 fmaz; \ 53 fmaz = (unsigned __int128)(x1) * (unsigned __int128)(y1) \ 54 + (unsigned __int128)(x2) * (unsigned __int128)(y2) \ 55 + (unsigned __int128)(v1) + (unsigned __int128)(v2); \ 56 (hi) = (uint64_t)(fmaz >> 64); \ 57 (lo) = (uint64_t)fmaz; \ 58 } while (0) 59 60 #elif BR_UMUL128 61 62 #include <intrin.h> 63 64 #define FMA1(hi, lo, x, y, v1, v2) do { \ 65 uint64_t fmahi, fmalo; \ 66 unsigned char fmacc; \ 67 fmalo = _umul128((x), (y), &fmahi); \ 68 fmacc = _addcarry_u64(0, fmalo, (v1), &fmalo); \ 69 _addcarry_u64(fmacc, fmahi, 0, &fmahi); \ 70 fmacc = _addcarry_u64(0, fmalo, (v2), &(lo)); \ 71 _addcarry_u64(fmacc, fmahi, 0, &(hi)); \ 72 } while (0) 73 74 /* 75 * Normally we should use _addcarry_u64() for FMA2 too, but it makes 76 * Visual Studio crash. Instead we use this version, which leverages 77 * the fact that the vx operands, and the products, are multiple of 4. 78 * This is unfortunately slower. 79 */ 80 #define FMA2(hi, lo, x1, y1, x2, y2, v1, v2) do { \ 81 uint64_t fma1hi, fma1lo; \ 82 uint64_t fma2hi, fma2lo; \ 83 uint64_t fmatt; \ 84 fma1lo = _umul128((x1), (y1), &fma1hi); \ 85 fma2lo = _umul128((x2), (y2), &fma2hi); \ 86 fmatt = (fma1lo >> 2) + (fma2lo >> 2) \ 87 + ((v1) >> 2) + ((v2) >> 2); \ 88 (lo) = fmatt << 2; \ 89 (hi) = fma1hi + fma2hi + (fmatt >> 62); \ 90 } while (0) 91 92 /* 93 * The FMA2 macro definition we would prefer to use, but it triggers 94 * an internal compiler error in Visual Studio 2015. 95 * 96 #define FMA2(hi, lo, x1, y1, x2, y2, v1, v2) do { \ 97 uint64_t fma1hi, fma1lo; \ 98 uint64_t fma2hi, fma2lo; \ 99 unsigned char fmacc; \ 100 fma1lo = _umul128((x1), (y1), &fma1hi); \ 101 fma2lo = _umul128((x2), (y2), &fma2hi); \ 102 fmacc = _addcarry_u64(0, fma1lo, (v1), &fma1lo); \ 103 _addcarry_u64(fmacc, fma1hi, 0, &fma1hi); \ 104 fmacc = _addcarry_u64(0, fma2lo, (v2), &fma2lo); \ 105 _addcarry_u64(fmacc, fma2hi, 0, &fma2hi); \ 106 fmacc = _addcarry_u64(0, fma1lo, fma2lo, &(lo)); \ 107 _addcarry_u64(fmacc, fma1hi, fma2hi, &(hi)); \ 108 } while (0) 109 */ 110 111 #endif 112 113 #define MASK62 ((uint64_t)0x3FFFFFFFFFFFFFFF) 114 #define MUL62_lo(x, y) (((uint64_t)(x) * (uint64_t)(y)) & MASK62) 115 116 /* 117 * Subtract b from a, and return the final carry. If 'ctl32' is 0, then 118 * a[] is kept unmodified, but the final carry is still computed and 119 * returned. 120 */ 121 static uint32_t 122 i62_sub(uint64_t *a, const uint64_t *b, size_t num, uint32_t ctl32) 123 { 124 uint64_t cc, mask; 125 size_t u; 126 127 cc = 0; 128 ctl32 = -ctl32; 129 mask = (uint64_t)ctl32 | ((uint64_t)ctl32 << 32); 130 for (u = 0; u < num; u ++) { 131 uint64_t aw, bw, dw; 132 133 aw = a[u]; 134 bw = b[u]; 135 dw = aw - bw - cc; 136 cc = dw >> 63; 137 dw &= MASK62; 138 a[u] = aw ^ (mask & (dw ^ aw)); 139 } 140 return (uint32_t)cc; 141 } 142 143 /* 144 * Montgomery multiplication, over arrays of 62-bit values. The 145 * destination array (d) must be distinct from the other operands 146 * (x, y and m). All arrays are in little-endian format (least 147 * significant word comes first) over 'num' words. 148 */ 149 static void 150 montymul(uint64_t *d, const uint64_t *x, const uint64_t *y, 151 const uint64_t *m, size_t num, uint64_t m0i) 152 { 153 uint64_t dh; 154 size_t u, num4; 155 156 num4 = 1 + ((num - 1) & ~(size_t)3); 157 memset(d, 0, num * sizeof *d); 158 dh = 0; 159 for (u = 0; u < num; u ++) { 160 size_t v; 161 uint64_t f, xu; 162 uint64_t r, zh; 163 uint64_t hi, lo; 164 165 xu = x[u] << 2; 166 f = MUL62_lo(d[0] + MUL62_lo(x[u], y[0]), m0i) << 2; 167 168 FMA2(hi, lo, xu, y[0], f, m[0], d[0] << 2, 0); 169 r = hi; 170 171 for (v = 1; v < num4; v += 4) { 172 FMA2(hi, lo, xu, y[v + 0], 173 f, m[v + 0], d[v + 0] << 2, r << 2); 174 r = hi + (r >> 62); 175 d[v - 1] = lo >> 2; 176 FMA2(hi, lo, xu, y[v + 1], 177 f, m[v + 1], d[v + 1] << 2, r << 2); 178 r = hi + (r >> 62); 179 d[v + 0] = lo >> 2; 180 FMA2(hi, lo, xu, y[v + 2], 181 f, m[v + 2], d[v + 2] << 2, r << 2); 182 r = hi + (r >> 62); 183 d[v + 1] = lo >> 2; 184 FMA2(hi, lo, xu, y[v + 3], 185 f, m[v + 3], d[v + 3] << 2, r << 2); 186 r = hi + (r >> 62); 187 d[v + 2] = lo >> 2; 188 } 189 for (; v < num; v ++) { 190 FMA2(hi, lo, xu, y[v], f, m[v], d[v] << 2, r << 2); 191 r = hi + (r >> 62); 192 d[v - 1] = lo >> 2; 193 } 194 195 zh = dh + r; 196 d[num - 1] = zh & MASK62; 197 dh = zh >> 62; 198 } 199 i62_sub(d, m, num, (uint32_t)dh | NOT(i62_sub(d, m, num, 0))); 200 } 201 202 /* 203 * Conversion back from Montgomery representation. 204 */ 205 static void 206 frommonty(uint64_t *x, const uint64_t *m, size_t num, uint64_t m0i) 207 { 208 size_t u, v; 209 210 for (u = 0; u < num; u ++) { 211 uint64_t f, cc; 212 213 f = MUL62_lo(x[0], m0i) << 2; 214 cc = 0; 215 for (v = 0; v < num; v ++) { 216 uint64_t hi, lo; 217 218 FMA1(hi, lo, f, m[v], x[v] << 2, cc); 219 cc = hi << 2; 220 if (v != 0) { 221 x[v - 1] = lo >> 2; 222 } 223 } 224 x[num - 1] = cc >> 2; 225 } 226 i62_sub(x, m, num, NOT(i62_sub(x, m, num, 0))); 227 } 228 229 /* see inner.h */ 230 uint32_t 231 br_i62_modpow_opt(uint32_t *x31, const unsigned char *e, size_t elen, 232 const uint32_t *m31, uint32_t m0i31, uint64_t *tmp, size_t twlen) 233 { 234 size_t u, mw31num, mw62num; 235 uint64_t *x, *m, *t1, *t2; 236 uint64_t m0i; 237 uint32_t acc; 238 int win_len, acc_len; 239 240 /* 241 * Get modulus size, in words. 242 */ 243 mw31num = (m31[0] + 31) >> 5; 244 mw62num = (mw31num + 1) >> 1; 245 246 /* 247 * In order to apply this function, we must have enough room to 248 * copy the operand and modulus into the temporary array, along 249 * with at least two temporaries. If there is not enough room, 250 * switch to br_i31_modpow(). We also use br_i31_modpow() if the 251 * modulus length is not at least four words (94 bits or more). 252 */ 253 if (mw31num < 4 || (mw62num << 2) > twlen) { 254 /* 255 * We assume here that we can split an aligned uint64_t 256 * into two properly aligned uint32_t. Since both types 257 * are supposed to have an exact width with no padding, 258 * then this property must hold. 259 */ 260 size_t txlen; 261 262 txlen = mw31num + 1; 263 if (twlen < txlen) { 264 return 0; 265 } 266 br_i31_modpow(x31, e, elen, m31, m0i31, 267 (uint32_t *)tmp, (uint32_t *)tmp + txlen); 268 return 1; 269 } 270 271 /* 272 * Convert x to Montgomery representation: this means that 273 * we replace x with x*2^z mod m, where z is the smallest multiple 274 * of the word size such that 2^z >= m. We want to reuse the 31-bit 275 * functions here (for constant-time operation), but we need z 276 * for a 62-bit word size. 277 */ 278 for (u = 0; u < mw62num; u ++) { 279 br_i31_muladd_small(x31, 0, m31); 280 br_i31_muladd_small(x31, 0, m31); 281 } 282 283 /* 284 * Assemble operands into arrays of 62-bit words. Note that 285 * all the arrays of 62-bit words that we will handle here 286 * are without any leading size word. 287 * 288 * We also adjust tmp and twlen to account for the words used 289 * for these extra arrays. 290 */ 291 m = tmp; 292 x = tmp + mw62num; 293 tmp += (mw62num << 1); 294 twlen -= (mw62num << 1); 295 for (u = 0; u < mw31num; u += 2) { 296 size_t v; 297 298 v = u >> 1; 299 if ((u + 1) == mw31num) { 300 m[v] = (uint64_t)m31[u + 1]; 301 x[v] = (uint64_t)x31[u + 1]; 302 } else { 303 m[v] = (uint64_t)m31[u + 1] 304 + ((uint64_t)m31[u + 2] << 31); 305 x[v] = (uint64_t)x31[u + 1] 306 + ((uint64_t)x31[u + 2] << 31); 307 } 308 } 309 310 /* 311 * Compute window size. We support windows up to 5 bits; for a 312 * window of size k bits, we need 2^k+1 temporaries (for k = 1, 313 * we use special code that uses only 2 temporaries). 314 */ 315 for (win_len = 5; win_len > 1; win_len --) { 316 if ((((uint32_t)1 << win_len) + 1) * mw62num <= twlen) { 317 break; 318 } 319 } 320 321 t1 = tmp; 322 t2 = tmp + mw62num; 323 324 /* 325 * Compute m0i, which is equal to -(1/m0) mod 2^62. We were 326 * provided with m0i31, which already fulfills this property 327 * modulo 2^31; the single expression below is then sufficient. 328 */ 329 m0i = (uint64_t)m0i31; 330 m0i = MUL62_lo(m0i, (uint64_t)2 + MUL62_lo(m0i, m[0])); 331 332 /* 333 * Compute window contents. If the window has size one bit only, 334 * then t2 is set to x; otherwise, t2[0] is left untouched, and 335 * t2[k] is set to x^k (for k >= 1). 336 */ 337 if (win_len == 1) { 338 memcpy(t2, x, mw62num * sizeof *x); 339 } else { 340 uint64_t *base; 341 342 memcpy(t2 + mw62num, x, mw62num * sizeof *x); 343 base = t2 + mw62num; 344 for (u = 2; u < ((unsigned)1 << win_len); u ++) { 345 montymul(base + mw62num, base, x, m, mw62num, m0i); 346 base += mw62num; 347 } 348 } 349 350 /* 351 * Set x to 1, in Montgomery representation. We again use the 352 * 31-bit code. 353 */ 354 br_i31_zero(x31, m31[0]); 355 x31[(m31[0] + 31) >> 5] = 1; 356 br_i31_muladd_small(x31, 0, m31); 357 if (mw31num & 1) { 358 br_i31_muladd_small(x31, 0, m31); 359 } 360 for (u = 0; u < mw31num; u += 2) { 361 size_t v; 362 363 v = u >> 1; 364 if ((u + 1) == mw31num) { 365 x[v] = (uint64_t)x31[u + 1]; 366 } else { 367 x[v] = (uint64_t)x31[u + 1] 368 + ((uint64_t)x31[u + 2] << 31); 369 } 370 } 371 372 /* 373 * We process bits from most to least significant. At each 374 * loop iteration, we have acc_len bits in acc. 375 */ 376 acc = 0; 377 acc_len = 0; 378 while (acc_len > 0 || elen > 0) { 379 int i, k; 380 uint32_t bits; 381 uint64_t mask1, mask2; 382 383 /* 384 * Get the next bits. 385 */ 386 k = win_len; 387 if (acc_len < win_len) { 388 if (elen > 0) { 389 acc = (acc << 8) | *e ++; 390 elen --; 391 acc_len += 8; 392 } else { 393 k = acc_len; 394 } 395 } 396 bits = (acc >> (acc_len - k)) & (((uint32_t)1 << k) - 1); 397 acc_len -= k; 398 399 /* 400 * We could get exactly k bits. Compute k squarings. 401 */ 402 for (i = 0; i < k; i ++) { 403 montymul(t1, x, x, m, mw62num, m0i); 404 memcpy(x, t1, mw62num * sizeof *x); 405 } 406 407 /* 408 * Window lookup: we want to set t2 to the window 409 * lookup value, assuming the bits are non-zero. If 410 * the window length is 1 bit only, then t2 is 411 * already set; otherwise, we do a constant-time lookup. 412 */ 413 if (win_len > 1) { 414 uint64_t *base; 415 416 memset(t2, 0, mw62num * sizeof *t2); 417 base = t2 + mw62num; 418 for (u = 1; u < ((uint32_t)1 << k); u ++) { 419 uint64_t mask; 420 size_t v; 421 422 mask = -(uint64_t)EQ(u, bits); 423 for (v = 0; v < mw62num; v ++) { 424 t2[v] |= mask & base[v]; 425 } 426 base += mw62num; 427 } 428 } 429 430 /* 431 * Multiply with the looked-up value. We keep the product 432 * only if the exponent bits are not all-zero. 433 */ 434 montymul(t1, x, t2, m, mw62num, m0i); 435 mask1 = -(uint64_t)EQ(bits, 0); 436 mask2 = ~mask1; 437 for (u = 0; u < mw62num; u ++) { 438 x[u] = (mask1 & x[u]) | (mask2 & t1[u]); 439 } 440 } 441 442 /* 443 * Convert back from Montgomery representation. 444 */ 445 frommonty(x, m, mw62num, m0i); 446 447 /* 448 * Convert result into 31-bit words. 449 */ 450 for (u = 0; u < mw31num; u += 2) { 451 uint64_t zw; 452 453 zw = x[u >> 1]; 454 x31[u + 1] = (uint32_t)zw & 0x7FFFFFFF; 455 if ((u + 1) < mw31num) { 456 x31[u + 2] = (uint32_t)(zw >> 31); 457 } 458 } 459 return 1; 460 } 461 462 #else 463 464 /* see inner.h */ 465 uint32_t 466 br_i62_modpow_opt(uint32_t *x31, const unsigned char *e, size_t elen, 467 const uint32_t *m31, uint32_t m0i31, uint64_t *tmp, size_t twlen) 468 { 469 size_t mwlen; 470 471 mwlen = (m31[0] + 63) >> 5; 472 if (twlen < mwlen) { 473 return 0; 474 } 475 return br_i31_modpow_opt(x31, e, elen, m31, m0i31, 476 (uint32_t *)tmp, twlen << 1); 477 } 478 479 #endif 480 481 /* see inner.h */ 482 uint32_t 483 br_i62_modpow_opt_as_i31(uint32_t *x31, const unsigned char *e, size_t elen, 484 const uint32_t *m31, uint32_t m0i31, uint32_t *tmp, size_t twlen) 485 { 486 /* 487 * As documented, this function expects the 'tmp' argument to be 488 * 64-bit aligned. This is OK since this function is internal (it 489 * is not part of BearSSL's public API). 490 */ 491 return br_i62_modpow_opt(x31, e, elen, m31, m0i31, 492 (uint64_t *)tmp, twlen >> 1); 493 }