diff --git a/src/private/utils.nim b/src/private/utils.nim index 964dd0f..55f0d4c 100644 --- a/src/private/utils.nim +++ b/src/private/utils.nim @@ -21,4 +21,33 @@ proc bit_length*[T: Natural](n: T): int {.noSideEffect.}= var y: T = n shr 1 while y > 0.T: y = y shr 1 - inc(result) \ No newline at end of file + inc(result) + +proc asUint*[T: MpUInt](n: T): auto {.noSideEffect, inline.}= + ## Casts a multiprecision integer to an uint of the same size + + when T.sizeof > 8: + raise newException("Unreachable. You are trying to cast a MpUint with more than 64-bit of precision") + elif T.sizeof == 8: + cast[uint64](n) + elif T.sizeof == 4: + cast[uint32](n) + elif T.sizeof == 2: + cast[uint16](n) + else: + raise newException("Unreachable. MpUInt must be 16-bit minimum and a power of 2") + +proc asUint*[T: SomeUnsignedInt](n: T): T {.noSideEffect, inline.}= + ## No-op overload of multi-precision int casting + n + +proc asDoubleUint*[T: BaseUint](n: T): auto {.noSideEffect, inline.} = + ## Convert an integer or MpUint to an uint with double the size + + type Double = ( + when T.sizeof == 4: uint64 + elif T.sizeof == 2: uint32 + else: uint16 + ) + + n.asUint.Double \ No newline at end of file diff --git a/src/uint_binary_ops.nim b/src/uint_binary_ops.nim index ac8d89d..50779ef 100644 --- a/src/uint_binary_ops.nim +++ b/src/uint_binary_ops.nim @@ -1,10 +1,11 @@ # Copyright (c) 2018 Status Research & Development GmbH # Distributed under the MIT License (license terms are at http://opensource.org/licenses/MIT). -import uint_type +import ./private/utils, + uint_type proc `+=`*[T: MpUint](a: var T, b: T) {.noSideEffect.}= - # In-place addition for multi-precision unsigned int + ## In-place addition for multi-precision unsigned int # # Optimized assembly should contain adc instruction (add with carry) # Clang on MacOS does with the -d:release switch and MpUint[uint32] (uint64) @@ -20,17 +21,74 @@ proc `+`*[T: MpUint](a, b: T): T {.noSideEffect, noInit, inline.}= result += b proc `-=`*[T: MpUint](a: var T, b: T) {.noSideEffect.}= - # In-place substraction for multi-precision unsigned int + ## In-place substraction for multi-precision unsigned int # - # Optimized assembly should contain sbc instruction (substract with carry) + # Optimized assembly should contain sbb instruction (substract with borrow) # Clang on MacOS does with the -d:release switch and MpUint[uint32] (uint64) - type Base = type a.lo + type MPBase = type a.lo let tmp = a.lo a.lo -= b.lo - a.hi -= (a.lo > tmp).Base + b.hi + a.hi -= (a.lo > tmp).MPBase + b.hi proc `-`*[T: MpUint](a, b: T): T {.noSideEffect, noInit, inline.}= # Substraction for multi-precision unsigned int result = a - result -= b \ No newline at end of file + result -= b + +proc karatsuba[T: BaseUint](a, b: T): MpUint[T] {.noSideEffect, noInit, inline.} + # Forward declaration + +proc `*`*[T: MpUint](a, b: T): T {.noSideEffect, noInit.}= + ## Multiplication for multi-precision unsigned uint + # + # We use a modified Karatsuba algorithm + # + # Karatsuba algorithm splits the operand into `hi * B + lo` + # For our representation, it is similar to school grade multiplication + # Consider hi and lo as if they were digits + # + # 12 + # X 15 + # ------ + # 10 lo*lo -> z0 + # 5 hi*lo -> z1 + # 2 lo*hi -> z1 + # 10 hi*hi -- z2 + # ------ + # 180 + # + # If T is a type + # For T * T --> T we don't need to compute z2 as it always overflow + # For T * T --> 2T (uint64 * uint64 --> uint128) we use the full precision Karatsuba algorithm + + result = karatsuba(a.lo, b.lo) + result.hi += (karatsuba(a.hi, b.lo) + karatsuba(a.lo, b.hi)).lo + +template karatsubaImpl[T: MpUint](x, y: T): MpUint[T] = + # https://en.wikipedia.org/wiki/Karatsuba_algorithm + let + z0 = karatsuba(x.lo, y.lo) + tmp = karatsuba(x.hi, y.lo) + + var z1 = tmp + z1 += karatsuba(x.hi, y.lo) + let z2 = (z1 < tmp).T + karatsuba(x.hi, y.hi) + + result.lo = z1.lo shl 32 + z0 + result.hi = z2 + z1.hi + +proc karatsuba[T: BaseUint](a, b: T): MpUint[T] {.noSideEffect, noInit, inline.}= + ## Karatsuba algorithm with full precision + + when T.sizeof in {1, 2, 4}: + # Use types twice bigger to do the multiplication + cast[type result](a.asDoubleUint * b.asDoubleUint) + + elif T.sizeof == 8: # uint64 or MpUint[uint32] + # We cannot double uint64 to uint128 + # We use the Karatsuba algorithm + karatsubaImpl(cast[MpUint[uint32]](a), cast[MpUint[uint32]](b)) + else: + # Case: at least uint128 * uint128 --> uint256 + karatsubaImpl(a, b) \ No newline at end of file diff --git a/tests/all_tests.nim b/tests/all_tests.nim index 1d7fb63..0a007da 100644 --- a/tests/all_tests.nim +++ b/tests/all_tests.nim @@ -2,4 +2,5 @@ # Distributed under the MIT License (license terms are at http://opensource.org/licenses/MIT).$ import test_endianness, - test_addsub \ No newline at end of file + test_addsub, + test_mul \ No newline at end of file diff --git a/tests/test_mul.nim b/tests/test_mul.nim new file mode 100644 index 0000000..64a222d --- /dev/null +++ b/tests/test_mul.nim @@ -0,0 +1,19 @@ +# Copyright (c) 2018 Status Research & Development GmbH +# Distributed under the MIT License (license terms are at http://opensource.org/licenses/MIT). + +import ../src/mpint, unittest + +suite "Testing multiplication implementation": + test "Multiplication with result fitting in low half": + + let a = initMpUint(10000, uint32) + let b = initMpUint(10000, uint32) + + check: cast[uint64](a*b) == 100_000_000'u64 # need 27-bits + + test "Multiplication with result overflowing low half": + + let a = initMpUint(1_000_000, uint32) + let b = initMpUint(1_000_000, uint32) + + check: cast[uint64](a*b) == 1_000_000_000_000'u64 # need 40 bits \ No newline at end of file