diff --git a/bigints.nimble b/bigints.nimble index 5de9609..17b8d86 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 ["tbigints.nim", "tbugs.nim", "trandom.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/random.nim b/src/bigints/random.nim new file mode 100644 index 0000000..e4c7058 --- /dev/null +++ b/src/bigints/random.nim @@ -0,0 +1,43 @@ +import ../bigints +import std/sequtils +import std/options +import std/random + +func rand*(r: var Rand, x: Slice[BigInt]): BigInt = + ## Return a random `BigInt`, within the given range, using the given state. + assert(x.a <= x.b, "invalid range") + let + spread = x.b - x.a + # number of bits *not* including leading bit + nbits = spread.fastLog2 + # number of limbs to generate completely randomly + nFullLimbs = max(nbits div 32 - 1, 0) + # highest possible value of the top two limbs. + hi64Max = (spread shr (nFullLimbs*32)).toInt[:uint64].get() + while true: + # these limbs can be generated completely arbitrarily + var limbs = newSeqWith(nFullLimbs, r.rand(uint32.low..uint32.high)) + # generate the top two limbs more carefully. This all but guarantees + # that the entire number is in the correct range + let hi64 = r.rand(uint64.low..hi64Max) + limbs.add(cast[uint32](hi64)) + limbs.add(cast[uint32](hi64 shr 32)) + result = initBigInt(limbs) + if result <= spread: + break + result += x.a + +func rand*(r: var Rand, max: BigInt): BigInt = + ## Return a random non-negative `BigInt`, up to `max`, using the given state. + rand(r, 0.initBigInt..max) + +# backwards compatibility with 1.4 +when not defined(randState): + var state = initRand(777) + proc randState(): var Rand = state + +proc rand*(x: Slice[BigInt]): BigInt = rand(randState(), x) + ## Return a random `BigInt`, within the given range. + +proc rand*(max: BigInt): BigInt = rand(randState(), max) + ## Return a random `BigInt`, up to `max`. diff --git a/tests/trandom.nim b/tests/trandom.nim new file mode 100644 index 0000000..c7f2b7e --- /dev/null +++ b/tests/trandom.nim @@ -0,0 +1,22 @@ +import bigints +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) +