x/exp/rand: new rand package

This is an approximately compatible new implementation of math/rand.
The differences are:

Source emits uint64s rather than positive int64s.
The method is now
	Uint64() uint64; not
	Int63() int64
There are corresponding new methods on Rand:
	func (r *Rand) Uint64() uint64
	func (r *Rand) Uint64n(n uint64) uint64
The default Source is now an exported type, PCGSource.

The default Source holds 128 bits of state for a 64-bit
result. This has good statistical properties but is slower,
largely because the multiplication step is inefficient.
That can be improved with assembler.

Thus the default Source has a two 64-bit words of state (in
math/rand it has 607 words). It is thus practical to have
millions of Sources in the address space, making it well
suited to lock-free simulations using random numbers.

Benchmarks:

benchmark                        old ns/op     new ns/op     delta
BenchmarkInt63Threadsafe-4       20.0          24.4          +22.00%
BenchmarkInt63Unthreadsafe-4     6.32          13.0          +105.70%
BenchmarkIntn1000-4              16.4          23.8          +45.12%
BenchmarkInt63n1000-4            25.5          23.8          -6.67%
BenchmarkInt31n1000-4            14.2          23.8          +67.61%
BenchmarkFloat32-4               11.8          21.0          +77.97%
BenchmarkFloat64-4               8.76          18.3          +108.90%
BenchmarkPerm3-4                 80.3          94.3          +17.43%
BenchmarkPerm30-4                627           814           +29.82%

The new generator is PCG XSL RR 128/64 (LCG) from
http://www.pcg-random.org/pdf/toms-oneill-pcg-family-v1.02.pdf
It has been tested against the C version and generates the same
output if initialized to the same value.
See also http://www.pcg-random.org/.

TODO: Improve performance, make initialization better.

Independently, this fixes a bug in the bias-prevention code that
appears, in this package, in Uint64n.

This also exports LockedSource:
Update golang/go#21393.

Change-Id: I48a410fade5d78b8ec993cc1210b96b7a9ab462f
Reviewed-on: https://go-review.googlesource.com/10161
Reviewed-by: Rob Pike <[email protected]>

rand/arith128_test.go

@@ -0,0 +1,126 @@
+// Copyright 2017 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 (
+	"math/big"
+	"math/rand"
+	"testing"
+)
+
+var bigMaxUint64 = big.NewInt(0).SetUint64(maxUint64)
+
+func bigInt(xHi, xLo uint64) *big.Int {
+	b := big.NewInt(0).SetUint64(xHi)
+	b.Lsh(b, 64)
+	b.Or(b, big.NewInt(0).SetUint64(xLo))
+	return b
+}
+
+func splitBigInt(b *big.Int) (outHi, outLo uint64) {
+	outHi = big.NewInt(0).Rsh(b, 64).Uint64()
+	outLo = big.NewInt(0).And(b, bigMaxUint64).Uint64()
+	return
+}
+
+func bigMulMod128bits(xHi, xLo, yHi, yLo uint64) (outHi, outLo uint64) {
+	bigX := bigInt(xHi, xLo)
+	bigY := bigInt(yHi, yLo)
+	return splitBigInt(bigX.Mul(bigX, bigY))
+}
+
+func bigAddMod128bits(xHi, xLo, yHi, yLo uint64) (outHi, outLo uint64) {
+	bigX := bigInt(xHi, xLo)
+	bigY := bigInt(yHi, yLo)
+	return splitBigInt(bigX.Add(bigX, bigY))
+}
+
+type arithTest struct {
+	xHi, xLo uint64
+}
+
+const (
+	iLo = increment & maxUint64
+	iHi = (increment >> 64) & maxUint64
+)
+
+var arithTests = []arithTest{
+	{0, 0},
+	{0, 1},
+	{1, 0},
+	{0, maxUint64},
+	{maxUint64, 0},
+	{maxUint64, maxUint64},
+	// Randomly generated 64-bit integers.
+	{3757956613005209672, 17983933746665545631},
+	{511324141977587414, 5626651684620191081},
+	{1534313104606153588, 2415006486399353367},
+	{6873586429837825902, 13854394671140464137},
+	{6617134480561088940, 18421520694158684312},
+}
+
+func TestPCGAdd(t *testing.T) {
+	for i, test := range arithTests {
+		p := &PCGSource{
+			low:  test.xLo,
+			high: test.xHi,
+		}
+		p.add()
+		expectHi, expectLo := bigAddMod128bits(test.xHi, test.xLo, iHi, iLo)
+		if p.low != expectLo || p.high != expectHi {
+			t.Errorf("%d: got hi=%d lo=%d; expect hi=%d lo=%d", i, p.high, p.low, expectHi, expectLo)
+		}
+	}
+}
+
+const (
+	mLo = multiplier & maxUint64
+	mHi = (multiplier >> 64) & maxUint64
+)
+
+func TestPCGMultiply(t *testing.T) {
+	for i, test := range arithTests {
+		p := &PCGSource{
+			low:  test.xLo,
+			high: test.xHi,
+		}
+		p.multiply()
+		expectHi, expectLo := bigMulMod128bits(test.xHi, test.xLo, mHi, mLo)
+		if p.low != expectLo || p.high != expectHi {
+			t.Errorf("%d: got hi=%d lo=%d; expect hi=%d lo=%d", i, p.high, p.low, expectHi, expectLo)
+		}
+	}
+}
+
+func TestPCGMultiplyLong(t *testing.T) {
+	if testing.Short() {
+		return
+	}
+	for i := 0; i < 1e6; i++ {
+		low := rand.Uint64()
+		high := rand.Uint64()
+		p := &PCGSource{
+			low:  low,
+			high: high,
+		}
+		p.multiply()
+		expectHi, expectLo := bigMulMod128bits(high, low, mHi, mLo)
+		if p.low != expectLo || p.high != expectHi {
+			t.Fatalf("%d: (%d,%d): got hi=%d lo=%d; expect hi=%d lo=%d", i, high, low, p.high, p.low, expectHi, expectLo)
+		}
+	}
+}
+
+func BenchmarkPCGMultiply(b *testing.B) {
+	low := rand.Uint64()
+	high := rand.Uint64()
+	p := &PCGSource{
+		low:  low,
+		high: high,
+	}
+	for i := 0; i < b.N; i++ {
+		p.multiply()
+	}
+}

rand/example_test.go

@@ -0,0 +1,102 @@
+// Copyright 2012 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_test
+
+import (
+	"fmt"
+	"os"
+	"text/tabwriter"
+
+	"golang.org/x/exp/rand"
+)
+
+// These tests serve as an example but also make sure we don't change
+// the output of the random number generator when given a fixed seed.
+
+func Example() {
+	rand.Seed(42) // Try changing this number!
+	answers := []string{
+		"It is certain",
+		"It is decidedly so",
+		"Without a doubt",
+		"Yes definitely",
+		"You may rely on it",
+		"As I see it yes",
+		"Most likely",
+		"Outlook good",
+		"Yes",
+		"Signs point to yes",
+		"Reply hazy try again",
+		"Ask again later",
+		"Better not tell you now",
+		"Cannot predict now",
+		"Concentrate and ask again",
+		"Don't count on it",
+		"My reply is no",
+		"My sources say no",
+		"Outlook not so good",
+		"Very doubtful",
+	}
+	fmt.Println("Magic 8-Ball says:", answers[rand.Intn(len(answers))])
+	// Output: Magic 8-Ball says: Most likely
+}
+
+// This example shows the use of each of the methods on a *Rand.
+// The use of the global functions is the same, without the receiver.
+func Example_rand() {
+	// Create and seed the generator.
+	// Typically a non-fixed seed should be used, such as time.Now().UnixNano().
+	// Using a fixed seed will produce the same output on every run.
+	r := rand.New(rand.NewSource(1234))
+
+	// The tabwriter here helps us generate aligned output.
+	w := tabwriter.NewWriter(os.Stdout, 1, 1, 1, ' ', 0)
+	defer w.Flush()
+	show := func(name string, v1, v2, v3 interface{}) {
+		fmt.Fprintf(w, "%s\t%v\t%v\t%v\n", name, v1, v2, v3)
+	}
+
+	// Float32 and Float64 values are in [0, 1).
+	show("Float32", r.Float32(), r.Float32(), r.Float32())
+	show("Float64", r.Float64(), r.Float64(), r.Float64())
+
+	// ExpFloat64 values have an average of 1 but decay exponentially.
+	show("ExpFloat64", r.ExpFloat64(), r.ExpFloat64(), r.ExpFloat64())
+
+	// NormFloat64 values have an average of 0 and a standard deviation of 1.
+	show("NormFloat64", r.NormFloat64(), r.NormFloat64(), r.NormFloat64())
+
+	// Int31, Int63, and Uint32 generate values of the given width.
+	// The Int method (not shown) is like either Int31 or Int63
+	// depending on the size of 'int'.
+	show("Int31", r.Int31(), r.Int31(), r.Int31())
+	show("Int63", r.Int63(), r.Int63(), r.Int63())
+	show("Uint32", r.Uint32(), r.Uint32(), r.Uint32())
+	show("Uint64", r.Uint64(), r.Uint64(), r.Uint64())
+
+	// Intn, Int31n, Int63n and Uint64n limit their output to be < n.
+	// They do so more carefully than using r.Int()%n.
+	show("Intn(10)", r.Intn(10), r.Intn(10), r.Intn(10))
+	show("Int31n(10)", r.Int31n(10), r.Int31n(10), r.Int31n(10))
+	show("Int63n(10)", r.Int63n(10), r.Int63n(10), r.Int63n(10))
+	show("Uint64n(10)", r.Uint64n(10), r.Uint64n(10), r.Uint64n(10))
+
+	// Perm generates a random permutation of the numbers [0, n).
+	show("Perm", r.Perm(5), r.Perm(5), r.Perm(5))
+	// Output:
+	// Float32     0.030719291          0.47512934           0.031019364
+	// Float64     0.6906635660087743   0.9898818576905045   0.2683634639782333
+	// ExpFloat64  1.24979080914592     0.3451975160045876   0.5456817760595064
+	// NormFloat64 0.879221333732727    -0.01508980368383761 -1.962250558270421
+	// Int31       2043816560           1870670250           1334960143
+	// Int63       7860766611810691572  1466711535823962239  3836585920276818709
+	// Uint32      2051241581           751073909            1353986074
+	// Uint64      10802154207635843641 14398820303406316826 11052107950969057042
+	// Intn(10)    3                    0                    1
+	// Int31n(10)  3                    8                    1
+	// Int63n(10)  4                    6                    0
+	// Uint64n(10) 2                    9                    4
+	// Perm        [1 3 4 0 2]          [2 4 0 3 1]          [3 2 0 4 1]
+}