rand: add full precision f32 and f64 random functions; fix f32/f64 multipliers (#16875)

Author: j-bm

vlib/math/bits/unsafe_bits.js.v

@@ -0,0 +1,59 @@
+module bits
+
+// f32_bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position.
+// f32_bits(f32_from_bits(x)) == x.
+pub fn f32_bits(f f32) u32 {
+	p := u32(0)
+	#let buffer = new ArrayBuffer(4)
+	#let floatArr = new Float32Array(buffer)
+	#floatArr[0] = f.val
+	#let uintArr = new Uint32Array(buffer)
+	#p.val = uintArr[0]
+
+	return p
+}
+
+// f32_from_bits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// f32_from_bits(f32_bits(x)) == x.
+pub fn f32_from_bits(b u32) f32 {
+	p := f32(0.0)
+	#let buffer = new ArrayBuffer(4)
+	#let floatArr = new Float32Array(buffer)
+	#let uintArr = new Uint32Array(buffer)
+	#uintArr[0] = Number(b.val)
+	#p.val = floatArr[0]
+
+	return p
+}
+
+// f64_bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position,
+// and f64_bits(f64_from_bits(x)) == x.
+pub fn f64_bits(f f64) u64 {
+	p := u64(0)
+	#let buffer = new ArrayBuffer(8)
+	#let floatArr = new Float64Array(buffer)
+	#floatArr[0] = f.val
+	#let uintArr = new BigUint64Array(buffer)
+	#p.val = uintArr[0]
+
+	return p
+}
+
+// f64_from_bits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// f64_from_bits(f64_bits(x)) == x.
+pub fn f64_from_bits(b u64) f64 {
+	p := 0.0
+	#let buffer = new ArrayBuffer(8)
+	#let floatArr = new Float64Array(buffer)
+	#let uintArr = new BigUint64Array(buffer)
+	#uintArr[0] = b.val
+	#p.val = floatArr[0]
+
+	return p
+}

vlib/math/bits/unsafe_bits.v

@@ -0,0 +1,42 @@
+// 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 bits
+
+// f32_bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position.
+// f32_bits(f32_from_bits(x)) == x.
+[inline]
+pub fn f32_bits(f f32) u32 {
+	p := *unsafe { &u32(&f) }
+	return p
+}
+
+// f32_from_bits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// f32_from_bits(f32_bits(x)) == x.
+[inline]
+pub fn f32_from_bits(b u32) f32 {
+	p := *unsafe { &f32(&b) }
+	return p
+}
+
+// f64_bits returns the IEEE 754 binary representation of f,
+// with the sign bit of f and the result in the same bit position,
+// and f64_bits(f64_from_bits(x)) == x.
+[inline]
+pub fn f64_bits(f f64) u64 {
+	p := *unsafe { &u64(&f) }
+	return p
+}
+
+// f64_from_bits returns the floating-point number corresponding
+// to the IEEE 754 binary representation b, with the sign bit of b
+// and the result in the same bit position.
+// f64_from_bits(f64_bits(x)) == x.
+[inline]
+pub fn f64_from_bits(b u64) f64 {
+	p := *unsafe { &f64(&b) }
+	return p
+}

vlib/rand/README.md

@@ -11,7 +11,7 @@ import rand
 ...
 
 // Optionally seed the default generator
-rand.seed([u32(3110), 50714])
+rand.seed([u32(3223878742), 1732001562])
 
 ...
 
@@ -61,21 +61,40 @@ Otherwise, there is feature parity between the generator functions and the top-l
 A PRNG is a Pseudo Random Number Generator.
 Computers cannot generate truly random numbers without an external source of noise or entropy.
 We can use algorithms to generate sequences of seemingly random numbers,
-but their outputs will always be deterministic.
-This is often useful for simulations that need the same starting seed.
+but their outputs will always be deterministic, according to the seed values.
+
+This is often useful for simulations that need the same starting seeds.
+You may be debugging a program and want to restart it with the same
+seeds, or you want to verify a working program is still
+operating identically after compiler or operating system updates.
 
 If you need truly random numbers that are going to be used for cryptography,
 use the `crypto.rand` module.
 
 # Seeding Functions
 
-All the generators are time-seeded.
+All the generators are initialized with time-based seeds.
 The helper functions publicly available in `rand.seed` module are:
 
 1. `time_seed_array()` - returns a `[]u32` that can be directly plugged into the `seed()` functions.
 2. `time_seed_32()` and `time_seed_64()` - 32-bit and 64-bit values respectively
    that are generated from the current time.
 
+When composing your own seeds, use "typical" u32 numbers, not small numbers.  This
+is especially important for PRNGs with large state, such as `mt19937`.  You can create
+random unsigned integers with openssl `rand` or with `v repl` as follows:
+
+```
+$ openssl rand -hex 4
+e3655862
+$ openssl rand -hex 4
+97c4b1db
+$ v repl
+>>> import rand
+>>> [rand.u32(),rand.u32()]
+[2132382944, 2443871665]
+```
+
 # Caveats
 
 Note that the `sys.SysRNG` struct (in the C backend) uses `C.srand()` which sets the seed globally.

vlib/rand/constants/constants.v

@@ -2,11 +2,16 @@ module constants
 
 // Commonly used constants across RNGs - some taken from "Numerical Recipes".
 pub const (
-	lower_mask     = u64(0x00000000FFFFFFFF)
-	max_u32        = u32(0xFFFFFFFF)
-	max_u64        = u64(0xFFFFFFFFFFFFFFFF)
-	max_u32_as_f32 = f32(max_u32) + 1
-	max_u64_as_f64 = f64(max_u64) + 1
-	u31_mask       = u32(0x7FFFFFFF)
-	u63_mask       = u64(0x7FFFFFFFFFFFFFFF)
+	lower_mask                = u64(0x00000000FFFFFFFF)
+	max_u32                   = u32(0xFFFFFFFF)
+	max_u64                   = u64(0xFFFFFFFFFFFFFFFF)
+	u31_mask                  = u32(0x7FFFFFFF)
+	u63_mask                  = u64(0x7FFFFFFFFFFFFFFF)
+	// 23 bits for f32
+	ieee754_mantissa_f32_mask = (u32(1) << 23) - 1
+	// 52 bits for f64
+	ieee754_mantissa_f64_mask = (u64(1) << 52) - 1
+	// smallest mantissa with exponent 0 (un normalized)
+	reciprocal_2_23rd         = 1.0 / f64(u32(1) << 23)
+	reciprocal_2_52nd         = 1.0 / f64(u64(1) << 52)
 )

vlib/rand/fp_test.v

@@ -0,0 +1,116 @@
+import rand
+
+struct Histo {
+mut:
+	lo f64
+	ct int
+}
+
+const (
+	// The sample size to be used; keep cpu time less than 5 seconds
+	count = 9100100
+	// Two sets of seeds
+	seeds = [[u32(2742798260), 2159764996], [u32(2135051596), 958016781]]
+)
+
+fn test_f32() {
+	mut histo := []Histo{}
+	mut candidate := f32(0.0)
+
+	histo << Histo{1.0e-9, 0}
+	histo << Histo{1.0e-6, 0}
+	histo << Histo{1.0e-4, 0}
+	histo << Histo{1.0e-2, 0}
+	histo << Histo{1.0e00, 0}
+
+	for seed in seeds {
+		rand.seed(seed)
+		for _ in 0 .. count {
+			candidate = rand.f32()
+			for mut p in histo {
+				if candidate < p.lo {
+					p.ct += 1
+				}
+			}
+		}
+	}
+	println(' f32  ')
+	println(histo)
+	assert histo[0].ct == 1
+	assert histo[1].ct == 16
+	assert histo[2].ct == 1802
+	assert histo[3].ct == 181963
+	assert histo[4].ct == 18200200
+	for mut p in histo {
+		p.ct = 0
+	}
+	for seed in seeds {
+		rand.seed(seed)
+		for _ in 0 .. count {
+			candidate = rand.f32cp()
+			for mut p in histo {
+				if candidate < p.lo {
+					p.ct += 1
+				}
+			}
+		}
+	}
+	println(' f32cp')
+	println(histo)
+	assert histo[0].ct == 0
+	assert histo[1].ct == 16
+	assert histo[2].ct == 1863
+	assert histo[3].ct == 142044
+	assert histo[4].ct == 18200200
+}
+
+fn test_f64() {
+	mut histo := []Histo{}
+	mut candidate := f64(0.0)
+
+	histo << Histo{1.0e-9, 0}
+	histo << Histo{1.0e-6, 0}
+	histo << Histo{1.0e-4, 0}
+	histo << Histo{1.0e-2, 0}
+	histo << Histo{1.0e00, 0}
+
+	for seed in seeds {
+		rand.seed(seed)
+		for _ in 0 .. count {
+			candidate = rand.f64()
+			for mut p in histo {
+				if candidate < p.lo {
+					p.ct += 1
+				}
+			}
+		}
+	}
+	println(' f64  ')
+	println(histo)
+	assert histo[0].ct == 0
+	assert histo[1].ct == 23
+	assert histo[2].ct == 1756
+	assert histo[3].ct == 182209
+	assert histo[4].ct == 18200200
+	for mut p in histo {
+		p.ct = 0
+	}
+	for seed in seeds {
+		rand.seed(seed)
+		for _ in 0 .. count {
+			candidate = rand.f64cp()
+			for mut p in histo {
+				if candidate < p.lo {
+					p.ct += 1
+				}
+			}
+		}
+	}
+	println(' f64cp')
+	println(histo)
+	assert histo[0].ct == 0
+	assert histo[1].ct == 17
+	assert histo[2].ct == 1878
+	assert histo[3].ct == 181754
+	assert histo[4].ct == 18200200
+}