Add and test initRandomBigInt for issue #101

bigints.nimble

@@ -16,7 +16,7 @@ task test, "Test bigints":
     echo "testing " & backend & " backend"
     for gc in ["refc", "arc", "orc"]:
       echo "  using " & gc & " GC"
-      for file in ["tbigints.nim", "tbugs.nim"]:
+      for file in ["trandom.nim", "tbigints.nim", "tbugs.nim"]:
         exec "nim r --hints:off --experimental:strictFuncs --backend:" & backend & " --gc:" & gc & " tests/" & file
       exec "nim doc --hints:off --backend:" & backend & " --gc:" & gc & " src/bigints.nim"
 

src/bigints.nim

@@ -1,6 +1,6 @@
 ## Arbitrary precision integers.
 
-import std/[algorithm, bitops, math, options]
+import std/[algorithm, bitops, math, options, random]
 
 type
   BigInt* = object
@@ -66,6 +66,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)
@@ -1198,3 +1242,4 @@ func powmod*(base, exponent, modulus: BigInt): BigInt =
         result = (result * basePow) mod modulus
       basePow = (basePow * basePow) mod modulus
       exponent = exponent shr 1
+

tests/trandom.nim

@@ -0,0 +1,67 @@
+import bigints
+import random
+
+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)
+
+main()