131 lines
3.1 KiB
V
131 lines
3.1 KiB
V
// Copyright (c) 2019-2022 Alexander Medvednikov. All rights reserved.
|
|
// Use of this source code is governed by an MIT license
|
|
// that can be found in the LICENSE file.
|
|
module rand
|
|
|
|
// NOTE: mini_math.v exists, so that we can avoid `import math`,
|
|
// just for the math.log and math.sqrt functions needed for the
|
|
// non uniform random number redistribution functions.
|
|
// Importing math is relatively heavy, both in terms of compilation
|
|
// speed (more source to process), and in terms of increases in the
|
|
// generated executable sizes (if the rest of the program does not use
|
|
// math already; many programs do not need math, for example the
|
|
// compiler itself does not, while needing random number generation.
|
|
|
|
const sqrt2 = 1.41421356237309504880168872420969807856967187537694807317667974
|
|
|
|
[inline]
|
|
fn msqrt(a f64) f64 {
|
|
if a == 0 {
|
|
return a
|
|
}
|
|
mut x := a
|
|
z, ex := frexp(x)
|
|
w := x
|
|
// approximate square root of number between 0.5 and 1
|
|
// relative error of approximation = 7.47e-3
|
|
x = 4.173075996388649989089e-1 + 5.9016206709064458299663e-1 * z // adjust for odd powers of 2
|
|
if (ex & 1) != 0 {
|
|
x *= rand.sqrt2
|
|
}
|
|
x = scalbn(x, ex >> 1)
|
|
// newton iterations
|
|
x = 0.5 * (x + w / x)
|
|
x = 0.5 * (x + w / x)
|
|
x = 0.5 * (x + w / x)
|
|
return x
|
|
}
|
|
|
|
// a simplified approximation (without the edge cases), see math.log
|
|
fn mlog(a f64) f64 {
|
|
ln2_lo := 1.90821492927058770002e-10
|
|
ln2_hi := 0.693147180369123816490
|
|
l1 := 0.6666666666666735130
|
|
l2 := 0.3999999999940941908
|
|
l3 := 0.2857142874366239149
|
|
l4 := 0.2222219843214978396
|
|
l5 := 0.1818357216161805012
|
|
l6 := 0.1531383769920937332
|
|
l7 := 0.1479819860511658591
|
|
x := a
|
|
mut f1, mut ki := frexp(x)
|
|
if f1 < rand.sqrt2 / 2 {
|
|
f1 *= 2
|
|
ki--
|
|
}
|
|
f := f1 - 1
|
|
k := f64(ki)
|
|
s := f / (2 + f)
|
|
s2 := s * s
|
|
s4 := s2 * s2
|
|
t1 := s2 * (l1 + s4 * (l3 + s4 * (l5 + s4 * l7)))
|
|
t2 := s4 * (l2 + s4 * (l4 + s4 * l6))
|
|
r := t1 + t2
|
|
hfsq := 0.5 * f * f
|
|
return k * ln2_hi - ((hfsq - (s * (hfsq + r) + k * ln2_lo)) - f)
|
|
}
|
|
|
|
fn frexp(x f64) (f64, int) {
|
|
mut y := f64_bits(x)
|
|
ee := int((y >> 52) & 0x7ff)
|
|
if ee == 0 {
|
|
if x != 0.0 {
|
|
x1p64 := f64_from_bits(u64(0x43f0000000000000))
|
|
z, e_ := frexp(x * x1p64)
|
|
return z, e_ - 64
|
|
}
|
|
return x, 0
|
|
} else if ee == 0x7ff {
|
|
return x, 0
|
|
}
|
|
e_ := ee - 0x3fe
|
|
y &= u64(0x800fffffffffffff)
|
|
y |= u64(0x3fe0000000000000)
|
|
return f64_from_bits(y), e_
|
|
}
|
|
|
|
fn scalbn(x f64, n_ int) f64 {
|
|
mut n := n_
|
|
x1p1023 := f64_from_bits(u64(0x7fe0000000000000))
|
|
x1p53 := f64_from_bits(u64(0x4340000000000000))
|
|
x1p_1022 := f64_from_bits(u64(0x0010000000000000))
|
|
|
|
mut y := x
|
|
if n > 1023 {
|
|
y *= x1p1023
|
|
n -= 1023
|
|
if n > 1023 {
|
|
y *= x1p1023
|
|
n -= 1023
|
|
if n > 1023 {
|
|
n = 1023
|
|
}
|
|
}
|
|
} else if n < -1022 {
|
|
/*
|
|
make sure final n < -53 to avoid double
|
|
rounding in the subnormal range
|
|
*/
|
|
y *= x1p_1022 * x1p53
|
|
n += 1022 - 53
|
|
if n < -1022 {
|
|
y *= x1p_1022 * x1p53
|
|
n += 1022 - 53
|
|
if n < -1022 {
|
|
n = -1022
|
|
}
|
|
}
|
|
}
|
|
return y * f64_from_bits(u64((0x3ff + n)) << 52)
|
|
}
|
|
|
|
[inline]
|
|
fn f64_from_bits(b u64) f64 {
|
|
return *unsafe { &f64(&b) }
|
|
}
|
|
|
|
[inline]
|
|
fn f64_bits(f f64) u64 {
|
|
return *unsafe { &u64(&f) }
|
|
}
|