Add divmod i.e. division + modulo done!

src/private/utils.nim

@@ -9,27 +9,34 @@ macro getSubType*(T: typedesc): untyped =
   ## MpUint[uint32] --> uint32
   getTypeInst(T)[1][1]
 
-
-proc bit_length*[T: BaseUint](n: T): int {.noSideEffect.}=
+proc bit_length*[T: SomeUnsignedInt](n: T): int {.noSideEffect.}=
   ## Calculates how many bits are necessary to represent the number
   result = 1
   var y: T = n shr 1
-  while y > 0.T:
+  while y != 0.T:
     y = y shr 1
     inc(result)
 
-
 proc bit_length*[T: Natural](n: T): int {.noSideEffect.}=
   ## Calculates how many bits are necessary to represent the number
   #
-  # For some reason using "BaseUint or Natural" directly makes Nim compiler
+  # For some reason using "SomeUnsignedInt or Natural" directly makes Nim compiler
   # throw a type mismatch
   result = 1
   var y: T = n shr 1
-  while y > 0.T:
+  while y != 0.T:
     y = y shr 1
     inc(result)
 
+proc bit_length*[T: MpUint](n: T): int {.noSideEffect.}=
+  ## Calculates how many bits are necessary to represent the number
+  const zero = T()
+  result = 1
+  var y: T = n shr 1
+  while y != zero:
+    y = y shr 1 # TODO: shr is slow!
+    inc(result)
+
 proc asUint*[T: MpUInt](n: T): auto {.noSideEffect, inline.}=
   ## Casts a multiprecision integer to an uint of the same size
 

src/uint_binary_ops.nim

@@ -106,3 +106,48 @@ proc naiveMul[T: BaseUint](x, y: T): MpUint[T] {.noSideEffect, noInit, inline.}=
     # Case: at least uint128 * uint128 --> uint256
     naiveMulImpl(x, y)
 
+
+proc divmod*[T: BaseUint](x, y: T): tuple[quot, rem: T] {.noSideEffect.}=
+  ## Division for multi-precision unsigned uint
+  ## Returns quotient + reminder in a (quot, rem) tuple
+  #
+  # Implementation through binary shift division
+  const zero = T()
+
+  when x.lo is MpUInt:
+    const one = T(lo: getSubType(T)(1))
+    const mpOne = one
+  else:
+    const one: getSubType(T) = 1
+    const mpOne = T(lo: getSubType(T)(1))
+
+  if y == zero:
+    raise newException(DivByZeroError, "You attempted to divide by zero")
+  elif y == mpOne:
+    result.quot = x
+    return
+
+  var
+    shift = x.bit_length - y.bit_length
+    d = y shl shift
+
+  result.rem  = x
+
+  while shift >= 0:
+    result.quot = result.quot shl 1
+    if result.rem >= d:
+      result.rem -= d
+      result.quot.lo = result.quot.lo or one
+
+    d = d shr 1
+    dec(shift)
+
+  # Performance note:
+  # The performance of this implementation is extremely dependant on shl and shr.
+  #
+  # Probably the most efficient algorithm that can benefit from MpUInt data structure is
+  # the recursive fast 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.
+  #  - Comparison of fast division algorithms fro large integers: http://bioinfo.ict.ac.cn/~dbu/AlgorithmCourses/Lectures/Hasselstrom2003.pdf

src/uint_bitwise_ops.nim

@@ -34,7 +34,7 @@ proc `shl`*[T: MpUint](x: T, y: SomeInteger): T {.noInit, noSideEffect.}=
   if y == 0:
     return x
 
-  let # cannot be const, compile-time sizeof only works for simple types
+  let
     size = T.sizeof * 8
     halfSize = size div 2
 
@@ -53,7 +53,7 @@ proc `shr`*[T: MpUint](x: T, y: SomeInteger): T {.noInit, noSideEffect.}=
   if y == 0:
     return x
 
-  let # cannot be const, compile-time sizeof only works for simple types
+  let
     size = T.sizeof * 8
     halfSize = size div 2
 
@@ -97,7 +97,7 @@ proc `shr`*[T: MpUint](x: T, y: SomeInteger): T {.noInit, noSideEffect.}=
 #   result.hi = (x.lo shl S) and not M2
 #   result.lo = (x.lo shl S) and M2
 #   result.hi = result.hi or ((
-#     x.hi shl S or (x.lo shr (size - S) and M1)
+#     x.hi shl S or (x.lo shr (halfSize - S) and M1)
 #   ) and M2)
 
 # proc `shr`*[T: MpUint](x: T, y: SomeInteger): T {.noInit, noSideEffect.}=
@@ -119,5 +119,5 @@ proc `shr`*[T: MpUint](x: T, y: SomeInteger): T {.noInit, noSideEffect.}=
 #   result.lo = (x.hi shr S) and not M2
 #   result.hi = (x.hi shr S) and M2
 #   result.lo = result.lo or ((
-#     x.lo shr S or (x.lo shl (size - S) and M1)
+#     x.lo shr S or (x.hi shl (halfSize - S) and M1)
 #   ) and M2)