diff --git a/benchmarks/benchBigints.nim b/benchmarks/benchBigints.nim new file mode 100644 index 0000000..f5b82da --- /dev/null +++ b/benchmarks/benchBigints.nim @@ -0,0 +1,61 @@ +import bigints +import benchy +import random + +block: # Bench random generation + randomize() + var n = 100 + timeIt "bench Random generation of bigints with 100 limbs": + keep initRandomBigint(n) + timeIt "bench Random generation of bigints with 1_000 limbs": + keep initRandomBigint(n) + n = 5_000 + timeIt "bench Random generation of bigints with 5_000 limbs": + keep initRandomBigint(n) + n = 10_000 + timeIt "bench Random generation of bigints with 10_000 limbs": + keep initRandomBigint(n) +block: # Bench multiplication + randomize() + var n = 100 + timeIt "bench Multiplication of bigints with 100 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a*b + n = 1_000 + timeIt "bench Multiplication of bigints with 1_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a*b + n = 5_000 + timeIt "bench Multiplication of bigints with 5_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a*b + n = 10_000 + timeIt "bench Multiplication of bigints with 10_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a*b +block: # Bench division + randomize() + var n = 100 + timeIt "bench Division of bigints with 100 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a div b + n = 1_000 + timeIt "bench Division of bigints with 1_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a div b + n = 5_000 + timeIt "bench Division of bigints with 5_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a div b + n = 10_000 + timeIt "bench Division of bigints with 10_000 limbs": + var a: Bigint = initRandomBigint(n) + var b: Bigint = initRandomBigint(n) + keep a div b diff --git a/benchmarks/benchExamples.nim b/benchmarks/benchExamples.nim new file mode 100644 index 0000000..dcff669 --- /dev/null +++ b/benchmarks/benchExamples.nim @@ -0,0 +1,35 @@ +import bigints +import benchy +from std/math import `^` +import pidigits +import rc_combperm + +block: # Binomial and permutations + timeIt "Permutation 1000, 969": + keep perm(1000, 969) + + timeIt "Binomial computation 1000, 969": + keep comb(1000, 969) + +block: # Power computation + timeIt "Power computation of 5^4^3^2": + keep 5.initBigInt.pow 4 ^ (3 ^ 2) + timeIt "Powers of two": + var power = 2.initBigInt + for _ in 1 .. 128000: + power = power * 2.initBigInt + +block: # Pidigits example + timeIt "Computation of 100 digits of Pi": + var i = 0 + while i < 100: + var d: int32 = findPiDigit() + inc(i) + eliminateDigit(d) + timeIt "Computation of 1000 digits of Pi": + var i = 0 + while i < 1000: + var d: int32 = findPiDigit() + inc(i) + eliminateDigit(d) + diff --git a/benchmarks/nim.cfg b/benchmarks/nim.cfg new file mode 100644 index 0000000..6feaed8 --- /dev/null +++ b/benchmarks/nim.cfg @@ -0,0 +1,2 @@ +path = "$projectPath/../src" +path = "$projectPath/../examples" diff --git a/bigints.nimble b/bigints.nimble index 17b8d86..d04d3a6 100644 --- a/bigints.nimble +++ b/bigints.nimble @@ -25,3 +25,20 @@ task checkExamples, "Check examples": for example in listFiles("examples"): if example.endsWith(".nim"): exec "nim check --hints:off " & example + +task benchAll, "Benchmark Library": + for backend in ["c", "cpp"]: + echo "Benchmarks with " & backend & " backend" + for gc in ["refc", "arc", "orc"]: + echo "Benchmark with " & gc & " garbage collector" + for benchmark in listFiles("benchmarks"): + if benchmark.endsWith(".nim"): + echo "Benchmark " & benchmark + exec "nim r --hints:off -d:danger --opt:speed -d:lto --backend:" & backend & " --gc:" & gc & " " & benchmark + +task benchOrcC, "Benchmark Library with orc garbage collector and C backend": + for benchmark in listFiles("benchmarks"): + if benchmark.endsWith(".nim"): + echo "Benchmark " & benchmark + exec "nim r --hints:off -d:danger --opt:speed -d:lto --backend:c --gc:orc " & benchmark + diff --git a/examples/elliptic.nim b/examples/elliptic.nim index 90e25a7..93380c5 100644 --- a/examples/elliptic.nim +++ b/examples/elliptic.nim @@ -35,7 +35,7 @@ proc ecDouble(a: tuple): (BigInt, BigInt) = lam = lam mod primeCurve result = (x, y) -proc ecMultiply(genPoint: tuple, scalarHex: BigInt): (BigInt, BigInt) = +proc ecMultiply*(genPoint: tuple, scalarHex: BigInt): (BigInt, BigInt) = if scalarHex == zero or scalarHex >= numberPoints: raise newException(Exception, "Invalid Scalar/Private Key") var @@ -47,7 +47,7 @@ proc ecMultiply(genPoint: tuple, scalarHex: BigInt): (BigInt, BigInt) = q = ecAdd(q, genPoint) result = q -proc main() = +when isMainModule: let publicKey = ecMultiply(Gpoint, privKey) echo "" @@ -65,5 +65,3 @@ proc main() = echo "the official Public Key - compressed:" echo if publicKey[1] mod two == one: "03" & publicKey[0].toString(base = 16).align(64, '0') else: "02" & publicKey[0].toString(base = 16).align(64, '0') - -main() diff --git a/examples/pidigits.nim b/examples/pidigits.nim index aa8448d..23ec4a5 100644 --- a/examples/pidigits.nim +++ b/examples/pidigits.nim @@ -34,7 +34,7 @@ proc extractDigit(): int32 = result = get(toInt[int32](tmp1 and mask)) -proc eliminateDigit(d: int32) = +proc eliminateDigit*(d: int32) = acc -= den * d.initBigInt acc *= ten num *= ten @@ -48,33 +48,34 @@ proc nextTerm() = den *= k2 num *= k -proc findPiDigit(): int32 = +proc findPiDigit*(): int32 = result = -1 while result < 0: nextTerm() result = extractDigit() -var i = 0 -if paramCount() == 0: - # prints an infinite amount of pi digits - while true: - var d: int32 = findPiDigit() - stdout.write chr(ord('0') + d) - inc i - if i == 40: - echo "" - i = 0 - eliminateDigit(d) +when isMainModule: + var i = 0 + if paramCount() == 0: + # prints an infinite amount of pi digits + while true: + var d: int32 = findPiDigit() + stdout.write chr(ord('0') + d) + inc i + if i == 40: + echo "" + i = 0 + eliminateDigit(d) -let n = parseInt(paramStr(1)) + let n = parseInt(paramStr(1)) -if n <= 0: - quit("The number you entered is negative. Please specify a strictly positive number") + if n <= 0: + quit("The number you entered is negative. Please specify a strictly positive number") -while i < n: - var d: int32 = findPiDigit() - stdout.write(chr(ord('0') + d)) - inc(i) - if i mod 40 == 0: - echo "\t:", i - eliminateDigit(d) + while i < n: + var d: int32 = findPiDigit() + stdout.write(chr(ord('0') + d)) + inc(i) + if i mod 40 == 0: + echo "\t:", i + eliminateDigit(d) diff --git a/examples/rc_combperm.nim b/examples/rc_combperm.nim index dc25dab..9549110 100644 --- a/examples/rc_combperm.nim +++ b/examples/rc_combperm.nim @@ -1,7 +1,7 @@ # Solution for https://rosettacode.org/wiki/Combinations_and_permutations import bigints -proc perm(n, k: int32): BigInt = +proc perm*(n, k: int32): BigInt = result = initBigInt(1) var k = initBigInt(n - k) @@ -10,12 +10,13 @@ proc perm(n, k: int32): BigInt = result *= n dec n -proc comb(n, k: int32): BigInt = +proc comb*(n, k: int32): BigInt = result = perm(n, k) var k = initBigInt(k) while k > 0.initBigInt: result = result div k dec k -echo "P(1000, 969) = ", perm(1000, 969) -echo "C(1000, 969) = ", comb(1000, 969) +when isMainModule: + echo "P(1000, 969) = ", perm(1000, 969) + echo "C(1000, 969) = ", comb(1000, 969) diff --git a/src/bigints.nim b/src/bigints.nim index a50e9b0..68ebc85 100644 --- a/src/bigints.nim +++ b/src/bigints.nim @@ -2,7 +2,7 @@ ## ## The bitwise operations behave as if negative numbers were represented in 2's complement. -import std/[algorithm, bitops, math, options] +import std/[algorithm, bitops, math, options, random] type BigInt* = object @@ -74,6 +74,50 @@ else: func initBigInt*(val: BigInt): BigInt = result = val +type + SizeDescriptor* = enum + Limbs, Bits + +proc initRandomBigInt*(number: Natural, unit: SizeDescriptor = Limbs): BigInt = + ## Initializes a standalone BigInt whose value is chosen randomly with exactly + ## `number` bits or limbs, depending on the value of `unit`. By default, the + ## BigInt is chosen with `number` limbs chosen randomly. + if unit == Limbs: + if number == 0: + raise newException(ValueError, "A Bigint must have at least one limb !") + result.limbs.setLen(number) + for i in 0 ..< result.limbs.len-1: + result.limbs[i] = rand(uint32) + var word = rand(uint32) + # Bigint's last limb can be zero, iff there is only one limb + # We can't normalize instead, since we need no less than number limbs + if number != 1: + while word == 0: # Very low probability + word = rand(uint32) + result.limbs[result.limbs.len-1] = word + + else: # unit == Bits + if number == 0: + return initBigInt(0) + let + remainder = number mod 32 + n_limbs = (if remainder == 0: number shr 5 else: number shr 5 + 1) + remainingBits = (if remainder == 0: 32 else: remainder) + result.limbs.setLen(n_limbs) + # mask ensures only remainingBits bits can be set to 1 + # Ensures the first bit is set to 1 + var + mask: uint32 = 0xFFFF_FFFF'u32 + mask2: uint32 = 0x8000_0000'u32 + if remainingBits != 32: + mask = 1'u32 shl remainingBits - 1 + mask2 = 1'u32 shl (remainingBits-1) + for i in 0 ..< result.limbs.len-1: + result.limbs[i] = rand(uint32) + let word = rand(uint32) + result.limbs[result.limbs.len-1] = word and mask or mask2 + + const zero = initBigInt(0) one = initBigInt(1) @@ -1304,3 +1348,4 @@ func powmod*(base, exponent, modulus: BigInt): BigInt = result = (result * basePow) mod modulus basePow = (basePow * basePow) mod modulus exponent = exponent shr 1 + diff --git a/tests/tbigints.nim b/tests/tbigints.nim index 379bc36..94cd7a5 100644 --- a/tests/tbigints.nim +++ b/tests/tbigints.nim @@ -1,5 +1,6 @@ import bigints import std/options +import random const zero = initBigInt(0) @@ -880,6 +881,11 @@ proc main() = doAssert pred(a, 3) == initBigInt(4) doAssert succ(a, 3) == initBigInt(10) + # block: + # randomize() + # let a: BigInt = initRandomBigInt(33) + # echo a -static: main() + +#static: main() main() diff --git a/tests/trandom.nim b/tests/trandom.nim index c7f2b7e..25a3f24 100644 --- a/tests/trandom.nim +++ b/tests/trandom.nim @@ -1,22 +1,87 @@ import bigints +import random import bigints/random import std/options -block: # check uniformity - let lo = pow(10.initBigInt, 90) - let hi = pow(10.initBigInt, 100) - var total = 0.initBigInt - let trials = 1000 - let nbuckets = 33 - var buckets = newSeq[int](nbuckets) - for x in 0..trials: - let r = rand(lo..hi) - doAssert(lo <= r) - doAssert(r <= hi) - total += r - let iBucket = (r-lo) div ((hi-lo) div initBigInt(nbuckets)) - buckets[iBucket.toInt[:int]().get()] += 1 - for x in buckets: - doAssert(trials/nbuckets*0.5 < float(x)) - doAssert(float(x) < trials/nbuckets*1.5) +type + MemSizeUnit = enum + o, Kio, Mio, Gio +const + zero = initBigInt(0) + one = initBigInt(1) + memSize = 2 # Max number of allocated memory for the tests + memSizeUnit = Mio # Unit in which memSize is expressed + +proc computeLimit(memSize: Natural, memSizeUnit: MemSizeUnit): Natural = + var factor = 1 + for _ in 1..ord(memSizeUnit): + factor *= 1024 + result = memSize * factor + +const + memLimit = computeLimit(memSize, memSizeUnit) # Number of octets + maxLimbs = memLimit div 8 + maxBits = 4*memLimit + +proc main() = + randomize() + + block: + let a: BigInt = initRandomBigInt(0, Bits) + doAssert a == zero + let b: BigInt = initRandomBigInt(1, Bits) + doAssert b == one + + block: + for nBits in [29, 32, 1037]: + for _ in 1 .. 5: # Repeat probabilistic tests + let a: BigInt = initRandomBigInt(nBits, Bits) + doAssert fastLog2(a) == (nBits - 1) + doAssert (toString(a, 2)).len == nBits + # For bigger bigints, remove the test with slow conversion to string + for nBits in [rand(1..maxBits), 32*rand(1..maxLimbs)]: + for _ in 1 .. 5: + let a: BigInt = initRandomBigInt(nBits, Bits) + doAssert fastLog2(a) == (nBits - 1) + + block: + for nLimbs in [1, 2, 3, 5, 10, 25, 100]: + for _ in 1 .. 5: + let a: BigInt = initRandomBigInt(nLimbs) + let n_bitsA = fastLog2(a) + 1 + doAssert n_bitsA <= 32*nlimbs + doAssert n_bitsA > 32*(nlimbs-1) + + block: # GCD properties but tested on random Bigints + let limitGCD = 100_000 # Special limit for the GCD, otherwise the tests run for hours + let (nBitsA, nBitsB, nBitsC) = (rand(1..limitGCD), rand(1..limitGCD), rand(1..limitGCD)) + let a = initRandomBigInt(nBitsA, Bits) + let b = initRandomBigInt(nBitsB, Bits) + let c = initRandomBigInt(nBitsC, Bits) + doAssert gcd(a, b) == gcd(b, a) + doAssert gcd(a, zero) == a + doAssert gcd(a, a) == a + doAssert gcd(c * a, c * b) == c * gcd(a,b) + doAssert gcd(a, gcd(b, c)) == gcd(gcd(a, b), c) + doAssert gcd(a, b) == gcd(b, a mod b) + + block: # check uniformity + let lo = pow(10.initBigInt, 90) + let hi = pow(10.initBigInt, 100) + var total = 0.initBigInt + let trials = 1000 + let nbuckets = 33 + var buckets = newSeq[int](nbuckets) + for x in 0..trials: + let r = rand(lo..hi) + doAssert(lo <= r) + doAssert(r <= hi) + total += r + let iBucket = (r-lo) div ((hi-lo) div initBigInt(nbuckets)) + buckets[iBucket.toInt[:int]().get()] += 1 + for x in buckets: + doAssert(trials/nbuckets*0.5 < float(x)) + doAssert(float(x) < trials/nbuckets*1.5) + +main()