random.nim 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714
  1. #
  2. #
  3. # Nim's Runtime Library
  4. # (c) Copyright 2017 Andreas Rumpf
  5. #
  6. # See the file "copying.txt", included in this
  7. # distribution, for details about the copyright.
  8. #
  9. ## Nim's standard random number generator (RNG).
  10. ##
  11. ## Its implementation is based on the `xoroshiro128+`
  12. ## (xor/rotate/shift/rotate) library.
  13. ## * More information: http://xoroshiro.di.unimi.it
  14. ## * C implementation: http://xoroshiro.di.unimi.it/xoroshiro128plus.c
  15. ##
  16. ## **Do not use this module for cryptographic purposes!**
  17. ##
  18. ## Basic usage
  19. ## ===========
  20. ##
  21. runnableExamples:
  22. # Call randomize() once to initialize the default random number generator.
  23. # If this is not called, the same results will occur every time these
  24. # examples are run.
  25. randomize()
  26. # Pick a number in 0..100.
  27. let num = rand(100)
  28. doAssert num in 0..100
  29. # Roll a six-sided die.
  30. let roll = rand(1..6)
  31. doAssert roll in 1..6
  32. # Pick a marble from a bag.
  33. let marbles = ["red", "blue", "green", "yellow", "purple"]
  34. let pick = sample(marbles)
  35. doAssert pick in marbles
  36. # Shuffle some cards.
  37. var cards = ["Ace", "King", "Queen", "Jack", "Ten"]
  38. shuffle(cards)
  39. doAssert cards.len == 5
  40. ## These examples all use the default RNG. The
  41. ## `Rand type <#Rand>`_ represents the state of an RNG.
  42. ## For convenience, this module contains a default Rand state that corresponds
  43. ## to the default RNG. Most procs in this module which do
  44. ## not take in a Rand parameter, including those called in the above examples,
  45. ## use the default generator. Those procs are **not** thread-safe.
  46. ##
  47. ## Note that the default generator always starts in the same state.
  48. ## The `randomize proc <#randomize>`_ can be called to initialize the default
  49. ## generator with a seed based on the current time, and it only needs to be
  50. ## called once before the first usage of procs from this module. If
  51. ## `randomize` is not called, the default generator will always produce
  52. ## the same results.
  53. ##
  54. ## RNGs that are independent of the default one can be created with the
  55. ## `initRand proc <#initRand,int64>`_.
  56. ##
  57. ## Again, it is important to remember that this module must **not** be used for
  58. ## cryptographic applications.
  59. ##
  60. ## See also
  61. ## ========
  62. ## * `std/sysrand module <sysrand.html>`_ for a cryptographically secure pseudorandom number generator
  63. ## * `mersenne module <mersenne.html>`_ for the Mersenne Twister random number generator
  64. ## * `math module <math.html>`_ for basic math routines
  65. ## * `stats module <stats.html>`_ for statistical analysis
  66. ## * `list of cryptographic and hashing modules <lib.html#pure-libraries-hashing>`_
  67. ## in the standard library
  68. import algorithm, math
  69. import std/private/since
  70. include system/inclrtl
  71. {.push debugger: off.}
  72. when defined(js):
  73. type Ui = uint32
  74. const randMax = 4_294_967_295u32
  75. else:
  76. type Ui = uint64
  77. const randMax = 18_446_744_073_709_551_615u64
  78. type
  79. Rand* = object ## State of a random number generator.
  80. ##
  81. ## Create a new Rand state using the `initRand proc <#initRand,int64>`_.
  82. ##
  83. ## The module contains a default Rand state for convenience.
  84. ## It corresponds to the default RNG's state.
  85. ## The default Rand state always starts with the same values, but the
  86. ## `randomize proc <#randomize>`_ can be used to seed the default generator
  87. ## with a value based on the current time.
  88. ##
  89. ## Many procs have two variations: one that takes in a Rand parameter and
  90. ## another that uses the default generator. The procs that use the default
  91. ## generator are **not** thread-safe!
  92. a0, a1: Ui
  93. when defined(js):
  94. var state = Rand(
  95. a0: 0x69B4C98Cu32,
  96. a1: 0xFED1DD30u32) # global for backwards compatibility
  97. else:
  98. const DefaultRandSeed = Rand(
  99. a0: 0x69B4C98CB8530805u64,
  100. a1: 0xFED1DD3004688D67CAu64)
  101. # racy for multi-threading but good enough for now:
  102. var state = DefaultRandSeed # global for backwards compatibility
  103. func isValid(r: Rand): bool {.inline.} =
  104. ## Check whether state of `r` is valid.
  105. ##
  106. ## In `xoroshiro128+`, if all bits of `a0` and `a1` are zero,
  107. ## they are always zero after calling `next(r: var Rand)`.
  108. not (r.a0 == 0 and r.a1 == 0)
  109. since (1, 5):
  110. template randState*(): untyped =
  111. ## Makes the default Rand state accessible from other modules.
  112. ## Useful for module authors.
  113. state
  114. proc rotl(x, k: Ui): Ui =
  115. result = (x shl k) or (x shr (Ui(64) - k))
  116. proc next*(r: var Rand): uint64 =
  117. ## Computes a random `uint64` number using the given state.
  118. ##
  119. ## **See also:**
  120. ## * `rand proc<#rand,Rand,Natural>`_ that returns an integer between zero and
  121. ## a given upper bound
  122. ## * `rand proc<#rand,Rand,range[]>`_ that returns a float
  123. ## * `rand proc<#rand,Rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  124. ## that accepts a slice
  125. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  126. ## * `skipRandomNumbers proc<#skipRandomNumbers,Rand>`_
  127. runnableExamples("-r:off"):
  128. var r = initRand(2019)
  129. assert r.next() == 13223559681708962501'u64 # implementation defined
  130. assert r.next() == 7229677234260823147'u64 # ditto
  131. let s0 = r.a0
  132. var s1 = r.a1
  133. result = s0 + s1
  134. s1 = s1 xor s0
  135. r.a0 = rotl(s0, 55) xor s1 xor (s1 shl 14) # a, b
  136. r.a1 = rotl(s1, 36) # c
  137. proc skipRandomNumbers*(s: var Rand) =
  138. ## The jump function for the generator.
  139. ##
  140. ## This proc is equivalent to `2^64` calls to `next <#next,Rand>`_, and it can
  141. ## be used to generate `2^64` non-overlapping subsequences for parallel
  142. ## computations.
  143. ##
  144. ## When multiple threads are generating random numbers, each thread must
  145. ## own the `Rand <#Rand>`_ state it is using so that the thread can safely
  146. ## obtain random numbers. However, if each thread creates its own Rand state,
  147. ## the subsequences of random numbers that each thread generates may overlap,
  148. ## even if the provided seeds are unique. This is more likely to happen as the
  149. ## number of threads and amount of random numbers generated increases.
  150. ##
  151. ## If many threads will generate random numbers concurrently, it is better to
  152. ## create a single Rand state and pass it to each thread. After passing the
  153. ## Rand state to a thread, call this proc before passing it to the next one.
  154. ## By using the Rand state this way, the subsequences of random numbers
  155. ## generated in each thread will never overlap as long as no thread generates
  156. ## more than `2^64` random numbers.
  157. ##
  158. ## **See also:**
  159. ## * `next proc<#next,Rand>`_
  160. runnableExamples("--threads:on"):
  161. import std/[random, threadpool]
  162. const spawns = 4
  163. const numbers = 100000
  164. proc randomSum(r: Rand): int =
  165. var r = r
  166. for i in 1..numbers:
  167. result += r.rand(0..10)
  168. var r = initRand(2019)
  169. var vals: array[spawns, FlowVar[int]]
  170. for val in vals.mitems:
  171. val = spawn randomSum(r)
  172. r.skipRandomNumbers()
  173. for val in vals:
  174. doAssert abs(^val - numbers * 5) / numbers < 0.1
  175. when defined(js):
  176. const helper = [0xbeac0467u32, 0xd86b048bu32]
  177. else:
  178. const helper = [0xbeac0467eba5facbu64, 0xd86b048b86aa9922u64]
  179. var
  180. s0 = Ui 0
  181. s1 = Ui 0
  182. for i in 0..high(helper):
  183. for b in 0 ..< 64:
  184. if (helper[i] and (Ui(1) shl Ui(b))) != 0:
  185. s0 = s0 xor s.a0
  186. s1 = s1 xor s.a1
  187. discard next(s)
  188. s.a0 = s0
  189. s.a1 = s1
  190. proc rand[T: uint | uint64](r: var Rand; max: T): T =
  191. # xxx export in future work
  192. if max == 0: return
  193. else:
  194. let max = uint64(max)
  195. when T.high.uint64 == uint64.high:
  196. if max == uint64.high: return T(next(r))
  197. while true:
  198. let x = next(r)
  199. # avoid `mod` bias
  200. if x <= randMax - (randMax mod max):
  201. return T(x mod (max + 1))
  202. proc rand*(r: var Rand; max: Natural): int {.benign.} =
  203. ## Returns a random integer in the range `0..max` using the given state.
  204. ##
  205. ## **See also:**
  206. ## * `rand proc<#rand,int>`_ that returns an integer using the default RNG
  207. ## * `rand proc<#rand,Rand,range[]>`_ that returns a float
  208. ## * `rand proc<#rand,Rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  209. ## that accepts a slice
  210. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  211. runnableExamples:
  212. var r = initRand(123)
  213. if false:
  214. assert r.rand(100) == 96 # implementation defined
  215. # bootstrap: can't use `runnableExamples("-r:off")`
  216. cast[int](rand(r, uint64(max)))
  217. # xxx toUnsigned pending https://github.com/nim-lang/Nim/pull/18445
  218. proc rand*(max: int): int {.benign.} =
  219. ## Returns a random integer in the range `0..max`.
  220. ##
  221. ## If `randomize <#randomize>`_ has not been called, the sequence of random
  222. ## numbers returned from this proc will always be the same.
  223. ##
  224. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  225. ##
  226. ## **See also:**
  227. ## * `rand proc<#rand,Rand,Natural>`_ that returns an integer using a
  228. ## provided state
  229. ## * `rand proc<#rand,float>`_ that returns a float
  230. ## * `rand proc<#rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  231. ## that accepts a slice
  232. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  233. runnableExamples("-r:off"):
  234. randomize(123)
  235. assert [rand(100), rand(100)] == [96, 63] # implementation defined
  236. rand(state, max)
  237. proc rand*(r: var Rand; max: range[0.0 .. high(float)]): float {.benign.} =
  238. ## Returns a random floating point number in the range `0.0..max`
  239. ## using the given state.
  240. ##
  241. ## **See also:**
  242. ## * `rand proc<#rand,float>`_ that returns a float using the default RNG
  243. ## * `rand proc<#rand,Rand,Natural>`_ that returns an integer
  244. ## * `rand proc<#rand,Rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  245. ## that accepts a slice
  246. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  247. runnableExamples:
  248. var r = initRand(234)
  249. let f = r.rand(1.0) # 8.717181376738381e-07
  250. let x = next(r)
  251. when defined(js):
  252. result = (float(x) / float(high(uint32))) * max
  253. else:
  254. let u = (0x3FFu64 shl 52u64) or (x shr 12u64)
  255. result = (cast[float](u) - 1.0) * max
  256. proc rand*(max: float): float {.benign.} =
  257. ## Returns a random floating point number in the range `0.0..max`.
  258. ##
  259. ## If `randomize <#randomize>`_ has not been called, the sequence of random
  260. ## numbers returned from this proc will always be the same.
  261. ##
  262. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  263. ##
  264. ## **See also:**
  265. ## * `rand proc<#rand,Rand,range[]>`_ that returns a float using a
  266. ## provided state
  267. ## * `rand proc<#rand,int>`_ that returns an integer
  268. ## * `rand proc<#rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  269. ## that accepts a slice
  270. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  271. runnableExamples:
  272. randomize(234)
  273. let f = rand(1.0) # 8.717181376738381e-07
  274. rand(state, max)
  275. proc rand*[T: Ordinal or SomeFloat](r: var Rand; x: HSlice[T, T]): T =
  276. ## For a slice `a..b`, returns a value in the range `a..b` using the given
  277. ## state.
  278. ##
  279. ## Allowed types for `T` are integers, floats, and enums without holes.
  280. ##
  281. ## **See also:**
  282. ## * `rand proc<#rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  283. ## that accepts a slice and uses the default RNG
  284. ## * `rand proc<#rand,Rand,Natural>`_ that returns an integer
  285. ## * `rand proc<#rand,Rand,range[]>`_ that returns a float
  286. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  287. runnableExamples:
  288. var r = initRand(345)
  289. assert r.rand(1..5) <= 5
  290. assert r.rand(-1.1 .. 1.2) >= -1.1
  291. assert x.a <= x.b
  292. when T is SomeFloat:
  293. result = rand(r, x.b - x.a) + x.a
  294. else: # Integers and Enum types
  295. when defined(js):
  296. result = cast[T](rand(r, cast[uint](x.b) - cast[uint](x.a)) + cast[uint](x.a))
  297. else:
  298. result = cast[T](rand(r, cast[uint64](x.b) - cast[uint64](x.a)) + cast[uint64](x.a))
  299. proc rand*[T: Ordinal or SomeFloat](x: HSlice[T, T]): T =
  300. ## For a slice `a..b`, returns a value in the range `a..b`.
  301. ##
  302. ## Allowed types for `T` are integers, floats, and enums without holes.
  303. ##
  304. ## If `randomize <#randomize>`_ has not been called, the sequence of random
  305. ## numbers returned from this proc will always be the same.
  306. ##
  307. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  308. ##
  309. ## **See also:**
  310. ## * `rand proc<#rand,Rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  311. ## that accepts a slice and uses a provided state
  312. ## * `rand proc<#rand,int>`_ that returns an integer
  313. ## * `rand proc<#rand,float>`_ that returns a floating point number
  314. ## * `rand proc<#rand,typedesc[T]>`_ that accepts an integer or range type
  315. runnableExamples:
  316. randomize(345)
  317. assert rand(1..6) <= 6
  318. result = rand(state, x)
  319. proc rand*[T: SomeInteger](t: typedesc[T]): T =
  320. ## Returns a random integer in the range `low(T)..high(T)`.
  321. ##
  322. ## If `randomize <#randomize>`_ has not been called, the sequence of random
  323. ## numbers returned from this proc will always be the same.
  324. ##
  325. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  326. ##
  327. ## **See also:**
  328. ## * `rand proc<#rand,int>`_ that returns an integer
  329. ## * `rand proc<#rand,float>`_ that returns a floating point number
  330. ## * `rand proc<#rand,HSlice[T: Ordinal or float or float32 or float64,T: Ordinal or float or float32 or float64]>`_
  331. ## that accepts a slice
  332. runnableExamples:
  333. randomize(567)
  334. if false: # implementation defined
  335. assert rand(int8) == -42
  336. assert rand(uint32) == 578980729'u32
  337. assert rand(range[1..16]) == 11
  338. # pending csources >= 1.4.0 or fixing https://github.com/timotheecour/Nim/issues/251#issuecomment-831599772,
  339. # use `runnableExamples("-r:off")` instead of `if false`
  340. when T is range:
  341. result = rand(state, low(T)..high(T))
  342. else:
  343. result = cast[T](state.next)
  344. proc sample*[T](r: var Rand; s: set[T]): T =
  345. ## Returns a random element from the set `s` using the given state.
  346. ##
  347. ## **See also:**
  348. ## * `sample proc<#sample,set[T]>`_ that uses the default RNG
  349. ## * `sample proc<#sample,Rand,openArray[T]>`_ for `openArray`s
  350. ## * `sample proc<#sample,Rand,openArray[T],openArray[U]>`_ that uses a
  351. ## cumulative distribution function
  352. runnableExamples:
  353. var r = initRand(987)
  354. let s = {1, 3, 5, 7, 9}
  355. assert r.sample(s) in s
  356. assert card(s) != 0
  357. var i = rand(r, card(s) - 1)
  358. for e in s:
  359. if i == 0: return e
  360. dec(i)
  361. proc sample*[T](s: set[T]): T =
  362. ## Returns a random element from the set `s`.
  363. ##
  364. ## If `randomize <#randomize>`_ has not been called, the order of outcomes
  365. ## from this proc will always be the same.
  366. ##
  367. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  368. ##
  369. ## **See also:**
  370. ## * `sample proc<#sample,Rand,set[T]>`_ that uses a provided state
  371. ## * `sample proc<#sample,openArray[T]>`_ for `openArray`s
  372. ## * `sample proc<#sample,openArray[T],openArray[U]>`_ that uses a
  373. ## cumulative distribution function
  374. runnableExamples:
  375. randomize(987)
  376. let s = {1, 3, 5, 7, 9}
  377. assert sample(s) in s
  378. sample(state, s)
  379. proc sample*[T](r: var Rand; a: openArray[T]): T =
  380. ## Returns a random element from `a` using the given state.
  381. ##
  382. ## **See also:**
  383. ## * `sample proc<#sample,openArray[T]>`_ that uses the default RNG
  384. ## * `sample proc<#sample,Rand,openArray[T],openArray[U]>`_ that uses a
  385. ## cumulative distribution function
  386. ## * `sample proc<#sample,Rand,set[T]>`_ for sets
  387. runnableExamples:
  388. let marbles = ["red", "blue", "green", "yellow", "purple"]
  389. var r = initRand(456)
  390. assert r.sample(marbles) in marbles
  391. result = a[r.rand(a.low..a.high)]
  392. proc sample*[T](a: openArray[T]): lent T =
  393. ## Returns a random element from `a`.
  394. ##
  395. ## If `randomize <#randomize>`_ has not been called, the order of outcomes
  396. ## from this proc will always be the same.
  397. ##
  398. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  399. ##
  400. ## **See also:**
  401. ## * `sample proc<#sample,Rand,openArray[T]>`_ that uses a provided state
  402. ## * `sample proc<#sample,openArray[T],openArray[U]>`_ that uses a
  403. ## cumulative distribution function
  404. ## * `sample proc<#sample,set[T]>`_ for sets
  405. runnableExamples:
  406. let marbles = ["red", "blue", "green", "yellow", "purple"]
  407. randomize(456)
  408. assert sample(marbles) in marbles
  409. result = a[rand(a.low..a.high)]
  410. proc sample*[T, U](r: var Rand; a: openArray[T]; cdf: openArray[U]): T =
  411. ## Returns an element from `a` using a cumulative distribution function
  412. ## (CDF) and the given state.
  413. ##
  414. ## The `cdf` argument does not have to be normalized, and it could contain
  415. ## any type of elements that can be converted to a `float`. It must be
  416. ## the same length as `a`. Each element in `cdf` should be greater than
  417. ## or equal to the previous element.
  418. ##
  419. ## The outcome of the `cumsum<math.html#cumsum,openArray[T]>`_ proc and the
  420. ## return value of the `cumsummed<math.html#cumsummed,openArray[T]>`_ proc,
  421. ## which are both in the math module, can be used as the `cdf` argument.
  422. ##
  423. ## **See also:**
  424. ## * `sample proc<#sample,openArray[T],openArray[U]>`_ that also utilizes
  425. ## a CDF but uses the default RNG
  426. ## * `sample proc<#sample,Rand,openArray[T]>`_ that does not use a CDF
  427. ## * `sample proc<#sample,Rand,set[T]>`_ for sets
  428. runnableExamples:
  429. from std/math import cumsummed
  430. let marbles = ["red", "blue", "green", "yellow", "purple"]
  431. let count = [1, 6, 8, 3, 4]
  432. let cdf = count.cumsummed
  433. var r = initRand(789)
  434. assert r.sample(marbles, cdf) in marbles
  435. assert(cdf.len == a.len) # Two basic sanity checks.
  436. assert(float(cdf[^1]) > 0.0)
  437. # While we could check cdf[i-1] <= cdf[i] for i in 1..cdf.len, that could get
  438. # awfully expensive even in debugging modes.
  439. let u = r.rand(float(cdf[^1]))
  440. a[cdf.upperBound(U(u))]
  441. proc sample*[T, U](a: openArray[T]; cdf: openArray[U]): T =
  442. ## Returns an element from `a` using a cumulative distribution function
  443. ## (CDF).
  444. ##
  445. ## This proc works similarly to
  446. ## `sample <#sample,Rand,openArray[T],openArray[U]>`_.
  447. ## See that proc's documentation for more details.
  448. ##
  449. ## If `randomize <#randomize>`_ has not been called, the order of outcomes
  450. ## from this proc will always be the same.
  451. ##
  452. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  453. ##
  454. ## **See also:**
  455. ## * `sample proc<#sample,Rand,openArray[T],openArray[U]>`_ that also utilizes
  456. ## a CDF but uses a provided state
  457. ## * `sample proc<#sample,openArray[T]>`_ that does not use a CDF
  458. ## * `sample proc<#sample,set[T]>`_ for sets
  459. runnableExamples:
  460. from std/math import cumsummed
  461. let marbles = ["red", "blue", "green", "yellow", "purple"]
  462. let count = [1, 6, 8, 3, 4]
  463. let cdf = count.cumsummed
  464. randomize(789)
  465. assert sample(marbles, cdf) in marbles
  466. state.sample(a, cdf)
  467. proc gauss*(r: var Rand; mu = 0.0; sigma = 1.0): float {.since: (1, 3).} =
  468. ## Returns a Gaussian random variate,
  469. ## with mean `mu` and standard deviation `sigma`
  470. ## using the given state.
  471. # Ratio of uniforms method for normal
  472. # https://www2.econ.osaka-u.ac.jp/~tanizaki/class/2013/econome3/13.pdf
  473. const K = sqrt(2 / E)
  474. var
  475. a = 0.0
  476. b = 0.0
  477. while true:
  478. a = rand(r, 1.0)
  479. b = (2.0 * rand(r, 1.0) - 1.0) * K
  480. if b * b <= -4.0 * a * a * ln(a): break
  481. result = mu + sigma * (b / a)
  482. proc gauss*(mu = 0.0, sigma = 1.0): float {.since: (1, 3).} =
  483. ## Returns a Gaussian random variate,
  484. ## with mean `mu` and standard deviation `sigma`.
  485. ##
  486. ## If `randomize <#randomize>`_ has not been called, the order of outcomes
  487. ## from this proc will always be the same.
  488. ##
  489. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  490. result = gauss(state, mu, sigma)
  491. proc initRand*(seed: int64): Rand =
  492. ## Initializes a new `Rand <#Rand>`_ state using the given seed.
  493. ##
  494. ## Providing a specific seed will produce the same results for that seed each time.
  495. ##
  496. ## The resulting state is independent of the default RNG's state. When `seed == 0`,
  497. ## we internally set the seed to an implementation defined non-zero value.
  498. ##
  499. ## **See also:**
  500. ## * `initRand proc<#initRand>`_ that uses the current time
  501. ## * `randomize proc<#randomize,int64>`_ that accepts a seed for the default RNG
  502. ## * `randomize proc<#randomize>`_ that initializes the default RNG using the current time
  503. runnableExamples:
  504. from std/times import getTime, toUnix, nanosecond
  505. var r1 = initRand(123)
  506. let now = getTime()
  507. var r2 = initRand(now.toUnix * 1_000_000_000 + now.nanosecond)
  508. const seedFallback0 = int32.high # arbitrary
  509. let seed = if seed != 0: seed else: seedFallback0 # because 0 is a fixed point
  510. result.a0 = Ui(seed shr 16)
  511. result.a1 = Ui(seed and 0xffff)
  512. when not defined(nimLegacyRandomInitRand):
  513. # calling `discard next(result)` (even a few times) would still produce
  514. # skewed numbers for the 1st call to `rand()`.
  515. skipRandomNumbers(result)
  516. discard next(result)
  517. proc randomize*(seed: int64) {.benign.} =
  518. ## Initializes the default random number generator with the given seed.
  519. ##
  520. ## Providing a specific seed will produce the same results for that seed each time.
  521. ##
  522. ## **See also:**
  523. ## * `initRand proc<#initRand,int64>`_ that initializes a Rand state
  524. ## with a given seed
  525. ## * `randomize proc<#randomize>`_ that uses the current time instead
  526. ## * `initRand proc<#initRand>`_ that initializes a Rand state using
  527. ## the current time
  528. runnableExamples:
  529. from std/times import getTime, toUnix, nanosecond
  530. randomize(123)
  531. let now = getTime()
  532. randomize(now.toUnix * 1_000_000_000 + now.nanosecond)
  533. state = initRand(seed)
  534. proc shuffle*[T](r: var Rand; x: var openArray[T]) =
  535. ## Shuffles a sequence of elements in-place using the given state.
  536. ##
  537. ## **See also:**
  538. ## * `shuffle proc<#shuffle,openArray[T]>`_ that uses the default RNG
  539. runnableExamples:
  540. var cards = ["Ace", "King", "Queen", "Jack", "Ten"]
  541. var r = initRand(678)
  542. r.shuffle(cards)
  543. import std/algorithm
  544. assert cards.sorted == @["Ace", "Jack", "King", "Queen", "Ten"]
  545. for i in countdown(x.high, 1):
  546. let j = r.rand(i)
  547. swap(x[i], x[j])
  548. proc shuffle*[T](x: var openArray[T]) =
  549. ## Shuffles a sequence of elements in-place.
  550. ##
  551. ## If `randomize <#randomize>`_ has not been called, the order of outcomes
  552. ## from this proc will always be the same.
  553. ##
  554. ## This proc uses the default RNG. Thus, it is **not** thread-safe.
  555. ##
  556. ## **See also:**
  557. ## * `shuffle proc<#shuffle,Rand,openArray[T]>`_ that uses a provided state
  558. runnableExamples:
  559. var cards = ["Ace", "King", "Queen", "Jack", "Ten"]
  560. randomize(678)
  561. shuffle(cards)
  562. import std/algorithm
  563. assert cards.sorted == @["Ace", "Jack", "King", "Queen", "Ten"]
  564. shuffle(state, x)
  565. when not defined(standalone):
  566. when defined(js):
  567. import std/times
  568. else:
  569. when defined(nimscript):
  570. import std/hashes
  571. else:
  572. import std/[hashes, os, sysrand, monotimes]
  573. when compileOption("threads"):
  574. import locks
  575. var baseSeedLock: Lock
  576. baseSeedLock.initLock
  577. var baseState: Rand
  578. proc initRand(): Rand =
  579. ## Initializes a new Rand state.
  580. ##
  581. ## The resulting state is independent of the default RNG's state.
  582. ##
  583. ## **Note:** Does not work for the compile-time VM.
  584. ##
  585. ## See also:
  586. ## * `initRand proc<#initRand,int64>`_ that accepts a seed for a new Rand state
  587. ## * `randomize proc<#randomize>`_ that initializes the default RNG using the current time
  588. ## * `randomize proc<#randomize,int64>`_ that accepts a seed for the default RNG
  589. when defined(js):
  590. let time = int64(times.epochTime() * 1000) and 0x7fff_ffff
  591. result = initRand(time)
  592. else:
  593. proc getRandomState(): Rand =
  594. when defined(nimscript):
  595. result = Rand(
  596. a0: CompileTime.hash.Ui,
  597. a1: CompileDate.hash.Ui)
  598. if not result.isValid:
  599. result = DefaultRandSeed
  600. else:
  601. var urand: array[sizeof(Rand), byte]
  602. for i in 0 .. 7:
  603. if sysrand.urandom(urand):
  604. copyMem(result.addr, urand[0].addr, sizeof(Rand))
  605. if result.isValid:
  606. break
  607. if not result.isValid:
  608. # Don't try to get alternative random values from other source like time or process/thread id,
  609. # because such code would be never tested and is a liability for security.
  610. quit("Failed to initializes baseState in random module as sysrand.urandom doesn't work.")
  611. when compileOption("threads"):
  612. baseSeedLock.withLock:
  613. if not baseState.isValid:
  614. baseState = getRandomState()
  615. result = baseState
  616. baseState.skipRandomNumbers
  617. else:
  618. if not baseState.isValid:
  619. baseState = getRandomState()
  620. result = baseState
  621. baseState.skipRandomNumbers
  622. since (1, 5, 1):
  623. export initRand
  624. proc randomize*() {.benign.} =
  625. ## Initializes the default random number generator with a seed based on
  626. ## random number source.
  627. ##
  628. ## This proc only needs to be called once, and it should be called before
  629. ## the first usage of procs from this module that use the default RNG.
  630. ##
  631. ## **Note:** Does not work for the compile-time VM.
  632. ##
  633. ## **See also:**
  634. ## * `randomize proc<#randomize,int64>`_ that accepts a seed
  635. ## * `initRand proc<#initRand>`_ that initializes a Rand state using
  636. ## the current time
  637. ## * `initRand proc<#initRand,int64>`_ that initializes a Rand state
  638. ## with a given seed
  639. state = initRand()
  640. {.pop.}