nim-stint/stint/private/uint_div.nim

283 lines
10 KiB
Nim
Raw Normal View History

# Stint
2023-06-08 10:53:38 +00:00
# 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.
2021-08-06 12:44:25 +00:00
import
# Status lib
stew/bitops2,
# Internal
./datatypes,
./uint_bitwise,
2022-01-12 17:25:55 +00:00
./primitives/[addcarry_subborrow, extended_precision]
2021-08-06 12:44:25 +00:00
# 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
2021-08-06 12:44:25 +00:00
# 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`
2023-06-14 10:41:02 +00:00
##
## Specialized for M being a multi-precision integer.
2021-08-06 12:44:25 +00:00
##
## With a word W = 2^WordBitWidth and a modulus M
## Does a <- a * W + c (mod M)
## and returns q = (a * W + c ) / M
2021-08-06 12:44:25 +00:00
##
## The modulus `M` most-significant bit at `mBits` MUST be set.
2023-06-14 10:41:02 +00:00
# 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`
2023-06-14 10:41:02 +00:00
##
## 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)
# (hi, lo) = a * 2^64 + c
2023-06-22 09:43:31 +00:00
if R == 0:
# We can delegate this R == 0 case to the
# shlAddMod_multi, with the same result.
# But isn't it faster to handle it here?
var q, r: Word
div2n1n(q, r, a[0], c, M[0])
a[0] = r
return q
else:
let clz = WordBitWidth-R
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):
2023-06-13 12:03:43 +00:00
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
2023-06-14 10:41:02 +00:00
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
2022-01-12 17:25:55 +00:00
# ######################################################################
# Division implementations
#
# Multi-precision division is a costly
# and also difficult to implement operation
2022-01-12 17:25:55 +00:00
# ##### 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
2023-06-14 10:41:02 +00:00
#
2022-01-12 17:25:55 +00:00
# - 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