aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--modexpng_fpga_model.py690
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")
+