src/math/rand/v2/rand_test.go

   1 // Copyright 2009 The Go Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package rand_test
   6
   7 import (
   8         "errors"
   9         "fmt"
  10         "internal/testenv"
  11         "math"
  12         . "math/rand/v2"
  13         "os"
  14         "runtime"
  15         "sync"
  16         "testing"
  17 )
  18
  19 const (
  20         numTestSamples = 10000
  21 )
  22
  23 var rn, kn, wn, fn = GetNormalDistributionParameters()
  24 var re, ke, we, fe = GetExponentialDistributionParameters()
  25
  26 type statsResults struct {
  27         mean        float64
  28         stddev      float64
  29         closeEnough float64
  30         maxError    float64
  31 }
  32
  33 func max(a, b float64) float64 {
  34         if a > b {
  35                 return a
  36         }
  37         return b
  38 }
  39
  40 func nearEqual(a, b, closeEnough, maxError float64) bool {
  41         absDiff := math.Abs(a - b)
  42         if absDiff < closeEnough { // Necessary when one value is zero and one value is close to zero.
  43                 return true
  44         }
  45         return absDiff/max(math.Abs(a), math.Abs(b)) < maxError
  46 }
  47
  48 var testSeeds = []int64{1, 1754801282, 1698661970, 1550503961}
  49
  50 // checkSimilarDistribution returns success if the mean and stddev of the
  51 // two statsResults are similar.
  52 func (this *statsResults) checkSimilarDistribution(expected *statsResults) error {
  53         if !nearEqual(this.mean, expected.mean, expected.closeEnough, expected.maxError) {
  54                 s := fmt.Sprintf("mean %v != %v (allowed error %v, %v)", this.mean, expected.mean, expected.closeEnough, expected.maxError)
  55                 fmt.Println(s)
  56                 return errors.New(s)
  57         }
  58         if !nearEqual(this.stddev, expected.stddev, expected.closeEnough, expected.maxError) {
  59                 s := fmt.Sprintf("stddev %v != %v (allowed error %v, %v)", this.stddev, expected.stddev, expected.closeEnough, expected.maxError)
  60                 fmt.Println(s)
  61                 return errors.New(s)
  62         }
  63         return nil
  64 }
  65
  66 func getStatsResults(samples []float64) *statsResults {
  67         res := new(statsResults)
  68         var sum, squaresum float64
  69         for _, s := range samples {
  70                 sum += s
  71                 squaresum += s * s
  72         }
  73         res.mean = sum / float64(len(samples))
  74         res.stddev = math.Sqrt(squaresum/float64(len(samples)) - res.mean*res.mean)
  75         return res
  76 }
  77
  78 func checkSampleDistribution(t *testing.T, samples []float64, expected *statsResults) {
  79         t.Helper()
  80         actual := getStatsResults(samples)
  81         err := actual.checkSimilarDistribution(expected)
  82         if err != nil {
  83                 t.Errorf(err.Error())
  84         }
  85 }
  86
  87 func checkSampleSliceDistributions(t *testing.T, samples []float64, nslices int, expected *statsResults) {
  88         t.Helper()
  89         chunk := len(samples) / nslices
  90         for i := 0; i < nslices; i++ {
  91                 low := i * chunk
  92                 var high int
  93                 if i == nslices-1 {
  94                         high = len(samples) - 1
  95                 } else {
  96                         high = (i + 1) * chunk
  97                 }
  98                 checkSampleDistribution(t, samples[low:high], expected)
  99         }
 100 }
 101
 102 //
 103 // Normal distribution tests
 104 //
 105
 106 func generateNormalSamples(nsamples int, mean, stddev float64, seed int64) []float64 {
 107         r := New(NewSource(seed))
 108         samples := make([]float64, nsamples)
 109         for i := range samples {
 110                 samples[i] = r.NormFloat64()*stddev + mean
 111         }
 112         return samples
 113 }
 114
 115 func testNormalDistribution(t *testing.T, nsamples int, mean, stddev float64, seed int64) {
 116         //fmt.Printf("testing nsamples=%v mean=%v stddev=%v seed=%v\n", nsamples, mean, stddev, seed);
 117
 118         samples := generateNormalSamples(nsamples, mean, stddev, seed)
 119         errorScale := max(1.0, stddev) // Error scales with stddev
 120         expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.08 * errorScale}
 121
 122         // Make sure that the entire set matches the expected distribution.
 123         checkSampleDistribution(t, samples, expected)
 124
 125         // Make sure that each half of the set matches the expected distribution.
 126         checkSampleSliceDistributions(t, samples, 2, expected)
 127
 128         // Make sure that each 7th of the set matches the expected distribution.
 129         checkSampleSliceDistributions(t, samples, 7, expected)
 130 }
 131
 132 // Actual tests
 133
 134 func TestStandardNormalValues(t *testing.T) {
 135         for _, seed := range testSeeds {
 136                 testNormalDistribution(t, numTestSamples, 0, 1, seed)
 137         }
 138 }
 139
 140 func TestNonStandardNormalValues(t *testing.T) {
 141         sdmax := 1000.0
 142         mmax := 1000.0
 143         if testing.Short() {
 144                 sdmax = 5
 145                 mmax = 5
 146         }
 147         for sd := 0.5; sd < sdmax; sd *= 2 {
 148                 for m := 0.5; m < mmax; m *= 2 {
 149                         for _, seed := range testSeeds {
 150                                 testNormalDistribution(t, numTestSamples, m, sd, seed)
 151                                 if testing.Short() {
 152                                         break
 153                                 }
 154                         }
 155                 }
 156         }
 157 }
 158
 159 //
 160 // Exponential distribution tests
 161 //
 162
 163 func generateExponentialSamples(nsamples int, rate float64, seed int64) []float64 {
 164         r := New(NewSource(seed))
 165         samples := make([]float64, nsamples)
 166         for i := range samples {
 167                 samples[i] = r.ExpFloat64() / rate
 168         }
 169         return samples
 170 }
 171
 172 func testExponentialDistribution(t *testing.T, nsamples int, rate float64, seed int64) {
 173         //fmt.Printf("testing nsamples=%v rate=%v seed=%v\n", nsamples, rate, seed);
 174
 175         mean := 1 / rate
 176         stddev := mean
 177
 178         samples := generateExponentialSamples(nsamples, rate, seed)
 179         errorScale := max(1.0, 1/rate) // Error scales with the inverse of the rate
 180         expected := &statsResults{mean, stddev, 0.10 * errorScale, 0.20 * errorScale}
 181
 182         // Make sure that the entire set matches the expected distribution.
 183         checkSampleDistribution(t, samples, expected)
 184
 185         // Make sure that each half of the set matches the expected distribution.
 186         checkSampleSliceDistributions(t, samples, 2, expected)
 187
 188         // Make sure that each 7th of the set matches the expected distribution.
 189         checkSampleSliceDistributions(t, samples, 7, expected)
 190 }
 191
 192 // Actual tests
 193
 194 func TestStandardExponentialValues(t *testing.T) {
 195         for _, seed := range testSeeds {
 196                 testExponentialDistribution(t, numTestSamples, 1, seed)
 197         }
 198 }
 199
 200 func TestNonStandardExponentialValues(t *testing.T) {
 201         for rate := 0.05; rate < 10; rate *= 2 {
 202                 for _, seed := range testSeeds {
 203                         testExponentialDistribution(t, numTestSamples, rate, seed)
 204                         if testing.Short() {
 205                                 break
 206                         }
 207                 }
 208         }
 209 }
 210
 211 //
 212 // Table generation tests
 213 //
 214
 215 func initNorm() (testKn []uint32, testWn, testFn []float32) {
 216         const m1 = 1 << 31
 217         var (
 218                 dn float64 = rn
 219                 tn         = dn
 220                 vn float64 = 9.91256303526217e-3
 221         )
 222
 223         testKn = make([]uint32, 128)
 224         testWn = make([]float32, 128)
 225         testFn = make([]float32, 128)
 226
 227         q := vn / math.Exp(-0.5*dn*dn)
 228         testKn[0] = uint32((dn / q) * m1)
 229         testKn[1] = 0
 230         testWn[0] = float32(q / m1)
 231         testWn[127] = float32(dn / m1)
 232         testFn[0] = 1.0
 233         testFn[127] = float32(math.Exp(-0.5 * dn * dn))
 234         for i := 126; i >= 1; i-- {
 235                 dn = math.Sqrt(-2.0 * math.Log(vn/dn+math.Exp(-0.5*dn*dn)))
 236                 testKn[i+1] = uint32((dn / tn) * m1)
 237                 tn = dn
 238                 testFn[i] = float32(math.Exp(-0.5 * dn * dn))
 239                 testWn[i] = float32(dn / m1)
 240         }
 241         return
 242 }
 243
 244 func initExp() (testKe []uint32, testWe, testFe []float32) {
 245         const m2 = 1 << 32
 246         var (
 247                 de float64 = re
 248                 te         = de
 249                 ve float64 = 3.9496598225815571993e-3
 250         )
 251
 252         testKe = make([]uint32, 256)
 253         testWe = make([]float32, 256)
 254         testFe = make([]float32, 256)
 255
 256         q := ve / math.Exp(-de)
 257         testKe[0] = uint32((de / q) * m2)
 258         testKe[1] = 0
 259         testWe[0] = float32(q / m2)
 260         testWe[255] = float32(de / m2)
 261         testFe[0] = 1.0
 262         testFe[255] = float32(math.Exp(-de))
 263         for i := 254; i >= 1; i-- {
 264                 de = -math.Log(ve/de + math.Exp(-de))
 265                 testKe[i+1] = uint32((de / te) * m2)
 266                 te = de
 267                 testFe[i] = float32(math.Exp(-de))
 268                 testWe[i] = float32(de / m2)
 269         }
 270         return
 271 }
 272
 273 // compareUint32Slices returns the first index where the two slices
 274 // disagree, or <0 if the lengths are the same and all elements
 275 // are identical.
 276 func compareUint32Slices(s1, s2 []uint32) int {
 277         if len(s1) != len(s2) {
 278                 if len(s1) > len(s2) {
 279                         return len(s2) + 1
 280                 }
 281                 return len(s1) + 1
 282         }
 283         for i := range s1 {
 284                 if s1[i] != s2[i] {
 285                         return i
 286                 }
 287         }
 288         return -1
 289 }
 290
 291 // compareFloat32Slices returns the first index where the two slices
 292 // disagree, or <0 if the lengths are the same and all elements
 293 // are identical.
 294 func compareFloat32Slices(s1, s2 []float32) int {
 295         if len(s1) != len(s2) {
 296                 if len(s1) > len(s2) {
 297                         return len(s2) + 1
 298                 }
 299                 return len(s1) + 1
 300         }
 301         for i := range s1 {
 302                 if !nearEqual(float64(s1[i]), float64(s2[i]), 0, 1e-7) {
 303                         return i
 304                 }
 305         }
 306         return -1
 307 }
 308
 309 func TestNormTables(t *testing.T) {
 310         testKn, testWn, testFn := initNorm()
 311         if i := compareUint32Slices(kn[0:], testKn); i >= 0 {
 312                 t.Errorf("kn disagrees at index %v; %v != %v", i, kn[i], testKn[i])
 313         }
 314         if i := compareFloat32Slices(wn[0:], testWn); i >= 0 {
 315                 t.Errorf("wn disagrees at index %v; %v != %v", i, wn[i], testWn[i])
 316         }
 317         if i := compareFloat32Slices(fn[0:], testFn); i >= 0 {
 318                 t.Errorf("fn disagrees at index %v; %v != %v", i, fn[i], testFn[i])
 319         }
 320 }
 321
 322 func TestExpTables(t *testing.T) {
 323         testKe, testWe, testFe := initExp()
 324         if i := compareUint32Slices(ke[0:], testKe); i >= 0 {
 325                 t.Errorf("ke disagrees at index %v; %v != %v", i, ke[i], testKe[i])
 326         }
 327         if i := compareFloat32Slices(we[0:], testWe); i >= 0 {
 328                 t.Errorf("we disagrees at index %v; %v != %v", i, we[i], testWe[i])
 329         }
 330         if i := compareFloat32Slices(fe[0:], testFe); i >= 0 {
 331                 t.Errorf("fe disagrees at index %v; %v != %v", i, fe[i], testFe[i])
 332         }
 333 }
 334
 335 func hasSlowFloatingPoint() bool {
 336         switch runtime.GOARCH {
 337         case "arm":
 338                 return os.Getenv("GOARM") == "5"
 339         case "mips", "mipsle", "mips64", "mips64le":
 340                 // Be conservative and assume that all mips boards
 341                 // have emulated floating point.
 342                 // TODO: detect what it actually has.
 343                 return true
 344         }
 345         return false
 346 }
 347
 348 func TestFloat32(t *testing.T) {
 349         // For issue 6721, the problem came after 7533753 calls, so check 10e6.
 350         num := int(10e6)
 351         // But do the full amount only on builders (not locally).
 352         // But ARM5 floating point emulation is slow (Issue 10749), so
 353         // do less for that builder:
 354         if testing.Short() && (testenv.Builder() == "" || hasSlowFloatingPoint()) {
 355                 num /= 100 // 1.72 seconds instead of 172 seconds
 356         }
 357
 358         r := New(NewSource(1))
 359         for ct := 0; ct < num; ct++ {
 360                 f := r.Float32()
 361                 if f >= 1 {
 362                         t.Fatal("Float32() should be in range [0,1). ct:", ct, "f:", f)
 363                 }
 364         }
 365 }
 366
 367 func TestShuffleSmall(t *testing.T) {
 368         // Check that Shuffle allows n=0 and n=1, but that swap is never called for them.
 369         r := New(NewSource(1))
 370         for n := 0; n <= 1; n++ {
 371                 r.Shuffle(n, func(i, j int) { t.Fatalf("swap called, n=%d i=%d j=%d", n, i, j) })
 372         }
 373 }
 374
 375 // encodePerm converts from a permuted slice of length n, such as Perm generates, to an int in [0, n!).
 376 // See https://en.wikipedia.org/wiki/Lehmer_code.
 377 // encodePerm modifies the input slice.
 378 func encodePerm(s []int) int {
 379         // Convert to Lehmer code.
 380         for i, x := range s {
 381                 r := s[i+1:]
 382                 for j, y := range r {
 383                         if y > x {
 384                                 r[j]--
 385                         }
 386                 }
 387         }
 388         // Convert to int in [0, n!).
 389         m := 0
 390         fact := 1
 391         for i := len(s) - 1; i >= 0; i-- {
 392                 m += s[i] * fact
 393                 fact *= len(s) - i
 394         }
 395         return m
 396 }
 397
 398 // TestUniformFactorial tests several ways of generating a uniform value in [0, n!).
 399 func TestUniformFactorial(t *testing.T) {
 400         r := New(NewSource(testSeeds[0]))
 401         top := 6
 402         if testing.Short() {
 403                 top = 3
 404         }
 405         for n := 3; n <= top; n++ {
 406                 t.Run(fmt.Sprintf("n=%d", n), func(t *testing.T) {
 407                         // Calculate n!.
 408                         nfact := 1
 409                         for i := 2; i <= n; i++ {
 410                                 nfact *= i
 411                         }
 412
 413                         // Test a few different ways to generate a uniform distribution.
 414                         p := make([]int, n) // re-usable slice for Shuffle generator
 415                         tests := [...]struct {
 416                                 name string
 417                                 fn   func() int
 418                         }{
 419                                 {name: "Int32N", fn: func() int { return int(r.Int32N(int32(nfact))) }},
 420                                 {name: "int31n", fn: func() int { return int(Int32NForTest(r, int32(nfact))) }},
 421                                 {name: "Perm", fn: func() int { return encodePerm(r.Perm(n)) }},
 422                                 {name: "Shuffle", fn: func() int {
 423                                         // Generate permutation using Shuffle.
 424                                         for i := range p {
 425                                                 p[i] = i
 426                                         }
 427                                         r.Shuffle(n, func(i, j int) { p[i], p[j] = p[j], p[i] })
 428                                         return encodePerm(p)
 429                                 }},
 430                         }
 431
 432                         for _, test := range tests {
 433                                 t.Run(test.name, func(t *testing.T) {
 434                                         // Gather chi-squared values and check that they follow
 435                                         // the expected normal distribution given n!-1 degrees of freedom.
 436                                         // See https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test and
 437                                         // https://www.johndcook.com/Beautiful_Testing_ch10.pdf.
 438                                         nsamples := 10 * nfact
 439                                         if nsamples < 200 {
 440                                                 nsamples = 200
 441                                         }
 442                                         samples := make([]float64, nsamples)
 443                                         for i := range samples {
 444                                                 // Generate some uniformly distributed values and count their occurrences.
 445                                                 const iters = 1000
 446                                                 counts := make([]int, nfact)
 447                                                 for i := 0; i < iters; i++ {
 448                                                         counts[test.fn()]++
 449                                                 }
 450                                                 // Calculate chi-squared and add to samples.
 451                                                 want := iters / float64(nfact)
 452                                                 var χ2 float64
 453                                                 for _, have := range counts {
 454                                                         err := float64(have) - want
 455                                                         χ2 += err * err
 456                                                 }
 457                                                 χ2 /= want
 458                                                 samples[i] = χ2
 459                                         }
 460
 461                                         // Check that our samples approximate the appropriate normal distribution.
 462                                         dof := float64(nfact - 1)
 463                                         expected := &statsResults{mean: dof, stddev: math.Sqrt(2 * dof)}
 464                                         errorScale := max(1.0, expected.stddev)
 465                                         expected.closeEnough = 0.10 * errorScale
 466                                         expected.maxError = 0.08 // TODO: What is the right value here? See issue 21211.
 467                                         checkSampleDistribution(t, samples, expected)
 468                                 })
 469                         }
 470                 })
 471         }
 472 }
 473
 474 // Benchmarks
 475
 476 func BenchmarkInt64Threadsafe(b *testing.B) {
 477         for n := b.N; n > 0; n-- {
 478                 Int64()
 479         }
 480 }
 481
 482 func BenchmarkInt64ThreadsafeParallel(b *testing.B) {
 483         b.RunParallel(func(pb *testing.PB) {
 484                 for pb.Next() {
 485                         Int64()
 486                 }
 487         })
 488 }
 489
 490 func BenchmarkInt64Unthreadsafe(b *testing.B) {
 491         r := New(NewSource(1))
 492         for n := b.N; n > 0; n-- {
 493                 r.Int64()
 494         }
 495 }
 496
 497 func BenchmarkIntN1000(b *testing.B) {
 498         r := New(NewSource(1))
 499         for n := b.N; n > 0; n-- {
 500                 r.IntN(1000)
 501         }
 502 }
 503
 504 func BenchmarkInt64N1000(b *testing.B) {
 505         r := New(NewSource(1))
 506         for n := b.N; n > 0; n-- {
 507                 r.Int64N(1000)
 508         }
 509 }
 510
 511 func BenchmarkInt32N1000(b *testing.B) {
 512         r := New(NewSource(1))
 513         for n := b.N; n > 0; n-- {
 514                 r.Int32N(1000)
 515         }
 516 }
 517
 518 func BenchmarkFloat32(b *testing.B) {
 519         r := New(NewSource(1))
 520         for n := b.N; n > 0; n-- {
 521                 r.Float32()
 522         }
 523 }
 524
 525 func BenchmarkFloat64(b *testing.B) {
 526         r := New(NewSource(1))
 527         for n := b.N; n > 0; n-- {
 528                 r.Float64()
 529         }
 530 }
 531
 532 func BenchmarkPerm3(b *testing.B) {
 533         r := New(NewSource(1))
 534         for n := b.N; n > 0; n-- {
 535                 r.Perm(3)
 536         }
 537 }
 538
 539 func BenchmarkPerm30(b *testing.B) {
 540         r := New(NewSource(1))
 541         for n := b.N; n > 0; n-- {
 542                 r.Perm(30)
 543         }
 544 }
 545
 546 func BenchmarkPerm30ViaShuffle(b *testing.B) {
 547         r := New(NewSource(1))
 548         for n := b.N; n > 0; n-- {
 549                 p := make([]int, 30)
 550                 for i := range p {
 551                         p[i] = i
 552                 }
 553                 r.Shuffle(30, func(i, j int) { p[i], p[j] = p[j], p[i] })
 554         }
 555 }
 556
 557 // BenchmarkShuffleOverhead uses a minimal swap function
 558 // to measure just the shuffling overhead.
 559 func BenchmarkShuffleOverhead(b *testing.B) {
 560         r := New(NewSource(1))
 561         for n := b.N; n > 0; n-- {
 562                 r.Shuffle(52, func(i, j int) {
 563                         if i < 0 || i >= 52 || j < 0 || j >= 52 {
 564                                 b.Fatalf("bad swap(%d, %d)", i, j)
 565                         }
 566                 })
 567         }
 568 }
 569
 570 func BenchmarkConcurrent(b *testing.B) {
 571         const goroutines = 4
 572         var wg sync.WaitGroup
 573         wg.Add(goroutines)
 574         for i := 0; i < goroutines; i++ {
 575                 go func() {
 576                         defer wg.Done()
 577                         for n := b.N; n > 0; n-- {
 578                                 Int64()
 579                         }
 580                 }()
 581         }
 582         wg.Wait()
 583 }