summary refs log tree commit diff stats
path: root/tests/stdlib/tmath.nim
blob: ff1f32d364ffbadd87c1a0cccfb08632cedb0c14 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
discard """
  action: run
  output: '''[Suite] random int

[Suite] random float

[Suite] cumsum

[Suite] random sample

[Suite] ^

'''
"""

import math, random, os
import unittest
import sets, tables

suite "random int":
  test "there might be some randomness":
    var set = initHashSet[int](128)

    for i in 1..1000:
      incl(set, rand(high(int)))
    check len(set) == 1000
  test "single number bounds work":

    var rand: int
    for i in 1..1000:
      rand = rand(1000)
      check rand < 1000
      check rand > -1
  test "slice bounds work":

    var rand: int
    for i in 1..1000:
      rand = rand(100..1000)
      check rand < 1000
      check rand >= 100
  test " again gives new numbers":

    var rand1 = rand(1000000)
    os.sleep(200)

    var rand2 = rand(1000000)
    check rand1 != rand2


suite "random float":
  test "there might be some randomness":
    var set = initHashSet[float](128)

    for i in 1..100:
      incl(set, rand(1.0))
    check len(set) == 100
  test "single number bounds work":

    var rand: float
    for i in 1..1000:
      rand = rand(1000.0)
      check rand < 1000.0
      check rand > -1.0
  test "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
  test " again gives new numbers":

    var rand1:float = rand(1000000.0)
    os.sleep(200)

    var rand2:float = rand(1000000.0)
    check rand1 != rand2

suite "cumsum":
  test "cumsum int seq return":
    let counts = [ 1, 2, 3, 4 ]
    check counts.cumsummed == [ 1, 3, 6, 10 ]

  test "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 ]

  test "cumsum int in-place":
    var counts = [ 1, 2, 3, 4 ]
    counts.cumsum
    check counts == [ 1, 3, 6, 10 ]

  test "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 ]

suite "random sample":
  test "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

  test "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

suite "^":
  test "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)

block:
  when not defined(js):
    # Check for no side effect annotation
    proc mySqrt(num: float): float {.noSideEffect.} =
      return sqrt(num)

    # check gamma function
    assert(gamma(5.0) == 24.0) # 4!
    assert(lgamma(1.0) == 0.0) # ln(1.0) == 0.0
    assert(erf(6.0) > erf(5.0))
    assert(erfc(6.0) < erfc(5.0))


    # Function for approximate comparison of floats
    proc `==~`(x, y: float): bool = (abs(x-y) < 1e-9)

    block: # prod
      doAssert prod([1, 2, 3, 4]) == 24
      doAssert prod([1.5, 3.4]) == 5.1
      let x: seq[float] = @[]
      doAssert prod(x) == 1.0

    block: # round() tests
      # 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

    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)
      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)
      doAssert(classify(trunc(f_nan.float32)) == fcNan)
      doAssert(classify(trunc(0.0'f32)) == fcZero)

    block: # sgn() tests
      assert sgn(1'i8) == 1
      assert sgn(1'i16) == 1
      assert sgn(1'i32) == 1
      assert sgn(1'i64) == 1
      assert sgn(1'u8) == 1
      assert sgn(1'u16) == 1
      assert sgn(1'u32) == 1
      assert sgn(1'u64) == 1
      assert sgn(-12342.8844'f32) == -1
      assert sgn(123.9834'f64) == 1
      assert sgn(0'i32) == 0
      assert sgn(0'f32) == 0
      assert sgn(NegInf) == -1
      assert sgn(Inf) == 1
      assert sgn(NaN) == 0

    block: # fac() tests
      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: # 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