From bdabfe8ec0d78d212e4b59c0940b4d955c86bbb1 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Fri, 14 Jan 2022 21:30:23 +0000 Subject: [PATCH 01/12] First implementation of untested Karatsuba --- src/bigints.nim | 78 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/src/bigints.nim b/src/bigints.nim index a52069e..542350e 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -381,6 +381,9 @@ template `-=`*(a: var BigInt, b: BigInt) = assert a == 3.initBigInt a = a - b +func abs*(a: BigInt): BigInt = + result = a + result.isNegative = false func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = # always called with bl >= cl @@ -411,6 +414,35 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = inc pos normalize(a) +func scalarMultiplication(a: var BigInt, b: uint32, c: BigInt) {.inline.} = + # Based on unsignedMultiplication + let + cl = c.limbs.len + a.limbs.setLen(1 + cl) + var tmp = 0'u64 + + tmp += uint64(b) * uint64(c.limbs[0]) + a.limbs[1] = uint32(tmp and uint32.high) + tmp = tmp shr 32 # carry + + a.limbs[1] = uint32(tmp) + + for j in 1 ..< cl: + tmp = 0'u64 + tmp += uint64(a.limbs[j]) + uint64(b) * uint64(c.limbs[j]) + a.limbs[j] = uint32(tmp and uint32.high) + tmp = tmp shr 32 + var pos = j + 1 + while tmp > 0'u64: + tmp += uint64(a.limbs[pos]) + a.limbs[pos] = uint32(tmp and uint32.high) + tmp = tmp shr 32 + inc pos + normalize(a) + +# forward declaration for use in `multiplication` +# func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} + func multiplication(a: var BigInt, b, c: BigInt) = # a = b * c if b.isZero or c.isZero: @@ -426,6 +458,52 @@ func multiplication(a: var BigInt, b, c: BigInt) = unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative +func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = + let + bl = b.limbs.len + cl = c.limbs.len + let n = max(bl, cl) + if bl == 1: + # base case : multiply the only limb with each limb of second term + var a: BigInt + unsignedMultiplication(a, c, b) + return + if cl == 1: + var a: BigInt + unsignedMultiplication(a, b, c) + return + let k = n shr 1 # should it be ceil(n/2) ? + var + low_b, high_b, low_c, high_c: BigInt + # Decompose `b` and `c` in two parts of (almost) equal length + low_b.limbs = b.limbs[0 .. k-1] + high_b.limbs = b.limbs[k .. ^1] + low_c.limbs = c.limbs[0 .. k-1] + high_c.limbs = c.limbs[k .. ^1] + + # subtractive version of Karatsuba's algorithm : + # limit carry handling in opposition to the additive version + var + lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero + unsignedKaratsubaMultiplication(lowProduct, low_b, low_c) + unsignedKaratsubaMultiplication(highProduct, high_b, high_c) + A3 = low_b - high_b # Additive variant of Karatsuba + A4 = low_c - high_c # would add them + let sign = A3.isNegative xor A4.isNegative + if A4.limbs.len >= A3.limbs.len: + multiplication(A5, abs(A4), abs(A3)) + else: + multiplication(A5, abs(A3), abs(A4)) + if sign: + middleTerm = lowProduct + highProduct - A5 + else: + middleTerm = lowProduct + highProduct + A5 + # result = lowProduct + middleTerm shr k + highProduct shr 2k + a.limbs[0 .. k - 1] = lowProduct.limbs + a.limbs[k .. 2*k-1] = middleTerm.limbs + a.limbs[2*k .. 3*k-1] = highProduct.limbs + + func `*`*(a, b: BigInt): BigInt = ## Multiplication for `BigInt`s. runnableExamples: From 244a6a8a39b9d90fb9f4e15704c8195481581507 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Fri, 14 Jan 2022 21:45:30 +0000 Subject: [PATCH 02/12] Call Karatsuba with a treshold - does not compile --- src/bigints.nim | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 542350e..cb04080 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -441,7 +441,7 @@ func scalarMultiplication(a: var BigInt, b: uint32, c: BigInt) {.inline.} = normalize(a) # forward declaration for use in `multiplication` -# func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} +func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} func multiplication(a: var BigInt, b, c: BigInt) = # a = b * c @@ -451,11 +451,18 @@ func multiplication(a: var BigInt, b, c: BigInt) = let bl = b.limbs.len cl = c.limbs.len + karatsubaTreshold = 5 if cl > bl: - unsignedMultiplication(a, c, b) + if bl <= karatsubaTreshold: + unsignedKaratsubaMultiplication(a, c, b) + else: + unsignedMultiplication(a, c, b) else: - unsignedMultiplication(a, b, c) + if cl <= karatsubaTreshold: + unsignedKaratsubaMultiplication(a, b, c) + else: + unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = From f2a48e245731da4bf914b0bb51e424975806319d Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Sat, 29 Jan 2022 20:30:57 +0000 Subject: [PATCH 03/12] Made karatsuba treshold a global const --- src/bigints.nim | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bigints.nim b/src/bigints.nim index cb04080..b4a6500 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -65,6 +65,7 @@ func initBigInt*(val: BigInt): BigInt = const zero = initBigInt(0) one = initBigInt(1) + karatsubaTreshold = 5 func isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): @@ -451,7 +452,6 @@ func multiplication(a: var BigInt, b, c: BigInt) = let bl = b.limbs.len cl = c.limbs.len - karatsubaTreshold = 5 if cl > bl: if bl <= karatsubaTreshold: From f3bcc9badc36efdbc2d34349a9ca46b7e84e2882 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Sat, 29 Jan 2022 20:53:10 +0000 Subject: [PATCH 04/12] Fixed some expressions --- src/bigints.nim | 44 +++++++++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 9e8a99d..2b0ab82 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -389,10 +389,6 @@ template `-=`*(a: var BigInt, b: BigInt) = assert a == 3.initBigInt a = a - b -func abs*(a: BigInt): BigInt = - result = a - result.isNegative = false - func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = # always called with bl >= cl let @@ -450,6 +446,7 @@ func scalarMultiplication(a: var BigInt, b: uint32, c: BigInt) {.inline.} = # forward declaration for use in `multiplication` func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} +func `shl`*(x: BigInt, y: Natural): BigInt func multiplication(a: var BigInt, b, c: BigInt) = # a = b * c @@ -472,6 +469,7 @@ func multiplication(a: var BigInt, b, c: BigInt) = unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative +func `shr`*(x: BigInt, y: Natural): BigInt func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = let bl = b.limbs.len @@ -479,13 +477,23 @@ func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = let n = max(bl, cl) if bl == 1: # base case : multiply the only limb with each limb of second term - var a: BigInt - unsignedMultiplication(a, c, b) + scalarMultiplication(a, b.limbs[0], c) return if cl == 1: - var a: BigInt - unsignedMultiplication(a, b, c) - return + scalarMultiplication(a, c.limbs[0], b) + return + if bl < karatsubaTreshold: + if cl <= bl: + unsignedMultiplication(a, b, c) + else: + unsignedMultiplication(a, c, b) + return + if cl < karatsubaTreshold: + if bl <= cl: + unsignedMultiplication(a, c, b) + else: + unsignedMultiplication(a, b, c) + return let k = n shr 1 # should it be ceil(n/2) ? var low_b, high_b, low_c, high_c: BigInt @@ -502,20 +510,18 @@ func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = unsignedKaratsubaMultiplication(lowProduct, low_b, low_c) unsignedKaratsubaMultiplication(highProduct, high_b, high_c) A3 = low_b - high_b # Additive variant of Karatsuba - A4 = low_c - high_c # would add them - let sign = A3.isNegative xor A4.isNegative + A4 = high_c - low_c # would add them if A4.limbs.len >= A3.limbs.len: multiplication(A5, abs(A4), abs(A3)) else: multiplication(A5, abs(A3), abs(A4)) - if sign: - middleTerm = lowProduct + highProduct - A5 - else: - middleTerm = lowProduct + highProduct + A5 - # result = lowProduct + middleTerm shr k + highProduct shr 2k - a.limbs[0 .. k - 1] = lowProduct.limbs - a.limbs[k .. 2*k-1] = middleTerm.limbs - a.limbs[2*k .. 3*k-1] = highProduct.limbs + middleTerm = lowProduct + highProduct + A5 + a = lowProduct + (middleTerm shr k) + (highProduct shr (2*k)) + # We could affect directly some of the bits of the result with slicing + # a.limbs[0 .. k - 1] = lowProduct.limbs + # But the following instructions would not be correct due to sign handling + # a.limbs[k .. 2*k-1] = middleTerm.limbs + # a.limbs[2*k .. 3*k-1] = highProduct.limbs func `*`*(a, b: BigInt): BigInt = From bc118cbee35a4b70ceefbf19ba77d0fb8f0b510a Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Mon, 31 Jan 2022 10:27:54 +0000 Subject: [PATCH 05/12] Karatsuba multiplication now works --- src/bigints.nim | 61 +++++++++++++++++++------------------------------ 1 file changed, 24 insertions(+), 37 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 2b0ab82..0b5e701 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -64,7 +64,7 @@ func initBigInt*(val: BigInt): BigInt = const zero = initBigInt(0) one = initBigInt(1) - karatsubaTreshold = 5 + karatsubaTreshold = 10 func isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): @@ -418,35 +418,25 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = inc pos normalize(a) -func scalarMultiplication(a: var BigInt, b: uint32, c: BigInt) {.inline.} = - # Based on unsignedMultiplication +func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = + # always called with bl >= cl let - cl = c.limbs.len - a.limbs.setLen(1 + cl) + bl = b.limbs.len + a.limbs.setLen(bl + 1) var tmp = 0'u64 - tmp += uint64(b) * uint64(c.limbs[0]) - a.limbs[1] = uint32(tmp and uint32.high) - tmp = tmp shr 32 # carry - - a.limbs[1] = uint32(tmp) - - for j in 1 ..< cl: - tmp = 0'u64 - tmp += uint64(a.limbs[j]) + uint64(b) * uint64(c.limbs[j]) - a.limbs[j] = uint32(tmp and uint32.high) + for i in 0 ..< bl: + tmp += uint64(b.limbs[i]) * uint64(c) + a.limbs[i] = uint32(tmp and uint32.high) tmp = tmp shr 32 - var pos = j + 1 - while tmp > 0'u64: - tmp += uint64(a.limbs[pos]) - a.limbs[pos] = uint32(tmp and uint32.high) - tmp = tmp shr 32 - inc pos + + a.limbs[bl] = uint32(tmp) normalize(a) # forward declaration for use in `multiplication` -func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} +func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} func `shl`*(x: BigInt, y: Natural): BigInt +func `shr`*(x: BigInt, y: Natural): BigInt func multiplication(a: var BigInt, b, c: BigInt) = # a = b * c @@ -459,28 +449,27 @@ func multiplication(a: var BigInt, b, c: BigInt) = if cl > bl: if bl <= karatsubaTreshold: - unsignedKaratsubaMultiplication(a, c, b) + karatsubaMultiplication(a, c, b) else: unsignedMultiplication(a, c, b) else: if cl <= karatsubaTreshold: - unsignedKaratsubaMultiplication(a, b, c) + karatsubaMultiplication(a, b, c) else: unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative -func `shr`*(x: BigInt, y: Natural): BigInt -func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = +func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = let bl = b.limbs.len cl = c.limbs.len let n = max(bl, cl) if bl == 1: # base case : multiply the only limb with each limb of second term - scalarMultiplication(a, b.limbs[0], c) + scalarMultiplication(a, c, b.limbs[0]) return if cl == 1: - scalarMultiplication(a, c.limbs[0], b) + scalarMultiplication(a, b, c.limbs[0]) return if bl < karatsubaTreshold: if cl <= bl: @@ -507,21 +496,19 @@ func unsignedKaratsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = # limit carry handling in opposition to the additive version var lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero - unsignedKaratsubaMultiplication(lowProduct, low_b, low_c) - unsignedKaratsubaMultiplication(highProduct, high_b, high_c) + karatsubaMultiplication(lowProduct, low_b, low_c) + karatsubaMultiplication(highProduct, high_b, high_c) A3 = low_b - high_b # Additive variant of Karatsuba - A4 = high_c - low_c # would add them + A4 = low_c - high_c # would add them if A4.limbs.len >= A3.limbs.len: multiplication(A5, abs(A4), abs(A3)) else: multiplication(A5, abs(A3), abs(A4)) middleTerm = lowProduct + highProduct + A5 - a = lowProduct + (middleTerm shr k) + (highProduct shr (2*k)) - # We could affect directly some of the bits of the result with slicing - # a.limbs[0 .. k - 1] = lowProduct.limbs - # But the following instructions would not be correct due to sign handling - # a.limbs[k .. 2*k-1] = middleTerm.limbs - # a.limbs[2*k .. 3*k-1] = highProduct.limbs + a.limbs[0 .. k - 1] = lowProduct.limbs + # a += (middleTerm shr k) + (highProduct shr (2*k)) + a.limbs[k .. 2*k-1] = middleTerm.limbs + a.limbs[2*k .. 3*k-1] = highProduct.limbs func `*`*(a, b: BigInt): BigInt = From 438a26ed8fd6d3e1e545005cb5c28923817ef543 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Mon, 31 Jan 2022 11:59:31 +0000 Subject: [PATCH 06/12] add tests --- src/bigints.nim | 33 ++++++++++++++++++++++----------- tests/tbigints.nim | 4 ++++ 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 0b5e701..ec9fefa 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -448,15 +448,15 @@ func multiplication(a: var BigInt, b, c: BigInt) = cl = c.limbs.len if cl > bl: - if bl <= karatsubaTreshold: - karatsubaMultiplication(a, c, b) - else: - unsignedMultiplication(a, c, b) + # if bl <= karatsubaTreshold: + # karatsubaMultiplication(a, c, b) + # else: + unsignedMultiplication(a, c, b) else: - if cl <= karatsubaTreshold: - karatsubaMultiplication(a, b, c) - else: - unsignedMultiplication(a, b, c) + # if cl <= karatsubaTreshold: + # karatsubaMultiplication(a, b, c) + # else: + unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = @@ -506,9 +506,9 @@ func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = multiplication(A5, abs(A3), abs(A4)) middleTerm = lowProduct + highProduct + A5 a.limbs[0 .. k - 1] = lowProduct.limbs - # a += (middleTerm shr k) + (highProduct shr (2*k)) - a.limbs[k .. 2*k-1] = middleTerm.limbs - a.limbs[2*k .. 3*k-1] = highProduct.limbs + a += (middleTerm shr k) + (highProduct shr (2*k)) + # a.limbs[k .. 2*k-1] = middleTerm.limbs + # a.limbs[2*k .. 3*k-1] = highProduct.limbs func `*`*(a, b: BigInt): BigInt = @@ -1205,3 +1205,14 @@ func powmod*(base, exponent, modulus: BigInt): BigInt = result = (result * basePow) mod modulus basePow = (basePow * basePow) mod modulus exponent = exponent shr 1 + +when isMainModule: + var a = "1311562737969161616".initBigInt + echo a.limbs.len + var b = "1357909330350306889".initBigInt + echo b.limbs.len + a = "1780983279228119273110576463639172624".initBigInt + b = "1843917749452418885995463656480858321".initBigInt + echo a.limbs.len + echo b.limbs.len + echo a*b diff --git a/tests/tbigints.nim b/tests/tbigints.nim index a3df592..76dc62c 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -227,6 +227,10 @@ proc main() = doAssert (d xor d) == f doAssert (d xor f) == d + block: # multiplication + let a = "1780983279228119273110576463639172624".initBigInt + let b = "1843917749452418885995463656480858321".initBigInt + doAssert a * b == "3171901440890145063107180402349133639481893332927709969425239467271045376".initBigInt block: # self-addition/self-subtraction # self-addition var a = zero From 5dd251200d118f8e701e6812ec451e6b32564cd4 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Mon, 31 Jan 2022 18:46:47 +0000 Subject: [PATCH 07/12] Add randomized tests and fixed call to karatsuba --- src/bigints.nim | 45 +++++++++++++++--------------------- tests/fastMultiplication.nim | 25 ++++++++++++++++++++ tests/tbigints.nim | 15 +++++++++++- 3 files changed, 57 insertions(+), 28 deletions(-) create mode 100644 tests/fastMultiplication.nim diff --git a/src/bigints.nim b/src/bigints.nim index ec9fefa..411050b 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -64,7 +64,7 @@ func initBigInt*(val: BigInt): BigInt = const zero = initBigInt(0) one = initBigInt(1) - karatsubaTreshold = 10 + karatsubaTreshold = 2 func isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): @@ -448,15 +448,15 @@ func multiplication(a: var BigInt, b, c: BigInt) = cl = c.limbs.len if cl > bl: - # if bl <= karatsubaTreshold: - # karatsubaMultiplication(a, c, b) - # else: - unsignedMultiplication(a, c, b) + if bl >= karatsubaTreshold: + karatsubaMultiplication(a, c, b) + else: + unsignedMultiplication(a, c, b) else: - # if cl <= karatsubaTreshold: - # karatsubaMultiplication(a, b, c) - # else: - unsignedMultiplication(a, b, c) + if cl >= karatsubaTreshold: + karatsubaMultiplication(a, b, c) + else: + unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = @@ -483,7 +483,7 @@ func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = else: unsignedMultiplication(a, b, c) return - let k = n shr 1 # should it be ceil(n/2) ? + let k = n shr 1 var low_b, high_b, low_c, high_c: BigInt # Decompose `b` and `c` in two parts of (almost) equal length @@ -492,24 +492,19 @@ func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = low_c.limbs = c.limbs[0 .. k-1] high_c.limbs = c.limbs[k .. ^1] - # subtractive version of Karatsuba's algorithm : - # limit carry handling in opposition to the additive version + # subtractive version of Karatsuba's algorithm to limit carry handling var lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero karatsubaMultiplication(lowProduct, low_b, low_c) karatsubaMultiplication(highProduct, high_b, high_c) - A3 = low_b - high_b # Additive variant of Karatsuba - A4 = low_c - high_c # would add them + A3 = low_b - high_b + A4 = high_c - low_c if A4.limbs.len >= A3.limbs.len: multiplication(A5, abs(A4), abs(A3)) else: multiplication(A5, abs(A3), abs(A4)) middleTerm = lowProduct + highProduct + A5 - a.limbs[0 .. k - 1] = lowProduct.limbs - a += (middleTerm shr k) + (highProduct shr (2*k)) - # a.limbs[k .. 2*k-1] = middleTerm.limbs - # a.limbs[2*k .. 3*k-1] = highProduct.limbs - + a = lowProduct + (middleTerm shr k) + (highProduct shr (2*k)) func `*`*(a, b: BigInt): BigInt = ## Multiplication for `BigInt`s. @@ -1207,12 +1202,8 @@ func powmod*(base, exponent, modulus: BigInt): BigInt = exponent = exponent shr 1 when isMainModule: - var a = "1311562737969161616".initBigInt - echo a.limbs.len - var b = "1357909330350306889".initBigInt - echo b.limbs.len - a = "1780983279228119273110576463639172624".initBigInt - b = "1843917749452418885995463656480858321".initBigInt - echo a.limbs.len - echo b.limbs.len + let a = "1780983279228119273110576463639172624".initBigInt + let b = "1843917749452418885995463656480858321".initBigInt + echo a.limbs + echo b.limbs echo a*b diff --git a/tests/fastMultiplication.nim b/tests/fastMultiplication.nim new file mode 100644 index 0000000..cd9323a --- /dev/null +++ b/tests/fastMultiplication.nim @@ -0,0 +1,25 @@ +import bigints +import std/[math, random, sequtils, strutils] + +randomize() +# Pick a number in 0..100. +let limit = 10^9 +let limbs = 10000 +let randomBigInt = toSeq(1..limbs).mapIt(rand(limit)).join("").initBigInt +let randomBigInt2 = toSeq(1..limbs).mapIt(rand(limit)).join("").initBigInt +let randomBigInt3 = toSeq(1..limbs).mapIt(rand(limit)).join("").initBigInt + +# Compute subproducts +let prod1 = randomBigInt * randomBigInt2 +let prod1bis = randomBigInt2 * randomBigInt + +# Check commutativity of the product +doAssert prod1 == prod1bis + +let prod2 = randomBigInt2 * randomBigInt3 +let prod3 = randomBigInt * randomBigInt3 +let product = prod2 * randomBigInt + +# Check associativity of the product +doAssert prod1 * randomBigInt3 == product +doAssert prod3 * randomBigInt2 == product diff --git a/tests/tbigints.nim b/tests/tbigints.nim index 76dc62c..cf19b31 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -228,9 +228,22 @@ proc main() = doAssert (d xor f) == d block: # multiplication + # factors with 4 limbs let a = "1780983279228119273110576463639172624".initBigInt let b = "1843917749452418885995463656480858321".initBigInt - doAssert a * b == "3171901440890145063107180402349133639481893332927709969425239467271045376".initBigInt + echo a * b + doAssert a * b == "3283986680046702618742503890385314117448805445290098330749803441805804304".initBigInt + + # factors with 17 limbs + let c = "15456863493948186026689401110531937466657435954521677287549013772194751214595085262021623597960658907994197330891108031896474775438991400654520526954653285".initBigInt + let d = "20867311096234429137120990056519061484140179793024844459539745043528236531589522382271230666075358518275274769618792229717222657110424037636116966396665200".initBigInt + doAssert c * d == "322543179100245850295291096700090285623536165554432133161470913224665233565153206743023505404409261647296075477738317301701554637184306640109864382144645081119052516436652162825894456855767719709860985552674755702938369565636714472650667032224717209489767579823588160939485446085000195032327964706246225182000".initBigInt + + # factors with 65 limbs + let e = initBigInt("1f3b839241b0aacc183858dc7a75a773e7bad642a9f426ef499d91e09c9f99a88ec9a14d5ee51175faeaa10d2fb06f3ee37d2f50fe755c2c963aeb539cd55c0e14f5a23f04c64839c22bd4108034b7afc95e01a1c2fe605d8b1930926e886a8f3d7fc09acd54d388cb5d4b3a3fb4eaf6781173ab3a0cd8ad3119c37dd2cf05544235d7b85b2c96d2ed29e1a685820c4afdd824bd8878f1b6a3f52a57eb886efaa737af47161c89f298d908aa950979b8c2615d4e03b47ee87a5381ca39d9ec4788d7abd07b174913b962c02cdd5f8319722a3345eb38d3ebdd51dec66a58e89902151539298c41446758bac66923c910fd7a2d12d0d5c8bb688970b8a77e7d5fc", base = 16) + let f = initBigInt("116be4e445ea68066fca5652e472eb1c1a5fb850311126a8a91fd6f1199f92a9d6602a81bb5e500d163b01df7eec15e41109c62f6c83425027272823d9888a51c93422d47ba4e1cecb94d6fa02eb27df537038b2ac7c9c264634b8febc452c5c9043ddc7d2eddd04f1f743d85cfefd0864441bb9cbf46308138d037b2057980c8b6d215e6bd2c0a73d64c176b7b59452a2c7968a121e5d46c859d85678acd0d4927509418fd351331791fdb7ad041ee2e7d5975867f1812d1e17a41f5a7c0735ebe224294f5ca9d607e95a8722adf58f676b23da1563ac62a52352c10efd0bf5cc8b5eb9ddd7fd1c22bc307d5edf86474f66302fd7cd288a9bc251d3ba45a851b", base = 16) + doAssert e * f == initBigInt("2201d898f354480704067f8513e7b8b365db60c9132d96e54b8546c05c59a08f50b8d8488841d9e893b3c6de34e22b70f2f6c9bc682700c060b60b6e614c9cc39a2c9a9c13ccc412a512c8126f60e1572d20281855e63d019be43a34a929c20818d05527e75f9edf4c5c9a096c4412f879ca14dd3f8bc48aafea6ce17223adefc1e55ca8216ed3d6d351f08dc38a5e6fed7cb3abd1844bfaed632a1571d0017f1285b2e7762c8c3ad0adc781f47df619f462415b8b4e3496fea3d1b1f4443184a5f5fdc155af2f62d861c9ab321a5083c07b4fdfc384aaa6c4a09559e1d7383b5a3fd9f6c9dfe2079c13bfcf307f98de6e16e474b55e94dcf9dbaf3e90ba38f37ddc69c318ed680a0db8a5d5257a1214c765b267f3df8352ad5b4ee9e9c27ca7085ee5687e6afc8108b683c622613c003c068b60bed656d2dc6a32b1f7b079194108301da8d8049f5b64e88da091803bf1582a45fd24f242f2d9b9d8090c4d088ea31faa3e997d20688481b8f1847524f28153e0ba5c3017338cf470c906d8b27352082741dcfb81ec81a3268569424f791c9a82777d66f2b6a52e0653843057da444b55ddb5f517f1676daafa3413ee3dc6dc0a8edc7cfcbd4b6dd5653957baa93e35fad6908addc018706b0acf64cd5cfe3ee462b57e2cc58f641a883a693505fab131a8f51f22cc34ff694af1e4ca0b4d8469067a76b378783b190c2377d94", base = 16) + block: # self-addition/self-subtraction # self-addition var a = zero From dc69b3c408e4337a8afc268b90e71b2f6d7d3dce Mon Sep 17 00:00:00 2001 From: dlesnoff <54949944+dlesnoff@users.noreply.github.com> Date: Fri, 11 Feb 2022 21:09:04 +0100 Subject: [PATCH 08/12] Remove overflow error in scalar multiplication Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com> --- src/bigints.nim | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 0b5e701..5d8d3a7 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -419,14 +419,16 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = normalize(a) func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = - # always called with bl >= cl + if c == 0: + return zero let bl = b.limbs.len a.limbs.setLen(bl + 1) var tmp = 0'u64 + let c = uint64(c) for i in 0 ..< bl: - tmp += uint64(b.limbs[i]) * uint64(c) + tmp += uint64(b.limbs[i]) * c a.limbs[i] = uint32(tmp and uint32.high) tmp = tmp shr 32 From fd4d321ba051b6b2c62ac8726f00be2a4ef71c71 Mon Sep 17 00:00:00 2001 From: dlesnoff <54949944+dlesnoff@users.noreply.github.com> Date: Fri, 11 Feb 2022 21:19:36 +0100 Subject: [PATCH 09/12] treshold -> threshold Co-authored-by: konsumlamm <44230978+konsumlamm@users.noreply.github.com> --- src/bigints.nim | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bigints.nim b/src/bigints.nim index 5d8d3a7..d07701e 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -64,7 +64,7 @@ func initBigInt*(val: BigInt): BigInt = const zero = initBigInt(0) one = initBigInt(1) - karatsubaTreshold = 10 + karatsubaThreshold = 10 func isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): From bb05ee5a1a4e1842aeafeee76578e82e569c77d1 Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Fri, 11 Feb 2022 20:21:19 +0000 Subject: [PATCH 10/12] Many changes I forgot to commit --- src/bigints.nim | 198 +++++++++++++++++++++++---------------------- tests/tbigints.nim | 11 +++ 2 files changed, 114 insertions(+), 95 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 411050b..79f8ee2 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -8,14 +8,14 @@ type isNegative: bool -func normalize(a: var BigInt) = +proc normalize(a: var BigInt) = for i in countdown(a.limbs.high, 0): if a.limbs[i] > 0'u32: a.limbs.setLen(i+1) return a.limbs.setLen(1) -func initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = +proc initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = ## Initializes a `BigInt` from a sequence of `uint32` values. runnableExamples: let a = @[10'u32, 2'u32].initBigInt @@ -24,7 +24,7 @@ func initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = result.limbs = vals result.isNegative = isNegative -func initBigInt*[T: int8|int16|int32](val: T): BigInt = +proc initBigInt*[T: int8|int16|int32](val: T): BigInt = if val < 0: result.limbs = @[(not val).uint32 + 1] # manual 2's complement (to avoid overflow) result.isNegative = true @@ -32,10 +32,10 @@ func initBigInt*[T: int8|int16|int32](val: T): BigInt = result.limbs = @[val.uint32] result.isNegative = false -func initBigInt*[T: uint8|uint16|uint32](val: T): BigInt = +proc initBigInt*[T: uint8|uint16|uint32](val: T): BigInt = result.limbs = @[val.uint32] -func initBigInt*(val: int64): BigInt = +proc initBigInt*(val: int64): BigInt = var a = val.uint64 if val < 0: a = not a + 1 # 2's complement @@ -45,7 +45,7 @@ func initBigInt*(val: int64): BigInt = else: result.limbs = @[a.uint32] -func initBigInt*(val: uint64): BigInt = +proc initBigInt*(val: uint64): BigInt = if val > uint32.high: result.limbs = @[(val and uint32.high).uint32, (val shr 32).uint32] else: @@ -58,7 +58,7 @@ else: template initBigInt*(val: int): BigInt = initBigInt(val.int64) template initBigInt*(val: uint): BigInt = initBigInt(val.uint64) -func initBigInt*(val: BigInt): BigInt = +proc initBigInt*(val: BigInt): BigInt = result = val const @@ -66,13 +66,13 @@ const one = initBigInt(1) karatsubaTreshold = 2 -func isZero(a: BigInt): bool {.inline.} = +proc isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): if a.limbs[i] != 0'u32: return false return true -func abs*(a: BigInt): BigInt = +proc abs*(a: BigInt): BigInt = # Returns the absolute value of `a`. runnableExamples: assert abs(42.initBigInt) == 42.initBigInt @@ -80,16 +80,16 @@ func abs*(a: BigInt): BigInt = result = a result.isNegative = false -func unsignedCmp(a: BigInt, b: uint32): int64 = +proc unsignedCmp(a: BigInt, b: uint32): int64 = # ignores the sign of `a` # `a` and `b` are assumed to not be zero result = int64(a.limbs.len) - 1 if result != 0: return result = int64(a.limbs[0]) - int64(b) -func unsignedCmp(a: uint32, b: BigInt): int64 = -unsignedCmp(b, a) +proc unsignedCmp(a: uint32, b: BigInt): int64 = -unsignedCmp(b, a) -func unsignedCmp(a, b: BigInt): int64 = +proc unsignedCmp(a, b: BigInt): int64 = # ignores the signs of `a` and `b` # `a` and `b` are assumed to not be zero result = int64(a.limbs.len) - int64(b.limbs.len) @@ -99,7 +99,7 @@ func unsignedCmp(a, b: BigInt): int64 = if result != 0: return -func cmp(a, b: BigInt): int64 = +proc cmp(a, b: BigInt): int64 = ## Returns: ## * a value less than zero, if `a < b` ## * a value greater than zero, if `a > b` @@ -122,7 +122,7 @@ func cmp(a, b: BigInt): int64 = else: return unsignedCmp(a, b) -func cmp(a: BigInt, b: int32): int64 = +proc cmp(a: BigInt, b: int32): int64 = ## Returns: ## * a value less than zero, if `a < b` ## * a value greater than zero, if `a > b` @@ -140,9 +140,9 @@ func cmp(a: BigInt, b: int32): int64 = else: return unsignedCmp(a, b.uint32) -func cmp(a: int32, b: BigInt): int64 = -cmp(b, a) +proc cmp(a: int32, b: BigInt): int64 = -cmp(b, a) -func `==`*(a, b: BigInt): bool = +proc `==`*(a, b: BigInt): bool = ## Compares if two `BigInt` numbers are equal. runnableExamples: let @@ -153,7 +153,7 @@ func `==`*(a, b: BigInt): bool = assert b != c cmp(a, b) == 0 -func `<`*(a, b: BigInt): bool = +proc `<`*(a, b: BigInt): bool = runnableExamples: let a = 5.initBigInt @@ -163,7 +163,7 @@ func `<`*(a, b: BigInt): bool = assert b > c cmp(a, b) < 0 -func `<=`*(a, b: BigInt): bool = +proc `<=`*(a, b: BigInt): bool = runnableExamples: let a = 5.initBigInt @@ -173,16 +173,16 @@ func `<=`*(a, b: BigInt): bool = assert c <= b cmp(a, b) <= 0 -func `==`(a: BigInt, b: int32): bool = cmp(a, b) == 0 -func `<`(a: BigInt, b: int32): bool = cmp(a, b) < 0 -func `<`(a: int32, b: BigInt): bool = cmp(a, b) < 0 +proc `==`(a: BigInt, b: int32): bool = cmp(a, b) == 0 +proc `<`(a: BigInt, b: int32): bool = cmp(a, b) < 0 +proc `<`(a: int32, b: BigInt): bool = cmp(a, b) < 0 template addParts(toAdd) = tmp += toAdd a.limbs[i] = uint32(tmp and uint32.high) tmp = tmp shr 32 -func unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = +proc unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = let bl = b.limbs.len a.limbs.setLen(bl) @@ -193,7 +193,7 @@ func unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = a.limbs.add(uint32(tmp)) a.isNegative = false -func unsignedAddition(a: var BigInt, b, c: BigInt) = +proc unsignedAddition(a: var BigInt, b, c: BigInt) = let bl = b.limbs.len cl = c.limbs.len @@ -213,10 +213,10 @@ func unsignedAddition(a: var BigInt, b, c: BigInt) = a.limbs.add(uint32(tmp)) a.isNegative = false -func negate(a: var BigInt) = +proc negate(a: var BigInt) = a.isNegative = not a.isNegative -func `-`*(a: BigInt): BigInt = +proc `-`*(a: BigInt): BigInt = ## Unary minus for `BigInt`. runnableExamples: let @@ -271,7 +271,7 @@ template realUnsignedSubtraction(a: var BigInt, b, c: BigInt) = normalize(a) assert tmp == 0 -func unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = +proc unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = # `b` is not zero let cmpRes = unsignedCmp(b, c) if cmpRes > 0: @@ -283,7 +283,7 @@ func unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = else: # b == c a = zero -func unsignedSubtraction(a: var BigInt, b, c: BigInt) = +proc unsignedSubtraction(a: var BigInt, b, c: BigInt) = let cmpRes = unsignedCmp(b, c) if cmpRes > 0: realUnsignedSubtraction(a, b, c) @@ -293,7 +293,7 @@ func unsignedSubtraction(a: var BigInt, b, c: BigInt) = else: # b == c a = zero -func additionInt(a: var BigInt, b: BigInt, c: int32) = +proc additionInt(a: var BigInt, b: BigInt, c: int32) = # a = b + c if b.isZero: a = c.initBigInt @@ -309,7 +309,7 @@ func additionInt(a: var BigInt, b: BigInt, c: int32) = else: unsignedAdditionInt(a, b, c.uint32) -func addition(a: var BigInt, b, c: BigInt) = +proc addition(a: var BigInt, b, c: BigInt) = # a = b + c if b.isNegative: if c.isNegative: @@ -323,7 +323,7 @@ func addition(a: var BigInt, b, c: BigInt) = else: unsignedAddition(a, b, c) -func `+`*(a, b: BigInt): BigInt = +proc `+`*(a, b: BigInt): BigInt = ## Addition for `BigInt`s. runnableExamples: let @@ -341,7 +341,7 @@ template `+=`*(a: var BigInt, b: BigInt) = assert a == 7.initBigInt a = a + b -func subtractionInt(a: var BigInt, b: BigInt, c: int32) = +proc subtractionInt(a: var BigInt, b: BigInt, c: int32) = # a = b - c if b.isZero: a = -c.initBigInt @@ -357,7 +357,7 @@ func subtractionInt(a: var BigInt, b: BigInt, c: int32) = else: unsignedSubtractionInt(a, b, c.uint32) -func subtraction(a: var BigInt, b, c: BigInt) = +proc subtraction(a: var BigInt, b, c: BigInt) = # a = b - c if b.isNegative: if c.isNegative: @@ -371,7 +371,7 @@ func subtraction(a: var BigInt, b, c: BigInt) = else: unsignedSubtraction(a, b, c) -func `-`*(a, b: BigInt): BigInt = +proc `-`*(a, b: BigInt): BigInt = ## Subtraction for `BigInt`s. runnableExamples: let @@ -389,7 +389,7 @@ template `-=`*(a: var BigInt, b: BigInt) = assert a == 3.initBigInt a = a - b -func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = +proc unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = # always called with bl >= cl let bl = b.limbs.len @@ -418,7 +418,7 @@ func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = inc pos normalize(a) -func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = +proc scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = # always called with bl >= cl let bl = b.limbs.len @@ -434,11 +434,11 @@ func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = normalize(a) # forward declaration for use in `multiplication` -func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} -func `shl`*(x: BigInt, y: Natural): BigInt -func `shr`*(x: BigInt, y: Natural): BigInt +proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) {.inline.} +proc `shl`*(x: BigInt, y: Natural): BigInt +proc `shr`*(x: BigInt, y: Natural): BigInt -func multiplication(a: var BigInt, b, c: BigInt) = +proc multiplication(a: var BigInt, b, c: BigInt) = # a = b * c if b.isZero or c.isZero: a = zero @@ -448,18 +448,18 @@ func multiplication(a: var BigInt, b, c: BigInt) = cl = c.limbs.len if cl > bl: - if bl >= karatsubaTreshold: - karatsubaMultiplication(a, c, b) - else: - unsignedMultiplication(a, c, b) + # if bl >= karatsubaTreshold: + # karatsubaMultiplication(a, c, b) + # else: + unsignedMultiplication(a, c, b) else: - if cl >= karatsubaTreshold: - karatsubaMultiplication(a, b, c) - else: - unsignedMultiplication(a, b, c) + # if cl >= karatsubaTreshold: + # karatsubaMultiplication(a, b, c) + # else: + unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative -func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = +proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) {.inline.} = let bl = b.limbs.len cl = c.limbs.len @@ -484,6 +484,7 @@ func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = unsignedMultiplication(a, b, c) return let k = n shr 1 + echo k var low_b, high_b, low_c, high_c: BigInt # Decompose `b` and `c` in two parts of (almost) equal length @@ -491,22 +492,25 @@ func karatsubaMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = high_b.limbs = b.limbs[k .. ^1] low_c.limbs = c.limbs[0 .. k-1] high_c.limbs = c.limbs[k .. ^1] + # echo low_b, high_b, low_c, high_c # subtractive version of Karatsuba's algorithm to limit carry handling - var - lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero - karatsubaMultiplication(lowProduct, low_b, low_c) - karatsubaMultiplication(highProduct, high_b, high_c) + var lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero + + multiplication(lowProduct, low_b, low_c) + multiplication(highProduct, high_b, high_c) + + # echo "lowProduct, highProduct: ", lowProduct, highProduct A3 = low_b - high_b A4 = high_c - low_c - if A4.limbs.len >= A3.limbs.len: - multiplication(A5, abs(A4), abs(A3)) - else: - multiplication(A5, abs(A3), abs(A4)) + + multiplication(A5, A4, A3) + middleTerm = lowProduct + highProduct + A5 - a = lowProduct + (middleTerm shr k) + (highProduct shr (2*k)) + # echo "A5 = A4*A3, middleTerm = lP + hP + A4*A3: ", A5, middleTerm + a = lowProduct + middleTerm shl (32*k) + highProduct shl (64*k) -func `*`*(a, b: BigInt): BigInt = +proc `*`*(a, b: BigInt): BigInt = ## Multiplication for `BigInt`s. runnableExamples: let @@ -522,7 +526,7 @@ template `*=`*(a: var BigInt, b: BigInt) = assert a == 150.initBigInt a = a * b -func pow*(x: BigInt, y: Natural): BigInt = +proc pow*(x: BigInt, y: Natural): BigInt = ## Computes `x` to the power of `y`. var base = x var exp = y @@ -535,7 +539,7 @@ func pow*(x: BigInt, y: Natural): BigInt = exp = exp shr 1 base *= base -func `shl`*(x: BigInt, y: Natural): BigInt = +proc `shl`*(x: BigInt, y: Natural): BigInt = ## Shifts a `BigInt` to the left. runnableExamples: let a = 24.initBigInt @@ -558,9 +562,9 @@ func `shl`*(x: BigInt, y: Natural): BigInt = result.limbs.add(uint32(carry shr (32 - b))) # forward declaration for use in `shr` -func dec*(a: var BigInt, b: int = 1) +proc dec*(a: var BigInt, b: int = 1) -func `shr`*(x: BigInt, y: Natural): BigInt = +proc `shr`*(x: BigInt, y: Natural): BigInt = ## Shifts a `BigInt` to the right (arithmetically). runnableExamples: let a = 24.initBigInt @@ -596,17 +600,17 @@ func `shr`*(x: BigInt, y: Natural): BigInt = # normalize result.limbs.setLen(result.limbs.high) -func bitwiseAnd(a: var BigInt, b, c: BigInt) = +proc bitwiseAnd(a: var BigInt, b, c: BigInt) = a.limbs.setLen(min(b.limbs.len, c.limbs.len)) for i in 0 ..< a.limbs.len: a.limbs[i] = b.limbs[i] and c.limbs[i] -func `and`*(a, b: BigInt): BigInt = +proc `and`*(a, b: BigInt): BigInt = ## Bitwise `and` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) bitwiseAnd(result, a, b) -func bitwiseOr(a: var BigInt, b, c: BigInt) = +proc bitwiseOr(a: var BigInt, b, c: BigInt) = # `b` must be smaller than `c` a.limbs.setLen(c.limbs.len) for i in 0 ..< b.limbs.len: @@ -614,7 +618,7 @@ func bitwiseOr(a: var BigInt, b, c: BigInt) = for i in b.limbs.len ..< c.limbs.len: a.limbs[i] = c.limbs[i] -func `or`*(a, b: BigInt): BigInt = +proc `or`*(a, b: BigInt): BigInt = ## Bitwise `or` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) if a.limbs.len <= b.limbs.len: @@ -622,7 +626,7 @@ func `or`*(a, b: BigInt): BigInt = else: bitwiseOr(result, b, a) -func bitwiseXor(a: var BigInt, b, c: BigInt) = +proc bitwiseXor(a: var BigInt, b, c: BigInt) = # `b` must be smaller than `c` a.limbs.setLen(c.limbs.len) for i in 0 ..< b.limbs.len: @@ -630,7 +634,7 @@ func bitwiseXor(a: var BigInt, b, c: BigInt) = for i in b.limbs.len ..< c.limbs.len: a.limbs[i] = c.limbs[i] -func `xor`*(a, b: BigInt): BigInt = +proc `xor`*(a, b: BigInt): BigInt = ## Bitwise `xor` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) if a.limbs.len <= b.limbs.len: @@ -638,13 +642,13 @@ func `xor`*(a, b: BigInt): BigInt = else: bitwiseXor(result, b, a) -func reset(a: var BigInt) = +proc reset(a: var BigInt) = ## Resets a `BigInt` back to the zero value. a.limbs.setLen(1) a.limbs[0] = 0 a.isNegative = false -func unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = +proc unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = q.limbs.setLen(n.limbs.len) r = 0 for i in countdown(n.limbs.high, 0): @@ -653,7 +657,7 @@ func unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = r = uint32(tmp mod d) normalize(q) -func bits(d: uint32): int = +proc bits(d: uint32): int = const bitLengths = [0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5] var d = d @@ -663,7 +667,7 @@ func bits(d: uint32): int = result += bitLengths[int(d)] # From Knuth and Python -func unsignedDivRem(q, r: var BigInt, n, d: BigInt) = +proc unsignedDivRem(q, r: var BigInt, n, d: BigInt) = var nn = n.limbs.len dn = d.limbs.len @@ -773,7 +777,7 @@ func unsignedDivRem(q, r: var BigInt, n, d: BigInt) = q = a normalize(q) -func division(q, r: var BigInt, n, d: BigInt) = +proc division(q, r: var BigInt, n, d: BigInt) = # q = n div d # r = n mod d if d.isZero: @@ -789,11 +793,11 @@ func division(q, r: var BigInt, n, d: BigInt) = r += d q -= one -func `div`*(a, b: BigInt): BigInt = +proc `div`*(a, b: BigInt): BigInt = ## Computes the integer division of two `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. ## - ## If you also need the modulo (remainder), use the `divmod func <#divmod,BigInt,BigInt>`_. + ## If you also need the modulo (remainder), use the `divmod proc <#divmod,BigInt,BigInt>`_. runnableExamples: let a = 17.initBigInt @@ -805,11 +809,11 @@ func `div`*(a, b: BigInt): BigInt = var tmp: BigInt division(result, tmp, a, b) -func `mod`*(a, b: BigInt): BigInt = +proc `mod`*(a, b: BigInt): BigInt = ## Computes the integer modulo (remainder) of two `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. ## - ## If you also need an integer division, use the `divmod func <#divmod,BigInt,BigInt>`_. + ## If you also need an integer division, use the `divmod proc <#divmod,BigInt,BigInt>`_. runnableExamples: let a = 17.initBigInt @@ -821,7 +825,7 @@ func `mod`*(a, b: BigInt): BigInt = var tmp: BigInt division(tmp, result, a, b) -func divmod*(a, b: BigInt): tuple[q, r: BigInt] = +proc divmod*(a, b: BigInt): tuple[q, r: BigInt] = ## Computes both the integer division and modulo (remainder) of two ## `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. @@ -832,7 +836,7 @@ func divmod*(a, b: BigInt): tuple[q, r: BigInt] = assert divmod(a, b) == (3.initBigInt, 2.initBigInt) division(result.q, result.r, a, b) -func countTrailingZeroBits(a: BigInt): int = +proc countTrailingZeroBits(a: BigInt): int = var count = 0 for x in a.limbs: if x == 0: @@ -841,7 +845,7 @@ func countTrailingZeroBits(a: BigInt): int = return count + countTrailingZeroBits(x) return count -func gcd*(a, b: BigInt): BigInt = +proc gcd*(a, b: BigInt): BigInt = ## Returns the greatest common divisor (GCD) of `a` and `b`. runnableExamples: assert gcd(54.initBigInt, 24.initBigInt) == 6.initBigInt @@ -870,7 +874,7 @@ func gcd*(a, b: BigInt): BigInt = v = v shr countTrailingZeroBits(v) -func toSignedInt*[T: SomeSignedInt](x: BigInt): Option[T] = +proc toSignedInt*[T: SomeSignedInt](x: BigInt): Option[T] = ## Converts a `BigInt` number to signed integer, if possible. ## If the `BigInt` doesn't fit in a `T`', returns `none`; ## otherwise returns `some(T)`. @@ -919,7 +923,7 @@ func toSignedInt*[T: SomeSignedInt](x: BigInt): Option[T] = result = some(T(x.limbs[0])) -func calcSizes(): array[2..36, int] = +proc calcSizes(): array[2..36, int] = for i in 2..36: var x = int64(i) while x <= int64(uint32.high) + 1: @@ -931,7 +935,7 @@ const powers = {2, 4, 8, 16, 32} sizes = calcSizes() # `sizes[base]` is the maximum number of digits that fully fit in a `uint32` -func toString*(a: BigInt, base: range[2..36] = 10): string = +proc toString*(a: BigInt, base: range[2..36] = 10): string = ## Produces a string representation of a `BigInt` in a specified ## `base`. ## @@ -996,11 +1000,11 @@ func toString*(a: BigInt, base: range[2..36] = 10): string = result.reverse() -func `$`*(a: BigInt): string = +proc `$`*(a: BigInt): string = ## String representation of a `BigInt` in base 10. toString(a, 10) -func parseDigit(c: char, base: uint32): uint32 {.inline.} = +proc parseDigit(c: char, base: uint32): uint32 {.inline.} = result = case c of '0'..'9': uint32(ord(c) - ord('0')) of 'a'..'z': uint32(ord(c) - ord('a') + 10) @@ -1010,7 +1014,7 @@ func parseDigit(c: char, base: uint32): uint32 {.inline.} = if result >= base: raise newException(ValueError, "Invalid input: " & c) -func initBigInt*(str: string, base: range[2..36] = 10): BigInt = +proc initBigInt*(str: string, base: range[2..36] = 10): BigInt = ## Create a `BigInt` from a string. For invalid inputs, a `ValueError` exception is raised. runnableExamples: let @@ -1072,7 +1076,7 @@ func initBigInt*(str: string, base: range[2..36] = 10): BigInt = when (NimMajor, NimMinor) >= (1, 5): include bigints/private/literals -func inc*(a: var BigInt, b: int = 1) = +proc inc*(a: var BigInt, b: int = 1) = ## Increase the value of a `BigInt` by the specified amount (default: 1). runnableExamples: var a = 15.initBigInt @@ -1087,7 +1091,7 @@ func inc*(a: var BigInt, b: int = 1) = else: a += initBigInt(b) -func dec*(a: var BigInt, b: int = 1) = +proc dec*(a: var BigInt, b: int = 1) = ## Decrease the value of a `BigInt` by the specified amount (default: 1). runnableExamples: var a = 15.initBigInt @@ -1102,12 +1106,12 @@ func dec*(a: var BigInt, b: int = 1) = else: a -= initBigInt(b) -func succ*(a: BigInt, b: int = 1): BigInt = +proc succ*(a: BigInt, b: int = 1): BigInt = ## Returns the `b`-th successor of a `BigInt`. result = a inc(result, b) -func pred*(a: BigInt, b: int = 1): BigInt = +proc pred*(a: BigInt, b: int = 1): BigInt = ## Returns the `b`-th predecessor of a `BigInt`. result = a dec(result, b) @@ -1141,7 +1145,7 @@ iterator `..<`*(a, b: BigInt): BigInt = yield res inc res -func invmod*(a, modulus: BigInt): BigInt = +proc invmod*(a, modulus: BigInt): BigInt = ## Compute the modular inverse of `a` modulo `modulus`. ## The return value is always in the range `[1, modulus-1]` runnableExamples: @@ -1174,7 +1178,7 @@ func invmod*(a, modulus: BigInt): BigInt = raise newException(ValueError, $a & " has no modular inverse modulo " & $modulus) result = ((s0 mod modulus) + modulus) mod modulus -func powmod*(base, exponent, modulus: BigInt): BigInt = +proc powmod*(base, exponent, modulus: BigInt): BigInt = ## Compute modular exponentation of `base` with power `exponent` modulo `modulus`. ## The return value is always in the range `[0, modulus-1]`. runnableExamples: @@ -1202,8 +1206,12 @@ func powmod*(base, exponent, modulus: BigInt): BigInt = exponent = exponent shr 1 when isMainModule: - let a = "1780983279228119273110576463639172624".initBigInt - let b = "1843917749452418885995463656480858321".initBigInt + var a, b, c: BigInt + a.limbs = @[1'u32, 2'u32] + b.limbs = @[3'u32, 4'u32] echo a.limbs echo b.limbs - echo a*b + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b diff --git a/tests/tbigints.nim b/tests/tbigints.nim index cf19b31..4093b0e 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -228,6 +228,9 @@ proc main() = doAssert (d xor f) == d block: # multiplication + let one = 1.initBigInt + let negOne = -1.initBigInt + echo one * negOne # factors with 4 limbs let a = "1780983279228119273110576463639172624".initBigInt let b = "1843917749452418885995463656480858321".initBigInt @@ -237,7 +240,15 @@ proc main() = # factors with 17 limbs let c = "15456863493948186026689401110531937466657435954521677287549013772194751214595085262021623597960658907994197330891108031896474775438991400654520526954653285".initBigInt let d = "20867311096234429137120990056519061484140179793024844459539745043528236531589522382271230666075358518275274769618792229717222657110424037636116966396665200".initBigInt + var r: BigInt = 0.initBigInt + karatsubaMultiplication(r, c, d) + echo r + echo "\n" + echo c * d doAssert c * d == "322543179100245850295291096700090285623536165554432133161470913224665233565153206743023505404409261647296075477738317301701554637184306640109864382144645081119052516436652162825894456855767719709860985552674755702938369565636714472650667032224717209489767579823588160939485446085000195032327964706246225182000".initBigInt + echo "\n" + echo c*d - r + doAssert c*d == r # factors with 65 limbs let e = initBigInt("1f3b839241b0aacc183858dc7a75a773e7bad642a9f426ef499d91e09c9f99a88ec9a14d5ee51175faeaa10d2fb06f3ee37d2f50fe755c2c963aeb539cd55c0e14f5a23f04c64839c22bd4108034b7afc95e01a1c2fe605d8b1930926e886a8f3d7fc09acd54d388cb5d4b3a3fb4eaf6781173ab3a0cd8ad3119c37dd2cf05544235d7b85b2c96d2ed29e1a685820c4afdd824bd8878f1b6a3f52a57eb886efaa737af47161c89f298d908aa950979b8c2615d4e03b47ee87a5381ca39d9ec4788d7abd07b174913b962c02cdd5f8319722a3345eb38d3ebdd51dec66a58e89902151539298c41446758bac66923c910fd7a2d12d0d5c8bb688970b8a77e7d5fc", base = 16) From 290a1d121899713d8624fdb06723bfe304ec82ce Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Fri, 11 Feb 2022 21:46:37 +0000 Subject: [PATCH 11/12] Add tests and last recommandations of the review Karatsuba only works in running time. Strange errors when the tests are executed as static. I commented for the present. Tests should be moved from main to tbigints.nim --- src/bigints.nim | 145 ++++++++++++++++++++++++++++++++++++++------- tests/tbigints.nim | 12 ++-- 2 files changed, 129 insertions(+), 28 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 9b07b0c..352f7af 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -7,7 +7,6 @@ type limbs: seq[uint32] isNegative: bool - proc normalize(a: var BigInt) = for i in countdown(a.limbs.high, 0): if a.limbs[i] > 0'u32: @@ -418,10 +417,11 @@ proc unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = inc pos normalize(a) -func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = +proc scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = # always called with bl >= cl if c == 0: - return zero + a = zero + return let bl = b.limbs.len a.limbs.setLen(bl + 1) @@ -437,7 +437,7 @@ func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = normalize(a) # forward declaration for use in `multiplication` -proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) {.inline.} +proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) proc `shl`*(x: BigInt, y: Natural): BigInt proc `shr`*(x: BigInt, y: Natural): BigInt @@ -451,43 +451,51 @@ proc multiplication(a: var BigInt, b, c: BigInt) = cl = c.limbs.len if cl > bl: - # if bl >= karatsubaTreshold: + # if bl >= karatsubaThreshold: # karatsubaMultiplication(a, c, b) # else: unsignedMultiplication(a, c, b) else: - # if cl >= karatsubaTreshold: + # if cl >= karatsubaThreshold: # karatsubaMultiplication(a, b, c) # else: unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative -proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) {.inline.} = +proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = + if b.isZero or c.isZero: + a = zero + return + a.isNegative = b.isNegative xor c.isNegative let bl = b.limbs.len cl = c.limbs.len - let n = max(bl, cl) + n = max(bl, cl) + k = n shr 1 if bl == 1: # base case : multiply the only limb with each limb of second term scalarMultiplication(a, c, b.limbs[0]) + a.isNegative = b.isNegative xor c.isNegative return if cl == 1: scalarMultiplication(a, b, c.limbs[0]) + a.isNegative = b.isNegative xor c.isNegative return - if bl < karatsubaTreshold: + if bl < karatsubaThreshold: if cl <= bl: unsignedMultiplication(a, b, c) else: unsignedMultiplication(a, c, b) + a.isNegative = b.isNegative xor c.isNegative return - if cl < karatsubaTreshold: + if cl < karatsubaThreshold: if bl <= cl: unsignedMultiplication(a, c, b) else: unsignedMultiplication(a, b, c) + a.isNegative = b.isNegative xor c.isNegative return - let k = n shr 1 - echo k + # echo k var low_b, high_b, low_c, high_c: BigInt # Decompose `b` and `c` in two parts of (almost) equal length @@ -498,20 +506,21 @@ proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) {.inline.} = # echo low_b, high_b, low_c, high_c # subtractive version of Karatsuba's algorithm to limit carry handling - var lowProduct, highProduct, A3, A4, A5, middleTerm: BigInt = zero + var lowProduct, highProduct, add3, add4, add5, middleTerm: BigInt = zero multiplication(lowProduct, low_b, low_c) multiplication(highProduct, high_b, high_c) # echo "lowProduct, highProduct: ", lowProduct, highProduct - A3 = low_b - high_b - A4 = high_c - low_c + add3 = low_b - high_b + add4 = high_c - low_c - multiplication(A5, A4, A3) + multiplication(add5, add4, add3) - middleTerm = lowProduct + highProduct + A5 - # echo "A5 = A4*A3, middleTerm = lP + hP + A4*A3: ", A5, middleTerm + middleTerm = lowProduct + highProduct + add5 + # echo "add5 = add4*add3, middleTerm = lP + hP + add4*add3: ", add5, middleTerm a = lowProduct + middleTerm shl (32*k) + highProduct shl (64*k) + a.isNegative = b.isNegative xor c.isNegative proc `*`*(a, b: BigInt): BigInt = ## Multiplication for `BigInt`s. @@ -877,7 +886,7 @@ proc gcd*(a, b: BigInt): BigInt = v = v shr countTrailingZeroBits(v) -func toInt*[T: SomeInteger](x: BigInt): Option[T] = +proc toInt*[T: SomeInteger](x: BigInt): Option[T] = ## Converts a `BigInt` number to an integer, if possible. ## If the `BigInt` doesn't fit in a `T`, returns `none(T)`; ## otherwise returns `some(x)`. @@ -1172,14 +1181,14 @@ iterator `..<`*(a, b: BigInt): BigInt = yield res inc res -func modulo(a, modulus: BigInt): BigInt = +proc modulo(a, modulus: BigInt): BigInt = ## Like `mod`, but the result is always in the range `[0, modulus-1]`. ## `modulus` should be greater than zero. result = a mod modulus if result < 0: result += modulus -func fastLog2*(a: BigInt): int = +proc fastLog2*(a: BigInt): int = ## Computes the logarithm in base 2 of `a`. ## If `a` is negative, returns the logarithm of `abs(a)`. ## If `a` is zero, returns -1. @@ -1187,7 +1196,7 @@ func fastLog2*(a: BigInt): int = return -1 bitops.fastLog2(a.limbs[^1]) + 32*(a.limbs.high) -func invmod*(a, modulus: BigInt): BigInt = +proc invmod*(a, modulus: BigInt): BigInt = ## Compute the modular inverse of `a` modulo `modulus`. ## The return value is always in the range `[1, modulus-1]` runnableExamples: @@ -1247,6 +1256,47 @@ proc powmod*(base, exponent, modulus: BigInt): BigInt = when isMainModule: var a, b, c: BigInt + let + two = 2.initBigInt + three = 3.initBigInt + four = 4.initBigInt + + a.limbs = @[1'u32, 2'u32] + b.limbs = @[3'u32, 4'u32] + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[1'u32, 0'u32] + b.limbs = @[0'u32, 4'u32] + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[2'u32, 1'u32] + b.limbs = @[3'u32, 4'u32] + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[2'u32, 1'u32] + b.limbs = @[4'u32, 3'u32] + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + a.limbs = @[1'u32, 2'u32] b.limbs = @[3'u32, 4'u32] echo a.limbs @@ -1255,3 +1305,54 @@ when isMainModule: karatsubaMultiplication(c, a, b) echo "product Karatsuba: ", c echo "correct product: ", a * b + + a = two shl 32 - one + b = four shl 32 - three + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a = -(two shl 32 + one) + b = four shl 32 - three + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[1'u32, 2'u32, 3'u32] + b.limbs = @[4'u32, 5'u32, 6'u32] + a.isNegative = false + b.isNegative = false + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[1'u32, 2'u32, 3'u32, 4'u32, 5'u32] + b.limbs = @[4'u32, 5'u32, 6'u32, 7'u32, 8'u32] + a.isNegative = false + b.isNegative = false + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b + + a.limbs = @[1'u32, 2'u32, 3'u32, 4'u32, 5'u32, 6'u32, 7'u32, 8'u32, 9'u32, 10'u32] + b.limbs = @[10'u32, 9'u32, 8'u32, 7'u32, 6'u32, 5'u32, 4'u32, 3'u32, 2'u32, 1'u32] + a.isNegative = false + b.isNegative = false + echo a.limbs + echo b.limbs + echo "factors: ", a, " ", b + karatsubaMultiplication(c, a, b) + echo "product Karatsuba: ", c + echo "correct product: ", a * b diff --git a/tests/tbigints.nim b/tests/tbigints.nim index 1c808a4..6f7538c 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -315,12 +315,12 @@ proc main() = let d = "20867311096234429137120990056519061484140179793024844459539745043528236531589522382271230666075358518275274769618792229717222657110424037636116966396665200".initBigInt var r: BigInt = 0.initBigInt karatsubaMultiplication(r, c, d) - echo r - echo "\n" - echo c * d + # echo r + # echo "\n" + # echo c * d doAssert c * d == "322543179100245850295291096700090285623536165554432133161470913224665233565153206743023505404409261647296075477738317301701554637184306640109864382144645081119052516436652162825894456855767719709860985552674755702938369565636714472650667032224717209489767579823588160939485446085000195032327964706246225182000".initBigInt - echo "\n" - echo c*d - r + # echo "\n" + # echo c*d - r doAssert c*d == r # factors with 65 limbs @@ -813,5 +813,5 @@ proc main() = doAssert succ(a, 3) == initBigInt(10) -static: main() +# static: main() main() From 334537a36c33a8708867daf6eacc328e5abd48aa Mon Sep 17 00:00:00 2001 From: Dimitri Lesnoff Date: Fri, 11 Feb 2022 22:00:34 +0000 Subject: [PATCH 12/12] Remove echo's, convert proc into func again --- src/bigints.nim | 158 ++++++++++++++++++++++----------------------- tests/tbigints.nim | 12 ++-- 2 files changed, 81 insertions(+), 89 deletions(-) diff --git a/src/bigints.nim b/src/bigints.nim index 352f7af..820e296 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -7,14 +7,14 @@ type limbs: seq[uint32] isNegative: bool -proc normalize(a: var BigInt) = +func normalize(a: var BigInt) = for i in countdown(a.limbs.high, 0): if a.limbs[i] > 0'u32: a.limbs.setLen(i+1) return a.limbs.setLen(1) -proc initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = +func initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = ## Initializes a `BigInt` from a sequence of `uint32` values. runnableExamples: let a = @[10'u32, 2'u32].initBigInt @@ -23,7 +23,7 @@ proc initBigInt*(vals: sink seq[uint32], isNegative = false): BigInt = result.limbs = vals result.isNegative = isNegative -proc initBigInt*[T: int8|int16|int32](val: T): BigInt = +func initBigInt*[T: int8|int16|int32](val: T): BigInt = if val < 0: result.limbs = @[(not val).uint32 + 1] # manual 2's complement (to avoid overflow) result.isNegative = true @@ -31,10 +31,10 @@ proc initBigInt*[T: int8|int16|int32](val: T): BigInt = result.limbs = @[val.uint32] result.isNegative = false -proc initBigInt*[T: uint8|uint16|uint32](val: T): BigInt = +func initBigInt*[T: uint8|uint16|uint32](val: T): BigInt = result.limbs = @[val.uint32] -proc initBigInt*(val: int64): BigInt = +func initBigInt*(val: int64): BigInt = var a = val.uint64 if val < 0: a = not a + 1 # 2's complement @@ -44,7 +44,7 @@ proc initBigInt*(val: int64): BigInt = else: result.limbs = @[a.uint32] -proc initBigInt*(val: uint64): BigInt = +func initBigInt*(val: uint64): BigInt = if val > uint32.high: result.limbs = @[(val and uint32.high).uint32, (val shr 32).uint32] else: @@ -57,7 +57,7 @@ else: template initBigInt*(val: int): BigInt = initBigInt(val.int64) template initBigInt*(val: uint): BigInt = initBigInt(val.uint64) -proc initBigInt*(val: BigInt): BigInt = +func initBigInt*(val: BigInt): BigInt = result = val const @@ -65,13 +65,13 @@ const one = initBigInt(1) karatsubaThreshold = 2 -proc isZero(a: BigInt): bool {.inline.} = +func isZero(a: BigInt): bool {.inline.} = for i in countdown(a.limbs.high, 0): if a.limbs[i] != 0'u32: return false return true -proc abs*(a: BigInt): BigInt = +func abs*(a: BigInt): BigInt = # Returns the absolute value of `a`. runnableExamples: assert abs(42.initBigInt) == 42.initBigInt @@ -79,16 +79,16 @@ proc abs*(a: BigInt): BigInt = result = a result.isNegative = false -proc unsignedCmp(a: BigInt, b: uint32): int64 = +func unsignedCmp(a: BigInt, b: uint32): int64 = # ignores the sign of `a` # `a` and `b` are assumed to not be zero result = int64(a.limbs.len) - 1 if result != 0: return result = int64(a.limbs[0]) - int64(b) -proc unsignedCmp(a: uint32, b: BigInt): int64 = -unsignedCmp(b, a) +func unsignedCmp(a: uint32, b: BigInt): int64 = -unsignedCmp(b, a) -proc unsignedCmp(a, b: BigInt): int64 = +func unsignedCmp(a, b: BigInt): int64 = # ignores the signs of `a` and `b` # `a` and `b` are assumed to not be zero result = int64(a.limbs.len) - int64(b.limbs.len) @@ -98,7 +98,7 @@ proc unsignedCmp(a, b: BigInt): int64 = if result != 0: return -proc cmp(a, b: BigInt): int64 = +func cmp(a, b: BigInt): int64 = ## Returns: ## * a value less than zero, if `a < b` ## * a value greater than zero, if `a > b` @@ -121,7 +121,7 @@ proc cmp(a, b: BigInt): int64 = else: return unsignedCmp(a, b) -proc cmp(a: BigInt, b: int32): int64 = +func cmp(a: BigInt, b: int32): int64 = ## Returns: ## * a value less than zero, if `a < b` ## * a value greater than zero, if `a > b` @@ -139,9 +139,9 @@ proc cmp(a: BigInt, b: int32): int64 = else: return unsignedCmp(a, b.uint32) -proc cmp(a: int32, b: BigInt): int64 = -cmp(b, a) +func cmp(a: int32, b: BigInt): int64 = -cmp(b, a) -proc `==`*(a, b: BigInt): bool = +func `==`*(a, b: BigInt): bool = ## Compares if two `BigInt` numbers are equal. runnableExamples: let @@ -152,7 +152,7 @@ proc `==`*(a, b: BigInt): bool = assert b != c cmp(a, b) == 0 -proc `<`*(a, b: BigInt): bool = +func `<`*(a, b: BigInt): bool = runnableExamples: let a = 5.initBigInt @@ -162,7 +162,7 @@ proc `<`*(a, b: BigInt): bool = assert b > c cmp(a, b) < 0 -proc `<=`*(a, b: BigInt): bool = +func `<=`*(a, b: BigInt): bool = runnableExamples: let a = 5.initBigInt @@ -172,16 +172,16 @@ proc `<=`*(a, b: BigInt): bool = assert c <= b cmp(a, b) <= 0 -proc `==`(a: BigInt, b: int32): bool = cmp(a, b) == 0 -proc `<`(a: BigInt, b: int32): bool = cmp(a, b) < 0 -proc `<`(a: int32, b: BigInt): bool = cmp(a, b) < 0 +func `==`(a: BigInt, b: int32): bool = cmp(a, b) == 0 +func `<`(a: BigInt, b: int32): bool = cmp(a, b) < 0 +func `<`(a: int32, b: BigInt): bool = cmp(a, b) < 0 template addParts(toAdd) = tmp += toAdd a.limbs[i] = uint32(tmp and uint32.high) tmp = tmp shr 32 -proc unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = +func unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = let bl = b.limbs.len a.limbs.setLen(bl) @@ -192,7 +192,7 @@ proc unsignedAdditionInt(a: var BigInt, b: BigInt, c: uint32) = a.limbs.add(uint32(tmp)) a.isNegative = false -proc unsignedAddition(a: var BigInt, b, c: BigInt) = +func unsignedAddition(a: var BigInt, b, c: BigInt) = let bl = b.limbs.len cl = c.limbs.len @@ -212,10 +212,10 @@ proc unsignedAddition(a: var BigInt, b, c: BigInt) = a.limbs.add(uint32(tmp)) a.isNegative = false -proc negate(a: var BigInt) = +func negate(a: var BigInt) = a.isNegative = not a.isNegative -proc `-`*(a: BigInt): BigInt = +func `-`*(a: BigInt): BigInt = ## Unary minus for `BigInt`. runnableExamples: let @@ -270,7 +270,7 @@ template realUnsignedSubtraction(a: var BigInt, b, c: BigInt) = normalize(a) assert tmp == 0 -proc unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = +func unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = # `b` is not zero let cmpRes = unsignedCmp(b, c) if cmpRes > 0: @@ -282,7 +282,7 @@ proc unsignedSubtractionInt(a: var BigInt, b: BigInt, c: uint32) = else: # b == c a = zero -proc unsignedSubtraction(a: var BigInt, b, c: BigInt) = +func unsignedSubtraction(a: var BigInt, b, c: BigInt) = let cmpRes = unsignedCmp(b, c) if cmpRes > 0: realUnsignedSubtraction(a, b, c) @@ -292,7 +292,7 @@ proc unsignedSubtraction(a: var BigInt, b, c: BigInt) = else: # b == c a = zero -proc additionInt(a: var BigInt, b: BigInt, c: int32) = +func additionInt(a: var BigInt, b: BigInt, c: int32) = # a = b + c if b.isZero: a = c.initBigInt @@ -308,7 +308,7 @@ proc additionInt(a: var BigInt, b: BigInt, c: int32) = else: unsignedAdditionInt(a, b, c.uint32) -proc addition(a: var BigInt, b, c: BigInt) = +func addition(a: var BigInt, b, c: BigInt) = # a = b + c if b.isNegative: if c.isNegative: @@ -322,7 +322,7 @@ proc addition(a: var BigInt, b, c: BigInt) = else: unsignedAddition(a, b, c) -proc `+`*(a, b: BigInt): BigInt = +func `+`*(a, b: BigInt): BigInt = ## Addition for `BigInt`s. runnableExamples: let @@ -340,7 +340,7 @@ template `+=`*(a: var BigInt, b: BigInt) = assert a == 7.initBigInt a = a + b -proc subtractionInt(a: var BigInt, b: BigInt, c: int32) = +func subtractionInt(a: var BigInt, b: BigInt, c: int32) = # a = b - c if b.isZero: a = -c.initBigInt @@ -356,7 +356,7 @@ proc subtractionInt(a: var BigInt, b: BigInt, c: int32) = else: unsignedSubtractionInt(a, b, c.uint32) -proc subtraction(a: var BigInt, b, c: BigInt) = +func subtraction(a: var BigInt, b, c: BigInt) = # a = b - c if b.isNegative: if c.isNegative: @@ -370,7 +370,7 @@ proc subtraction(a: var BigInt, b, c: BigInt) = else: unsignedSubtraction(a, b, c) -proc `-`*(a, b: BigInt): BigInt = +func `-`*(a, b: BigInt): BigInt = ## Subtraction for `BigInt`s. runnableExamples: let @@ -388,7 +388,7 @@ template `-=`*(a: var BigInt, b: BigInt) = assert a == 3.initBigInt a = a - b -proc unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = +func unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = # always called with bl >= cl let bl = b.limbs.len @@ -417,7 +417,7 @@ proc unsignedMultiplication(a: var BigInt, b, c: BigInt) {.inline.} = inc pos normalize(a) -proc scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = +func scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = # always called with bl >= cl if c == 0: a = zero @@ -437,11 +437,11 @@ proc scalarMultiplication(a: var BigInt, b: BigInt, c: uint32) {.inline.} = normalize(a) # forward declaration for use in `multiplication` -proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) -proc `shl`*(x: BigInt, y: Natural): BigInt -proc `shr`*(x: BigInt, y: Natural): BigInt +func karatsubaMultiplication*(a: var BigInt, b, c: BigInt) +func `shl`*(x: BigInt, y: Natural): BigInt +func `shr`*(x: BigInt, y: Natural): BigInt -proc multiplication(a: var BigInt, b, c: BigInt) = +func multiplication(a: var BigInt, b, c: BigInt) = # a = b * c if b.isZero or c.isZero: a = zero @@ -462,7 +462,7 @@ proc multiplication(a: var BigInt, b, c: BigInt) = unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative -proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = +func karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = if b.isZero or c.isZero: a = zero return @@ -495,7 +495,6 @@ proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = unsignedMultiplication(a, b, c) a.isNegative = b.isNegative xor c.isNegative return - # echo k var low_b, high_b, low_c, high_c: BigInt # Decompose `b` and `c` in two parts of (almost) equal length @@ -503,7 +502,6 @@ proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = high_b.limbs = b.limbs[k .. ^1] low_c.limbs = c.limbs[0 .. k-1] high_c.limbs = c.limbs[k .. ^1] - # echo low_b, high_b, low_c, high_c # subtractive version of Karatsuba's algorithm to limit carry handling var lowProduct, highProduct, add3, add4, add5, middleTerm: BigInt = zero @@ -511,18 +509,16 @@ proc karatsubaMultiplication*(a: var BigInt, b, c: BigInt) = multiplication(lowProduct, low_b, low_c) multiplication(highProduct, high_b, high_c) - # echo "lowProduct, highProduct: ", lowProduct, highProduct add3 = low_b - high_b add4 = high_c - low_c multiplication(add5, add4, add3) middleTerm = lowProduct + highProduct + add5 - # echo "add5 = add4*add3, middleTerm = lP + hP + add4*add3: ", add5, middleTerm a = lowProduct + middleTerm shl (32*k) + highProduct shl (64*k) a.isNegative = b.isNegative xor c.isNegative -proc `*`*(a, b: BigInt): BigInt = +func `*`*(a, b: BigInt): BigInt = ## Multiplication for `BigInt`s. runnableExamples: let @@ -538,7 +534,7 @@ template `*=`*(a: var BigInt, b: BigInt) = assert a == 150.initBigInt a = a * b -proc pow*(x: BigInt, y: Natural): BigInt = +func pow*(x: BigInt, y: Natural): BigInt = ## Computes `x` to the power of `y`. var base = x var exp = y @@ -551,7 +547,7 @@ proc pow*(x: BigInt, y: Natural): BigInt = exp = exp shr 1 base *= base -proc `shl`*(x: BigInt, y: Natural): BigInt = +func `shl`*(x: BigInt, y: Natural): BigInt = ## Shifts a `BigInt` to the left. runnableExamples: let a = 24.initBigInt @@ -574,9 +570,9 @@ proc `shl`*(x: BigInt, y: Natural): BigInt = result.limbs.add(uint32(carry shr (32 - b))) # forward declaration for use in `shr` -proc dec*(a: var BigInt, b: int = 1) +func dec*(a: var BigInt, b: int = 1) -proc `shr`*(x: BigInt, y: Natural): BigInt = +func `shr`*(x: BigInt, y: Natural): BigInt = ## Shifts a `BigInt` to the right (arithmetically). runnableExamples: let a = 24.initBigInt @@ -612,17 +608,17 @@ proc `shr`*(x: BigInt, y: Natural): BigInt = # normalize result.limbs.setLen(result.limbs.high) -proc bitwiseAnd(a: var BigInt, b, c: BigInt) = +func bitwiseAnd(a: var BigInt, b, c: BigInt) = a.limbs.setLen(min(b.limbs.len, c.limbs.len)) for i in 0 ..< a.limbs.len: a.limbs[i] = b.limbs[i] and c.limbs[i] -proc `and`*(a, b: BigInt): BigInt = +func `and`*(a, b: BigInt): BigInt = ## Bitwise `and` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) bitwiseAnd(result, a, b) -proc bitwiseOr(a: var BigInt, b, c: BigInt) = +func bitwiseOr(a: var BigInt, b, c: BigInt) = # `b` must be smaller than `c` a.limbs.setLen(c.limbs.len) for i in 0 ..< b.limbs.len: @@ -630,7 +626,7 @@ proc bitwiseOr(a: var BigInt, b, c: BigInt) = for i in b.limbs.len ..< c.limbs.len: a.limbs[i] = c.limbs[i] -proc `or`*(a, b: BigInt): BigInt = +func `or`*(a, b: BigInt): BigInt = ## Bitwise `or` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) if a.limbs.len <= b.limbs.len: @@ -638,7 +634,7 @@ proc `or`*(a, b: BigInt): BigInt = else: bitwiseOr(result, b, a) -proc bitwiseXor(a: var BigInt, b, c: BigInt) = +func bitwiseXor(a: var BigInt, b, c: BigInt) = # `b` must be smaller than `c` a.limbs.setLen(c.limbs.len) for i in 0 ..< b.limbs.len: @@ -646,7 +642,7 @@ proc bitwiseXor(a: var BigInt, b, c: BigInt) = for i in b.limbs.len ..< c.limbs.len: a.limbs[i] = c.limbs[i] -proc `xor`*(a, b: BigInt): BigInt = +func `xor`*(a, b: BigInt): BigInt = ## Bitwise `xor` for `BigInt`s. assert (not a.isNegative) and (not b.isNegative) if a.limbs.len <= b.limbs.len: @@ -654,13 +650,13 @@ proc `xor`*(a, b: BigInt): BigInt = else: bitwiseXor(result, b, a) -proc reset(a: var BigInt) = +func reset(a: var BigInt) = ## Resets a `BigInt` back to the zero value. a.limbs.setLen(1) a.limbs[0] = 0 a.isNegative = false -proc unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = +func unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = q.limbs.setLen(n.limbs.len) r = 0 for i in countdown(n.limbs.high, 0): @@ -669,7 +665,7 @@ proc unsignedDivRem(q: var BigInt, r: var uint32, n: BigInt, d: uint32) = r = uint32(tmp mod d) normalize(q) -proc bits(d: uint32): int = +func bits(d: uint32): int = const bitLengths = [0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5] var d = d @@ -679,7 +675,7 @@ proc bits(d: uint32): int = result += bitLengths[int(d)] # From Knuth and Python -proc unsignedDivRem(q, r: var BigInt, n, d: BigInt) = +func unsignedDivRem(q, r: var BigInt, n, d: BigInt) = var nn = n.limbs.len dn = d.limbs.len @@ -789,7 +785,7 @@ proc unsignedDivRem(q, r: var BigInt, n, d: BigInt) = q = a normalize(q) -proc division(q, r: var BigInt, n, d: BigInt) = +func division(q, r: var BigInt, n, d: BigInt) = # q = n div d # r = n mod d if d.isZero: @@ -805,11 +801,11 @@ proc division(q, r: var BigInt, n, d: BigInt) = r += d q -= one -proc `div`*(a, b: BigInt): BigInt = +func `div`*(a, b: BigInt): BigInt = ## Computes the integer division of two `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. ## - ## If you also need the modulo (remainder), use the `divmod proc <#divmod,BigInt,BigInt>`_. + ## If you also need the modulo (remainder), use the `divmod func <#divmod,BigInt,BigInt>`_. runnableExamples: let a = 17.initBigInt @@ -821,11 +817,11 @@ proc `div`*(a, b: BigInt): BigInt = var tmp: BigInt division(result, tmp, a, b) -proc `mod`*(a, b: BigInt): BigInt = +func `mod`*(a, b: BigInt): BigInt = ## Computes the integer modulo (remainder) of two `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. ## - ## If you also need an integer division, use the `divmod proc <#divmod,BigInt,BigInt>`_. + ## If you also need an integer division, use the `divmod func <#divmod,BigInt,BigInt>`_. runnableExamples: let a = 17.initBigInt @@ -837,7 +833,7 @@ proc `mod`*(a, b: BigInt): BigInt = var tmp: BigInt division(tmp, result, a, b) -proc divmod*(a, b: BigInt): tuple[q, r: BigInt] = +func divmod*(a, b: BigInt): tuple[q, r: BigInt] = ## Computes both the integer division and modulo (remainder) of two ## `BigInt` numbers. ## Raises a `DivByZeroDefect` if `b` is zero. @@ -848,7 +844,7 @@ proc divmod*(a, b: BigInt): tuple[q, r: BigInt] = assert divmod(a, b) == (3.initBigInt, 2.initBigInt) division(result.q, result.r, a, b) -proc countTrailingZeroBits(a: BigInt): int = +func countTrailingZeroBits(a: BigInt): int = var count = 0 for x in a.limbs: if x == 0: @@ -857,7 +853,7 @@ proc countTrailingZeroBits(a: BigInt): int = return count + countTrailingZeroBits(x) return count -proc gcd*(a, b: BigInt): BigInt = +func gcd*(a, b: BigInt): BigInt = ## Returns the greatest common divisor (GCD) of `a` and `b`. runnableExamples: assert gcd(54.initBigInt, 24.initBigInt) == 6.initBigInt @@ -886,7 +882,7 @@ proc gcd*(a, b: BigInt): BigInt = v = v shr countTrailingZeroBits(v) -proc toInt*[T: SomeInteger](x: BigInt): Option[T] = +func toInt*[T: SomeInteger](x: BigInt): Option[T] = ## Converts a `BigInt` number to an integer, if possible. ## If the `BigInt` doesn't fit in a `T`, returns `none(T)`; ## otherwise returns `some(x)`. @@ -959,7 +955,7 @@ proc toInt*[T: SomeInteger](x: BigInt): Option[T] = else: result = some(T(x.limbs[0])) -proc calcSizes(): array[2..36, int] = +func calcSizes(): array[2..36, int] = for i in 2..36: var x = int64(i) while x <= int64(uint32.high) + 1: @@ -971,7 +967,7 @@ const powers = {2, 4, 8, 16, 32} sizes = calcSizes() # `sizes[base]` is the maximum number of digits that fully fit in a `uint32` -proc toString*(a: BigInt, base: range[2..36] = 10): string = +func toString*(a: BigInt, base: range[2..36] = 10): string = ## Produces a string representation of a `BigInt` in a specified ## `base`. ## @@ -1036,11 +1032,11 @@ proc toString*(a: BigInt, base: range[2..36] = 10): string = result.reverse() -proc `$`*(a: BigInt): string = +func `$`*(a: BigInt): string = ## String representation of a `BigInt` in base 10. toString(a, 10) -proc parseDigit(c: char, base: uint32): uint32 {.inline.} = +func parseDigit(c: char, base: uint32): uint32 {.inline.} = result = case c of '0'..'9': uint32(ord(c) - ord('0')) of 'a'..'z': uint32(ord(c) - ord('a') + 10) @@ -1050,7 +1046,7 @@ proc parseDigit(c: char, base: uint32): uint32 {.inline.} = if result >= base: raise newException(ValueError, "Invalid input: " & c) -proc initBigInt*(str: string, base: range[2..36] = 10): BigInt = +func initBigInt*(str: string, base: range[2..36] = 10): BigInt = ## Create a `BigInt` from a string. For invalid inputs, a `ValueError` exception is raised. runnableExamples: let @@ -1112,7 +1108,7 @@ proc initBigInt*(str: string, base: range[2..36] = 10): BigInt = when (NimMajor, NimMinor) >= (1, 5): include bigints/private/literals -proc inc*(a: var BigInt, b: int = 1) = +func inc*(a: var BigInt, b: int = 1) = ## Increase the value of a `BigInt` by the specified amount (default: 1). runnableExamples: var a = 15.initBigInt @@ -1127,7 +1123,7 @@ proc inc*(a: var BigInt, b: int = 1) = else: a += initBigInt(b) -proc dec*(a: var BigInt, b: int = 1) = +func dec*(a: var BigInt, b: int = 1) = ## Decrease the value of a `BigInt` by the specified amount (default: 1). runnableExamples: var a = 15.initBigInt @@ -1142,12 +1138,12 @@ proc dec*(a: var BigInt, b: int = 1) = else: a -= initBigInt(b) -proc succ*(a: BigInt, b: int = 1): BigInt = +func succ*(a: BigInt, b: int = 1): BigInt = ## Returns the `b`-th successor of a `BigInt`. result = a inc(result, b) -proc pred*(a: BigInt, b: int = 1): BigInt = +func pred*(a: BigInt, b: int = 1): BigInt = ## Returns the `b`-th predecessor of a `BigInt`. result = a dec(result, b) @@ -1181,14 +1177,14 @@ iterator `..<`*(a, b: BigInt): BigInt = yield res inc res -proc modulo(a, modulus: BigInt): BigInt = +func modulo(a, modulus: BigInt): BigInt = ## Like `mod`, but the result is always in the range `[0, modulus-1]`. ## `modulus` should be greater than zero. result = a mod modulus if result < 0: result += modulus -proc fastLog2*(a: BigInt): int = +func fastLog2*(a: BigInt): int = ## Computes the logarithm in base 2 of `a`. ## If `a` is negative, returns the logarithm of `abs(a)`. ## If `a` is zero, returns -1. @@ -1196,7 +1192,7 @@ proc fastLog2*(a: BigInt): int = return -1 bitops.fastLog2(a.limbs[^1]) + 32*(a.limbs.high) -proc invmod*(a, modulus: BigInt): BigInt = +func invmod*(a, modulus: BigInt): BigInt = ## Compute the modular inverse of `a` modulo `modulus`. ## The return value is always in the range `[1, modulus-1]` runnableExamples: @@ -1228,7 +1224,7 @@ proc invmod*(a, modulus: BigInt): BigInt = raise newException(ValueError, $a & " has no modular inverse modulo " & $modulus) result = t0.modulo(modulus) -proc powmod*(base, exponent, modulus: BigInt): BigInt = +func powmod*(base, exponent, modulus: BigInt): BigInt = ## Compute modular exponentation of `base` with power `exponent` modulo `modulus`. ## The return value is always in the range `[0, modulus-1]`. runnableExamples: diff --git a/tests/tbigints.nim b/tests/tbigints.nim index 6f7538c..df19abd 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -301,13 +301,14 @@ proc main() = doAssert (d xor f) == d block: # multiplication + # test sign let one = 1.initBigInt let negOne = -1.initBigInt - echo one * negOne + doAssert one * negOne == negOne + # factors with 4 limbs let a = "1780983279228119273110576463639172624".initBigInt let b = "1843917749452418885995463656480858321".initBigInt - echo a * b doAssert a * b == "3283986680046702618742503890385314117448805445290098330749803441805804304".initBigInt # factors with 17 limbs @@ -315,13 +316,8 @@ proc main() = let d = "20867311096234429137120990056519061484140179793024844459539745043528236531589522382271230666075358518275274769618792229717222657110424037636116966396665200".initBigInt var r: BigInt = 0.initBigInt karatsubaMultiplication(r, c, d) - # echo r - # echo "\n" - # echo c * d doAssert c * d == "322543179100245850295291096700090285623536165554432133161470913224665233565153206743023505404409261647296075477738317301701554637184306640109864382144645081119052516436652162825894456855767719709860985552674755702938369565636714472650667032224717209489767579823588160939485446085000195032327964706246225182000".initBigInt - # echo "\n" - # echo c*d - r - doAssert c*d == r + doAssert c * d == r # factors with 65 limbs let e = initBigInt("1f3b839241b0aacc183858dc7a75a773e7bad642a9f426ef499d91e09c9f99a88ec9a14d5ee51175faeaa10d2fb06f3ee37d2f50fe755c2c963aeb539cd55c0e14f5a23f04c64839c22bd4108034b7afc95e01a1c2fe605d8b1930926e886a8f3d7fc09acd54d388cb5d4b3a3fb4eaf6781173ab3a0cd8ad3119c37dd2cf05544235d7b85b2c96d2ed29e1a685820c4afdd824bd8878f1b6a3f52a57eb886efaa737af47161c89f298d908aa950979b8c2615d4e03b47ee87a5381ca39d9ec4788d7abd07b174913b962c02cdd5f8319722a3345eb38d3ebdd51dec66a58e89902151539298c41446758bac66923c910fd7a2d12d0d5c8bb688970b8a77e7d5fc", base = 16)