• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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.
mod_u16(uint32_t n,uint16_t d,uint32_t p,uint32_t m)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).
shift_and_add_mod_u16(uint16_t r,uint32_t a,uint16_t d,uint32_t p,uint32_t m)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 
bn_mod_u16_consttime(const BIGNUM * bn,uint16_t d)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 }
88