274 lines
9.8 KiB
Nim
274 lines
9.8 KiB
Nim
# Stint
|
|
# Copyright 2018-2023 Status Research & Development GmbH
|
|
# Licensed under either of
|
|
#
|
|
# * Apache License, version 2.0, ([LICENSE-APACHE](LICENSE-APACHE) or http://www.apache.org/licenses/LICENSE-2.0)
|
|
# * MIT license ([LICENSE-MIT](LICENSE-MIT) or http://opensource.org/licenses/MIT)
|
|
#
|
|
# at your option. This file may not be copied, modified, or distributed except according to those terms.
|
|
|
|
import
|
|
# Status lib
|
|
stew/bitops2,
|
|
# Internal
|
|
./datatypes,
|
|
./uint_bitwise,
|
|
./primitives/[addcarry_subborrow, extended_precision]
|
|
|
|
# Division
|
|
# --------------------------------------------------------
|
|
|
|
func shortDiv*(a: var Limbs, k: Word): Word =
|
|
## Divide `a` by k in-place and return the remainder
|
|
result = Word(0)
|
|
|
|
let clz = leadingZeros(k)
|
|
let normK = k shl clz
|
|
|
|
for i in countdown(a.len-1, 0):
|
|
# dividend = 2^64 * remainder + a[i]
|
|
var hi = result
|
|
var lo = a[i]
|
|
if hi == 0:
|
|
if lo < k:
|
|
a[i] = 0
|
|
elif lo == k:
|
|
a[i] = 1
|
|
result = 0
|
|
continue
|
|
# Normalize, shifting the remainder by clz(k) cannot overflow.
|
|
hi = (hi shl clz) or (lo shr (WordBitWidth - clz))
|
|
lo = lo shl clz
|
|
div2n1n(a[i], result, hi, lo, normK)
|
|
# Undo normalization
|
|
result = result shr clz
|
|
|
|
func shlAddMod_multi(a: var openArray[Word], c: Word,
|
|
M: openArray[Word], mBits: int): Word =
|
|
## Fused modular left-shift + add
|
|
## Shift input `a` by a word and add `c` modulo `M`
|
|
##
|
|
## Specialized for M being a multi-precision integer.
|
|
##
|
|
## With a word W = 2^WordBitWidth and a modulus M
|
|
## Does a <- a * W + c (mod M)
|
|
## and returns q = (a * W + c ) / M
|
|
##
|
|
## The modulus `M` most-significant bit at `mBits` MUST be set.
|
|
|
|
# Assuming 64-bit words
|
|
let hi = a[^1] # Save the high word to detect carries
|
|
let R = mBits and (WordBitWidth - 1) # R = mBits mod 64
|
|
|
|
var a0, a1, m0: Word
|
|
if R == 0: # If the number of mBits is a multiple of 64
|
|
a0 = a[^1] #
|
|
copyWords(a, 1, a, 0, a.len-1) # we can just shift words
|
|
a[0] = c # and replace the first one by c
|
|
a1 = a[^1]
|
|
m0 = M[^1]
|
|
else: # Else: need to deal with partial word shifts at the edge.
|
|
let clz = WordBitWidth-R
|
|
a0 = (a[^1] shl clz) or (a[^2] shr R)
|
|
copyWords(a, 1, a, 0, a.len-1)
|
|
a[0] = c
|
|
a1 = (a[^1] shl clz) or (a[^2] shr R)
|
|
m0 = (M[^1] shl clz) or (M[^2] shr R)
|
|
|
|
# m0 has its high bit set. (a0, a1)/m0 fits in a limb.
|
|
# Get a quotient q, at most we will be 2 iterations off
|
|
# from the true quotient
|
|
var q: Word # Estimate quotient
|
|
if a0 == m0: # if a_hi == divisor
|
|
q = high(Word) # quotient = MaxWord (0b1111...1111)
|
|
elif a0 == 0 and a1 < m0: # elif q == 0, true quotient = 0
|
|
q = 0
|
|
return q
|
|
else:
|
|
var r: Word
|
|
div2n1n(q, r, a0, a1, m0) # else instead of being of by 0, 1 or 2
|
|
q -= 1 # we return q-1 to be off by -1, 0 or 1
|
|
|
|
# Now substract a*2^64 - q*m
|
|
var carry = Word(0)
|
|
var overM = true # Track if quotient greater than the modulus
|
|
|
|
for i in 0 ..< M.len:
|
|
var qm_lo: Word
|
|
block: # q*m
|
|
# q * p + carry (doubleword) carry from previous limb
|
|
muladd1(carry, qm_lo, q, M[i], carry)
|
|
|
|
block: # a*2^64 - q*m
|
|
var borrow: Borrow
|
|
subB(borrow, a[i], a[i], qm_lo, Borrow(0))
|
|
carry += Word(borrow) # Adjust if borrow
|
|
|
|
if a[i] != M[i]:
|
|
overM = a[i] > M[i]
|
|
|
|
# Fix quotient, the true quotient is either q-1, q or q+1
|
|
#
|
|
# if carry < q or carry == q and overM we must do "a -= M"
|
|
# if carry > hi (negative result) we must do "a += M"
|
|
if carry > hi:
|
|
var c = Carry(0)
|
|
for i in 0 ..< a.len:
|
|
addC(c, a[i], a[i], M[i], c)
|
|
q -= 1
|
|
elif overM or (carry < hi):
|
|
var b = Borrow(0)
|
|
for i in 0 ..< a.len:
|
|
subB(b, a[i], a[i], M[i], b)
|
|
q += 1
|
|
|
|
return q
|
|
|
|
func shlAddMod(a: var openArray[Word], c: Word,
|
|
M: openArray[Word], mBits: int): Word {.inline.}=
|
|
## Fused modular left-shift + add
|
|
## Shift input `a` by a word and add `c` modulo `M`
|
|
##
|
|
## With a word W = 2^WordBitWidth and a modulus M
|
|
## Does a <- a * W + c (mod M)
|
|
## and returns q = (a * W + c ) / M
|
|
##
|
|
## The modulus `M` most-significant bit at `mBits` MUST be set.
|
|
if mBits <= WordBitWidth:
|
|
# If M fits in a single limb
|
|
|
|
# We normalize M with clz so that the MSB is set
|
|
# And normalize (a * 2^64 + c) by R as well to maintain the result
|
|
# This ensures that (a0, a1)/p0 fits in a limb.
|
|
let R = mBits and (WordBitWidth - 1)
|
|
let clz = WordBitWidth-R
|
|
|
|
# (hi, lo) = a * 2^64 + c
|
|
let hi = (a[0] shl clz) or (c shr R)
|
|
let lo = c shl clz
|
|
let m0 = M[0] shl clz
|
|
|
|
var q, r: Word
|
|
div2n1n(q, r, hi, lo, m0)
|
|
a[0] = r shr clz
|
|
return q
|
|
else:
|
|
return shlAddMod_multi(a, c, M, mBits)
|
|
|
|
func divRem*(
|
|
q, r: var openArray[Word],
|
|
a, b: openArray[Word]
|
|
) =
|
|
let (aBits, aLen) = usedBitsAndWords(a)
|
|
let (bBits, bLen) = usedBitsAndWords(b)
|
|
let rLen = bLen
|
|
|
|
if unlikely(bBits == 0):
|
|
raise newException(DivByZeroDefect, "You attempted to divide by zero")
|
|
|
|
if aBits < bBits:
|
|
# if a uses less bits than b,
|
|
# a < b, so q = 0 and r = a
|
|
copyWords(r, 0, a, 0, aLen)
|
|
for i in aLen ..< r.len: # r.len >= rLen
|
|
r[i] = 0
|
|
for i in 0 ..< q.len:
|
|
q[i] = 0
|
|
else:
|
|
# The length of a is at least the divisor
|
|
# We can copy bLen-1 words
|
|
# and modular shift-lef-add the rest
|
|
let aOffset = aLen - bLen
|
|
copyWords(r, 0, a, aOffset+1, bLen-1)
|
|
r[rLen-1] = 0
|
|
# Now shift-left the copied words while adding the new word mod b
|
|
|
|
when nimvm:
|
|
# workaround nim bug #22095
|
|
var rr = @(r.toOpenArray(0, rLen-1))
|
|
var bb = @(b.toOpenArray(0, bLen-1))
|
|
for i in countdown(aOffset, 0):
|
|
q[i] = shlAddMod(
|
|
rr,
|
|
a[i],
|
|
bb,
|
|
bBits
|
|
)
|
|
for i in 0..rLen-1:
|
|
r[i] = rr[i]
|
|
else:
|
|
for i in countdown(aOffset, 0):
|
|
q[i] = shlAddMod(
|
|
r.toOpenArray(0, rLen-1),
|
|
a[i],
|
|
b.toOpenArray(0, bLen-1),
|
|
bBits
|
|
)
|
|
|
|
# Clean up extra words
|
|
for i in aOffset+1 ..< q.len:
|
|
q[i] = 0
|
|
for i in rLen ..< r.len:
|
|
r[i] = 0
|
|
|
|
# ######################################################################
|
|
# Division implementations
|
|
#
|
|
# Multi-precision division is a costly
|
|
# and also difficult to implement operation
|
|
|
|
# ##### Research #####
|
|
|
|
# Overview of division algorithms:
|
|
# - https://gmplib.org/manual/Division-Algorithms.html#Division-Algorithms
|
|
# - https://gmplib.org/~tege/division-paper.pdf
|
|
# - Comparison of fast division algorithms for large integers: http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Lec5-Fast-Division-Hasselstrom2003.pdf
|
|
|
|
# Schoolbook / Knuth Division (Algorithm D)
|
|
# - https://skanthak.homepage.t-online.de/division.html
|
|
# Review of implementation flaws
|
|
# - Hacker's Delight https://github.com/hcs0/Hackers-Delight/blob/master/divmnu64.c.txt
|
|
# - LLVM: https://github.com/llvm-mirror/llvm/blob/2c4ca68/lib/Support/APInt.cpp#L1289-L1451
|
|
# - ctbignum: https://github.com/niekbouman/ctbignum/blob/v0.5/include/ctbignum/division.hpp
|
|
# - Modern Computer Arithmetic - https://members.loria.fr/PZimmermann/mca/mca-cup-0.5.9.pdf
|
|
# p14 - 1.4.1 Naive Division
|
|
# - Handbook of Applied Cryptography - https://cacr.uwaterloo.ca/hac/about/chap14.pdf
|
|
# Chapter 14 algorithm 14.2.5
|
|
|
|
# Smith Method (and derivatives)
|
|
# This method improves Knuth algorithm by ~3x by removing regular normalization
|
|
# - A Multiple-Precision Division Algorithm, David M Smith
|
|
# American mathematical Society, 1996
|
|
# https://www.ams.org/journals/mcom/1996-65-213/S0025-5718-96-00688-6/S0025-5718-96-00688-6.pdf
|
|
#
|
|
# - An Efficient Multiple-Precision Division Algorithm,
|
|
# Liusheng Huang, Hong Zhong, Hong Shen, Yonglong Luo, 2005
|
|
# https://ieeexplore.ieee.org/document/1579076
|
|
#
|
|
# - Efficient multiple-precision integer division algorithm
|
|
# Debapriyay Mukhopadhyaya, Subhas C.Nandy, 2014
|
|
# https://www.sciencedirect.com/science/article/abs/pii/S0020019013002627
|
|
|
|
# Recursive division by Burnikel and Ziegler (http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz):
|
|
# - Python implementation: https://bugs.python.org/file11060/fast_div.py and discussion https://bugs.python.org/issue3451
|
|
# - C++ implementation: https://github.com/linbox-team/givaro/blob/master/src/kernel/recint/rudiv.h
|
|
# - The Handbook of Elliptic and Hyperelliptic Cryptography Algorithm 10.35 on page 188 has a more explicit version of the div2NxN algorithm. This algorithm is directly recursive and avoids the mutual recursion of the original paper's calls between div2NxN and div3Nx2N.
|
|
# - Modern Computer Arithmetic - https://members.loria.fr/PZimmermann/mca/mca-cup-0.5.9.pdf
|
|
# p18 - 1.4.3 Divide and Conquer Division
|
|
|
|
# Newton Raphson Iterations
|
|
# - Putty (constant-time): https://github.com/github/putty/blob/0.74/mpint.c#L1818-L2112
|
|
# - Modern Computer Arithmetic - https://members.loria.fr/PZimmermann/mca/mca-cup-0.5.9.pdf
|
|
# p18 - 1.4.3 Divide and Conquer Division
|
|
|
|
# Other libraries that can be used as reference for alternative (?) implementations:
|
|
# - TTMath: https://github.com/status-im/nim-ttmath/blob/8f6ff2e57b65a350479c4012a53699e262b19975/src/headers/ttmathuint.h#L1530-L2383
|
|
# - LibTomMath: https://github.com/libtom/libtommath
|
|
# - Google Abseil for uint128: https://github.com/abseil/abseil-cpp/tree/master/absl/numeric
|
|
# Note: GMP/MPFR are GPL. The papers can be used but not their code.
|
|
|
|
# Related research
|
|
# - Efficient divide-and-conquer multiprecision integer division
|
|
# William Hart, IEEE 2015
|
|
# https://github.com/wbhart/bsdnt
|
|
# https://ieeexplore.ieee.org/document/7203801 |