diff --git a/stint.nim b/stint.nim index 2631d54..b6c02e7 100644 --- a/stint.nim +++ b/stint.nim @@ -10,8 +10,8 @@ # import stint/[bitops2, endians2, intops, io, modular_arithmetic, literals_stint] # export bitops2, endians2, intops, io, modular_arithmetic, literals_stint -import stint/[io, uintops, bitops2] -export io, uintops, bitops2 +import stint/[io, uintops] +export io, uintops type # Int128* = Stint[128] diff --git a/stint/private/uint_div.nim b/stint/private/uint_div.nim index fb62b73..e3bfaed 100644 --- a/stint/private/uint_div.nim +++ b/stint/private/uint_div.nim @@ -13,7 +13,8 @@ import # Internal ./datatypes, ./uint_bitwise, - ./uint_shift + ./uint_shift, + ./primitives/[addcarry_subborrow, extended_precision] # Division # -------------------------------------------------------- @@ -36,31 +37,31 @@ func shortDiv*(a: var Limbs, k: Word): Word = # Undo normalization result = result shr clz -func binaryShiftDiv[qLen, rLen, uLen, vLen: static int]( - q: var Limbs[qLen], - r: var Limbs[rLen], - u: Limbs[uLen], - v: Limbs[vLen]) = - ## Division for multi-precision unsigned uint - ## Implementation through binary shift division - doAssert y.isZero.not() # This should be checked on release mode in the divmod caller proc +# func binaryShiftDiv[qLen, rLen, uLen, vLen: static int]( +# q: var Limbs[qLen], +# r: var Limbs[rLen], +# u: Limbs[uLen], +# v: Limbs[vLen]) = +# ## Division for multi-precision unsigned uint +# ## Implementation through binary shift division +# doAssert y.isZero.not() # This should be checked on release mode in the divmod caller proc - type SubTy = type x.lo +# type SubTy = type x.lo - var - shift = y.leadingZeros - x.leadingZeros - d = y shl shift +# var +# shift = y.leadingZeros - x.leadingZeros +# d = y shl shift - r = x +# r = x - while shift >= 0: - q += q - if r >= d: - r -= d - q.lo = q.lo or one(SubTy) +# while shift >= 0: +# q += q +# if r >= d: +# r -= d +# q.lo = q.lo or one(SubTy) - d = d shr 1 - dec(shift) +# d = d shr 1 +# dec(shift) func knuthDivLE[qLen, rLen, uLen, vLen: static int]( q: var Limbs[qLen], @@ -75,11 +76,9 @@ func knuthDivLE[qLen, rLen, uLen, vLen: static int]( ## - r must be of size vLen (assuming v uses all words) ## - uLen >= vLen ## - ## Knuth Division - ## - Knuth's "Algorithm D", The Art of Computer Programming, 1998 - ## - Warren, Hacker's Delight, 2013 - ## ## For now only LittleEndian is implemented + # + # Resources at the bottom of the file # Find the most significant word with actual set bits # and get the leading zero count there @@ -92,9 +91,10 @@ func knuthDivLE[qLen, rLen, uLen, vLen: static int]( else: divisorLen -= 1 - doAssert msw != 0, "Division by zero. Abandon ship!" + doAssert divisorLen != 0, "Division by zero. Abandon ship!" - if mswLen == 1: + # Divisor is a single word. + if divisorLen == 1: q.copyFrom(u) r.leastSignificantWord() = q.shortDiv(v.leastSignificantWord()) # zero all but the least significant word @@ -114,23 +114,64 @@ func knuthDivLE[qLen, rLen, uLen, vLen: static int]( un.shlSmallOverflowing(u, clz) vn.shlSmall(v, clz) - static: doAssert cpuEndian == littleEndian, "As it is the division algorithm requires little endian ordering of the limbs". + static: doAssert cpuEndian == littleEndian, "Currently the division algorithm requires little endian ordering of the limbs" # TODO: is it worth it to have the uint be the exact same extended precision representation # as a wide int (say uint128 or uint256)? # in big-endian, the following loop must go the other way and the -1 must be +1 + + let vhi = vn[divisorLen-1] + let vlo = vn[divisorLen-2] + for j in countdown(uLen - divisorLen, 0, 1): # Compute qhat estimate of q[j] (off by 0, 1 and rarely 2) var qhat, rhat: Word - let hi = un[j+divisorLen] - let lo = un[j+divisorLen-1] - div2n1n(qhat, rhat, hi, lo, vn[divisorLen-1]) + let uhi = un[j+divisorLen] + let ulo = un[j+divisorLen-1] + div2n1n(qhat, rhat, uhi, ulo, vhi) + var mhi, mlo: Word + var rhi, rlo: Word + mul(mhi, mlo, qhat, vlo) + rhi = rhat + rlo = ulo + # if r < m, adjust approximation, up to twice + while rhi < mhi or (rhi == mhi and rlo < mlo): + qhat -= 1 + rhi += vhi + + # Found the quotient + q[j] = qhat + + # un -= qhat * v + var borrow = Borrow(0) + var qvhi, qvlo: Word + for i in 0 ..< divisorLen-1: + mul(qvhi, qvlo, qhat, v[i]) + subB(borrow, un[j+i], un[j+i], qvlo, borrow) + subB(borrow, un[j+i+1], un[j+i+1], qvhi, borrow) + # Last step + mul(qvhi, qvlo, qhat, v[divisorLen-1]) + subB(borrow, un[j+divisorLen-1], un[j+divisorLen-1], qvlo, borrow) + qvhi += Word(borrow) + let isNeg = un[j+divisorLen] < qvhi + un[j+divisorLen] -= qvhi + + if isNeg: + # oops, too big by one, add back + q[j] -= 1 + var carry = Carry(0) + for i in 0 ..< divisorLen: + addC(carry, u[j+i], u[j+i], v[i], carry) + + # Quotient is found, if remainder is needed we need to un-normalize un + if needRemainder: + r.shrSmall(un, clz) const BinaryShiftThreshold = 8 # If the difference in bit-length is below 8 # binary shift is probably faster -func divmod*[T](x, y: UintImpl[T]): tuple[quot, rem: UintImpl[T]]= - +func divmod(q, r: var Stuint, + x: Limbs[xLen], y: Limbs[yLen], needRemainder: bool) = let x_clz = x.leadingZeros() let y_clz = y.leadingZeros() @@ -139,30 +180,99 @@ func divmod*[T](x, y: UintImpl[T]): tuple[quot, rem: UintImpl[T]]= raise newException(DivByZeroDefect, "You attempted to divide by zero") elif y_clz == (bitsof(y) - 1): # y is one - result.quot = x - elif (x.hi or y.hi).isZero: - # If computing just on the low part is enough - (result.quot.lo, result.rem.lo) = divmod(x.lo, y.lo) - elif (y and (y - one(type y))).isZero: - # y is a power of 2. (this also matches 0 but it was eliminated earlier) - # TODO. Would it be faster to use countTrailingZero (ctz) + clz == size(y) - 1? - # Especially because we shift by ctz after. - let y_ctz = bitsof(y) - y_clz - 1 - result.quot = x shr y_ctz - result.rem = x and (y - one(type y)) + q = x + # elif (x.hi or y.hi).isZero: + # # If computing just on the low part is enough + # (result.quot.lo, result.rem.lo) = divmod(x.lo, y.lo, needRemainder) + # elif (y and (y - one(type y))).isZero: + # # y is a power of 2. (this also matches 0 but it was eliminated earlier) + # # TODO. Would it be faster to use countTrailingZero (ctz) + clz == size(y) - 1? + # # Especially because we shift by ctz after. + # let y_ctz = bitsof(y) - y_clz - 1 + # result.quot = x shr y_ctz + # if needRemainder: + # result.rem = x and (y - one(type y)) elif x == y: - result.quot.lo = one(T) + q.setOne() elif x < y: - result.rem = x - elif (y_clz - x_clz) < BinaryShiftThreshold: - binaryShiftDiv(x, y, result.quot, result.rem) + r = x + # elif (y_clz - x_clz) < BinaryShiftThreshold: + # binaryShiftDiv(x, y, result.quot, result.rem) else: - divmodBZ(x, y, result.quot, result.rem) + knuthDivLE(q, r, x, y, needRemainder) -func `div`*(x, y: UintImpl): UintImpl {.inline.} = +func `div`*(x, y: Stuint): Stuint {.inline.} = ## Division operation for multi-precision unsigned uint - divmod(x,y).quot + var tmp{.noInit.}: Stuint + divmod(result, tmp, x, y, needRemainder = false) -func `mod`*(x, y: UintImpl): UintImpl {.inline.} = - ## Division operation for multi-precision unsigned uint - divmod(x,y).rem +func `mod`*(x, y: Stuint): Stuint {.inline.} = + ## Remainder operation for multi-precision unsigned uint + var tmp{.noInit.}: Stuint + divmod(tmp, result, x,y, needRemainder = true) + +func divmod*(x, y: Stuint): tuple[quot, rem: Stuint] = + ## Division and remainder operations for multi-precision unsigned uint + divmod(result.quot, result.rem, x, y, needRemainder = true) + +# ###################################################################### +# 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 \ No newline at end of file diff --git a/stint/uintops.nim b/stint/uintops.nim index e89987c..738a1b9 100644 --- a/stint/uintops.nim +++ b/stint/uintops.nim @@ -13,6 +13,8 @@ import ./private/uint_bitwise, ./private/uint_shift, ./private/uint_addsub, + ./private/uint_mul, + ./private/uint_div, ./private/primitives/addcarry_subborrow export StUint @@ -171,7 +173,6 @@ export `+=` # - It's implemented at the limb-level so that # in the future Stuint[254] and Stuint256] share a common codepath -import ./private/uint_mul {.push raises: [], inline, noInit, gcsafe.} func `*`*(a, b: Stuint): Stuint = @@ -227,3 +228,5 @@ func pow*[aBits, eBits](a: Stuint[aBits], e: Stuint[eBits]): Stuint[aBits] = # Division & Modulo # -------------------------------------------------------- + +export uint_div \ No newline at end of file