diff --git a/internal/rand/rand.go b/internal/rand/rand.go deleted file mode 100644 index 521d8f40..00000000 --- a/internal/rand/rand.go +++ /dev/null @@ -1,393 +0,0 @@ -// Copyright ©2024 The Gonum Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -// Copyright 2009 The Go Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -// Package rand implements pseudo-random number generators. -// -// Random numbers are generated by a Source. Top-level functions, such as -// Float64 and Int, use a default shared Source that produces a deterministic -// sequence of values each time a program is run. Use the Seed function to -// initialize the default Source if different behavior is required for each run. -// The default Source, a LockedSource, is safe for concurrent use by multiple -// goroutines, but Sources created by NewSource are not. However, Sources are small -// and it is reasonable to have a separate Source for each goroutine, seeded -// differently, to avoid locking. -// -// For random numbers suitable for security-sensitive work, see the crypto/rand -// package. -package rand - -import ( - "math/rand/v2" - "sync" -) - -// A Source represents a source of uniformly-distributed -// pseudo-random int64 values in the range [0, 1<<64). -type Source interface { - Uint64() uint64 - Seed(seed uint64) -} - -// NewSource returns a new pseudo-random Source seeded with the given value. -func NewSource(seed uint64) Source { - return &pcgShim{rand.NewPCG(seed, seed)} -} - -type pcgShim struct { - *rand.PCG -} - -func (p *pcgShim) Seed(seed uint64) { - p.PCG.Seed(seed, seed) -} - -// A Rand is a source of random numbers. -type Rand struct { - src Source - - // readVal contains remainder of 64-bit integer used for bytes - // generation during most recent Read call. - // It is saved so next Read call can start where the previous - // one finished. - readVal uint64 - // readPos indicates the number of low-order bytes of readVal - // that are still valid. - readPos int8 -} - -// New returns a new Rand that uses random values from src -// to generate other random values. -func New(src Source) *Rand { - return &Rand{src: src} -} - -func (r *Rand) NormFloat64() float64 { - return rand.New(r.src).NormFloat64() -} - -func (r *Rand) ExpFloat64() float64 { - return rand.New(r.src).ExpFloat64() -} - -// Seed uses the provided seed value to initialize the generator to a deterministic state. -// Seed should not be called concurrently with any other Rand method. -func (r *Rand) Seed(seed uint64) { - if lk, ok := r.src.(*LockedSource); ok { - lk.seedPos(seed, &r.readPos) - return - } - - r.src.Seed(seed) - r.readPos = 0 -} - -// Uint64 returns a pseudo-random 64-bit integer as a uint64. -func (r *Rand) Uint64() uint64 { return r.src.Uint64() } - -// Int63 returns a non-negative pseudo-random 63-bit integer as an int64. -func (r *Rand) Int63() int64 { return int64(r.src.Uint64() &^ (1 << 63)) } - -// Uint32 returns a pseudo-random 32-bit value as a uint32. -func (r *Rand) Uint32() uint32 { return uint32(r.Uint64() >> 32) } - -// Int31 returns a non-negative pseudo-random 31-bit integer as an int32. -func (r *Rand) Int31() int32 { return int32(r.Uint64() >> 33) } - -// Int returns a non-negative pseudo-random int. -func (r *Rand) Int() int { - u := uint(r.Uint64()) - return int(u << 1 >> 1) // clear sign bit. -} - -const maxUint64 = (1 << 64) - 1 - -// Uint64n returns, as a uint64, a pseudo-random number in [0,n). -// It is guaranteed more uniform than taking a Source value mod n -// for any n that is not a power of 2. -func (r *Rand) Uint64n(n uint64) uint64 { - if n&(n-1) == 0 { // n is power of two, can mask - if n == 0 { - panic("invalid argument to Uint64n") - } - return r.Uint64() & (n - 1) - } - // If n does not divide v, to avoid bias we must not use - // a v that is within maxUint64%n of the top of the range. - v := r.Uint64() - if v > maxUint64-n { // Fast check. - ceiling := maxUint64 - maxUint64%n - for v >= ceiling { - v = r.Uint64() - } - } - - return v % n -} - -// Int63n returns, as an int64, a non-negative pseudo-random number in [0,n). -// It panics if n <= 0. -func (r *Rand) Int63n(n int64) int64 { - if n <= 0 { - panic("invalid argument to Int63n") - } - return int64(r.Uint64n(uint64(n))) -} - -// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n). -// It panics if n <= 0. -func (r *Rand) Int31n(n int32) int32 { - if n <= 0 { - panic("invalid argument to Int31n") - } - // TODO: Avoid some 64-bit ops to make it more efficient on 32-bit machines. - return int32(r.Uint64n(uint64(n))) -} - -// Intn returns, as an int, a non-negative pseudo-random number in [0,n). -// It panics if n <= 0. -func (r *Rand) Intn(n int) int { - if n <= 0 { - panic("invalid argument to Intn") - } - // TODO: Avoid some 64-bit ops to make it more efficient on 32-bit machines. - return int(r.Uint64n(uint64(n))) -} - -// Float64 returns, as a float64, a pseudo-random number in [0.0,1.0). -func (r *Rand) Float64() float64 { - // There is one bug in the value stream: r.Int63() may be so close - // to 1<<63 that the division rounds up to 1.0, and we've guaranteed - // that the result is always less than 1.0. - // - // We tried to fix this by mapping 1.0 back to 0.0, but since float64 - // values near 0 are much denser than near 1, mapping 1 to 0 caused - // a theoretically significant overshoot in the probability of returning 0. - // Instead of that, if we round up to 1, just try again. - // Getting 1 only happens 1/2⁵³ of the time, so most clients - // will not observe it anyway. -again: - f := float64(r.Uint64n(1<<53)) / (1 << 53) - if f == 1.0 { - goto again // resample; this branch is taken O(never) - } - return f -} - -// Float32 returns, as a float32, a pseudo-random number in [0.0,1.0). -func (r *Rand) Float32() float32 { - // We do not want to return 1.0. - // This only happens 1/2²⁴ of the time (plus the 1/2⁵³ of the time in Float64). -again: - f := float32(r.Float64()) - if f == 1 { - goto again // resample; this branch is taken O(very rarely) - } - return f -} - -// Perm returns, as a slice of n ints, a pseudo-random permutation of the integers [0,n). -func (r *Rand) Perm(n int) []int { - m := make([]int, n) - // In the following loop, the iteration when i=0 always swaps m[0] with m[0]. - // A change to remove this useless iteration is to assign 1 to i in the init - // statement. But Perm also effects r. Making this change will affect - // the final state of r. So this change can't be made for compatibility - // reasons for Go 1. - for i := 0; i < n; i++ { - j := r.Intn(i + 1) - m[i] = m[j] - m[j] = i - } - return m -} - -// Shuffle pseudo-randomizes the order of elements. -// n is the number of elements. Shuffle panics if n < 0. -// swap swaps the elements with indexes i and j. -func (r *Rand) Shuffle(n int, swap func(i, j int)) { - if n < 0 { - panic("invalid argument to Shuffle") - } - - // Fisher-Yates shuffle: https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle - // Shuffle really ought not be called with n that doesn't fit in 32 bits. - // Not only will it take a very long time, but with 2³¹! possible permutations, - // there's no way that any PRNG can have a big enough internal state to - // generate even a minuscule percentage of the possible permutations. - // Nevertheless, the right API signature accepts an int n, so handle it as best we can. - i := n - 1 - for ; i > 1<<31-1-1; i-- { - j := int(r.Int63n(int64(i + 1))) - swap(i, j) - } - for ; i > 0; i-- { - j := int(r.Int31n(int32(i + 1))) - swap(i, j) - } -} - -// Read generates len(p) random bytes and writes them into p. It -// always returns len(p) and a nil error. -// Read should not be called concurrently with any other Rand method unless -// the underlying source is a LockedSource. -func (r *Rand) Read(p []byte) (n int, err error) { - if lk, ok := r.src.(*LockedSource); ok { - return lk.Read(p, &r.readVal, &r.readPos) - } - return read(p, r.src, &r.readVal, &r.readPos) -} - -func read(p []byte, src Source, readVal *uint64, readPos *int8) (n int, err error) { - pos := *readPos - val := *readVal - rng, _ := src.(*pcgShim) - for n = 0; n < len(p); n++ { - if pos == 0 { - if rng != nil { - val = rng.Uint64() - } else { - val = src.Uint64() - } - pos = 8 - } - p[n] = byte(val) - val >>= 8 - pos-- - } - *readPos = pos - *readVal = val - return -} - -/* - * Top-level convenience functions - */ - -var globalRand = New(&LockedSource{src: *NewSource(1).(*pcgShim)}) - -// Type assert that globalRand's source is a LockedSource whose src is a PCGSource. -var _ pcgShim = globalRand.src.(*LockedSource).src - -// Seed uses the provided seed value to initialize the default Source to a -// deterministic state. If Seed is not called, the generator behaves as -// if seeded by Seed(1). -// Seed, unlike the Rand.Seed method, is safe for concurrent use. -func Seed(seed uint64) { globalRand.Seed(seed) } - -// Int63 returns a non-negative pseudo-random 63-bit integer as an int64 -// from the default Source. -func Int63() int64 { return globalRand.Int63() } - -// Uint32 returns a pseudo-random 32-bit value as a uint32 -// from the default Source. -func Uint32() uint32 { return globalRand.Uint32() } - -// Uint64 returns a pseudo-random 64-bit value as a uint64 -// from the default Source. -func Uint64() uint64 { return globalRand.Uint64() } - -// Int31 returns a non-negative pseudo-random 31-bit integer as an int32 -// from the default Source. -func Int31() int32 { return globalRand.Int31() } - -// Int returns a non-negative pseudo-random int from the default Source. -func Int() int { return globalRand.Int() } - -// Int63n returns, as an int64, a non-negative pseudo-random number in [0,n) -// from the default Source. -// It panics if n <= 0. -func Int63n(n int64) int64 { return globalRand.Int63n(n) } - -// Int31n returns, as an int32, a non-negative pseudo-random number in [0,n) -// from the default Source. -// It panics if n <= 0. -func Int31n(n int32) int32 { return globalRand.Int31n(n) } - -// Intn returns, as an int, a non-negative pseudo-random number in [0,n) -// from the default Source. -// It panics if n <= 0. -func Intn(n int) int { return globalRand.Intn(n) } - -// Float64 returns, as a float64, a pseudo-random number in [0.0,1.0) -// from the default Source. -func Float64() float64 { return globalRand.Float64() } - -// Float32 returns, as a float32, a pseudo-random number in [0.0,1.0) -// from the default Source. -func Float32() float32 { return globalRand.Float32() } - -// Perm returns, as a slice of n ints, a pseudo-random permutation of the integers [0,n) -// from the default Source. -func Perm(n int) []int { return globalRand.Perm(n) } - -// Shuffle pseudo-randomizes the order of elements using the default Source. -// n is the number of elements. Shuffle panics if n < 0. -// swap swaps the elements with indexes i and j. -func Shuffle(n int, swap func(i, j int)) { globalRand.Shuffle(n, swap) } - -// Read generates len(p) random bytes from the default Source and -// writes them into p. It always returns len(p) and a nil error. -// Read, unlike the Rand.Read method, is safe for concurrent use. -func Read(p []byte) (n int, err error) { return globalRand.Read(p) } - -// NormFloat64 returns a normally distributed float64 in the range -// [-math.MaxFloat64, +math.MaxFloat64] with -// standard normal distribution (mean = 0, stddev = 1) -// from the default Source. -// To produce a different normal distribution, callers can -// adjust the output using: -// -// sample = NormFloat64() * desiredStdDev + desiredMean -func NormFloat64() float64 { return globalRand.NormFloat64() } - -// ExpFloat64 returns an exponentially distributed float64 in the range -// (0, +math.MaxFloat64] with an exponential distribution whose rate parameter -// (lambda) is 1 and whose mean is 1/lambda (1) from the default Source. -// To produce a distribution with a different rate parameter, -// callers can adjust the output using: -// -// sample = ExpFloat64() / desiredRateParameter -func ExpFloat64() float64 { return globalRand.ExpFloat64() } - -// LockedSource is an implementation of Source that is concurrency-safe. -// A Rand using a LockedSource is safe for concurrent use. -// -// The zero value of LockedSource is valid, but should be seeded before use. -type LockedSource struct { - lk sync.Mutex - src pcgShim -} - -func (s *LockedSource) Uint64() (n uint64) { - s.lk.Lock() - n = s.src.Uint64() - s.lk.Unlock() - return -} - -func (s *LockedSource) Seed(seed uint64) { - s.lk.Lock() - s.src.Seed(seed) - s.lk.Unlock() -} - -// seedPos implements Seed for a LockedSource without a race condiiton. -func (s *LockedSource) seedPos(seed uint64, readPos *int8) { - s.lk.Lock() - s.src.Seed(seed) - *readPos = 0 - s.lk.Unlock() -} - -// Read implements Read for a LockedSource. -func (s *LockedSource) Read(p []byte, readVal *uint64, readPos *int8) (n int, err error) { - s.lk.Lock() - n, err = read(p, &s.src, readVal, readPos) - s.lk.Unlock() - return -} diff --git a/internal/rand/rand_test.go b/internal/rand/rand_test.go deleted file mode 100644 index a528fc0e..00000000 --- a/internal/rand/rand_test.go +++ /dev/null @@ -1,368 +0,0 @@ -// Copyright ©2024 The Gonum Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -// Copyright 2009 The Go Authors. All rights reserved. -// Use of this source code is governed by a BSD-style -// license that can be found in the LICENSE file. - -package rand - -import ( - "bytes" - "errors" - "fmt" - "io" - "math" - "os" - "runtime" - "strings" - "testing" - "testing/iotest" -) - -const ( - numTestSamples = 10000 -) - -type statsResults struct { - mean float64 - stddev float64 - closeEnough float64 - maxError float64 -} - -func nearEqual(a, b, closeEnough, maxError float64) bool { - absDiff := math.Abs(a - b) - if absDiff < closeEnough { // Necessary when one value is zero and one value is close to zero. - return true - } - return absDiff/max(math.Abs(a), math.Abs(b)) < maxError -} - -var testSeeds = []uint64{1, 1754801282, 1698661970, 1550503961} - -// checkSimilarDistribution returns success if the mean and stddev of the -// two statsResults are similar. -func (sr *statsResults) checkSimilarDistribution(expected *statsResults) error { - if !nearEqual(sr.mean, expected.mean, expected.closeEnough, expected.maxError) { - s := fmt.Sprintf("mean %v != %v (allowed error %v, %v)", sr.mean, expected.mean, expected.closeEnough, expected.maxError) - fmt.Println(s) - return errors.New(s) - } - if !nearEqual(sr.stddev, expected.stddev, expected.closeEnough, expected.maxError) { - s := fmt.Sprintf("stddev %v != %v (allowed error %v, %v)", sr.stddev, expected.stddev, expected.closeEnough, expected.maxError) - fmt.Println(s) - return errors.New(s) - } - return nil -} - -func getStatsResults(samples []float64) *statsResults { - res := new(statsResults) - var sum, squaresum float64 - for _, s := range samples { - sum += s - squaresum += s * s - } - res.mean = sum / float64(len(samples)) - res.stddev = math.Sqrt(squaresum/float64(len(samples)) - res.mean*res.mean) - return res -} - -func checkSampleDistribution(t *testing.T, samples []float64, expected *statsResults) { - t.Helper() - actual := getStatsResults(samples) - err := actual.checkSimilarDistribution(expected) - if err != nil { - t.Error(err) - } -} - -func checkSampleSliceDistributions(t *testing.T, samples []float64, nslices int, expected *statsResults) { - t.Helper() - chunk := len(samples) / nslices - for i := 0; i < nslices; i++ { - low := i * chunk - var high int - if i == nslices-1 { - high = len(samples) - 1 - } else { - high = (i + 1) * chunk - } - checkSampleDistribution(t, samples[low:high], expected) - } -} - -// -// Normal distribution tests -// - -func generateNormalSamples(nsamples int, mean, stddev float64, seed uint64) []float64 { - r := New(NewSource(seed)) - samples := make([]float64, nsamples) - for i := range samples { - samples[i] = r.NormFloat64()*stddev + mean - } - return samples -} - -func testNormalDistribution(t *testing.T, nsamples int, mean, stddev float64, seed uint64) { - //fmt.Printf("testing nsamples=%v mean=%v stddev=%v seed=%v\n", nsamples, mean, stddev, seed); - - samples := generateNormalSamples(nsamples, mean, stddev, seed) - errorScale := max(1.0, stddev) // Error scales with stddev - expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.08 * errorScale} - - // Make sure that the entire set matches the expected distribution. - checkSampleDistribution(t, samples, expected) - - // Make sure that each half of the set matches the expected distribution. - checkSampleSliceDistributions(t, samples, 2, expected) - - // Make sure that each 7th of the set matches the expected distribution. - checkSampleSliceDistributions(t, samples, 7, expected) -} - -// Actual tests - -func TestStandardNormalValues(t *testing.T) { - for _, seed := range testSeeds { - testNormalDistribution(t, numTestSamples, 0, 1, seed) - } -} - -func TestNonStandardNormalValues(t *testing.T) { - sdmax := 1000.0 - mmax := 1000.0 - if testing.Short() { - sdmax = 5 - mmax = 5 - } - for sd := 0.5; sd < sdmax; sd *= 2 { - for m := 0.5; m < mmax; m *= 2 { - for _, seed := range testSeeds { - testNormalDistribution(t, numTestSamples, m, sd, seed) - if testing.Short() { - break - } - } - } - } -} - -// -// Exponential distribution tests -// - -func generateExponentialSamples(nsamples int, rate float64, seed uint64) []float64 { - r := New(NewSource(seed)) - samples := make([]float64, nsamples) - for i := range samples { - samples[i] = r.ExpFloat64() / rate - } - return samples -} - -func testExponentialDistribution(t *testing.T, nsamples int, rate float64, seed uint64) { - //fmt.Printf("testing nsamples=%v rate=%v seed=%v\n", nsamples, rate, seed) - - mean := 1 / rate - stddev := mean - - samples := generateExponentialSamples(nsamples, rate, seed) - errorScale := max(1.0, 1/rate) // Error scales with the inverse of the rate - expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.20 * errorScale} - - // Make sure that the entire set matches the expected distribution. - checkSampleDistribution(t, samples, expected) - - // Make sure that each half of the set matches the expected distribution. - checkSampleSliceDistributions(t, samples, 2, expected) - - // Make sure that each 7th of the set matches the expected distribution. - checkSampleSliceDistributions(t, samples, 7, expected) -} - -// Actual tests - -func TestStandardExponentialValues(t *testing.T) { - for _, seed := range testSeeds { - testExponentialDistribution(t, numTestSamples, 1, seed) - } -} - -func TestNonStandardExponentialValues(t *testing.T) { - for rate := 0.05; rate < 10; rate *= 2 { - for _, seed := range testSeeds { - testExponentialDistribution(t, numTestSamples, rate, seed) - if testing.Short() { - break - } - } - } -} - -// compareUint32Slices returns the first index where the two slices -// disagree, or <0 if the lengths are the same and all elements -// are identical. -func compareUint32Slices(s1, s2 []uint32) int { - if len(s1) != len(s2) { - if len(s1) > len(s2) { - return len(s2) + 1 - } - return len(s1) + 1 - } - for i := range s1 { - if s1[i] != s2[i] { - return i - } - } - return -1 -} - -// compareFloat32Slices returns the first index where the two slices -// disagree, or <0 if the lengths are the same and all elements -// are identical. -func compareFloat32Slices(s1, s2 []float32) int { - if len(s1) != len(s2) { - if len(s1) > len(s2) { - return len(s2) + 1 - } - return len(s1) + 1 - } - for i := range s1 { - if !nearEqual(float64(s1[i]), float64(s2[i]), 0, 1e-7) { - return i - } - } - return -1 -} - -func hasSlowFloatingPoint() bool { - switch runtime.GOARCH { - case "arm": - return os.Getenv("GOARM") == "5" || strings.HasSuffix(os.Getenv("GOARM"), ",softfloat") - case "mips", "mipsle", "mips64", "mips64le": - // Be conservative and assume that all mips boards - // have emulated floating point. - // TODO: detect what it actually has. - return true - } - return false -} - -func TestFloat32(t *testing.T) { - // For issue 6721, the problem came after 7533753 calls, so check 10e6. - num := int(10e6) - // But do the full amount only on builders (not locally). - // But ARM5 floating point emulation is slow (Issue 10749), so - // do less for that builder: - if testing.Short() && hasSlowFloatingPoint() { // TODO: (testenv.Builder() == "" || hasSlowFloatingPoint()) - num /= 100 // 1.72 seconds instead of 172 seconds - } - - r := New(NewSource(1)) - for ct := 0; ct < num; ct++ { - f := r.Float32() - if f >= 1 { - t.Fatal("Float32() should be in range [0,1). ct:", ct, "f:", f) - } - } -} - -func testReadUniformity(t *testing.T, n int, seed uint64) { - r := New(NewSource(seed)) - buf := make([]byte, n) - nRead, err := r.Read(buf) - if err != nil { - t.Errorf("Read err %v", err) - } - if nRead != n { - t.Errorf("Read returned unexpected n; %d != %d", nRead, n) - } - - // Expect a uniform distribution of byte values, which lie in [0, 255]. - var ( - mean = 255.0 / 2 - stddev = 256.0 / math.Sqrt(12.0) - errorScale = stddev / math.Sqrt(float64(n)) - ) - - expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.08 * errorScale} - - // Cast bytes as floats to use the common distribution-validity checks. - samples := make([]float64, n) - for i, val := range buf { - samples[i] = float64(val) - } - // Make sure that the entire set matches the expected distribution. - checkSampleDistribution(t, samples, expected) -} - -func TestReadUniformity(t *testing.T) { - testBufferSizes := []int{ - 2, 4, 7, 64, 1024, 1 << 16, 1 << 20, - } - for _, seed := range testSeeds { - for _, n := range testBufferSizes { - testReadUniformity(t, n, seed) - } - } -} - -func TestReadEmpty(t *testing.T) { - r := New(NewSource(1)) - buf := make([]byte, 0) - n, err := r.Read(buf) - if err != nil { - t.Errorf("Read err into empty buffer; %v", err) - } - if n != 0 { - t.Errorf("Read into empty buffer returned unexpected n of %d", n) - } -} - -func TestReadByOneByte(t *testing.T) { - r := New(NewSource(1)) - b1 := make([]byte, 100) - _, err := io.ReadFull(iotest.OneByteReader(r), b1) - if err != nil { - t.Errorf("read by one byte: %v", err) - } - r = New(NewSource(1)) - b2 := make([]byte, 100) - _, err = r.Read(b2) - if err != nil { - t.Errorf("read: %v", err) - } - if !bytes.Equal(b1, b2) { - t.Errorf("read by one byte vs single read:\n%x\n%x", b1, b2) - } -} - -func TestReadSeedReset(t *testing.T) { - r := New(NewSource(42)) - b1 := make([]byte, 128) - _, err := r.Read(b1) - if err != nil { - t.Errorf("read: %v", err) - } - r.Seed(42) - b2 := make([]byte, 128) - _, err = r.Read(b2) - if err != nil { - t.Errorf("read: %v", err) - } - if !bytes.Equal(b1, b2) { - t.Errorf("mismatch after re-seed:\n%x\n%x", b1, b2) - } -} - -func TestShuffleSmall(t *testing.T) { - // Check that Shuffle allows n=0 and n=1, but that swap is never called for them. - r := New(NewSource(1)) - for n := 0; n <= 1; n++ { - r.Shuffle(n, func(i, j int) { t.Fatalf("swap called, n=%d i=%d j=%d", n, i, j) }) - } -} diff --git a/internal/rand/staticcheck.conf b/internal/rand/staticcheck.conf deleted file mode 100644 index d533ce21..00000000 --- a/internal/rand/staticcheck.conf +++ /dev/null @@ -1 +0,0 @@ -checks = ["inherit", "-U1000"] \ No newline at end of file