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

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