ModExpNG ("Next Generation") math model.

1 files changed, 690 insertions, 0 deletions

diff --git a/modexpng_fpga_model.py b/modexpng_fpga_model.py
new file mode 100644
index 0000000..1152bdf
--- /dev/null
+++ b/modexpng_fpga_model.py
@@ -0,0 +1,690 @@
+#!/usr/bin/python3
+#
+#
+# ModExpNG core math model.
+#
+#
+# Copyright (c) 2019, NORDUnet A/S
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+# - Redistributions of source code must retain the above copyright notice,
+#   this list of conditions and the following disclaimer.
+#
+# - Redistributions in binary form must reproduce the above copyright
+#   notice, this list of conditions and the following disclaimer in the
+#   documentation and/or other materials provided with the distribution.
+#
+# - Neither the name of the NORDUnet nor the names of its contributors may
+#   be used to endorse or promote products derived from this software
+#   without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+# TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+# TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+#
+
+
+# -------
+# Imports
+#--------
+
+import sys
+import importlib
+
+
+# --------------
+# Model Settings
+# --------------
+
+# length of public key
+KEY_LENGTH = 1024
+
+# how many parallel multipliers to use
+NUM_MULTS  = 8
+
+
+# ---------------
+# Internal Values
+# ---------------
+
+# half of key length
+_KEY_LENGTH_HALF = KEY_LENGTH // 2
+
+# width of internal math pipeline
+_WORD_WIDTH = 16
+
+# folder with test vector scripts
+_VECTOR_PATH = "/vector"
+
+# name of test vector class
+_VECTOR_CLASS = "Vector"
+
+
+#
+# Multi-Precision Integer
+#
+class ModExpNG_Operand():
+
+    def __init__(self, number, length, words = None):
+
+        if words is None:
+
+            # length must be divisible by word width
+            if (length % _WORD_WIDTH) > 0:
+                raise Exception("Bad number length!")
+
+            self._init_from_number(number, length)
+
+        else:
+
+            # length must match words count
+            if len(words) != length:
+                raise Exception("Bad words count!")
+
+            self._init_from_words(words, length)
+
+
+    def _init_from_words(self, words, count):
+
+        for i in range(count):
+
+            # word must not exceed 17 bits
+            if words[i] >= (2 ** (_WORD_WIDTH + 1)):
+                raise Exception("Word is too large!")
+
+        self.words = words
+
+    def _init_from_number(self, number, length):
+
+        num_hexchars_per_word = _WORD_WIDTH // 4
+        num_hexchars_total = length // num_hexchars_per_word
+
+        value_hex = format(number, 'x')
+
+        # value must not be larger than specified, but it can be smaller, so
+        # we may need to prepend it with zeroes
+        if len(value_hex) > num_hexchars_total:
+            raise Exception("Number is too large!")
+        else:
+            while len(value_hex) < num_hexchars_total:
+                value_hex = "0" + value_hex
+
+        # create empty list
+        self.words = list()
+
+        # fill in words
+        while len(value_hex) > 0:
+            value_hex_part = value_hex[-num_hexchars_per_word:]
+            value_hex = value_hex[:-num_hexchars_per_word]
+            self.words.append(int(value_hex_part, 16))
+
+    def number(self):
+        ret = 0
+        shift = 0
+        for word in self.words:
+            ret += word << shift
+            shift += _WORD_WIDTH
+        return ret
+
+
+#
+# Test Vector
+#
+class ModExpNG_TestVector():
+
+    def __init__(self):
+
+        # format target filename
+        filename = "vector_" + str(KEY_LENGTH) + "_randomized"
+
+        # add ./vector to import search path
+        sys.path.insert(1, sys.path[0] + _VECTOR_PATH)
+
+        # import from filename
+        vector_module = importlib.import_module(filename)
+
+        # get vector class
+        vector_class = getattr(vector_module, _VECTOR_CLASS)
+
+        # instantiate vector class
+        vector_inst = vector_class()
+
+        # obtain parts of vector
+        self.m        = ModExpNG_Operand(vector_inst.m,         KEY_LENGTH)
+        self.n        = ModExpNG_Operand(vector_inst.n,         KEY_LENGTH)
+        self.d        = ModExpNG_Operand(vector_inst.d,         KEY_LENGTH)
+        self.p        = ModExpNG_Operand(vector_inst.p,        _KEY_LENGTH_HALF)
+        self.q        = ModExpNG_Operand(vector_inst.q,        _KEY_LENGTH_HALF)
+        self.dp       = ModExpNG_Operand(vector_inst.dp,       _KEY_LENGTH_HALF)
+        self.dq       = ModExpNG_Operand(vector_inst.dq,       _KEY_LENGTH_HALF)
+        self.qinv     = ModExpNG_Operand(vector_inst.qinv,     _KEY_LENGTH_HALF)
+        self.n_factor = ModExpNG_Operand(vector_inst.n_factor,  KEY_LENGTH)
+        self.p_factor = ModExpNG_Operand(vector_inst.p_factor, _KEY_LENGTH_HALF)
+        self.q_factor = ModExpNG_Operand(vector_inst.q_factor, _KEY_LENGTH_HALF)
+        self.n_coeff  = ModExpNG_Operand(vector_inst.n_coeff,   KEY_LENGTH      + _WORD_WIDTH)
+        self.p_coeff  = ModExpNG_Operand(vector_inst.p_coeff,  _KEY_LENGTH_HALF + _WORD_WIDTH)
+        self.q_coeff  = ModExpNG_Operand(vector_inst.q_coeff,  _KEY_LENGTH_HALF + _WORD_WIDTH)
+
+
+class ModExpNG_PartRecombinator():
+
+    def _bit_select(self, x, msb, lsb):
+        y = 0
+        for pos in range(lsb, msb+1):
+            y |= (x & (1 << pos)) >> lsb
+        return y
+
+    def _flush_pipeline(self):
+        self.z0, self.y0, self.x0 = 0, 0, 0
+
+    def _push_pipeline(self, part):
+
+        # split next part into 16-bit words
+        z = self._bit_select(part, 47, 32)
+        y = self._bit_select(part, 31, 16)
+        x = self._bit_select(part, 15,  0)
+
+        # shift to the right
+        z1 = z
+        y1 = y + self.z0
+        x1 = x + self.y0 + (self.x0 >> 16) # IMPORTANT: This carry can be up to two bits wide!!
+
+        # save lower 16 bits of the rightmost cell
+        t = self.x0 & 0xffff
+
+        # update internal latches
+        self.z0, self.y0, self.x0 = z1, y1, x1
+
+        # done
+        return t
+
+    def recombine_square(self, parts, ab_num_words):
+
+        # empty result so far
+        words = list()
+
+        # flush recombinator pipeline
+        self._flush_pipeline()
+
+        # the first tick produces null result, the last part produces
+        # two words, so we need (2*n - 1) + 2 = 2*n + 1 ticks total
+        # and should only save the result word during the last 2 * n ticks
+        for i in range(2 * ab_num_words + 1):
+
+            next_part = parts[i] if i < (2 * ab_num_words - 1) else 0
+            next_word = self._push_pipeline(next_part)
+
+            if i > 0:
+                words.append(next_word)
+
+        return words
+
+    def recombine_triangle(self, parts, ab_num_words):
+
+        # empty result so far
+        words = list()
+
+        # flush recombinator pipeline
+        self._flush_pipeline()
+
+        # the first tick produces null result, so we need n + 1 + 1 = n + 2
+        # ticks total and should only save the result word during the last n ticks
+        for i in range(ab_num_words + 2):
+
+            next_part = parts[i] if i < (ab_num_words + 1) else 0
+            next_word = self._push_pipeline(next_part)
+
+            if i > 0:
+                words.append(next_word)
+
+        return words
+
+    def recombine_rectangle(self, parts, ab_num_words):
+
+        # empty result so far
+        words = list()
+
+        # flush recombinator pipeline
+        self._flush_pipeline()
+
+        # the first tick produces null result, the last part produces
+        # two words, so we need 2 * n + 2 ticks total and should only save
+        # the result word during the last 2 * n + 1 ticks
+        for i in range(2 * ab_num_words + 2):
+
+            next_part = parts[i] if i < (2 * ab_num_words) else 0
+            next_word = self._push_pipeline(next_part)
+
+            if i > 0:
+                words.append(next_word)
+
+        return words
+
+
+class ModExpNG_WordMultiplier():
+
+    def __init__(self):
+
+        self._macs = list()
+        self._indices = list()
+
+        self._mac_aux = list()
+        self._index_aux = list()
+
+        for x in range(NUM_MULTS):
+            self._macs.append(0)
+            self._indices.append(0)
+
+        self._mac_aux.append(0)
+        self._index_aux.append(0)
+
+    def _clear_all_macs(self):
+        for x in range(NUM_MULTS):
+            self._macs[x] = 0
+
+    def _clear_one_mac(self, x):
+        self._macs[x] = 0
+
+    def _clear_mac_aux(self):
+        self._mac_aux[0] = 0
+
+    def _update_one_mac(self, x, value):
+        self._macs[x] += value
+
+    def _update_mac_aux(self, value):
+        self._mac_aux[0] += value
+
+    def _preset_indices(self, col):
+        for x in range(len(self._indices)):
+            self._indices[x] = col * len(self._indices) + x
+
+    def _preset_index_aux(self, num_cols):
+        self._index_aux[0] = num_cols * len(self._indices)
+
+    def _rotate_indices(self, num_words):
+        for x in range(len(self._indices)):
+            self._indices[x] -= 1
+            if self._indices[x] < 0:
+                self._indices[x] += num_words
+
+    def _rotate_index_aux(self):
+        self._index_aux[0] -= 1
+
+    def multiply_square(self, a_wide, b_narrow, ab_num_words):
+
+        num_cols = ab_num_words // NUM_MULTS
+
+        parts = list()
+        for i in range(2 * ab_num_words - 1):
+            parts.append(0)
+
+        for col in range(num_cols):
+
+            self._clear_all_macs()
+            self._preset_indices(col)
+
+            for t in range(ab_num_words):
+
+                # current b-word
+                bt = b_narrow.words[t]
+
+                # multiply by a-words
+                for x in range(NUM_MULTS):
+                    ax = a_wide.words[self._indices[x]]
+                    self._update_one_mac(x, ax * bt)
+
+                    if t == (col * NUM_MULTS + x):
+                        parts[t] = self._macs[x]
+                        self._clear_one_mac(x)
+
+                # save the uppers part of product at end of column,
+                # for the last column don't save the very last part
+                if t == (ab_num_words - 1):
+                    for x in range(NUM_MULTS):
+                        if not (col == (num_cols - 1) and x == (NUM_MULTS - 1)):
+                            parts[ab_num_words + col * NUM_MULTS + x] = self._macs[x]
+
+                self._rotate_indices(ab_num_words)
+
+        return parts
+
+    def multiply_triangle(self, a_wide, b_narrow, ab_num_words):
+
+        num_cols = ab_num_words // NUM_MULTS
+
+        parts = list()
+        for i in range(ab_num_words + 1):
+            parts.append(0)
+
+        for col in range(num_cols):
+
+            last_col = col == (num_cols - 1)
+
+            self._clear_all_macs()
+            self._preset_indices(col)
+
+            if last_col:
+                self._clear_mac_aux()
+                self._preset_index_aux(num_cols)
+
+            for t in range(ab_num_words + 1):
+
+                # current b-word
+                bt = b_narrow.words[t]
+
+                # multiply by a-words
+                for x in range(NUM_MULTS):
+                    ax = a_wide.words[self._indices[x]]
+                    self._update_one_mac(x, ax * bt)
+
+                    if t == (col * NUM_MULTS + x):
+                        parts[t] = self._macs[x]
+
+                # aux multiplier
+                if last_col:
+                    ax = a_wide.words[self._index_aux[0]]
+                    self._update_mac_aux(ax * bt)
+
+                    if t == ab_num_words:
+                        parts[t] = self._mac_aux[0]
+
+                # shortcut
+                if not last_col:
+                    if t == (NUM_MULTS * (col + 1) - 1): break
+
+                # advance indices
+                self._rotate_indices(ab_num_words)
+                if last_col:
+                    self._rotate_index_aux()
+
+        return parts
+
+    def multiply_rectangle(self, a_wide, b_narrow, ab_num_words):
+
+        num_cols = ab_num_words // NUM_MULTS
+
+        parts = list()
+        for i in range(2 * ab_num_words):
+            parts.append(0)
+
+        for col in range(num_cols):
+
+            self._clear_all_macs()
+            self._preset_indices(col)
+
+            for t in range(ab_num_words+1):
+
+                # current b-word
+                bt = b_narrow.words[t]
+
+                # multiply by a-words
+                for x in range(NUM_MULTS):
+                    ax = a_wide.words[self._indices[x]]
+                    self._update_one_mac(x, ax * bt)
+
+                    # don't save one value for the very last time instant per column
+                    if t < ab_num_words and t == (col * NUM_MULTS + x):
+                        parts[t] = self._macs[x]
+                        self._clear_one_mac(x)
+
+                # save the uppers part of product at end of column
+                if t == ab_num_words:
+                    for x in range(NUM_MULTS):
+                        parts[ab_num_words + col * NUM_MULTS + x] = self._macs[x]
+
+                self._rotate_indices(ab_num_words)
+
+        return parts
+
+
+class ModExpNG_LowlevelOperator():
+
+    def __init__(self):
+        self._word_mask = 0
+        for x in range(_WORD_WIDTH):
+            self._word_mask |= (1 << x)
+
+    def _check_word(self, a):
+        if a < 0 or a >= (2 ** _WORD_WIDTH):
+            raise Exception("Word out of range!")
+
+    def _check_carry_borrow(self, cb):
+        if cb < 0 or cb > 1:
+            raise Exception("Carry or borrow out of range!")
+
+    def add_words(self, a, b, c_in):
+
+        self._check_word(a)
+        self._check_word(b)
+        self._check_carry_borrow(c_in)
+
+        sum = a + b + c_in
+
+        sum_s = sum & self._word_mask
+        sum_c = (sum >> _WORD_WIDTH) & 1
+
+        return (sum_c, sum_s)
+
+    def sub_words(self, a, b, b_in):
+        self._check_word(a)
+        self._check_word(b)
+        self._check_carry_borrow(b_in)
+
+        dif = a - b - b_in
+
+        if dif < 0:
+            dif_b = 1
+            dif_d = dif + 2 ** _WORD_WIDTH
+        else:
+            dif_b = 0
+            dif_d = dif
+
+        return (dif_b, dif_d)
+
+
+class ModExpNG_Worker():
+
+    def __init__(self):
+        self.recombinator = ModExpNG_PartRecombinator()
+        self.multiplier   = ModExpNG_WordMultiplier()
+        self.lowlevel     = ModExpNG_LowlevelOperator()
+
+    def exponentiate(self, iz, bz, e, n, n_factor, n_coeff, num_words):
+
+        # working variables
+        t1, t2 = iz, bz
+
+        # length-1, length-2, length-3, ..., 1, 0 (left-to-right)
+        for bit in range(_WORD_WIDTH * num_words - 1, -1, -1):
+
+            if e.number() & (1 << bit):
+                p1 = self.multiply(t1, t2, n, n_coeff, num_words)
+                p2 = self.multiply(t2, t2, n, n_coeff, num_words)
+            else:
+                p1 = self.multiply(t1, t1, n, n_coeff, num_words)
+                p2 = self.multiply(t2, t1, n, n_coeff, num_words)
+
+            t1, t2 = p1, p2
+
+            if (bit % 8) == 0:
+                pct = float((_WORD_WIDTH * num_words - bit) / (_WORD_WIDTH * num_words)) * 100.0
+                print("\rpct: %5.1f%%" % pct, end='')
+
+        print("")
+
+        return t1
+
+    def subtract(self, a, b, n, ab_num_words):
+
+        c_in = 0
+        b_in = 0
+
+        ab = list()
+        ab_n = list()
+
+        for x in range(ab_num_words):
+
+            a_word = a.words[x]
+            b_word = b.words[x]
+
+            (b_out, d_out) = self.lowlevel.sub_words(a_word, b_word, b_in)
+            (c_out, s_out) = self.lowlevel.add_words(d_out, n.words[x], c_in)
+
+            ab.append(d_out)
+            ab_n.append(s_out)
+
+            (c_in, b_in) = (c_out, b_out)
+
+        d = ab if not b_out else ab_n
+
+        return ModExpNG_Operand(None, ab_num_words, d)
+
+    def add(self, a, b, ab_num_words):
+
+        c_in = 0
+
+        ab = list()
+
+        for x in range(2 * ab_num_words):
+
+            a_word = a.words[x] if x < ab_num_words else 0
+            b_word = b.words[x]
+
+            (c_out, s_out) = self.lowlevel.add_words(a_word, b_word, c_in)
+
+            ab.append(s_out)
+
+            c_in = c_out
+
+        return ModExpNG_Operand(None, 2*ab_num_words, ab)
+
+    def multiply(self, a, b, n, n_coeff, ab_num_words, reduce_only=False, multiply_only=False):
+
+        # 1.
+        if reduce_only:
+            ab = a
+        else:
+            ab_parts = self.multiplier.multiply_square(a, b, ab_num_words)
+            ab_words = self.recombinator.recombine_square(ab_parts, ab_num_words)
+            ab = ModExpNG_Operand(None, 2 * ab_num_words, ab_words)
+
+        if multiply_only:
+            return ModExpNG_Operand(None, 2*ab_num_words, ab_words)
+
+        # 2.
+        q_parts = self.multiplier.multiply_triangle(ab, n_coeff, ab_num_words)
+        q_words = self.recombinator.recombine_triangle(q_parts, ab_num_words)
+        q = ModExpNG_Operand(None, ab_num_words + 1, q_words)
+
+        # 3.
+        m_parts = self.multiplier.multiply_rectangle(n, q, ab_num_words)
+        m_words = self.recombinator.recombine_rectangle(m_parts, ab_num_words)
+        m = ModExpNG_Operand(None, 2 * ab_num_words + 1, m_words)
+
+        # 4.
+        r_xwords = list()
+        for i in range(2*ab_num_words):
+            r_xwords.append(ab.words[i] + m.words[i])
+
+        r_xwords.append(m.words[2 * ab_num_words])
+
+        cy = 0
+        for i in range(ab_num_words+1):
+            s = r_xwords[i] + cy
+            cy = s >> 16
+
+        R = list()
+        for i in range(ab_num_words):
+            R.append(0)
+
+        R[0] += cy # !!!
+
+        for i in range(ab_num_words):
+            R[i] += r_xwords[ab_num_words + i + 1]
+
+        return ModExpNG_Operand(None, ab_num_words, R)
+
+
+if __name__ == "__main__":
+
+    # load test vector
+    vector = ModExpNG_TestVector()
+
+    # create worker
+    worker = ModExpNG_Worker()
+
+    # number of words
+    pq_num_words = _KEY_LENGTH_HALF // _WORD_WIDTH
+
+    # obtain known good reference values with built-in math
+    s_known  = pow(vector.m.number(), vector.d.number(),  vector.n.number())
+    sp_known = pow(vector.m.number(), vector.dp.number(), vector.p.number())
+    sq_known = pow(vector.m.number(), vector.dq.number(), vector.q.number())
+
+    # first reduce message, this glues 2**-r to the message as a side effect
+    mpa = worker.multiply(vector.m, None, vector.p, vector.p_coeff, pq_num_words, reduce_only=True)
+    mqa = worker.multiply(vector.m, None, vector.q, vector.q_coeff, pq_num_words, reduce_only=True)
+
+    # unglue 2**-r from message by gluing 2**r to it to compensate
+    mp = worker.multiply(mpa, vector.p_factor, vector.p, vector.p_coeff, pq_num_words)
+    mq = worker.multiply(mqa, vector.q_factor, vector.q, vector.q_coeff, pq_num_words)
+
+    # one
+    i = ModExpNG_Operand(1, _KEY_LENGTH_HALF)
+
+    # bring one into Montgomery domain (glue 2**r to one)
+    ipz = worker.multiply(i, vector.p_factor, vector.p, vector.p_coeff, pq_num_words)
+    iqz = worker.multiply(i, vector.q_factor, vector.q, vector.q_coeff, pq_num_words)
+
+    # bring message into Montgomery domain (glue 2**r to message)
+    mpz = worker.multiply(mp, vector.p_factor, vector.p, vector.p_coeff, pq_num_words)
+    mqz = worker.multiply(mq, vector.q_factor, vector.q, vector.q_coeff, pq_num_words)
+
+    # do "easier" exponentiations
+    spz = worker.exponentiate(ipz, mpz, vector.dp, vector.p, vector.p_factor, vector.p_coeff, pq_num_words)
+    sqz = worker.exponentiate(iqz, mqz, vector.dq, vector.q, vector.q_factor, vector.q_coeff, pq_num_words)
+
+    # return "easier" parts from Montgomery domain (unglue 2**r from result)
+    sp = worker.multiply(i, spz, vector.p, vector.p_coeff, pq_num_words)
+    sq = worker.multiply(i, sqz, vector.q, vector.q_coeff, pq_num_words)
+
+    # check "easier" results
+    if sp.number() == sp_known: print("sp is OK")
+    else:              print("sp is WRONG!")
+
+    if sq.number() == sq_known: print("sq is OK")
+    else:              print("sq is WRONG!")
+
+
+    # do the "Garner's formula" part
+
+    # 1. r = sp - sq mod p
+    sr = worker.subtract(sp, sq, vector.p, pq_num_words)
+
+    # 2. sr_qinv = sr * qinv mod p
+    sr_qinv_a = worker.multiply(sr, vector.qinv, vector.p, vector.p_coeff, pq_num_words)
+    sr_qinv = worker.multiply(sr_qinv_a, vector.p_factor, vector.p, vector.p_coeff, pq_num_words)
+
+    # 3. q_sr_qinv = q * sr_qinv
+    q_sr_qinv = worker.multiply(vector.q, sr_qinv, None, None, pq_num_words, multiply_only=True)
+
+    # 4. s_crt = sq + q_sr_qinv
+    s_crt = worker.add(sq, q_sr_qinv, pq_num_words)
+
+    # check
+    if s_crt.number() != s_known:
+        print("ERROR: s_crt != s_known!")
+    else:
+        print("s is OK")
+