This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org

Python and Viv · Random numbers

Seeded random numbers

What to notice

Open the Viv in the workbench The atlas of the Viv py2viv translates this Python; the translation prints the same output.
Python seeded.py · 40 lines
Viv seeded.viv · 46 lines
py2viv translated from seeded.py
What py2viv adds before the program
(nothing)
(nothing)

1390 lines: a header and the guard functions its translation calls

1###
2Translated from seeded.py by py2viv 0.2.0 (numeric Python -> Viv).
3Python source sha256: e2f027bb7a2b45ceaa1f97d5e08abe2f0a04e1bf84e518108f4f9e523ed5b635
4Each print(...) is one print(...) line; fmt(x) writes a float as Python's repr does.
5A `raise`, a failed assert and every guard `stop` with Python's exception line.
6###
7
8#| 0, 1, ..., n - 1.
9def py_iota(n) {
10 a = zeros(n)
11 for i to n {
12 a[i] = i
13 }
14 return a
15}
16
17#| Every entry of a real array, row-major, as a vector (a copy).
18def py_flat(a: Arr) {
19 return take(a, py_iota(numel(a)))
20}
21
22#| np.arange(start, stop, step), as NumPy fills it: ceil((stop - start) / step) entries,
23#| entry k = start + k * delta with delta = (start + step) - start.
24def py_arange(start, stop_, step) {
25 if step == 0 {
26 stop("ZeroDivisionError: division by zero")
27 }
28 if not (abs((stop_ - start) / step) < 1.0 / 0.0) {
29 stop("py2viv: np.arange of a length that is not finite (NumPy raises ValueError)")
30 }
31 n = max(ceil((stop_ - start) / step), 0)
32 a = zeros(n)
33 delta = (start + step) - start
34 for k to n {
35 a[k] = start + k * delta
36 }
37 return a
38}
39
40#| np.mean(a): NumPy's pairwise sum over numel(a).
41def py_mean(a: Arr) {
42 return py_npsum(a) / numel(a)
43}
44
45#| np.var(a): the pairwise sum of the squared deviations over n, as NumPy's _var computes it.
46def py_var(a: Arr) {
47 d = a - py_mean(a)
48 return py_npsum(d * d) / numel(a)
49}
50
51#| np.std(a): divides by n, as NumPy does (lib/stats `std`).
52def py_std(a: Arr) {
53 return sqrt(py_var(a))
54}
55
56#| An int result, as Python's int: a -0.0 becomes 0 (an int has no sign on zero), and the run
57#| stops at 2^53, past which a double no longer holds every integer (Python's int does).
58def py_int(x) {
59 if abs(x) >= 9007199254740992 {
60 stop("py2viv: an integer reached 2^53; Python's int is exact there and a double is not")
61 }
62 return x + 0.0
63}
64
65#| a / b of two ints: ZeroDivisionError on a zero divisor.
66def py_div_i(a, b) {
67 if b == 0 {
68 stop("ZeroDivisionError: division by zero")
69 }
70 return a / b
71}
72
73#| NumPy's pairwise sum of f[lo], ..., f[lo + n - 1] (numpy/core/src/umath/loops_utils.h):
74#| below 8 terms left to right from 0; to 128, eight running sums combined as a tree, then
75#| the rest; above, the two halves (the first a multiple of 8 long) summed apart.
76def py_pw(f: Arr, lo, n) {
77 if n < 8 {
78 s = 0.0
79 for i to n {
80 s += f[lo + i]
81 }
82 return s
83 }
84 if n <= 128 {
85 r = zeros(8)
86 for j to 8 {
87 r[j] = f[lo + j]
88 }
89 i = 8
90 while i < n - n % 8 {
91 for j to 8 {
92 r[j] += f[lo + i + j]
93 }
94 i += 8
95 }
96 s = ((r[0] + r[1]) + (r[2] + r[3])) + ((r[4] + r[5]) + (r[6] + r[7]))
97 while i < n {
98 s += f[lo + i]
99 i += 1
100 }
101 return s
102 }
103 n2 = floor(n / 2)
104 n2 = n2 - n2 % 8
105 return py_pw(f, lo, n2) + py_pw(f, lo + n2, n - n2)
106}
107
108#| NumPy's sum of f[lo], ..., f[lo + n - 1] along one axis: from 0.0, a pairwise sum of each
109#| block of NumPy's buffer size (8192 entries), added in order.
110def py_blk(f: Arr, lo, n) {
111 bufsize = 8192
112 s = 0.0
113 b = 0
114 while b < n {
115 s = s + py_pw(f, lo + b, min(bufsize, n - b))
116 b += bufsize
117 }
118 return s
119}
120
121#| np.sum(a), a.sum() of a float array: NumPy's sum over the entries in row-major order (py_blk).
122def py_npsum(a: Arr) {
123 f = py_flat(a)
124 return py_blk(f, 0, len(f))
125}
126
127#| The sum of an integer array or list, exactly: every partial sum checked below 2^53.
128def py_isum(a: Arr) {
129 f = py_flat(a)
130 s = 0
131 for i to len(f) {
132 s = py_int(s + f[i])
133 }
134 return s
135}
136
137#| An array size, as NumPy takes it: ValueError for a negative one (Viv would make an empty
138#| array).
139def py_size(n) {
140 if n < 0 {
141 stop("ValueError: negative dimensions are not allowed")
142 }
143 return n
144}
145
146#| The shape NumPy broadcasts two shapes da and db to (aligned at the last axis, a length 1
147#| stretched); the run stops where NumPy raises ValueError.
148def py_bdims(da: Arr, db: Arr) {
149 n = len(da)
150 if len(db) > n {
151 n = len(db)
152 }
153 e = zeros(n)
154 for i to n {
155 x = 1
156 y = 1
157 if i >= n - len(da) {
158 x = da[i - n + len(da)]
159 }
160 if i >= n - len(db) {
161 y = db[i - n + len(db)]
162 }
163 if x == y or y == 1 {
164 e[i] = x
165 }
166 else {
167 if x == 1 {
168 e[i] = y
169 }
170 else {
171 stop("ValueError: operands could not be broadcast together with shapes")
172 }
173 }
174 }
175 return e
176}
177
178#| The flat indices, in row-major order of the shape e, of the entries of an array of shape d
179#| that broadcasting to e reads (a length 1 of d is read again along that axis of e).
180def py_bidx(d: Arr, e: Arr) {
181 n = len(e)
182 m = len(d)
183 st = zeros(n)
184 s = 1
185 for k to m {
186 j = m - 1 - k
187 if d[j] != 1 {
188 st[j + n - m] = s
189 }
190 s = s * d[j]
191 }
192 idx = zeros(1)
193 for i to n {
194 w = e[i]
195 nxt = zeros(len(idx) * w)
196 for p to len(idx) {
197 for c to w {
198 nxt[p * w + c] = idx[p] + c * st[i]
199 }
200 }
201 idx = nxt
202 }
203 return idx
204}
205
206#| a broadcast to the shape e (a shape py_bdims gave for a and another array): a itself when
207#| it has that shape, else a new array whose entries are copies of a's.
208def py_bto(a: Arr, e: Arr) {
209 d = shape(a)
210 if len(d) == len(e) {
211 if all(d == e) {
212 return a
213 }
214 }
215 return py_shaped(take(a, py_bidx(d, e)), e)
216}
217
218#| a + b of two real arrays, broadcast as NumPy broadcasts them.
219def py_badd(a: Arr, b: Arr) {
220 e = py_bdims(shape(a), shape(b))
221 return py_bto(a, e) + py_bto(b, e)
222}
223
224#| The flat indices of the nonzero entries of m, in order (a NaN is nonzero, -0.0 is zero).
225def py_nzflat(m: Arr) {
226 f = py_flat(m)
227 c = 0
228 for i to len(f) {
229 if f[i] != 0 {
230 c += 1
231 }
232 }
233 r = zeros(c)
234 c = 0
235 for i to len(f) {
236 if f[i] != 0 {
237 r[c] = i
238 c += 1
239 }
240 }
241 return r
242}
243
244#| The entries of f, row-major, as an array of shape e (a vector of 1 to 4 lengths).
245def py_shaped(f: Arr, e: Arr) {
246 if len(e) == 1 {
247 return f
248 }
249 if len(e) == 2 {
250 z = zeros(e[0], e[1])
251 }
252 else {
253 if len(e) == 3 {
254 z = zeros(e[0], e[1], e[2])
255 }
256 else {
257 z = zeros(e[0], e[1], e[2], e[3])
258 }
259 }
260 put(z, py_iota(numel(z)), f)
261 return z
262}
263
264NPR_TWO32 = 4294967296.0
265NPR_TWO53 = 9007199254740992.0
266NPR_U32_MAX = 4294967295
267NPR_TWO16 = 65536.0
268
269#| `a * b mod 2^32` for words `a` and `b`: a = hi * 2^16 + lo, and each partial
270#| product stays below 2^48.
271def npr_mul32(a, b) {
272 a_hi = floor(a / NPR_TWO16)
273 a_lo = a - a_hi * NPR_TWO16
274 return mod(a_lo * b + mod(a_hi * b, NPR_TWO16) * NPR_TWO16, NPR_TWO32)
275}
276
277#| `(a + b) mod 2^32`, and `(a - b) mod 2^32`, for words.
278def npr_add32(a, b) {
279 return mod(a + b, NPR_TWO32)
280}
281
282def npr_sub32(a, b) {
283 return mod(a - b + NPR_TWO32, NPR_TWO32)
284}
285
286#| The smallest `2^k - 1` at or above `r` (NumPy's gen_mask), for 0 <= r < 2^63.
287def npr_mask(r: Int) -> Int {
288 m: Int = 0
289 while m < r {
290 m = m * 2 + 1
291 }
292 return m
293}
294
295
298
299NPR_N = 624
300NPR_M = 397
301# 0x9908b0df, 0x80000000, 0x7fffffff
302NPR_MATRIX_A = 2567483615
303NPR_UPPER_MASK = 2147483648
304NPR_LOWER_MASK = 2147483647
305# the tempering masks 0x9d2c5680 and 0xefc60000, and the shifts
306NPR_TEMPER_B = 2636928640
307NPR_TEMPER_C = 4022730752
308NPR_TEMPER_U = 11
309NPR_TEMPER_S = 7
310NPR_TEMPER_T = 15
311NPR_TEMPER_L = 18
312# Knuth's multiplier (mt19937_seed, init_genrand), init_by_array's seed and multipliers
313NPR_SEED_MULT = 1812433253
314NPR_SEED_SHIFT = 30
315NPR_INIT_SEED = 19650218
316NPR_INIT_MULT1 = 1664525
317NPR_INIT_MULT2 = 1566083941
318# the double: a 27-bit and a 26-bit part
319NPR_DOUBLE_SHIFT_A = 5
320NPR_DOUBLE_SHIFT_B = 6
321NPR_TWO26 = 67108864.0
322# where the position and the seeded flag live in the state
323NPR_MT_POS = NPR_N
324NPR_MT_SEEDED = NPR_N + 1
325NPR_MT_LEN = NPR_N + 2
326
327#| A state nothing has seeded yet.
328def npr_mt_new() {
329 return zeros(NPR_MT_LEN)
330}
331
332#| mt19937_seed: np.random.seed(s) of an int s in [0, 2^32 - 1].
333def npr_mt_seed(st: Arr, seed) {
334 s = seed
335 for pos to NPR_N {
336 st[pos] = s
337 s = npr_add32(npr_mul32(NPR_SEED_MULT, bxor(s, shr(s, NPR_SEED_SHIFT))), pos + 1)
338 }
339 st[NPR_MT_POS] = NPR_N
340 st[NPR_MT_SEEDED] = 1
341}
342
343#| init_genrand: the start of init_by_array.
344def npr_mt_init_genrand(st: Arr, s) {
345 st[0] = s
346 for i from 1 to NPR_N {
347 st[i] = npr_add32(npr_mul32(NPR_SEED_MULT, bxor(st[i - 1], shr(st[i - 1], NPR_SEED_SHIFT))), i)
348 }
349 st[NPR_MT_POS] = NPR_N
350}
351
352#| mt19937_init_by_array: np.random.seed of an array of words.
353def npr_mt_seed_array(st: Arr, key: Arr) {
354 npr_mt_init_genrand(st, NPR_INIT_SEED)
355 key_length = len(key)
356 i = 1
357 j = 0
358 k = max(NPR_N, key_length)
359 while k > 0 {
360 prev = st[i - 1]
361 mixed = npr_mul32(bxor(prev, shr(prev, NPR_SEED_SHIFT)), NPR_INIT_MULT1)
362 st[i] = npr_add32(npr_add32(bxor(st[i], mixed), key[j]), j)
363 i = i + 1
364 j = j + 1
365 if i >= NPR_N {
366 st[0] = st[NPR_N - 1]
367 i = 1
368 }
369 if j >= key_length {
370 j = 0
371 }
372 k = k - 1
373 }
374 k = NPR_N - 1
375 while k > 0 {
376 prev = st[i - 1]
377 mixed = npr_mul32(bxor(prev, shr(prev, NPR_SEED_SHIFT)), NPR_INIT_MULT2)
378 st[i] = npr_sub32(bxor(st[i], mixed), i)
379 i = i + 1
380 if i >= NPR_N {
381 st[0] = st[NPR_N - 1]
382 i = 1
383 }
384 k = k - 1
385 }
386 st[0] = NPR_UPPER_MASK
387 st[NPR_MT_POS] = NPR_N
388 st[NPR_MT_SEEDED] = 1
389}
390
391#| mt19937_gen: all 624 words at once, in the reference's three runs. Word i
392#| takes the upper bit of key[i] and the lower 31 of key[i + 1], shifted right
393#| once, XORed with MATRIX_A when its lowest bit is set, and with key[i + M].
394def npr_mt_gen(st: Arr) {
395 for i to NPR_N - NPR_M {
396 y = bor(band(st[i], NPR_UPPER_MASK), band(st[i + 1], NPR_LOWER_MASK))
397 v = bxor(st[i + NPR_M], shr(y, 1))
398 if band(y, 1) == 1 {
399 v = bxor(v, NPR_MATRIX_A)
400 }
401 st[i] = v
402 }
403 for i from NPR_N - NPR_M to NPR_N - 1 {
404 y = bor(band(st[i], NPR_UPPER_MASK), band(st[i + 1], NPR_LOWER_MASK))
405 v = bxor(st[i + NPR_M - NPR_N], shr(y, 1))
406 if band(y, 1) == 1 {
407 v = bxor(v, NPR_MATRIX_A)
408 }
409 st[i] = v
410 }
411 y = bor(band(st[NPR_N - 1], NPR_UPPER_MASK), band(st[0], NPR_LOWER_MASK))
412 v = bxor(st[NPR_M - 1], shr(y, 1))
413 if band(y, 1) == 1 {
414 v = bxor(v, NPR_MATRIX_A)
415 }
416 st[NPR_N - 1] = v
417 st[NPR_MT_POS] = 0
418}
419
420#| mt19937_next: the next word, tempered.
421def npr_mt_next32(st: Arr) {
422 if st[NPR_MT_POS] == NPR_N {
423 npr_mt_gen(st)
424 }
425 pos = st[NPR_MT_POS]
426 y = st[pos]
427 st[NPR_MT_POS] = pos + 1
428 y = bxor(y, shr(y, NPR_TEMPER_U))
429 y = bxor(y, band(shl(y, NPR_TEMPER_S), NPR_TEMPER_B))
430 y = bxor(y, band(shl(y, NPR_TEMPER_T), NPR_TEMPER_C))
431 return bxor(y, shr(y, NPR_TEMPER_L))
432}
433
434#| mt19937_next_double: a double in [0, 1) from two words.
435def npr_mt_double(st: Arr) {
436 a = shr(npr_mt_next32(st), NPR_DOUBLE_SHIFT_A)
437 b = shr(npr_mt_next32(st), NPR_DOUBLE_SHIFT_B)
438 return (a * NPR_TWO26 + b) / NPR_TWO53
439}
440
441#| mt19937_next64 ANDed with a mask: the first word drawn is the high one.
442def npr_mt_masked64(st: Arr, mask: Int) -> Int {
443 hi: Int = int(npr_mt_next32(st))
444 lo: Int = int(npr_mt_next32(st))
445 return band(bor(shl(hi, 32), lo), mask)
446}
447
450
451#| random_uniform: low + rng * a double in [0, 1).
452def npr_lg_uniform(st: Arr, low, rng) {
453 return low + rng * npr_mt_double(st)
454}
455
456#| A masked 32-bit draw at or below rng (rng < 2^32 - 1, mask = gen_mask(rng)).
457def npr_lg_masked32(st: Arr, rng, mask: Int) {
458 v = band(int(npr_mt_next32(st)), mask)
459 while v > rng {
460 v = band(int(npr_mt_next32(st)), mask)
461 }
462 return float(v)
463}
464
465#| A masked 64-bit draw at or below rng (2^32 <= rng < 2^53).
466def npr_lg_masked64(st: Arr, rng, mask: Int) {
467 v = npr_mt_masked64(st, mask)
468 while v > rng {
469 v = npr_mt_masked64(st, mask)
470 }
471 return float(v)
472}
473
474#| random_bounded_uint64 with use_masked: an int in [0, rng], 0 <= rng < 2^53.
475def npr_lg_bounded(st: Arr, rng) {
476 if rng == 0 {
477 return 0
478 }
479 if rng == NPR_U32_MAX {
480 return npr_mt_next32(st)
481 }
482 if rng < NPR_U32_MAX {
483 return npr_lg_masked32(st, rng, npr_mask(int(rng)))
484 }
485 return npr_lg_masked64(st, rng, npr_mask(int(rng)))
486}
487
488#| random_interval: shuffle's draw of an int in [0, mx].
489def npr_lg_interval(st: Arr, mx) {
490 if mx == 0 {
491 return 0
492 }
493 if mx <= NPR_U32_MAX {
494 return npr_lg_masked32(st, mx, npr_mask(int(mx)))
495 }
496 return npr_lg_masked64(st, mx, npr_mask(int(mx)))
497}
498
499#| RandomState.shuffle of a vector, in place.
500def npr_lg_shuffle(st: Arr, a: Arr) {
501 n = len(a)
502 for k to n - 1 {
503 i = n - 1 - k
504 j = npr_lg_interval(st, i)
505 t = a[i]
506 a[i] = a[j]
507 a[j] = t
508 }
509}
510
511#| RandomState.shuffle of a matrix: its rows, in place.
512def npr_lg_shuffle_rows(st: Arr, A: Arr) {
513 n = shape(A)[0]
514 m = shape(A)[1]
515 for k to n - 1 {
516 i = n - 1 - k
517 j = npr_lg_interval(st, i)
518 if i != j {
519 for c to m {
520 t = A[i, c]
521 A[i, c] = A[j, c]
522 A[j, c] = t
523 }
524 }
525 }
526}
527
528#| kahan_sum (numpy/random/_common.pyx): the sum choice checks p with.
529def npr_kahan_sum(p: Arr) {
530 n = len(p)
531 if n == 0 {
532 return 0.0
533 }
534 s = p[0]
535 c = 0.0
536 for i from 1 to n {
537 y = p[i] - c
538 t = s + y
539 c = (t - s) - y
540 s = t
541 }
542 return s
543}
544
545#| choice's cdf: p.cumsum() (in order), then each entry divided by the last.
546def npr_lg_cdf(p: Arr) {
547 cdf = cumsum(p)
548 return cdf / cdf[len(p) - 1]
549}
550
551#| searchsorted(u, side='right') in a non-decreasing cdf: how many entries are <= u.
552def npr_lg_search(cdf: Arr, u) {
553 lo = 0
554 hi = len(cdf)
555 while lo < hi {
556 mid = floor((lo + hi) / 2)
557 if cdf[mid] <= u {
558 lo = mid + 1
559 }
560 else {
561 hi = mid
562 }
563 }
564 return lo
565}
566
569
570NPR_SS_POOL_SIZE = 4
571NPR_SS_XSHIFT = 16
572# 0x43b0d7e5, 0x931e8875, 0x8b51f9dd, 0x58f38ded, 0xca01f9dd, 0x4973f715
573NPR_SS_INIT_A = 1135663077
574NPR_SS_MULT_A = 2468251765
575NPR_SS_INIT_B = 2337405405
576NPR_SS_MULT_B = 1492356589
577NPR_SS_MIX_MULT_L = 3389127133
578NPR_SS_MIX_MULT_R = 1232336661
579
580#| The 32-bit words of n (0 <= n < 2^53), lowest first.
581def npr_ss_words(n) {
582 if n < NPR_TWO32 {
583 return [n]
584 }
585 return [mod(n, NPR_TWO32), floor(n / NPR_TWO32)]
586}
587
588#| hashmix: value ^ hash_const, times the next hash_const, xor-shifted. The hash
589#| constant lives in hc[0] and moves on with each call.
590def npr_ss_hashmix(value, hc: Arr) {
591 v = bxor(value, hc[0])
592 hc[0] = npr_mul32(hc[0], NPR_SS_MULT_A)
593 v = npr_mul32(v, hc[0])
594 return bxor(v, shr(v, NPR_SS_XSHIFT))
595}
596
597#| mix: MIX_MULT_L * x - MIX_MULT_R * y, xor-shifted.
598def npr_ss_mix(x, y) {
599 r = npr_sub32(npr_mul32(NPR_SS_MIX_MULT_L, x), npr_mul32(NPR_SS_MIX_MULT_R, y))
600 return bxor(r, shr(r, NPR_SS_XSHIFT))
601}
602
603#| mix_entropy: the pool of four words from the entropy words.
604def npr_ss_pool(entropy: Arr) {
605 pool = zeros(NPR_SS_POOL_SIZE)
606 hc = [NPR_SS_INIT_A]
607 n = len(entropy)
608 for i to NPR_SS_POOL_SIZE {
609 if i < n {
610 pool[i] = npr_ss_hashmix(entropy[i], hc)
611 }
612 else {
613 pool[i] = npr_ss_hashmix(0, hc)
614 }
615 }
616 for i_src to NPR_SS_POOL_SIZE {
617 for i_dst to NPR_SS_POOL_SIZE {
618 if i_src != i_dst {
619 pool[i_dst] = npr_ss_mix(pool[i_dst], npr_ss_hashmix(pool[i_src], hc))
620 }
621 }
622 }
623 for i_src from NPR_SS_POOL_SIZE to max(n, NPR_SS_POOL_SIZE) {
624 for i_dst to NPR_SS_POOL_SIZE {
625 pool[i_dst] = npr_ss_mix(pool[i_dst], npr_ss_hashmix(entropy[i_src], hc))
626 }
627 }
628 return pool
629}
630
631#| generate_state: n_words 32-bit words from the pool, cycling through it.
632def npr_ss_generate(pool: Arr, n_words) {
633 out = zeros(n_words)
634 hc = NPR_SS_INIT_B
635 for i to n_words {
636 v = bxor(pool[mod(i, NPR_SS_POOL_SIZE)], hc)
637 hc = npr_mul32(hc, NPR_SS_MULT_B)
638 v = npr_mul32(v, hc)
639 out[i] = bxor(v, shr(v, NPR_SS_XSHIFT))
640 }
641 return out
642}
643
646
647# PCG_DEFAULT_MULTIPLIER_128 = 0x2360ed051fc65da44385df649fccf645, its halves
648NPR_PCG_MULT_HI: Int = 2549297995355413924
649NPR_PCG_MULT_LO: Int = 4865540595714422341
650# the output rotation is the state's top 6 bits
651NPR_PCG_ROT_SHIFT: Int = 58
652NPR_BITS64: Int = 64
653NPR_HALF64: Int = 32
654NPR_LOW32: Int = 4294967295
655# uint64_to_double: the top 53 bits of the output over 2^53
656NPR_PCG_DOUBLE_SHIFT: Int = 11
657NPR_PCG_SEED_WORDS = 8
658# the state array: the state, the increment, the kept half, the last output, the seeded flag
659NPR_PCG_LO = 0
660NPR_PCG_HI = 1
661NPR_PCG_INC = 2
662NPR_PCG_HAS32 = 4
663NPR_PCG_KEPT32 = 5
664NPR_PCG_OUT = 6
665NPR_PCG_SEEDED = 7
666NPR_PCG_LEN = 8
667
668#| A generator nothing has seeded.
669def npr_pcg_new() -> [Int] {
670 return int_zeros(NPR_PCG_LEN)
671}
672
673#| state = state * MULT + inc (mod 2^128), in 64-bit halves: the low half of
674#| lo * MULT_LO, then its high half plus the cross products, then the
675#| increment with the carry out of the low half.
676def npr_pcg_step(st: [Int]) {
677 lo = st[NPR_PCG_LO]
678 hi = st[NPR_PCG_HI]
679 p_lo = wrap_mul(lo, NPR_PCG_MULT_LO)
680 p_hi = wrap_add(umulhi(lo, NPR_PCG_MULT_LO), wrap_add(wrap_mul(hi, NPR_PCG_MULT_LO), wrap_mul(lo, NPR_PCG_MULT_HI)))
681 n_lo = wrap_add(p_lo, st[NPR_PCG_INC])
682 carry: Int = 0
683 if ult(n_lo, p_lo) {
684 carry = 1
685 }
686 st[NPR_PCG_HI] = wrap_add(wrap_add(p_hi, st[NPR_PCG_INC + 1]), carry)
687 st[NPR_PCG_LO] = n_lo
688}
689
690#| pcg64_next64: one step, then the XSL-RR output: (hi ^ lo) rotated right by
691#| the state's top 6 bits (a shift by 64 is 0, as the rotation by 0 needs).
692def npr_pcg_next64(st: [Int]) -> Int {
693 npr_pcg_step(st)
694 x = bxor(st[NPR_PCG_HI], st[NPR_PCG_LO])
695 r = shr(st[NPR_PCG_HI], NPR_PCG_ROT_SHIFT)
696 out = bor(shr(x, r), shl(x, NPR_BITS64 - r))
697 st[NPR_PCG_OUT] = out
698 return out
699}
700
701#| pcg64_next32: the low half of a fresh output; the high half is kept for the next call.
702def npr_pcg_next32(st: [Int]) -> Int {
703 if st[NPR_PCG_HAS32] == 1 {
704 st[NPR_PCG_HAS32] = 0
705 return st[NPR_PCG_KEPT32]
706 }
707 out = npr_pcg_next64(st)
708 st[NPR_PCG_HAS32] = 1
709 st[NPR_PCG_KEPT32] = shr(out, NPR_HALF64)
710 return band(out, NPR_LOW32)
711}
712
713#| uint64_to_double(pcg64_next64): the output's top 53 bits over 2^53.
714def npr_pcg_double(st: [Int]) {
715 return float(shr(npr_pcg_next64(st), NPR_PCG_DOUBLE_SHIFT)) / NPR_TWO53
716}
717
718#| A 64-bit Int from two 32-bit words held in numbers: lo + 2^32 hi.
719#@ 199eb438
720def npr_word64(lo, hi) -> Int {
721 l: Int = int(lo)
722 h: Int = int(hi)
723 return bor(l, shl(h, NPR_HALF64))
724}
725
726#| PCG64.__init__: seed from generate_state(4, np.uint64), given as its eight
727#| 32-bit words w (the uint64 k is w[2k] + 2^32 w[2k + 1]): the state is
728#| u64[0] * 2^64 + u64[1] and the sequence u64[2] * 2^64 + u64[3]. Then
729#| pcg_setseq_128_srandom_r: inc = (seq << 1) | 1, a step, state += init, a step.
730def npr_pcg_seed(st: [Int], w: Arr) {
731 for i to NPR_PCG_LEN {
732 st[i] = 0
733 }
734 seq_lo = npr_word64(w[6], w[7])
735 seq_hi = npr_word64(w[4], w[5])
736 st[NPR_PCG_INC] = bor(shl(seq_lo, 1), 1)
737 st[NPR_PCG_INC + 1] = bor(shl(seq_hi, 1), shr(seq_lo, NPR_BITS64 - 1))
738 npr_pcg_step(st)
739 init_lo = npr_word64(w[2], w[3])
740 init_hi = npr_word64(w[0], w[1])
741 s_lo = wrap_add(st[NPR_PCG_LO], init_lo)
742 carry: Int = 0
743 if ult(s_lo, init_lo) {
744 carry = 1
745 }
746 st[NPR_PCG_HI] = wrap_add(wrap_add(st[NPR_PCG_HI], init_hi), carry)
747 st[NPR_PCG_LO] = s_lo
748 npr_pcg_step(st)
749 st[NPR_PCG_SEEDED] = 1
750}
751
752#| default_rng(seed) for the entropy words of the seed (SeedSequence(seed)).
753def npr_pcg_from_entropy(entropy: Arr) -> [Int] {
754 st = npr_pcg_new()
755 npr_pcg_seed(st, npr_ss_generate(npr_ss_pool(entropy), NPR_PCG_SEED_WORDS))
756 return st
757}
758
761
762NPR_ZIG_R = 3.654152885361009
763NPR_ZIG_INV_R = 0.2736612373297583
764
765npr_zig_ki = [
766 4208095142473578, 0, 3387314423973544, 3838760076542274,
767 4030768804392682, 4136731738896254, 4203757248105145, 4249917568205994,
768 4283617341590296, 4309289223136604, 4329489775174550, 4345795907393188,
769 4359232558744730, 4370494503737299, 4380069246215646, 4388308869042394,
770 4395473957549321, 4401761481783924, 4407323076021240, 4412277362218204,
771 4416718463613199, 4420722014516422, 4424349484777079, 4427651345409294,
772 4430669422005229, 4433438668975191, 4435988524278344, 4438343955930065,
773 4440526279077425, 4442553800234660, 4444442329865861, 4446205593658138,
774 4447855565093316, 4449402736340121, 4450856340408624, 4452224534496486,
775 4453514552210512, 4454732830656798, 4455885117109368, 4456976558985043,
776 4458011780094444, 4458994945550386, 4459929817254120, 4460819801517196,
777 4461667990089170, 4462477195632268, 4463249982500384, 4463988693531856,
778 4464695473445501, 4465372289331869, 4466020948651920, 4466643115089764,
779 4467240322552142, 4467813987562542, 4468365420260672, 4468895834186994,
780 4469406355006040, 4469898028300364, 4470371826548633, 4470828655385770,
781 4471269359229841, 4471694726349190, 4472105493433674, 4472502349725738,
782 4472885940759935, 4473256871753524, 4473615710685532, 4473962991097124,
783 4474299214642296, 4474624853414418, 4474940352071305, 4475246129778808,
784 4475542581990776, 4475830082081194, 4476108982842610, 4476379617863426,
785 4476642302795321, 4476897336520866, 4477145002230339, 4477385568415884,
786 4477619289790266, 4477846408136804, 4478067153096380, 4478281742896886,
787 4478490385029917, 4478693276879082, 4478890606303906, 4479082552182886,
788 4479269284918997, 4479450966910588, 4479627752990372, 4479799790834988,
789 4479967221347354, 4480130179013872, 4480288792238368, 4480443183654460,
790 4480593470417939, 4480739764480586, 4480882172846772, 4481020797814010,
791 4481155737198612, 4481287084547452, 4481414929336784, 4481539357158974,
792 4481660449897960, 4481778285894165, 4481892940099539, 4482004484223382,
793 4482112986869492, 4482218513665204, 4482321127382802, 4482420888053758,
794 4482517853076245, 4482612077316275, 4482703613202871, 4482792510817576,
795 4482878817978627, 4482962580320076, 4483043841366126, 4483122642600925,
796 4483199023534056, 4483273021761922, 4483344673025224, 4483414011262724,
797 4483481068661428, 4483545875703378, 4483608461209170, 4483668852378323,
798 4483727074826624, 4483783152620564, 4483837108308932, 4483888962951686,
799 4483938736146144, 4483986446050596, 4484032109405372, 4484075741551420,
800 4484117356446452, 4484156966678662, 4484194583478081, 4484230216725550,
801 4484263874959345, 4484295565379450, 4484325293849474, 4484353064896186,
802 4484378881706674, 4484402746123075, 4484424658634833, 4484444618368474,
803 4484462623074794, 4484478669113436, 4484492751434740, 4484504863558830,
804 4484514997551788, 4484523143998833, 4484529291974394, 4484533429008906,
805 4484535541052219, 4484535612433424, 4484533625816926, 4484529562154580,
806 4484523400633636, 4484515118620291, 4484504691598554, 4484492093104164,
807 4484477294653230, 4484460265665252, 4484440973380154, 4484419382768918,
808 4484395456437370, 4484369154522621, 4484340434581640, 4484309251471359,
809 4484275557219678, 4484239300886654, 4484200428415112, 4484158882469814,
810 4484114602264271, 4484067523374160, 4484017577536216, 4483964692431365,
811 4483908791450714, 4483849793442887, 4483787612441036, 4483722157367660,
812 4483653331715198, 4483581033200083, 4483505153387764, 4483425577285833,
813 4483342182902157, 4483254840764470, 4483163413397547, 4483067754753536,
814 4482967709590562, 4482863112794072, 4482753788634692, 4482639549955636,
815 4482520197281720, 4482395517841076, 4482265284489409, 4482129254525304,
816 4481987168383486, 4481838748191074, 4481683696169781, 4481521692864464,
817 4481352395175570, 4481175434169564, 4480990412637506, 4480796902367134,
818 4480594441088331, 4480382529045225, 4480160625140311, 4479928142586662,
819 4479684443993061, 4479428835793398, 4479160561915451, 4478878796564388,
820 4478582635972392, 4478271088936406, 4477943065929958, 4477597366530538,
821 4477232664848704, 4476847492576192, 4476440219183781, 4476009028690434,
822 4475551892286424, 4475066535915646, 4474550401693506, 4474000601739904,
823 4473413862618200, 4472786458058295, 4472114126959004, 4471391972746494,
824 4470614338917719, 4469774653883156, 4468865235838896, 4467877045039530,
825 4466799366045354, 4465619395558397, 4464321701199635, 4462887501169282,
826 4461293691124341, 4459511507635972, 4457504658253067, 4455226650325010,
827 4452616884242348, 4449594783440798, 4446050695647666, 4441831266659618,
828 4436714892174061, 4430368316897338, 4422264825074740, 4411517007702132,
829 4396496531309976, 4373832704204284, 4335125104963628, 4251099761679434
830]
831
832npr_zig_wi = [
833 8.683627060801306e-16, 4.779330175727737e-17, 6.354352417405262e-17,
834 7.454870481247696e-17, 8.3293668157931e-17, 9.068060405059482e-17,
835 9.714860076567762e-17, 1.0294750314241019e-16, 1.0823430288447684e-16,
836 1.131147019610903e-16, 1.176635945702292e-16, 1.2193617278714363e-16,
837 1.2597439914637093e-16, 1.2981099886264032e-16, 1.3347203736824123e-16,
838 1.3697864842571203e-16, 1.4034823001242382e-16, 1.4359529452056943e-16,
839 1.4673208742364422e-16, 1.4976904668391037e-16, 1.5271515003596198e-16,
840 1.5557818169460764e-16, 1.5836494009290885e-16, 1.6108140175274928e-16,
841 1.6373285203969853e-16, 1.6632399058420835e-16, 1.6885901708676596e-16,
842 1.713417017655966e-16, 1.737754436586486e-16, 1.7616331923000996e-16,
843 1.7850812316976727e-16, 1.8081240285799152e-16, 1.830784876482675e-16,
844 1.853085138861802e-16, 1.8750444639373882e-16, 1.896680970077476e-16,
845 1.918011406483862e-16, 1.9390512930625104e-16, 1.9598150426628824e-16,
846 1.9803160683128174e-16, 2.000566877627333e-16, 2.0205791562071654e-16,
847 2.0403638415480212e-16, 2.0599311887403706e-16, 2.079290829041402e-16,
848 2.0984518222370352e-16, 2.1174227035760342e-16, 2.1362115259449868e-16,
849 2.1548258978581458e-16, 2.1732730177564367e-16, 2.191559705042727e-16,
850 2.2096924282235318e-16, 2.2276773304789553e-16, 2.2455202529414355e-16,
851 2.263226755928568e-16, 2.280802138345017e-16, 2.2982514554424684e-16,
852 2.3155795351040804e-16, 2.3327909928004356e-16, 2.3498902453470955e-16,
853 2.3668815235791604e-16, 2.3837688840454243e-16, 2.4005562198135063e-16,
854 2.4172472704675025e-16, 2.433845631371103e-16, 2.4503547622614954e-16,
855 2.466777995232705e-16, 2.4831185421610877e-16, 2.4993795016204524e-16,
856 2.515563865329658e-16, 2.5316745241713583e-16, 2.547714273816944e-16,
857 2.563685819989397e-16, 2.579591783392867e-16, 2.5954347043351707e-16,
858 2.6112170470670194e-16, 2.6269412038597256e-16, 2.6426094988411895e-16,
859 2.658224191608307e-16, 2.6737874806323633e-16, 2.689301506472616e-16,
860 2.704768354811995e-16, 2.720190059327732e-16, 2.735568604408679e-16,
861 2.7509059277301666e-16, 2.7662039226963903e-16, 2.781464440759544e-16,
862 2.79668929362423e-16, 2.8118802553450207e-16, 2.827039064324479e-16,
863 2.842167425218406e-16, 2.8572670107546015e-16, 2.87233946347098e-16,
864 2.887386397378482e-16, 2.9024093995538423e-16, 2.9174100316669455e-16,
865 2.9323898314471816e-16, 2.947350314092935e-16, 2.9622929736280665e-16,
866 2.977219284209029e-16, 2.992130701386013e-16, 3.007028663321331e-16,
867 3.0219145919680615e-16, 3.036789894211802e-16, 3.051655962978219e-16,
868 3.0665141783089545e-16, 3.081365908408297e-16, 3.0962125106629225e-16,
869 3.111055332636893e-16, 3.125895713043999e-16, 3.140734982699446e-16,
870 3.1555744654528006e-16, 3.1704154791040285e-16, 3.1852593363044065e-16,
871 3.2001073454440114e-16, 3.214960811527447e-16, 3.2298210370394156e-16,
872 3.244689322801698e-16, 3.2595669688230784e-16, 3.2744552751437067e-16,
873 3.2893555426753697e-16, 3.3042690740391284e-16, 3.3191971744017523e-16,
874 3.3341411523123725e-16, 3.3491023205407785e-16, 3.364081996918765e-16,
875 3.37908150518595e-16, 3.394102175841489e-16, 3.409145347003126e-16,
876 3.424212365275018e-16, 3.4393045866258313e-16, 3.454423377278584e-16,
877 3.4695701146137835e-16, 3.4847461880874137e-16, 3.499953000165381e-16,
878 3.5151919672760744e-16, 3.53046452078274e-16, 3.5457721079774357e-16,
879 3.5611161930983884e-16, 3.5764982583726505e-16, 3.59191980508603e-16,
880 3.6073823546823514e-16, 3.6228874498941915e-16, 3.6384366559073444e-16,
881 3.65403156156137e-16, 3.669673780588701e-16, 3.685364952894914e-16,
882 3.7011067458828983e-16, 3.716900855823823e-16, 3.7327490092779435e-16,
883 3.7486529645684887e-16, 3.7646145133120287e-16, 3.7806354820089604e-16,
884 3.7967177336979443e-16, 3.8128631696783774e-16, 3.829073731305243e-16,
885 3.8453514018609596e-16, 3.8616982085091493e-16, 3.878116224335587e-16,
886 3.894607570481926e-16, 3.9111744183782054e-16, 3.9278189920805415e-16,
887 3.944543570720877e-16, 3.9613504910761354e-16, 3.9782421502646826e-16,
888 3.995221008578565e-16, 4.012289592460629e-16, 4.029450497636328e-16,
889 4.04670639241075e-16, 4.0640600211422504e-16, 4.0815142079049387e-16,
890 4.0990718603532664e-16, 4.1167359738030257e-16, 4.134509635544236e-16,
891 4.1523960294026883e-16, 4.170398440568316e-16, 4.1885202607101123e-16,
892 4.206764993399015e-16, 4.2251362598620494e-16, 4.243637805093078e-16,
893 4.262273504347798e-16, 4.2810473700531167e-16, 4.2999635591638323e-16,
894 4.3190263810026294e-16, 4.338240305622791e-16, 4.357609972736849e-16,
895 4.3771402012585875e-16, 4.3968359995105214e-16, 4.4167025761542035e-16,
896 4.4367453519065673e-16, 4.456969972112043e-16, 4.477382320247534e-16,
897 4.49798853244555e-16, 4.518795013130059e-16, 4.539808451870034e-16,
898 4.561035841567422e-16, 4.582484498109567e-16, 4.604162081631153e-16,
899 4.626076619547846e-16, 4.648236531543207e-16, 4.670650656712631e-16,
900 4.693328283093329e-16, 4.716279179838351e-16, 4.739513632325867e-16,
901 4.763042480533137e-16, 4.786877161048723e-16, 4.811029753147417e-16,
902 4.835513029411525e-16, 4.860340511450812e-16, 4.885526531353603e-16,
903 4.91108629959527e-16, 4.937035980240335e-16, 4.963392774403987e-16,
904 4.990175013091822e-16, 5.017402260718089e-16, 5.045095430818727e-16,
905 5.073276915733542e-16, 5.101970732341562e-16, 5.131202686306784e-16,
906 5.161000557743228e-16, 5.191394311757699e-16, 5.222416338000234e-16,
907 5.254101724177597e-16, 5.286488569504945e-16, 5.3196183453384e-16,
908 5.353536311816497e-16, 5.388292001334053e-16, 5.423939782201712e-16,
909 5.46053951907478e-16, 5.498157350892814e-16, 5.536866612467876e-16,
910 5.576748932926576e-16, 5.617895553555417e-16, 5.660408920082422e-16,
911 5.704404621291389e-16, 5.750013768919895e-16, 5.797385945724594e-16,
912 5.846692893455479e-16, 5.898133176477899e-16, 5.951938149641444e-16,
913 6.008379696271908e-16, 6.067780409333449e-16, 6.130527208725282e-16,
914 6.197089894581626e-16, 6.268046963301284e-16, 6.344122407127506e-16,
915 6.426239659548055e-16, 6.515603317344994e-16, 6.613827885097664e-16,
916 6.723150462505587e-16, 6.846803417564259e-16, 6.98971833638762e-16,
917 7.159994934830664e-16, 7.372424301798799e-16, 7.658936370805573e-16,
918 8.113849337656484e-16
919]
920
921npr_zig_fi = [
922 1.0, 0.9771017012676716, 0.9598790918001067,
923 0.9451989534422996, 0.9320600759592305, 0.919991505039347,
924 0.9087264400521309, 0.8980959218983434, 0.8879846607558334,
925 0.8783096558089174, 0.869008688036857, 0.8600336211963315,
926 0.851346258458678, 0.8429156531122042, 0.8347162929868834,
927 0.8267268339462214, 0.8189291916037024, 0.8113078743126563,
928 0.8038494831709643, 0.796542330422959, 0.7893761435660246,
929 0.7823418326548025, 0.7754313049811872, 0.7686373157984863,
930 0.7619533468367954, 0.7553735065070961, 0.7488924472191568,
931 0.742505296340151, 0.7362075981268627, 0.7299952645614762,
932 0.7238645334686302, 0.717811932630722, 0.7118342488782484,
933 0.7059285013327543, 0.7000919181365116, 0.6943219161261167,
934 0.6886160830046718, 0.6829721616449949, 0.6773880362187735,
935 0.6718617198970821, 0.6663913439087501, 0.6609751477766631,
936 0.6556114705796973, 0.6502987431108167, 0.6450354808208223,
937 0.6398202774530566, 0.6346517992876236, 0.6295287799248367,
938 0.6244500155470265, 0.6194143606058343, 0.6144207238889139,
939 0.6094680649257734, 0.6045553906974678, 0.5996817526191253,
940 0.5948462437679874, 0.590047996332826, 0.5852861792633715,
941 0.5805599961007909, 0.5758686829723537, 0.5712115067352532,
942 0.5665877632561644, 0.5619967758145243, 0.557437893618766,
943 0.5529104904258323, 0.5484139632552658, 0.5439477311900263,
944 0.5395112342569521, 0.5351039323804576, 0.5307253044036621,
945 0.5263748471716845, 0.5220520746723218, 0.5177565172297564,
946 0.513487720747327, 0.5092452459957479, 0.5050286679434681,
947 0.5008375751261487, 0.4966715690524897, 0.49253026364386854,
948 0.48841328470545803, 0.4843202694266833, 0.48025086590904675,
949 0.47620473271950586, 0.4721815384677302, 0.4681809614056936,
950 0.46420268904817436, 0.46024641781284287, 0.45631185267871643,
951 0.4523987068618485, 0.44850670150720306, 0.4446355653957394,
952 0.440785034665804, 0.43695485254798555, 0.43314476911265226,
953 0.4293545410294414, 0.42558393133802197, 0.4218327092294959,
954 0.4181006498378482, 0.4143875340408911, 0.41069314827018816,
955 0.40701728432947337, 0.4033597392211145, 0.3997203149801972,
956 0.39609881851583245, 0.3924950614593156, 0.3889088600187887,
957 0.3853400348400773, 0.38178841087339366, 0.3782538172456192,
958 0.37473608713789114, 0.3712350576682395, 0.3677505697790326,
959 0.36428246812900406, 0.36083060098964803, 0.3573948201457805,
960 0.3539749808000768, 0.3505709414814061, 0.34718256395679364,
961 0.3438097131468507, 0.34045225704452187, 0.33711006663700605,
962 0.33378301583071845, 0.3304709813791636, 0.3271738428136014,
963 0.3238914823763911, 0.32062378495690536, 0.3173706380299136,
964 0.3141319315963372, 0.3109075581262865, 0.30769741250429206,
965 0.30450139197665, 0.30131939610080305, 0.2981513266966855,
966 0.2949970877999618, 0.2918565856170952, 0.2887297284821829,
967 0.28561642681550176, 0.2825165930837076, 0.27943014176163794,
968 0.2763569892956683, 0.27329705406857707, 0.27025025636587546,
969 0.26721651834356147, 0.2641957639972612, 0.2611879191327212,
970 0.25819291133761924, 0.25521066995466196, 0.2522411260559422,
971 0.24928421241852852, 0.24633986350126383, 0.2434080154227503,
972 0.2404886059405006, 0.2375815744312381, 0.23468686187233,
973 0.23180441082433872, 0.22893416541468034, 0.22607607132238028,
974 0.22323007576391748, 0.220396127480152, 0.21757417672433113,
975 0.21476417525117358, 0.21196607630703018, 0.20917983462112508,
976 0.2064054063978808, 0.2036427493103349, 0.2008918224946566,
977 0.19815258654577514, 0.1954250035141343, 0.19270903690358918,
978 0.19000465167046499, 0.1873118142238003, 0.18463049242679927,
979 0.18196065559952251, 0.17930227452284758, 0.17665532144373486,
980 0.17401977008183855, 0.17139559563750575, 0.1687827748012113,
981 0.1661812857644819, 0.16359110823236558, 0.161012223437511,
982 0.15844461415592428, 0.1558882647244792, 0.15334316106026286,
983 0.15080929068184568, 0.14828664273257455, 0.14577520800599403,
984 0.14327497897351346, 0.1407859498144447, 0.13830811644855073,
985 0.13584147657125376, 0.13338602969166916, 0.13094177717364436,
986 0.12850872227999957, 0.1260868702201859, 0.12367622820159657,
987 0.1212768054847903, 0.11888861344291006, 0.11651166562561087,
988 0.11414597782783849, 0.11179156816383809, 0.1094484571468118,
989 0.1071166677746838, 0.10479622562248707, 0.10248715894193525,
990 0.10018949876881002, 0.09790327903886246, 0.095628536713009,
991 0.09336531191269101, 0.09111364806637376, 0.08887359206827589,
992 0.08664519445055807, 0.08442850957035347, 0.0822235958132029,
993 0.08003051581466307, 0.07784933670209612, 0.07568013035892718,
994 0.07352297371398132, 0.0713779490588904, 0.06924514439700676,
995 0.0671246538277885, 0.0650165779712429, 0.06292102443775814,
996 0.06083810834953988, 0.05876795292093374, 0.0567106901062029,
997 0.05466646132488892, 0.05263541827679219, 0.05061772386094778,
998 0.04861355321586854, 0.04662309490193038, 0.044646552251294463,
999 0.04268414491647446, 0.04073611065594094, 0.03880270740452615,
1000 0.036884215688567305, 0.034980941461716125, 0.03309321945857858,
1001 0.0312214171919203, 0.02936593975813336, 0.027527235669603113,
1002 0.02570580400854891, 0.02390220330579588, 0.02211706270730885,
1003 0.02035109623004451, 0.018605121275724622, 0.016880083152543142,
1004 0.01517708830793531, 0.013497450601739867, 0.011842757857907879,
1005 0.010214971439701459, 0.008616582769398726, 0.007050875471373222,
1006 0.0055224032992509916, 0.0040379725933630236, 0.0026090727461021593,
1007 0.001260285930498598
1008]
1009
1010# a strip index is the low 8 bits of the 64-bit draw, the sign the next bit,
1011# and |x| the 52 bits above that
1012NPR_ZIG_IDX_MASK: Int = 255
1013NPR_ZIG_SIGN_SHIFT: Int = 8
1014NPR_ZIG_ABS_SHIFT: Int = 9
1015NPR_ZIG_ABS_MASK: Int = 4503599627370495
1016# the wedge test's margin, in units of 2^-52 of exp(-x^2/2)
1017NPR_ZIG_WEDGE_ULPS = 8
1018NPR_EPS = 2.220446049250313e-16
1019# choice(replace=False): the partial shuffle's thresholds
1020NPR_CHOICE_POP_MIN = 10000
1021NPR_CHOICE_CUTOFF_SHUFFLE = 50
1022NPR_CHOICE_CUTOFF_NOSHUFFLE = 20
1023# 2^32, and 2^62 (2^64 mod r is 4 (2^62 mod r) mod r, with every value an Int)
1024NPR_TWO32_INT: Int = 4294967296
1025NPR_TWO62: Int = 4611686018427387904
1026NPR_FOUR: Int = 4
1027
1028#| random_uniform: low + rng * a double in [0, 1).
1029def npr_gn_uniform(st: [Int], low, rng) {
1030 return low + rng * npr_pcg_double(st)
1031}
1032
1033#| buffered_bounded_lemire_uint32 for rng < 2^32 - 1: m = next32 * (rng + 1), a
1034#| 64-bit product, its high word the draw, its low word checked for bias.
1035def npr_gn_lemire32(st: [Int], rng) {
1036 r: Int = int(rng) + 1
1037 m = wrap_mul(npr_pcg_next32(st), r)
1038 leftover = band(m, NPR_LOW32)
1039 if leftover < r {
1040 threshold = (NPR_TWO32_INT - r) % r
1041 while leftover < threshold {
1042 m = wrap_mul(npr_pcg_next32(st), r)
1043 leftover = band(m, NPR_LOW32)
1044 }
1045 }
1046 return float(shr(m, NPR_HALF64))
1047}
1048
1049#| bounded_lemire_uint64 for 2^32 - 1 < rng < 2^53: the 128-bit product of the
1050#| 64-bit draw and rng + 1; its high half is the draw, its low half (as an
1051#| unsigned number) is checked against 2^64 mod (rng + 1).
1052def npr_gn_lemire64(st: [Int], rng) {
1053 r: Int = int(rng) + 1
1054 x = npr_pcg_next64(st)
1055 leftover = wrap_mul(x, r)
1056 if ult(leftover, r) {
1057 threshold = (NPR_TWO62 % r) * NPR_FOUR % r
1058 while ult(leftover, threshold) {
1059 x = npr_pcg_next64(st)
1060 leftover = wrap_mul(x, r)
1061 }
1062 }
1063 return float(umulhi(x, r))
1064}
1065
1066#| random_bounded_uint64 without the mask: an int in [0, rng], 0 <= rng < 2^53.
1067def npr_gn_bounded(st: [Int], rng) {
1068 if rng == 0 {
1069 return 0
1070 }
1071 if rng == NPR_U32_MAX {
1072 return float(npr_pcg_next32(st))
1073 }
1074 if rng < NPR_U32_MAX {
1075 return npr_gn_lemire32(st, rng)
1076 }
1077 return npr_gn_lemire64(st, rng)
1078}
1079
1080#| random_interval: an int in [0, mx] by masked rejection (shuffle's draw).
1081def npr_gn_interval(st: [Int], mx) {
1082 if mx == 0 {
1083 return 0
1084 }
1085 mask = npr_mask(int(mx))
1086 if mx <= NPR_U32_MAX {
1087 v = band(npr_pcg_next32(st), mask)
1088 while v > mx {
1089 v = band(npr_pcg_next32(st), mask)
1090 }
1091 return float(v)
1092 }
1093 v = band(npr_pcg_next64(st), mask)
1094 while v > mx {
1095 v = band(npr_pcg_next64(st), mask)
1096 }
1097 return float(v)
1098}
1099
1100#| Generator.shuffle of a vector, in place.
1101def npr_gn_shuffle(st: [Int], a: Arr) {
1102 n = len(a)
1103 for k to n - 1 {
1104 i = n - 1 - k
1105 j = npr_gn_interval(st, i)
1106 t = a[i]
1107 a[i] = a[j]
1108 a[j] = t
1109 }
1110}
1111
1112#| Generator.shuffle of a matrix (axis 0): its rows, in place.
1113def npr_gn_shuffle_rows(st: [Int], A: Arr) {
1114 n = shape(A)[0]
1115 m = shape(A)[1]
1116 for k to n - 1 {
1117 i = n - 1 - k
1118 j = npr_gn_interval(st, i)
1119 if i != j {
1120 for c to m {
1121 t = A[i, c]
1122 A[i, c] = A[j, c]
1123 A[j, c] = t
1124 }
1125 }
1126 }
1127}
1128
1129#| _shuffle_int: swap a[i] with a[random_bounded_uint64(0, i)] for i = n - 1 down to first.
1130def npr_gn_shuffle_int(st: [Int], a: Arr, n, first) {
1131 i = n - 1
1132 while i >= first {
1133 j = npr_gn_bounded(st, i)
1134 t = a[j]
1135 a[j] = a[i]
1136 a[i] = t
1137 i = i - 1
1138 }
1139}
1140
1141#| Generator.choice(pop, size, replace=False) without p: `size` distinct ints in [0, pop).
1142def npr_gn_choice_noreplace(st: [Int], pop, size, do_shuffle) {
1143 cutoff = NPR_CHOICE_CUTOFF_NOSHUFFLE
1144 if do_shuffle == 1 {
1145 cutoff = NPR_CHOICE_CUTOFF_SHUFFLE
1146 }
1147 if pop > NPR_CHOICE_POP_MIN and size > floor(pop / cutoff) {
1148 idx = zeros(pop)
1149 for i to pop {
1150 idx[i] = i
1151 }
1152 npr_gn_shuffle_int(st, idx, pop, max(pop - size, 1))
1153 out = zeros(size)
1154 for i to size {
1155 out[i] = idx[pop - size + i]
1156 }
1157 return out
1158 }
1159 out = zeros(size)
1160 taken = zeros(pop)
1161 for j from pop - size to pop {
1162 val = npr_gn_bounded(st, j)
1163 if taken[val] == 0 {
1164 taken[val] = 1
1165 out[j - pop + size] = val
1166 }
1167 else {
1168 taken[j] = 1
1169 out[j - pop + size] = j
1170 }
1171 }
1172 if do_shuffle == 1 {
1173 npr_gn_shuffle_int(st, out, size, 1)
1174 }
1175 return out
1176}
1177
1178#| random_standard_normal: the ziggurat; NaN where NumPy's value rests on the
1179#| platform's libm (see the header).
1180def npr_gn_normal(st: [Int]) {
1181 while true {
1182 r = npr_pcg_next64(st)
1183 idx = band(r, NPR_ZIG_IDX_MASK)
1184 sign = band(shr(r, NPR_ZIG_SIGN_SHIFT), 1)
1185 rabs = float(band(shr(r, NPR_ZIG_ABS_SHIFT), NPR_ZIG_ABS_MASK))
1186 x = rabs * npr_zig_wi[idx]
1187 if sign == 1 {
1188 x = -x
1189 }
1190 if rabs < npr_zig_ki[idx] {
1191 return x
1192 }
1193 if idx == 0 {
1194 # NumPy draws the tail with log1p: not reproducible here
1195 return 0.0 / 0.0
1196 }
1197 lhs = (npr_zig_fi[idx - 1] - npr_zig_fi[idx]) * npr_pcg_double(st) + npr_zig_fi[idx]
1198 e = exp(-0.5 * x * x)
1199 if abs(lhs - e) <= NPR_ZIG_WEDGE_ULPS * NPR_EPS * e {
1200 return 0.0 / 0.0
1201 }
1202 if lhs < e {
1203 return x
1204 }
1205 }
1206 return 0.0
1207}
1208
1211
1212#| The sign bit of x (a -0.0 has it), false for NaN: NumPy's check of a non-negative parameter.
1213def py_npr_signbit(x) {
1214 return x < 0.0 or x == 0.0 and 1.0 / x < 0.0
1215}
1216
1217#| np.random.default_rng(seed) of an int: SeedSequence(seed), then PCG64.
1218def py_npg_new(seed) {
1219 if seed < 0 {
1220 stop("ValueError: expected non-negative integer")
1221 }
1222 return npr_pcg_from_entropy(npr_ss_words(seed))
1223}
1224
1225#| default_rng of an array of ints: each entry's words in turn (_coerce_to_uint32_array).
1226def py_npg_new_arr(seeds: Arr) {
1227 words = []
1228 for i to len(seeds) {
1229 if seeds[i] < 0 {
1230 stop("ValueError: expected non-negative integer")
1231 }
1232 w = npr_ss_words(seeds[i])
1233 for k to len(w) {
1234 words = append(words, w[k])
1235 }
1236 }
1237 return npr_pcg_from_entropy(words)
1238}
1239
1240#| Generator.random(): one double; _v and _m fill a vector or matrix row by row.
1241def py_npg_random_v(st: [Int], n) {
1242 a = zeros(py_size(n))
1243 for i to n {
1244 a[i] = npr_pcg_double(st)
1245 }
1246 return a
1247}
1248
1249def py_npg_random_m(st: [Int], rows, cols) {
1250 A = zeros(py_size(rows), py_size(cols))
1251 for i to rows {
1252 for j to cols {
1253 A[i, j] = npr_pcg_double(st)
1254 }
1255 }
1256 return A
1257}
1258
1259#| Generator.standard_normal() and normal(loc, scale) = loc + scale * z. The ziggurat returns
1260#| NaN where NumPy's value depends on the platform's libm (the tail's log1p, a near-tie of the
1261#| wedge test's exp); the run stops there. scale must not carry a sign bit (NumPy's "scale <
1262#| 0", checked before the size).
1263def py_npg_z(st: [Int]) {
1264 v = npr_gn_normal(st)
1265 if v != v {
1266 stop("py2viv: a normal draw reached the ziggurat's tail (or a near-tie in its wedge test), where NumPy's value depends on the platform's log1p (or exp), which Viv's libm does not reproduce bit for bit")
1267 }
1268 return v
1269}
1270
1271def py_npg_scale(scale) {
1272 if py_npr_signbit(scale) {
1273 stop("ValueError: scale < 0")
1274 }
1275 return scale
1276}
1277
1278def py_npg_normal(st: [Int], loc, scale) {
1279 s = py_npg_scale(scale)
1280 return loc + s * py_npg_z(st)
1281}
1282
1283def py_npg_normal_v(st: [Int], loc, scale, n) {
1284 s = py_npg_scale(scale)
1285 a = zeros(py_size(n))
1286 for i to n {
1287 a[i] = loc + s * py_npg_z(st)
1288 }
1289 return a
1290}
1291
1292def py_npg_normal_m(st: [Int], loc, scale, rows, cols) {
1293 s = py_npg_scale(scale)
1294 A = zeros(py_size(rows), py_size(cols))
1295 for i to rows {
1296 for j to cols {
1297 A[i, j] = loc + s * py_npg_z(st)
1298 }
1299 }
1300 return A
1301}
1302
1303#| Generator.integers(low, high, endpoint): low + a Lemire draw in [0, high - low] (high - 1
1304#| without the endpoint). An empty size gives an empty array before the bounds are checked.
1305def py_npg_irange(low, high, closed) {
1306 hi = high
1307 if closed == 0 {
1308 hi = high - 1
1309 }
1310 if low > hi {
1311 if closed == 1 {
1312 if low == 0 {
1313 stop("ValueError: high < 0")
1314 }
1315 stop("ValueError: low > high")
1316 }
1317 if low == 0 {
1318 stop("ValueError: high <= 0")
1319 }
1320 stop("ValueError: low >= high")
1321 }
1322 r = hi - low
1323 if r >= NPR_TWO53 {
1324 stop("py2viv: a random int range of 2^53 or more; a double holds the ints exactly only below it")
1325 }
1326 return r
1327}
1328
1329def py_npg_integers(st: [Int], low, high, closed) {
1330 r = py_npg_irange(low, high, closed)
1331 return low + npr_gn_bounded(st, r)
1332}
1333
1334def py_npg_integers_v(st: [Int], low, high, n, closed) {
1335 if n == 0 {
1336 return zeros(0)
1337 }
1338 r = py_npg_irange(low, high, closed)
1339 a = zeros(py_size(n))
1340 for i to n {
1341 a[i] = low + npr_gn_bounded(st, r)
1342 }
1343 return a
1344}
1345
1346def py_npg_integers_m(st: [Int], low, high, rows, cols, closed) {
1347 if rows * cols == 0 {
1348 return zeros(py_size(rows), py_size(cols))
1349 }
1350 r = py_npg_irange(low, high, closed)
1351 A = zeros(py_size(rows), py_size(cols))
1352 for i to rows {
1353 for j to cols {
1354 A[i, j] = low + npr_gn_bounded(st, r)
1355 }
1356 }
1357 return A
1358}
1359
1360#| Generator.shuffle(x) of a vector (or a Python list), and of a matrix's rows, in place; and
1361#| Generator.permutation of an int, a vector and a matrix.
1362def py_npg_shuffle(st: [Int], a: Arr) {
1363 if len(a) > 1 {
1364 npr_gn_shuffle(st, a)
1365 }
1366 return 0
1367}
1368
1369def py_npg_shuffle_rows(st: [Int], A: Arr) {
1370 if numel(A) > 0 and shape(A)[0] > 1 {
1371 npr_gn_shuffle_rows(st, A)
1372 }
1373 return 0
1374}
1375
1376def py_npg_perm_n(st: [Int], n) {
1377 a = zeros(max(n, 0))
1378 for i to len(a) {
1379 a[i] = i
1380 }
1381 py_npg_shuffle(st, a)
1382 return a
1383}
1384
1385def py_npg_perm_v(st: [Int], x: Arr) {
1386 a = x * 1.0
1387 py_npg_shuffle(st, a)
1388 return a
1389}
1390
1391def py_npg_perm_m(st: [Int], X: Arr) {
1392 rows = shape(X)[0]
1393 cols = shape(X)[1]
1394 idx = py_npg_perm_n(st, rows)
1395 A = zeros(rows, cols)
1396 for i to rows {
1397 for j to cols {
1398 A[i, j] = X[idx[i], j]
1399 }
1400 }
1401 return A
1402}
1 imports
4import numpy as np
2from nprandom.generator import npr_gn_bounded, npr_gn_normal, npr_gn_shuffle
3from nprandom.pcg64 import npr_pcg_double, npr_pcg_from_entropy
(nothing)
2 a generator from a seed
7SEED = 42
8rng = np.random.default_rng(SEED)
6SEED = 42
7rng = npr_pcg_from_entropy([SEED])
1408SEED = 42
1409rng = py_npg_new(SEED)
3 uniform draws in [0, 1)
11FIRST = 4
12u = rng.random(FIRST)
13print("first draws", " ".join(f"{x:.6f}" for x in u))
10FIRST = 4
11print("first draws", " ".join(f"{npr_pcg_double(rng):.6f}" for i to FIRST))
1411FIRST = 4
1412u = py_npg_random_v(rng, FIRST)
1413print("first draws", join(" ", [fmt(py_x3, ".6f") for py_x3 in u]))
4 estimate pi from random points in the unit square
16N = 20000
17x = rng.random(N)
18y = rng.random(N)
19inside = np.count_nonzero(x * x + y * y < 1.0)
20print(f"pi ~ {4 * inside / N:.4f} from {N} points")
14N = 20000
15x = [npr_pcg_double(rng) for i to N]
16y = [npr_pcg_double(rng) for i to N]
17inside = sum(x * x + y * y < 1.0)
18print(f"pi ~ {4 * inside / N:.4f} from {N:d} points")
1415N = 20000
1416x = py_npg_random_v(rng, N)
1417y = py_npg_random_v(rng, N)
1418inside = len(py_nzflat(py_badd(x * x, y * y) < 1.0))
1419print("pi ~", fmt(py_div_i(py_int(4 * inside), N), ".4f"), "from", fmt(N + 0.0, "d"), "points")
5 dice
23ROLLS = 12
24rolls = rng.integers(1, 7, size=ROLLS)
25print("dice", rolls.tolist())
21ROLLS = 12
22SIDES = 6
23rolls = [1 + npr_gn_bounded(rng, SIDES - 1) for i to ROLLS]
24print("dice", rolls)
1421ROLLS = 12
1422rolls = py_npg_integers_v(rng, 1, 7, ROLLS, 0)
1423print("dice", fmt(rolls, "d"))
6 shuffle a deck
28CARDS = 10
29deck = np.arange(CARDS)
30rng.shuffle(deck)
31print("deck", deck.tolist())
27CARDS = 10
28deck = [i for i to CARDS]
29npr_gn_shuffle(rng, deck)
30print("deck", deck)
1425CARDS = 10
1426deck = py_arange(0, CARDS, 1)
1427py_npg_shuffle(rng, deck)
1428print("deck", fmt(deck, "d"))
7 normal draws
34M = 1000
35z = rng.standard_normal(M)
36print(f"mean {z.mean():.4f}, std {z.std():.4f}")
33M = 1000
34z = [npr_gn_normal(rng) for i to M]
35mean = sum(z) / M
36std = sqrt(sum((z - mean) ^ 2) / M)
37print(f"mean {mean:.4f}, std {std:.4f}")
1430M = 1000
1431z = py_npg_normal_v(rng, 0.0, 1.0, M)
1432print("mean ", fmt(py_mean(z), ".4f"), ", std ", fmt(py_std(z), ".4f"), sep="")
8 checks
39assert inside == np.sum(x * x + y * y < 1.0)
40assert rolls.min() >= 1 and rolls.max() <= 6
40test "draws" {
41 assert_true(min(rolls) >= 1 and max(rolls) <= SIDES)
42 assert_true(all(z == z))
43}
1434if not (inside == py_isum(py_badd(x * x, y * y) < 1.0)) {
1435 stop("AssertionError")
1436}
1437if not (min(rolls) >= 1 and max(rolls) <= 6) {
1438 stop("AssertionError")
1439}

What they print

identical — the two programs print the same text.

Dimmed lines are left out of the comparison. Each Viv test block prints one PASS line; the Python checks the same facts with assert, which prints nothing.

$ python3 seeded.py

first draws 0.773956 0.438878 0.858598 0.697368
pi ~ 3.1432 from 20000 points
dice [3, 6, 1, 4, 2, 2, 3, 4, 5, 4, 1, 2]
deck [1, 4, 7, 6, 2, 9, 3, 8, 0, 5]
mean -0.0596, std 0.9912

$ viv run seeded.viv

first draws 0.773956 0.438878 0.858598 0.697368
pi ~ 3.1432 from 20000 points
dice [3, 6, 1, 4, 2, 2, 3, 4, 5, 4, 1, 2]
deck [1, 4, 7, 6, 2, 9, 3, 8, 0, 5]
mean -0.0596, std 0.9912
PASS draws

$ viv run seeded.py (py2viv’s translation)

The translation prints exactly what the Python prints. 

first draws 0.773956 0.438878 0.858598 0.697368
pi ~ 3.1432 from 20000 points
dice [3, 6, 1, 4, 2, 2, 3, 4, 5, 4, 1, 2]
deck [1, 4, 7, 6, 2, 9, 3, 8, 0, 5]
mean -0.0596, std 0.9912