diff --git a/bigints.nimble b/bigints.nimble
index 5de9609..b2659c7 100644
--- a/bigints.nimble
+++ b/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"
 
diff --git a/src/bigints.nim b/src/bigints.nim
index 9ef5abb..68f2e7c 100644
--- a/src/bigints.nim
+++ b/src/bigints.nim
@@ -66,6 +66,7 @@ else:
 func initBigInt*(val: BigInt): BigInt =
   result = val
 
+
 const
   zero = initBigInt(0)
   one = initBigInt(1)
@@ -1198,3 +1199,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/src/bigints/utilities.nim b/src/bigints/utilities.nim
new file mode 100644
index 0000000..0b9dde0
--- /dev/null
+++ b/src/bigints/utilities.nim
@@ -0,0 +1,67 @@
+import std/random
+import ../bigints
+
+const zero = initBigInt(0)
+type
+  RandomMode* = enum
+    Limbs, Bits
+
+proc randomizeBigInt(container: var seq[uint32], number: Natural, mode: RandomMode = Limbs) =
+  case mode
+  of Limbs:
+    if number == 0:
+      raise newException(ValueError, "A Bigint must have at least one limb !")
+    # result.limbs.setLen(number)
+    for i in 0 ..< number-1:
+      container[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)
+    container[number-1] = word
+
+  of Bits: # unit == Bits
+    if number == 0:
+      container = @[]
+    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
+    # mask2 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 ..< container.high:
+      container[i] = rand(uint32)
+    let word = rand(uint32)
+    container[container.high] = word and mask or mask2
+
+proc initRandomBigInt*(number: Natural, mode: RandomMode = Limbs): BigInt =
+  ## Initializes a `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.
+  ## Generates only positive bigints.
+  var limbs: seq[uint32]
+  let
+    remainder = number mod 32
+    n_limbs = (if remainder == 0: number shr 5 else: number shr 5 + 1)
+  case mode
+  of Limbs:
+    limbs.setLen(number)
+  of Bits:
+    if number == 0:
+      return zero
+    let
+      remainder = number mod 32
+      len_limbs = (if remainder == 0: number shr 5 else: number shr 5 + 1)
+    limbs.setLen(len_limbs)
+  randomizeBigInt(limbs, number, mode)
+  result = initBigInt(limbs, false)
+
diff --git a/tests/trandom.nim b/tests/trandom.nim
new file mode 100644
index 0000000..23f8c19
--- /dev/null
+++ b/tests/trandom.nim
@@ -0,0 +1,67 @@
+import bigints
+import ../src/bigints/utilities
+import std/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 =
+  result = memSize
+  for _ in 1..ord(memSizeUnit):
+    result *= 1024
+
+const
+  memLimit = computeLimit(memSize, memSizeUnit) # Number of bytes
+  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()