Adam Vartanian | bfcf3a7 | 2018-08-10 14:55:24 +0100 | [diff] [blame] | 1 | /* Copyright (c) 2018, Google Inc. |
| 2 | * |
| 3 | * Permission to use, copy, modify, and/or distribute this software for any |
| 4 | * purpose with or without fee is hereby granted, provided that the above |
| 5 | * copyright notice and this permission notice appear in all copies. |
| 6 | * |
| 7 | * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
| 8 | * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
| 9 | * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY |
| 10 | * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
| 11 | * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION |
| 12 | * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN |
| 13 | * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ |
| 14 | |
| 15 | #include <openssl/bn.h> |
| 16 | |
| 17 | #include <assert.h> |
| 18 | |
| 19 | #include "internal.h" |
| 20 | |
| 21 | |
| 22 | // The following functions use a Barrett reduction variant to avoid leaking the |
| 23 | // numerator. See http://ridiculousfish.com/blog/posts/labor-of-division-episode-i.html |
| 24 | // |
| 25 | // We use 32-bit numerator and 16-bit divisor for simplicity. This allows |
| 26 | // computing |m| and |q| without architecture-specific code. |
| 27 | |
| 28 | // mod_u16 returns |n| mod |d|. |p| and |m| are the "magic numbers" for |d| (see |
| 29 | // reference). For proof of correctness in Coq, see |
| 30 | // https://github.com/davidben/fiat-crypto/blob/barrett/src/Arithmetic/BarrettReduction/RidiculousFish.v |
| 31 | // Note the Coq version of |mod_u16| additionally includes the computation of |
| 32 | // |p| and |m| from |bn_mod_u16_consttime| below. |
| 33 | static uint16_t mod_u16(uint32_t n, uint16_t d, uint32_t p, uint32_t m) { |
| 34 | // Compute floor(n/d) per steps 3 through 5. |
| 35 | uint32_t q = ((uint64_t)m * n) >> 32; |
| 36 | // Note there is a typo in the reference. We right-shift by one, not two. |
| 37 | uint32_t t = ((n - q) >> 1) + q; |
| 38 | t = t >> (p - 1); |
| 39 | |
| 40 | // Multiply and subtract to get the remainder. |
| 41 | n -= d * t; |
| 42 | assert(n < d); |
| 43 | return n; |
| 44 | } |
| 45 | |
| 46 | // shift_and_add_mod_u16 returns |r| * 2^32 + |a| mod |d|. |p| and |m| are the |
| 47 | // "magic numbers" for |d| (see reference). |
| 48 | static uint16_t shift_and_add_mod_u16(uint16_t r, uint32_t a, uint16_t d, |
| 49 | uint32_t p, uint32_t m) { |
| 50 | // Incorporate |a| in two 16-bit chunks. |
| 51 | uint32_t t = r; |
| 52 | t <<= 16; |
| 53 | t |= a >> 16; |
| 54 | t = mod_u16(t, d, p, m); |
| 55 | |
| 56 | t <<= 16; |
| 57 | t |= a & 0xffff; |
| 58 | t = mod_u16(t, d, p, m); |
| 59 | return t; |
| 60 | } |
| 61 | |
| 62 | uint16_t bn_mod_u16_consttime(const BIGNUM *bn, uint16_t d) { |
| 63 | if (d <= 1) { |
| 64 | return 0; |
| 65 | } |
| 66 | |
| 67 | // Compute the "magic numbers" for |d|. See steps 1 and 2. |
| 68 | // This computes p = ceil(log_2(d)). |
| 69 | uint32_t p = BN_num_bits_word(d - 1); |
| 70 | // This operation is not constant-time, but |p| and |d| are public values. |
| 71 | // Note that |p| is at most 16, so the computation fits in |uint64_t|. |
| 72 | assert(p <= 16); |
| 73 | uint32_t m = ((UINT64_C(1) << (32 + p)) + d - 1) / d; |
| 74 | |
| 75 | uint16_t ret = 0; |
| 76 | for (int i = bn->width - 1; i >= 0; i--) { |
| 77 | #if BN_BITS2 == 32 |
| 78 | ret = shift_and_add_mod_u16(ret, bn->d[i], d, p, m); |
| 79 | #elif BN_BITS2 == 64 |
| 80 | ret = shift_and_add_mod_u16(ret, bn->d[i] >> 32, d, p, m); |
| 81 | ret = shift_and_add_mod_u16(ret, bn->d[i] & 0xffffffff, d, p, m); |
| 82 | #else |
| 83 | #error "Unknown BN_ULONG size" |
| 84 | #endif |
| 85 | } |
| 86 | return ret; |
| 87 | } |