pgregory.net/rand@v1.0.3-0.20230808192358-a0b8ce02f4da/bias_test.go (about) 1 // Copyright 2022 Gregory Petrosyan <gregory.petrosyan@gmail.com> 2 // 3 // This Source Code Form is subject to the terms of the Mozilla Public 4 // License, v. 2.0. If a copy of the MPL was not distributed with this 5 // file, You can obtain one at https://mozilla.org/MPL/2.0/. 6 7 package rand_test 8 9 import ( 10 "flag" 11 "fmt" 12 "math" 13 "math/bits" 14 "pgregory.net/rand" 15 "testing" 16 ) 17 18 var ( 19 biasFlag = flag.Bool("bias", false, "run bias-detection tests") 20 ) 21 22 func uint16nModulo(r *rand.Rand, n uint16) uint16 { 23 return uint16(r.Uint32()) % n // biased 24 } 25 26 func uint16nFixedPoint(r *rand.Rand, n uint16, m int) uint16 { 27 res, frac := bits.Mul32(uint32(n), r.Uint32()&(1<<(16+m)-1)) 28 return (uint16(res) << (16 - m)) | (uint16(frac>>16) >> m) // biased with probability 2^-m 29 } 30 31 func uint16nFixedPoint0(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 0) } 32 func uint16nFixedPoint2(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 2) } 33 func uint16nFixedPoint4(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 4) } 34 func uint16nFixedPoint8(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 8) } 35 func uint16nFixedPoint10(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 10) } 36 func uint16nFixedPoint12(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 12) } 37 func uint16nFixedPoint14(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 14) } 38 func uint16nFixedPoint16(r *rand.Rand, n uint16) uint16 { return uint16nFixedPoint(r, n, 16) } 39 40 func uint16nCanon(r *rand.Rand, n uint16) uint16 { 41 res, frac := bits.Mul64(uint64(n), r.Uint64()) 42 hi, _ := bits.Mul64(uint64(n), r.Uint64()) 43 _, carry := bits.Add64(frac, hi, 0) 44 return uint16(res + carry) // biased with probability 2^-64 45 } 46 47 func uint16nLemire(r *rand.Rand, n uint16) uint16 { 48 v := uint16(r.Uint32()) 49 prod := uint32(v) * uint32(n) 50 low := uint16(prod) 51 if low < n { 52 thresh := -n % n 53 for low < thresh { 54 v = uint16(r.Uint32()) 55 prod = uint32(v) * uint32(n) 56 low = uint16(prod) 57 } 58 } 59 return uint16(prod >> 16) // unbiased 60 } 61 62 func TestRand_Uint16nBias_Modulo(t *testing.T) { testRandUint16nBias(t, uint16nModulo) } 63 func TestRand_Uint16nBias_FixedPoint0(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint0) } 64 func TestRand_Uint16nBias_FixedPoint2(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint2) } 65 func TestRand_Uint16nBias_FixedPoint4(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint4) } 66 func TestRand_Uint16nBias_FixedPoint8(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint8) } 67 func TestRand_Uint16nBias_FixedPoint10(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint10) } 68 func TestRand_Uint16nBias_FixedPoint12(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint12) } 69 func TestRand_Uint16nBias_FixedPoint14(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint14) } 70 func TestRand_Uint16nBias_FixedPoint16(t *testing.T) { testRandUint16nBias(t, uint16nFixedPoint16) } 71 func TestRand_Uint16nBias_Canon(t *testing.T) { testRandUint16nBias(t, uint16nCanon) } 72 func TestRand_Uint16nBias_Lemire(t *testing.T) { testRandUint16nBias(t, uint16nLemire) } 73 74 func testRandUint16nBias(t *testing.T, gen func(*rand.Rand, uint16) uint16) { 75 t.Helper() 76 77 if !*biasFlag { 78 t.Skip("specify -bias flag to run bias tests") 79 } 80 81 const bound = math.MaxUint16 / 4 * 3 82 for pow := 10; pow <= 50; pow++ { 83 t.Run(fmt.Sprintf("%d/%dbit", bound, pow), func(t *testing.T) { 84 r := rand.New() 85 data := make([]uint64, bound) 86 attempts := 1 << int64(pow) 87 for i := 0; i < attempts; i++ { 88 ix := gen(r, bound) 89 data[ix]++ 90 } 91 92 var chiSq float64 93 expected := float64(attempts) / float64(bound) 94 for _, n := range data { 95 obs := float64(n) 96 chiSq += (obs - expected) * (obs - expected) 97 } 98 chiSq /= expected 99 100 df := float64(bound - 1) 101 t.Logf("χ2 = %.1f, DoF = %v (%v attempts, delta = %.1f%%)", chiSq, df, attempts, (chiSq-df)/df*100) 102 }) 103 } 104 }