Improve math module (#17019)

* Improve documentation for math

Support empty input for cumsummed
Use runnableExamples
Move some examples to tests
Add more tests

* Update tests/stdlib/tmath.nim

Move some tests to trandom.nim
Move tests into main template where possible
Add test for #17017

* Add more tests for gamma & lgamma

Remove gamma(-1.0) example
Small fixes/changes

* Move more tests into template main()

* Fix typos

* Add edge case examples for copySign

2 files changed, 358 insertions, 311 deletions

diff --git a/tests/stdlib/tmath.nim b/tests/stdlib/tmath.nim
index 7af9eb73e..894c79369 100644
--- a/tests/stdlib/tmath.nim
+++ b/tests/stdlib/tmath.nim
@@ -6,296 +6,197 @@ discard """
 # xxx: there should be a test with `-d:nimTmathCase2 -d:danger --passc:-ffast-math`,
 # but it requires disabling certain lines with `when not defined(nimTmathCase2)`
 
-import std/[math, random, os]
-import std/[unittest]
-import std/[sets, tables]
-
+import std/math
 
 # Function for approximate comparison of floats
-proc `==~`(x, y: float): bool = (abs(x-y) < 1e-9)
-
-block: # random int
-  block: # there might be some randomness
-    var set = initHashSet[int](128)
-
-    for i in 1..1000:
-      incl(set, rand(high(int)))
-    check len(set) == 1000
-
-  block: # single number bounds work
-    var rand: int
-    for i in 1..1000:
-      rand = rand(1000)
-      check rand < 1000
-      check rand > -1
-
-  block: # slice bounds work
-    var rand: int
-    for i in 1..1000:
-      rand = rand(100..1000)
-      when defined(js): # xxx bug: otherwise fails
-        check rand <= 1000
-      else:
-        check rand < 1000
-      check rand >= 100
-
-  block: # again gives new numbers
-    var rand1 = rand(1000000)
-    when not defined(js):
-      os.sleep(200)
-
-    var rand2 = rand(1000000)
-    check rand1 != rand2
-
-
-block: # random float
-  block: # there might be some randomness
-    var set = initHashSet[float](128)
-
-    for i in 1..100:
-      incl(set, rand(1.0))
-    check len(set) == 100
-
-  block: # single number bounds work
-    var rand: float
-    for i in 1..1000:
-      rand = rand(1000.0)
-      check rand < 1000.0
-      check rand > -1.0
-
-  block: # slice bounds work
-    var rand: float
-    for i in 1..1000:
-      rand = rand(100.0..1000.0)
-      check rand < 1000.0
-      check rand >= 100.0
-
-  block: # again gives new numbers
-
-    var rand1:float = rand(1000000.0)
-    when not defined(js):
-      os.sleep(200)
-
-    var rand2:float = rand(1000000.0)
-    check rand1 != rand2
-
-block: # cumsum
-  block: # cumsum int seq return
-    let counts = [ 1, 2, 3, 4 ]
-    check counts.cumsummed == [ 1, 3, 6, 10 ]
-
-  block: # cumsum float seq return
-    let counts = [ 1.0, 2.0, 3.0, 4.0 ]
-    check counts.cumsummed == [ 1.0, 3.0, 6.0, 10.0 ]
-
-  block: # cumsum int in-place
-    var counts = [ 1, 2, 3, 4 ]
-    counts.cumsum
-    check counts == [ 1, 3, 6, 10 ]
-
-  block: # cumsum float in-place
-    var counts = [ 1.0, 2.0, 3.0, 4.0 ]
-    counts.cumsum
-    check counts == [ 1.0, 3.0, 6.0, 10.0 ]
-
-block: # random sample
-  block: # "non-uniform array sample unnormalized int CDF
-    let values = [ 10, 20, 30, 40, 50 ] # values
-    let counts = [ 4, 3, 2, 1, 0 ]      # weights aka unnormalized probabilities
-    var histo = initCountTable[int]()
-    let cdf = counts.cumsummed          # unnormalized CDF
-    for i in 0 ..< 5000:
-      histo.inc(sample(values, cdf))
-    check histo.len == 4                # number of non-zero in `counts`
-    # Any one bin is a binomial random var for n samples, each with prob p of
-    # adding a count to k; E[k]=p*n, Var k=p*(1-p)*n, approximately Normal for
-    # big n.  So, P(abs(k - p*n)/sqrt(p*(1-p)*n))>3.0) =~ 0.0027, while
-    # P(wholeTestFails) =~ 1 - P(binPasses)^4 =~ 1 - (1-0.0027)^4 =~ 0.01.
-    for i, c in counts:
-      if c == 0:
-        check values[i] notin histo
-        continue
-      let p = float(c) / float(cdf[^1])
-      let n = 5000.0
-      let expected = p * n
-      let stdDev = sqrt(n * p * (1.0 - p))
-      check abs(float(histo[values[i]]) - expected) <= 3.0 * stdDev
-
-  block: # non-uniform array sample normalized float CDF
-    let values = [ 10, 20, 30, 40, 50 ]     # values
-    let counts = [ 0.4, 0.3, 0.2, 0.1, 0 ]  # probabilities
-    var histo = initCountTable[int]()
-    let cdf = counts.cumsummed              # normalized CDF
-    for i in 0 ..< 5000:
-      histo.inc(sample(values, cdf))
-    check histo.len == 4                    # number of non-zero in ``counts``
-    for i, c in counts:
-      if c == 0:
-        check values[i] notin histo
-        continue
-      let p = float(c) / float(cdf[^1])
-      let n = 5000.0
-      let expected = p * n
-      let stdDev = sqrt(n * p * (1.0 - p))
-      # NOTE: like unnormalized int CDF test, P(wholeTestFails) =~ 0.01.
-      check abs(float(histo[values[i]]) - expected) <= 3.0 * stdDev
-
-block: # ^
-  block: # compiles for valid types
-    check: compiles(5 ^ 2)
-    check: compiles(5.5 ^ 2)
-    check: compiles(5.5 ^ 2.int8)
-    check: compiles(5.5 ^ 2.uint)
-    check: compiles(5.5 ^ 2.uint8)
-    check: not compiles(5.5 ^ 2.2)
+proc `==~`(x, y: float): bool = abs(x - y) < 1e-9
 
 block:
   when not defined(js):
-    # Check for no side effect annotation
+    # check for no side effect annotation
     proc mySqrt(num: float): float {.noSideEffect.} =
       # xxx unused
-      return sqrt(num)
+      sqrt(num)
 
     # check gamma function
-    doAssert(gamma(5.0) == 24.0) # 4!
-    doAssert(lgamma(1.0) == 0.0) # ln(1.0) == 0.0
-    doAssert(erf(6.0) > erf(5.0))
-    doAssert(erfc(6.0) < erfc(5.0))
-
-    block: # prod
-      doAssert prod([1, 2, 3, 4]) == 24
-      doAssert prod([1.5, 3.4]).almostEqual 5.1
-      let x: seq[float] = @[]
-      doAssert prod(x) == 1.0
-
-    block: # splitDecimal() tests
-      doAssert splitDecimal(54.674).intpart == 54.0
-      doAssert splitDecimal(54.674).floatpart ==~ 0.674
-      doAssert splitDecimal(-693.4356).intpart == -693.0
-      doAssert splitDecimal(-693.4356).floatpart ==~ -0.4356
-      doAssert splitDecimal(0.0).intpart == 0.0
-      doAssert splitDecimal(0.0).floatpart == 0.0
-
-    block: # trunc tests for vcc
-      doAssert(trunc(-1.1) == -1)
-      doAssert(trunc(1.1) == 1)
-      doAssert(trunc(-0.1) == -0)
-      doAssert(trunc(0.1) == 0)
-
-      #special case
-      doAssert(classify(trunc(1e1000000)) == fcInf)
-      doAssert(classify(trunc(-1e1000000)) == fcNegInf)
-      when not defined(nimTmathCase2):
-        doAssert(classify(trunc(0.0/0.0)) == fcNan)
-      doAssert(classify(trunc(0.0)) == fcZero)
-
-      #trick the compiler to produce signed zero
-      let
-        f_neg_one = -1.0
-        f_zero = 0.0
-        f_nan = f_zero / f_zero
-
-      doAssert(classify(trunc(f_neg_one*f_zero)) == fcNegZero)
-
-      doAssert(trunc(-1.1'f32) == -1)
-      doAssert(trunc(1.1'f32) == 1)
-      doAssert(trunc(-0.1'f32) == -0)
-      doAssert(trunc(0.1'f32) == 0)
-      doAssert(classify(trunc(1e1000000'f32)) == fcInf)
-      doAssert(classify(trunc(-1e1000000'f32)) == fcNegInf)
-      when not defined(nimTmathCase2):
-        doAssert(classify(trunc(f_nan.float32)) == fcNan)
-      doAssert(classify(trunc(0.0'f32)) == fcZero)
-
-    block: # sgn() tests
-      doAssert sgn(1'i8) == 1
-      doAssert sgn(1'i16) == 1
-      doAssert sgn(1'i32) == 1
-      doAssert sgn(1'i64) == 1
-      doAssert sgn(1'u8) == 1
-      doAssert sgn(1'u16) == 1
-      doAssert sgn(1'u32) == 1
-      doAssert sgn(1'u64) == 1
-      doAssert sgn(-12342.8844'f32) == -1
-      doAssert sgn(123.9834'f64) == 1
-      doAssert sgn(0'i32) == 0
-      doAssert sgn(0'f32) == 0
-      doAssert sgn(NegInf) == -1
-      doAssert sgn(Inf) == 1
-      doAssert sgn(NaN) == 0
-
-    block: # fac() tests
+    doAssert gamma(5.0) == 24.0 # 4!
+    doAssert almostEqual(gamma(0.5), sqrt(PI))
+    doAssert almostEqual(gamma(-0.5), -2 * sqrt(PI))
+    doAssert lgamma(1.0) == 0.0 # ln(1.0) == 0.0
+    doAssert almostEqual(lgamma(0.5), 0.5 * ln(PI))
+    doAssert erf(6.0) > erf(5.0)
+    doAssert erfc(6.0) < erfc(5.0)
+
+when not defined(js) and not defined(windows): # xxx pending bug #17017
+  doAssert gamma(-1.0).isNaN
+
+template main() =
+  block: # sgn() tests
+    doAssert sgn(1'i8) == 1
+    doAssert sgn(1'i16) == 1
+    doAssert sgn(1'i32) == 1
+    doAssert sgn(1'i64) == 1
+    doAssert sgn(1'u8) == 1
+    doAssert sgn(1'u16) == 1
+    doAssert sgn(1'u32) == 1
+    doAssert sgn(1'u64) == 1
+    doAssert sgn(-12342.8844'f32) == -1
+    doAssert sgn(123.9834'f64) == 1
+    doAssert sgn(0'i32) == 0
+    doAssert sgn(0'f32) == 0
+    doAssert sgn(NegInf) == -1
+    doAssert sgn(Inf) == 1
+    doAssert sgn(NaN) == 0
+
+  block: # fac() tests
+    when nimvm: discard
+    else:
       try:
         discard fac(-1)
       except AssertionDefect:
         discard
 
-      doAssert fac(0) == 1
-      doAssert fac(1) == 1
-      doAssert fac(2) == 2
-      doAssert fac(3) == 6
-      doAssert fac(4) == 24
-
-    block: # floorMod/floorDiv
-      doAssert floorDiv(8, 3) == 2
-      doAssert floorMod(8, 3) == 2
-
-      doAssert floorDiv(8, -3) == -3
-      doAssert floorMod(8, -3) == -1
-
-      doAssert floorDiv(-8, 3) == -3
-      doAssert floorMod(-8, 3) == 1
-
-      doAssert floorDiv(-8, -3) == 2
-      doAssert floorMod(-8, -3) == -2
-
-      doAssert floorMod(8.0, -3.0) == -1.0
-      doAssert floorMod(-8.5, 3.0) == 0.5
-
-    block: # euclDiv/euclMod
-      doAssert euclDiv(8, 3) == 2
-      doAssert euclMod(8, 3) == 2
-
-      doAssert euclDiv(8, -3) == -2
-      doAssert euclMod(8, -3) == 2
-
-      doAssert euclDiv(-8, 3) == -3
-      doAssert euclMod(-8, 3) == 1
+    doAssert fac(0) == 1
+    doAssert fac(1) == 1
+    doAssert fac(2) == 2
+    doAssert fac(3) == 6
+    doAssert fac(4) == 24
+    doAssert fac(5) == 120
+
+  block: # floorMod/floorDiv
+    doAssert floorDiv(8, 3) == 2
+    doAssert floorMod(8, 3) == 2
+
+    doAssert floorDiv(8, -3) == -3
+    doAssert floorMod(8, -3) == -1
+
+    doAssert floorDiv(-8, 3) == -3
+    doAssert floorMod(-8, 3) == 1
+
+    doAssert floorDiv(-8, -3) == 2
+    doAssert floorMod(-8, -3) == -2
+
+    doAssert floorMod(8.0, -3.0) == -1.0
+    doAssert floorMod(-8.5, 3.0) == 0.5
+
+  block: # euclDiv/euclMod
+    doAssert euclDiv(8, 3) == 2
+    doAssert euclMod(8, 3) == 2
+
+    doAssert euclDiv(8, -3) == -2
+    doAssert euclMod(8, -3) == 2
+
+    doAssert euclDiv(-8, 3) == -3
+    doAssert euclMod(-8, 3) == 1
+
+    doAssert euclDiv(-8, -3) == 3
+    doAssert euclMod(-8, -3) == 1
+
+    doAssert euclMod(8.0, -3.0) == 2.0
+    doAssert euclMod(-8.5, 3.0) == 0.5
+
+    doAssert euclDiv(9, 3) == 3
+    doAssert euclMod(9, 3) == 0
+
+    doAssert euclDiv(9, -3) == -3
+    doAssert euclMod(9, -3) == 0
+
+    doAssert euclDiv(-9, 3) == -3
+    doAssert euclMod(-9, 3) == 0
+
+    doAssert euclDiv(-9, -3) == 3
+    doAssert euclMod(-9, -3) == 0
+
+  block: # splitDecimal() tests
+    doAssert splitDecimal(54.674).intpart == 54.0
+    doAssert splitDecimal(54.674).floatpart ==~ 0.674
+    doAssert splitDecimal(-693.4356).intpart == -693.0
+    doAssert splitDecimal(-693.4356).floatpart ==~ -0.4356
+    doAssert splitDecimal(0.0).intpart == 0.0
+    doAssert splitDecimal(0.0).floatpart == 0.0
+
+  block: # trunc tests for vcc
+    doAssert trunc(-1.1) == -1
+    doAssert trunc(1.1) == 1
+    doAssert trunc(-0.1) == -0
+    doAssert trunc(0.1) == 0
+
+    # special case
+    doAssert classify(trunc(1e1000000)) == fcInf
+    doAssert classify(trunc(-1e1000000)) == fcNegInf
+    when not defined(nimTmathCase2):
+      doAssert classify(trunc(0.0/0.0)) == fcNan
+    doAssert classify(trunc(0.0)) == fcZero
+
+    # trick the compiler to produce signed zero
+    let
+      f_neg_one = -1.0
+      f_zero = 0.0
+      f_nan = f_zero / f_zero
+
+    doAssert classify(trunc(f_neg_one*f_zero)) == fcNegZero
+
+    doAssert trunc(-1.1'f32) == -1
+    doAssert trunc(1.1'f32) == 1
+    doAssert trunc(-0.1'f32) == -0
+    doAssert trunc(0.1'f32) == 0
+    doAssert classify(trunc(1e1000000'f32)) == fcInf
+    doAssert classify(trunc(-1e1000000'f32)) == fcNegInf
+    when not defined(nimTmathCase2):
+      doAssert classify(trunc(f_nan.float32)) == fcNan
+    doAssert classify(trunc(0.0'f32)) == fcZero
+
+  block: # log
+    doAssert log(4.0, 3.0) ==~ ln(4.0) / ln(3.0)
+    doAssert log2(8.0'f64) == 3.0'f64
+    doAssert log2(4.0'f64) == 2.0'f64
+    doAssert log2(2.0'f64) == 1.0'f64
+    doAssert log2(1.0'f64) == 0.0'f64
+    doAssert classify(log2(0.0'f64)) == fcNegInf
+
+    doAssert log2(8.0'f32) == 3.0'f32
+    doAssert log2(4.0'f32) == 2.0'f32
+    doAssert log2(2.0'f32) == 1.0'f32
+    doAssert log2(1.0'f32) == 0.0'f32
+    doAssert classify(log2(0.0'f32)) == fcNegInf
+
+  block: # cumsum
+    block: # cumsum int seq return
+      let counts = [1, 2, 3, 4]
+      doAssert counts.cumsummed == @[1, 3, 6, 10]
+      let empty: seq[int] = @[]
+      doAssert empty.cumsummed == @[]
+
+    block: # cumsum float seq return
+      let counts = [1.0, 2.0, 3.0, 4.0]
+      doAssert counts.cumsummed == @[1.0, 3.0, 6.0, 10.0]
+      let empty: seq[float] = @[]
+      doAssert empty.cumsummed == @[]
+
+    block: # cumsum int in-place
+      var counts = [1, 2, 3, 4]
+      counts.cumsum
+      doAssert counts == [1, 3, 6, 10]
+      var empty: seq[int] = @[]
+      empty.cumsum
+      doAssert empty == @[]
+
+    block: # cumsum float in-place
+      var counts = [1.0, 2.0, 3.0, 4.0]
+      counts.cumsum
+      doAssert counts == [1.0, 3.0, 6.0, 10.0]
+      var empty: seq[float] = @[]
+      empty.cumsum
+      doAssert empty == @[]
+
+  block: # ^ compiles for valid types
+    doAssert: compiles(5 ^ 2)
+    doAssert: compiles(5.5 ^ 2)
+    doAssert: compiles(5.5 ^ 2.int8)
+    doAssert: compiles(5.5 ^ 2.uint)
+    doAssert: compiles(5.5 ^ 2.uint8)
+    doAssert: not compiles(5.5 ^ 2.2)
 
-      doAssert euclDiv(-8, -3) == 3
-      doAssert euclMod(-8, -3) == 1
-
-      doAssert euclMod(8.0, -3.0) == 2.0
-      doAssert euclMod(-8.5, 3.0) == 0.5
-
-      doAssert euclDiv(9, 3) == 3
-      doAssert euclMod(9, 3) == 0
-
-      doAssert euclDiv(9, -3) == -3
-      doAssert euclMod(9, -3) == 0
-
-      doAssert euclDiv(-9, 3) == -3
-      doAssert euclMod(-9, 3) == 0
-
-      doAssert euclDiv(-9, -3) == 3
-      doAssert euclMod(-9, -3) == 0
-
-    block: # log
-      doAssert log(4.0, 3.0) ==~ ln(4.0) / ln(3.0)
-      doAssert log2(8.0'f64) == 3.0'f64
-      doAssert log2(4.0'f64) == 2.0'f64
-      doAssert log2(2.0'f64) == 1.0'f64
-      doAssert log2(1.0'f64) == 0.0'f64
-      doAssert classify(log2(0.0'f64)) == fcNegInf
-
-      doAssert log2(8.0'f32) == 3.0'f32
-      doAssert log2(4.0'f32) == 2.0'f32
-      doAssert log2(2.0'f32) == 1.0'f32
-      doAssert log2(1.0'f32) == 0.0'f32
-      doAssert classify(log2(0.0'f32)) == fcNegInf
-
-template main =
-  # xxx wrap all under `main` so it also gets tested in vm.
   block: # signbit
     doAssert not signbit(0.0)
     doAssert signbit(-0.0)
@@ -310,6 +211,8 @@ template main =
     doAssert NaN.isNaN
     doAssert not Inf.isNaN
     doAssert isNaN(Inf - Inf)
+    doAssert not isNaN(3.1415926)
+    doAssert not isNaN(0'f32)
 
   block: # signbit
     let x1 = NaN
@@ -324,16 +227,19 @@ template main =
     doAssert signbit(x3)
 
   block: # copySign
+    doAssert copySign(10.0, 1.0) == 10.0
     doAssert copySign(10.0, -1.0) == -10.0
     doAssert copySign(-10.0, -1.0) == -10.0
     doAssert copySign(-10.0, 1.0) == 10.0
     doAssert copySign(float(10), -1.0) == -10.0
 
+    doAssert copySign(10.0'f64, 1.0) == 10.0
     doAssert copySign(10.0'f64, -1.0) == -10.0
     doAssert copySign(-10.0'f64, -1.0) == -10.0
     doAssert copySign(-10.0'f64, 1.0) == 10.0
     doAssert copySign(10'f64, -1.0) == -10.0
 
+    doAssert copySign(10.0'f32, 1.0) == 10.0
     doAssert copySign(10.0'f32, -1.0) == -10.0
     doAssert copySign(-10.0'f32, -1.0) == -10.0
     doAssert copySign(-10.0'f32, 1.0) == 10.0
@@ -363,37 +269,47 @@ template main =
     doAssert copySign(-NaN, 0.0).isNaN
     doAssert copySign(-NaN, -0.0).isNaN
 
-    block: # round() tests
-      block: # Round to 0 decimal places
-        doAssert round(54.652) == 55.0
-        doAssert round(54.352) == 54.0
-        doAssert round(-54.652) == -55.0
-        doAssert round(-54.352) == -54.0
-        doAssert round(0.0) == 0.0
-        doAssert 1 / round(0.0) == Inf
-        doAssert 1 / round(-0.0) == -Inf
-        doAssert round(Inf) == Inf
-        doAssert round(-Inf) == -Inf
-        doAssert round(NaN).isNaN
-        doAssert round(-NaN).isNaN
-        doAssert round(-0.5) == -1.0
-        doAssert round(0.5) == 1.0
-        doAssert round(-1.5) == -2.0
-        doAssert round(1.5) == 2.0
-        doAssert round(-2.5) == -3.0
-        doAssert round(2.5) == 3.0
-        doAssert round(2.5'f32) == 3.0'f32
-        doAssert round(2.5'f64) == 3.0'f64
-      block: # func round*[T: float32|float64](x: T, places: int): T
-        doAssert round(54.345, 0) == 54.0
-        template fn(x) =
-          doAssert round(x, 2).almostEqual 54.35
-          doAssert round(x, 2).almostEqual 54.35
-          doAssert round(x, -1).almostEqual 50.0
-          doAssert round(x, -2).almostEqual 100.0
-          doAssert round(x, -3).almostEqual 0.0
-        fn(54.346)
-        fn(54.346'f32)
+  block: # almostEqual
+    doAssert almostEqual(3.141592653589793, 3.1415926535897936)
+    doAssert almostEqual(1.6777215e7'f32, 1.6777216e7'f32)
+    doAssert almostEqual(Inf, Inf)
+    doAssert almostEqual(-Inf, -Inf)
+    doAssert not almostEqual(Inf, -Inf)
+    doAssert not almostEqual(-Inf, Inf)
+    doAssert not almostEqual(Inf, NaN)
+    doAssert not almostEqual(NaN, NaN)
+
+  block: # round() tests
+    block: # Round to 0 decimal places
+      doAssert round(54.652) == 55.0
+      doAssert round(54.352) == 54.0
+      doAssert round(-54.652) == -55.0
+      doAssert round(-54.352) == -54.0
+      doAssert round(0.0) == 0.0
+      doAssert 1 / round(0.0) == Inf
+      doAssert 1 / round(-0.0) == -Inf
+      doAssert round(Inf) == Inf
+      doAssert round(-Inf) == -Inf
+      doAssert round(NaN).isNaN
+      doAssert round(-NaN).isNaN
+      doAssert round(-0.5) == -1.0
+      doAssert round(0.5) == 1.0
+      doAssert round(-1.5) == -2.0
+      doAssert round(1.5) == 2.0
+      doAssert round(-2.5) == -3.0
+      doAssert round(2.5) == 3.0
+      doAssert round(2.5'f32) == 3.0'f32
+      doAssert round(2.5'f64) == 3.0'f64
+    block: # func round*[T: float32|float64](x: T, places: int): T
+      doAssert round(54.345, 0) == 54.0
+      template fn(x) =
+        doAssert round(x, 2).almostEqual 54.35
+        doAssert round(x, 2).almostEqual 54.35
+        doAssert round(x, -1).almostEqual 50.0
+        doAssert round(x, -2).almostEqual 100.0
+        doAssert round(x, -3).almostEqual 0.0
+      fn(54.346)
+      fn(54.346'f32)
 
     when nimvm:
       discard
@@ -416,5 +332,36 @@ template main =
     doAssert abs(NaN).isNaN
     doAssert abs(-NaN).isNaN
 
+  block: # classify
+    doAssert classify(0.3) == fcNormal
+    doAssert classify(-0.3) == fcNormal
+    doAssert classify(5.0e-324) == fcSubnormal
+    doAssert classify(-5.0e-324) == fcSubnormal
+    doAssert classify(0.0) == fcZero
+    doAssert classify(-0.0) == fcNegZero
+    doAssert classify(NaN) == fcNan
+    doAssert classify(0.3 / 0.0) == fcInf
+    doAssert classify(Inf) == fcInf
+    doAssert classify(-0.3 / 0.0) == fcNegInf
+    doAssert classify(-Inf) == fcNegInf
+
+  block: # sum
+    let empty: seq[int] = @[]
+    doAssert sum(empty) == 0
+    doAssert sum([1, 2, 3, 4]) == 10
+    doAssert sum([-4, 3, 5]) == 4
+
+  block: # prod
+    let empty: seq[int] = @[]
+    doAssert prod(empty) == 1
+    doAssert prod([1, 2, 3, 4]) == 24
+    doAssert prod([-4, 3, 5]) == -60
+    doAssert almostEqual(prod([1.5, 3.4]), 5.1)
+    let x: seq[float] = @[]
+    doAssert prod(x) == 1.0
+
+  when not defined(windows): # xxx pending bug #17017
+    doAssert sqrt(-1.0).isNaN
+
 static: main()
 main()
diff --git a/tests/stdlib/trandom.nim b/tests/stdlib/trandom.nim
index f8797944e..9fc68fca1 100644
--- a/tests/stdlib/trandom.nim
+++ b/tests/stdlib/trandom.nim
@@ -3,11 +3,11 @@ discard """
   targets: "c js"
 """
 
-import std/[random, stats]
+import std/[random, math, os, stats, sets, tables]
 
 randomize(233)
 
-proc main =
+proc main() =
   var occur: array[1000, int]
 
   for i in 0..100_000:
@@ -33,11 +33,8 @@ proc main =
   # don't use causes integer overflow
   doAssert compiles(rand[int](low(int) .. high(int)))
 
-
 main()
 
-import math
-
 block:
   when not defined(js):
     doAssert almostEqual(rand(12.5), 4.012897747078944)
@@ -64,3 +61,106 @@ block:
     if rand(5.DiceRoll) < 3:
       flag = true
   doAssert flag # because of: rand(max: int): int
+
+
+block: # random int
+  block: # there might be some randomness
+    var set = initHashSet[int](128)
+
+    for i in 1..1000:
+      incl(set, rand(high(int)))
+    doAssert len(set) == 1000
+
+  block: # single number bounds work
+    var rand: int
+    for i in 1..1000:
+      rand = rand(1000)
+      doAssert rand <= 1000
+      doAssert rand >= 0
+
+  block: # slice bounds work
+    var rand: int
+    for i in 1..1000:
+      rand = rand(100..1000)
+      doAssert rand <= 1000
+      doAssert rand >= 100
+
+  block: # again gives new numbers
+    var rand1 = rand(1000000)
+    when not defined(js):
+      os.sleep(200)
+
+    var rand2 = rand(1000000)
+    doAssert rand1 != rand2
+
+block: # random float
+  block: # there might be some randomness
+    var set = initHashSet[float](128)
+
+    for i in 1..100:
+      incl(set, rand(1.0))
+    doAssert len(set) == 100
+
+  block: # single number bounds work
+    var rand: float
+    for i in 1..1000:
+      rand = rand(1000.0)
+      doAssert rand <= 1000.0
+      doAssert rand >= 0.0
+
+  block: # slice bounds work
+    var rand: float
+    for i in 1..1000:
+      rand = rand(100.0..1000.0)
+      doAssert rand <= 1000.0
+      doAssert rand >= 100.0
+
+  block: # again gives new numbers
+    var rand1: float = rand(1000000.0)
+    when not defined(js):
+      os.sleep(200)
+
+    var rand2: float = rand(1000000.0)
+    doAssert rand1 != rand2
+
+block: # random sample
+  block: # "non-uniform array sample unnormalized int CDF
+    let values = [10, 20, 30, 40, 50] # values
+    let counts = [4, 3, 2, 1, 0]      # weights aka unnormalized probabilities
+    var histo = initCountTable[int]()
+    let cdf = counts.cumsummed        # unnormalized CDF
+    for i in 0 ..< 5000:
+      histo.inc(sample(values, cdf))
+    doAssert histo.len == 4              # number of non-zero in `counts`
+    # Any one bin is a binomial random var for n samples, each with prob p of
+    # adding a count to k; E[k]=p*n, Var k=p*(1-p)*n, approximately Normal for
+    # big n.  So, P(abs(k - p*n)/sqrt(p*(1-p)*n))>3.0) =~ 0.0027, while
+    # P(wholeTestFails) =~ 1 - P(binPasses)^4 =~ 1 - (1-0.0027)^4 =~ 0.01.
+    for i, c in counts:
+      if c == 0:
+        doAssert values[i] notin histo
+        continue
+      let p = float(c) / float(cdf[^1])
+      let n = 5000.0
+      let expected = p * n
+      let stdDev = sqrt(n * p * (1.0 - p))
+      doAssert abs(float(histo[values[i]]) - expected) <= 3.0 * stdDev
+
+  block: # non-uniform array sample normalized float CDF
+    let values = [10, 20, 30, 40, 50]     # values
+    let counts = [0.4, 0.3, 0.2, 0.1, 0]  # probabilities
+    var histo = initCountTable[int]()
+    let cdf = counts.cumsummed            # normalized CDF
+    for i in 0 ..< 5000:
+      histo.inc(sample(values, cdf))
+    doAssert histo.len == 4                  # number of non-zero in ``counts``
+    for i, c in counts:
+      if c == 0:
+        doAssert values[i] notin histo
+        continue
+      let p = float(c) / float(cdf[^1])
+      let n = 5000.0
+      let expected = p * n
+      let stdDev = sqrt(n * p * (1.0 - p))
+      # NOTE: like unnormalized int CDF test, P(wholeTestFails) =~ 0.01.
+      doAssert abs(float(histo[values[i]]) - expected) <= 3.0 * stdDev