Stash constant-time decimal conversion for constantine

README.md

Algorithms

Continued fractions

Continued fractions are used to convert

size_in_bits <=> size_in_decimal

for constant-time buffer preallocation when converting to decimal.

• https://en.wikipedia.org/wiki/Continued_fraction

Division by constants

References for division by constants as used to convert to decimal.

• Barret Reduction, Barrett, 1986 Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor https://en.wikipedia.org/wiki/Barrett_reduction https://link.springer.com/content/pdf/10.1007%2F3-540-47721-7_24.pdf

• Division by Invariant Integers using Multiplication Granlund, Montgomery, 1994, https://gmplib.org/~tege/divcnst-pldi94.pdf

• Reciprocal Multiplication, a tutorial Douglas W. Jones, 1999 https://homepage.divms.uiowa.edu/~jones/bcd/divide.html

• Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries Daniel Lemire, Owen Kaser, Nathan Kurz, 2019 https://arxiv.org/abs/1902.01961

Modular inverses mod 2ⁿ

We use "Dumas iterations" Explanation p11 "Dumas iterations" based on Newton-Raphson:

• Cetin Kaya Koc (2017), https://eprint.iacr.org/2017/411
• Jean-Guillaume Dumas (2012), https://arxiv.org/pdf/1209.6626v2.pdf
• Colin Plumb (1994), http://groups.google.com/groups?selm=1994Apr6.093116.27805%40mnemosyne.cs.du.edu Other sources:
• https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
• https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
• http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html

stash_constant_time_decimal_conv.nim

# Conversion to decimal
# ----------------------------------------------------------------
#
# The first problem to solve is precomputing the final size
# We also use continued fractions to approximate log₁₀(2)
#
# decimal_length = log₁₀(2) * binary_length
#
# sage: [(frac, numerical_approx(frac - log(2,10))) for frac in continued_fraction(log(2,10)).convergents()[0:10]]
# [(0, -0.301029995663981),
#  (1/3, 0.0323033376693522),
#  (3/10, -0.00102999566398115),
#  (28/93, 0.0000452731532231687),
#  (59/196, -9.58750071583525e-6),
#  (146/485, 9.32171070389121e-7),
#  (643/2136, -3.31171646772432e-8),
#  (4004/13301, 2.08054934391910e-9),
#  (8651/28738, -5.35579747218407e-10),
#  (12655/42039, 2.92154800352051e-10)]

const log10_2_Num = 12655
const log10_2_Denom = 42039

# Then the naive way to serialize is to repeatedly do
# const intToCharMap = "0123456789"
# rest = number
# while rest != 0:
#   digitToPrint = rest mod 10
#   result.add intToCharMap[digitToPrint]
#   rest /= 10
#
# For constant-time we:
# 1. can't compare with 0 as a stopping condition
# 2. repeatedly add to a buffer
# 3. can't use naive indexing (cache timing attacks though very unlikely given the small size)
# 4. need (fast) constant-time division
#
# 1 and 2 is solved by precomputing the length and make the number of add be fixed.
# 3 is easily solved by doing "digitToPrint + ord('0')" instead
#
# For 4 we note that a/10 = (a * 2ⁿ/10) / 2ⁿ
# So we can precompute 2ⁿ/10 at compile-time
# and do multiplication + shift at runtime.
#
# See references in the README, for proof that the approximation
# leads to correct result despite rounding to integer.

# No exceptions for the in-place API
{.push raises: [].}

func invmod_bitlength(r: var BigInt, a: BigInt) =
  ## This algorithm is needed to find 1/10 mod 2ⁿ
  ## It requires a to be odd so we will actually use it with 1/5.
  ##
  ## The algorithm does O(log(log a)) iterations and is constant-time
  ## Note that within the loop it does use multiplication which is O(n²)
  ##
  # Modular inverse algorithm:
  # Explanation p11 "Dumas iterations" based on Newton-Raphson:
  # - Cetin Kaya Koc (2017), https://eprint.iacr.org/2017/411
  # - Jean-Guillaume Dumas (2012), https://arxiv.org/pdf/1209.6626v2.pdf
  # - Colin Plumb (1994), http://groups.google.com/groups?selm=1994Apr6.093116.27805%40mnemosyne.cs.du.edu
  # Other sources:
  # - https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
  # - https://mumble.net/~campbell/2015/01/21/inverse-mod-power-of-two
  # - http://marc-b-reynolds.github.io/math/2017/09/18/ModInverse.html
  assert bool(BaseType(a.limbs[0]) and 1), "Internal Error: the input must be odd."

  let k = nextPowerOf2(BigInt.bits.uint64)
  let a = a # Prevent aliasing

  r.setOne()        # Start from an inverse of a modulo 2,
                    # a is odd so -1 is its inverse
  for _ in 0 ..< k: # at each iteration we get the inverse mod(2^2k)
    var t {.noInit.}, t2 {.noInit.}: BigInt
    t.prod(a, r)
    t2.setUint(2'u32)
    t2 -= t
    r.prod(r, t2)   # x' = x(2 - ax)

  # We overshoot, we have nextPowerOf2(k) instead of k
  # Hence for 2^381 bit number we have computed inverse mod 2^512
  # Since we inverse on the bitlength,
  # we only need to fix the last word, bits 381 to 383
  const WordMask = WordBitWidth - 1
  const DigitMask = block:
    const shift = BigInt.bits and WordMask
    if shift == 0:
      high(BaseType)
    else:
      (BaseType(1) shl shift) - 1
  # cast because "Error: -1 can't be converted to uint64"
  # due to the VM representing everything at compile-time as a signed int :/
  r.limbs[^1] = r.limbs[^1] and cast[SecretWord](Digitmask)

func barretQuotRem[nBits, dBits: static int](
       q: var BigInt[nBits], r: var BigInt[dBits],
       num: BigInt[nBits], d: BigInt[dBits], k: int) =
  ## Compute:
  ## q = num / (d * 2^k)
  ## r = num mod (d * 2^k)
  ##
  ## d must be odd,
  ## dividing by 10 should be d = 5 and k = 1
  ##
  ## dbits must be able to hold the real `d` bitlength
  ## plus the shift k
  ## For example for 10, dBits should be 4 even if you pass 5 (3 bits) and k = 1
  assert bool(BaseType(d.limbs[0]) and 1), "Internal Error: the input must be odd."

  # We shift by (2ⁿ)^w with n = 32 or 64 (word size)
  # and w chosen to have the next multiple of the word size
  const D = wordsRequired(dBits)
  const N = (D + wordsRequired(nBits)) * WordBitWidth
  debugecho "shift+dBits: ", N

  var m, shifted: BigInt[N] # zero-init required
  # shifted = num << 2^nw
  for i in 0 ..< num.limbs.len:
    shifted.limbs[i+D] = num.limbs[i]

  debugecho "num[", num.bits, "]:     ", num.tohex()
  debugecho "shifted[", shifted.bits, "]: ", shifted.tohex()

  # m = n << 2^nw / 10
  for i in 0 ..< d.limbs.len:
    m.limbs[i] = d.limbs[i]

  debugecho "m:        ", m.toHex()
  m.invmod_bitlength(m)
  debugecho "m:        ", m.toHex()
  let m2 = m
  m.prod(m2, shifted)
  debugecho "m:        ", m.toHex()
  m.shiftRight(k)

  debugecho "m:        ", m.toHex()
  debugecho "expected: ", "0x299b4fdd28cca42a11c5d9239edf7af23a5878d4b8d4eacbd84e1dce578189d3644666644eeccccc5ccb33333332aaab333333333333333"

  # q = (num * 2^nw/10) / 2^nw
  q.prod_high_words(num, m, D)

  # r = num - q*d
  var qd {.noInit.}: BigInt[nBits]
  qd.prod(q, d)

  var t {.noInit.}: BigInt[nBits]
  discard t.diff(num, qd)

  for i in 0 ..< r.limbs.len:
    r.limbs[i] = t.limbs[i]

var q: BigInt[381]
var r: BigInt[4]

let num = BigInt[381].fromHex("0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab")
let five = BigInt[4].fromUint(5'u32)

barretQuotRem(q, r, num, five, 1)

echo "q: ", q.toHex()
echo "r: ", r.toHex()