v/vlib/rand/mini_math.v

131 lines
3.1 KiB
V
Raw Normal View History

// 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) }
}