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  }