diff --git a/rand/arith128_test.go b/rand/arith128_test.go
new file mode 100644
index 0000000..eed655b
--- /dev/null
+++ b/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()
+	}
+}
diff --git a/rand/example_test.go b/rand/example_test.go
new file mode 100644
index 0000000..3b119fc
--- /dev/null
+++ b/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]
+}
diff --git a/rand/exp.go b/rand/exp.go
new file mode 100644
index 0000000..4bc110f
--- /dev/null
+++ b/rand/exp.go
@@ -0,0 +1,222 @@
+// 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 (
+	"math"
+)
+
+/*
+ * Exponential distribution
+ *
+ * See "The Ziggurat Method for Generating Random Variables"
+ * (Marsaglia & Tsang, 2000)
+ * http://www.jstatsoft.org/v05/i08/paper [pdf]
+ */
+
+const (
+	re = 7.69711747013104972
+)
+
+// 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).
+// To produce a distribution with a different rate parameter,
+// callers can adjust the output using:
+//
+//  sample = ExpFloat64() / desiredRateParameter
+//
+func (r *Rand) ExpFloat64() float64 {
+	for {
+		j := r.Uint32()
+		i := j & 0xFF
+		x := float64(j) * float64(we[i])
+		if j < ke[i] {
+			return x
+		}
+		if i == 0 {
+			return re - math.Log(r.Float64())
+		}
+		if fe[i]+float32(r.Float64())*(fe[i-1]-fe[i]) < float32(math.Exp(-x)) {
+			return x
+		}
+	}
+}
+
+var ke = [256]uint32{
+	0xe290a139, 0x0, 0x9beadebc, 0xc377ac71, 0xd4ddb990,
+	0xde893fb8, 0xe4a8e87c, 0xe8dff16a, 0xebf2deab, 0xee49a6e8,
+	0xf0204efd, 0xf19bdb8e, 0xf2d458bb, 0xf3da104b, 0xf4b86d78,
+	0xf577ad8a, 0xf61de83d, 0xf6afb784, 0xf730a573, 0xf7a37651,
+	0xf80a5bb6, 0xf867189d, 0xf8bb1b4f, 0xf9079062, 0xf94d70ca,
+	0xf98d8c7d, 0xf9c8928a, 0xf9ff175b, 0xfa319996, 0xfa6085f8,
+	0xfa8c3a62, 0xfab5084e, 0xfadb36c8, 0xfaff0410, 0xfb20a6ea,
+	0xfb404fb4, 0xfb5e2951, 0xfb7a59e9, 0xfb95038c, 0xfbae44ba,
+	0xfbc638d8, 0xfbdcf892, 0xfbf29a30, 0xfc0731df, 0xfc1ad1ed,
+	0xfc2d8b02, 0xfc3f6c4d, 0xfc5083ac, 0xfc60ddd1, 0xfc708662,
+	0xfc7f8810, 0xfc8decb4, 0xfc9bbd62, 0xfca9027c, 0xfcb5c3c3,
+	0xfcc20864, 0xfccdd70a, 0xfcd935e3, 0xfce42ab0, 0xfceebace,
+	0xfcf8eb3b, 0xfd02c0a0, 0xfd0c3f59, 0xfd156b7b, 0xfd1e48d6,
+	0xfd26daff, 0xfd2f2552, 0xfd372af7, 0xfd3eeee5, 0xfd4673e7,
+	0xfd4dbc9e, 0xfd54cb85, 0xfd5ba2f2, 0xfd62451b, 0xfd68b415,
+	0xfd6ef1da, 0xfd750047, 0xfd7ae120, 0xfd809612, 0xfd8620b4,
+	0xfd8b8285, 0xfd90bcf5, 0xfd95d15e, 0xfd9ac10b, 0xfd9f8d36,
+	0xfda43708, 0xfda8bf9e, 0xfdad2806, 0xfdb17141, 0xfdb59c46,
+	0xfdb9a9fd, 0xfdbd9b46, 0xfdc170f6, 0xfdc52bd8, 0xfdc8ccac,
+	0xfdcc542d, 0xfdcfc30b, 0xfdd319ef, 0xfdd6597a, 0xfdd98245,
+	0xfddc94e5, 0xfddf91e6, 0xfde279ce, 0xfde54d1f, 0xfde80c52,
+	0xfdeab7de, 0xfded5034, 0xfdefd5be, 0xfdf248e3, 0xfdf4aa06,
+	0xfdf6f984, 0xfdf937b6, 0xfdfb64f4, 0xfdfd818d, 0xfdff8dd0,
+	0xfe018a08, 0xfe03767a, 0xfe05536c, 0xfe07211c, 0xfe08dfc9,
+	0xfe0a8fab, 0xfe0c30fb, 0xfe0dc3ec, 0xfe0f48b1, 0xfe10bf76,
+	0xfe122869, 0xfe1383b4, 0xfe14d17c, 0xfe1611e7, 0xfe174516,
+	0xfe186b2a, 0xfe19843e, 0xfe1a9070, 0xfe1b8fd6, 0xfe1c8289,
+	0xfe1d689b, 0xfe1e4220, 0xfe1f0f26, 0xfe1fcfbc, 0xfe2083ed,
+	0xfe212bc3, 0xfe21c745, 0xfe225678, 0xfe22d95f, 0xfe234ffb,
+	0xfe23ba4a, 0xfe241849, 0xfe2469f2, 0xfe24af3c, 0xfe24e81e,
+	0xfe25148b, 0xfe253474, 0xfe2547c7, 0xfe254e70, 0xfe25485a,
+	0xfe25356a, 0xfe251586, 0xfe24e88f, 0xfe24ae64, 0xfe2466e1,
+	0xfe2411df, 0xfe23af34, 0xfe233eb4, 0xfe22c02c, 0xfe22336b,
+	0xfe219838, 0xfe20ee58, 0xfe20358c, 0xfe1f6d92, 0xfe1e9621,
+	0xfe1daef0, 0xfe1cb7ac, 0xfe1bb002, 0xfe1a9798, 0xfe196e0d,
+	0xfe1832fd, 0xfe16e5fe, 0xfe15869d, 0xfe141464, 0xfe128ed3,
+	0xfe10f565, 0xfe0f478c, 0xfe0d84b1, 0xfe0bac36, 0xfe09bd73,
+	0xfe07b7b5, 0xfe059a40, 0xfe03644c, 0xfe011504, 0xfdfeab88,
+	0xfdfc26e9, 0xfdf98629, 0xfdf6c83b, 0xfdf3ec01, 0xfdf0f04a,
+	0xfdedd3d1, 0xfdea953d, 0xfde7331e, 0xfde3abe9, 0xfddffdfb,
+	0xfddc2791, 0xfdd826cd, 0xfdd3f9a8, 0xfdcf9dfc, 0xfdcb1176,
+	0xfdc65198, 0xfdc15bb3, 0xfdbc2ce2, 0xfdb6c206, 0xfdb117be,
+	0xfdab2a63, 0xfda4f5fd, 0xfd9e7640, 0xfd97a67a, 0xfd908192,
+	0xfd8901f2, 0xfd812182, 0xfd78d98e, 0xfd7022bb, 0xfd66f4ed,
+	0xfd5d4732, 0xfd530f9c, 0xfd48432b, 0xfd3cd59a, 0xfd30b936,
+	0xfd23dea4, 0xfd16349e, 0xfd07a7a3, 0xfcf8219b, 0xfce7895b,
+	0xfcd5c220, 0xfcc2aadb, 0xfcae1d5e, 0xfc97ed4e, 0xfc7fe6d4,
+	0xfc65ccf3, 0xfc495762, 0xfc2a2fc8, 0xfc07ee19, 0xfbe213c1,
+	0xfbb8051a, 0xfb890078, 0xfb5411a5, 0xfb180005, 0xfad33482,
+	0xfa839276, 0xfa263b32, 0xf9b72d1c, 0xf930a1a2, 0xf889f023,
+	0xf7b577d2, 0xf69c650c, 0xf51530f0, 0xf2cb0e3c, 0xeeefb15d,
+	0xe6da6ecf,
+}
+var we = [256]float32{
+	2.0249555e-09, 1.486674e-11, 2.4409617e-11, 3.1968806e-11,
+	3.844677e-11, 4.4228204e-11, 4.9516443e-11, 5.443359e-11,
+	5.905944e-11, 6.344942e-11, 6.7643814e-11, 7.1672945e-11,
+	7.556032e-11, 7.932458e-11, 8.298079e-11, 8.654132e-11,
+	9.0016515e-11, 9.3415074e-11, 9.674443e-11, 1.0001099e-10,
+	1.03220314e-10, 1.06377254e-10, 1.09486115e-10, 1.1255068e-10,
+	1.1557435e-10, 1.1856015e-10, 1.2151083e-10, 1.2442886e-10,
+	1.2731648e-10, 1.3017575e-10, 1.3300853e-10, 1.3581657e-10,
+	1.3860142e-10, 1.4136457e-10, 1.4410738e-10, 1.4683108e-10,
+	1.4953687e-10, 1.5222583e-10, 1.54899e-10, 1.5755733e-10,
+	1.6020171e-10, 1.6283301e-10, 1.6545203e-10, 1.6805951e-10,
+	1.7065617e-10, 1.732427e-10, 1.7581973e-10, 1.7838787e-10,
+	1.8094774e-10, 1.8349985e-10, 1.8604476e-10, 1.8858298e-10,
+	1.9111498e-10, 1.9364126e-10, 1.9616223e-10, 1.9867835e-10,
+	2.0119004e-10, 2.0369768e-10, 2.0620168e-10, 2.087024e-10,
+	2.1120022e-10, 2.136955e-10, 2.1618855e-10, 2.1867974e-10,
+	2.2116936e-10, 2.2365775e-10, 2.261452e-10, 2.2863202e-10,
+	2.311185e-10, 2.3360494e-10, 2.360916e-10, 2.3857874e-10,
+	2.4106667e-10, 2.4355562e-10, 2.4604588e-10, 2.485377e-10,
+	2.5103128e-10, 2.5352695e-10, 2.560249e-10, 2.585254e-10,
+	2.6102867e-10, 2.6353494e-10, 2.6604446e-10, 2.6855745e-10,
+	2.7107416e-10, 2.7359479e-10, 2.761196e-10, 2.7864877e-10,
+	2.8118255e-10, 2.8372119e-10, 2.8626485e-10, 2.888138e-10,
+	2.9136826e-10, 2.939284e-10, 2.9649452e-10, 2.9906677e-10,
+	3.016454e-10, 3.0423064e-10, 3.0682268e-10, 3.0942177e-10,
+	3.1202813e-10, 3.1464195e-10, 3.1726352e-10, 3.19893e-10,
+	3.2253064e-10, 3.251767e-10, 3.2783135e-10, 3.3049485e-10,
+	3.3316744e-10, 3.3584938e-10, 3.3854083e-10, 3.4124212e-10,
+	3.4395342e-10, 3.46675e-10, 3.4940711e-10, 3.5215003e-10,
+	3.5490397e-10, 3.5766917e-10, 3.6044595e-10, 3.6323455e-10,
+	3.660352e-10, 3.6884823e-10, 3.7167386e-10, 3.745124e-10,
+	3.773641e-10, 3.802293e-10, 3.8310827e-10, 3.860013e-10,
+	3.8890866e-10, 3.918307e-10, 3.9476775e-10, 3.9772008e-10,
+	4.0068804e-10, 4.0367196e-10, 4.0667217e-10, 4.09689e-10,
+	4.1272286e-10, 4.1577405e-10, 4.1884296e-10, 4.2192994e-10,
+	4.250354e-10, 4.281597e-10, 4.313033e-10, 4.3446652e-10,
+	4.3764986e-10, 4.408537e-10, 4.4407847e-10, 4.4732465e-10,
+	4.5059267e-10, 4.5388301e-10, 4.571962e-10, 4.6053267e-10,
+	4.6389292e-10, 4.6727755e-10, 4.70687e-10, 4.741219e-10,
+	4.7758275e-10, 4.810702e-10, 4.845848e-10, 4.8812715e-10,
+	4.9169796e-10, 4.9529775e-10, 4.989273e-10, 5.0258725e-10,
+	5.0627835e-10, 5.100013e-10, 5.1375687e-10, 5.1754584e-10,
+	5.21369e-10, 5.2522725e-10, 5.2912136e-10, 5.330522e-10,
+	5.370208e-10, 5.4102806e-10, 5.45075e-10, 5.491625e-10,
+	5.532918e-10, 5.5746385e-10, 5.616799e-10, 5.6594107e-10,
+	5.7024857e-10, 5.746037e-10, 5.7900773e-10, 5.834621e-10,
+	5.8796823e-10, 5.925276e-10, 5.971417e-10, 6.018122e-10,
+	6.065408e-10, 6.113292e-10, 6.1617933e-10, 6.2109295e-10,
+	6.260722e-10, 6.3111916e-10, 6.3623595e-10, 6.4142497e-10,
+	6.4668854e-10, 6.5202926e-10, 6.5744976e-10, 6.6295286e-10,
+	6.6854156e-10, 6.742188e-10, 6.79988e-10, 6.858526e-10,
+	6.9181616e-10, 6.978826e-10, 7.04056e-10, 7.103407e-10,
+	7.167412e-10, 7.2326256e-10, 7.2990985e-10, 7.366886e-10,
+	7.4360473e-10, 7.5066453e-10, 7.5787476e-10, 7.6524265e-10,
+	7.7277595e-10, 7.80483e-10, 7.883728e-10, 7.9645507e-10,
+	8.047402e-10, 8.1323964e-10, 8.219657e-10, 8.309319e-10,
+	8.401528e-10, 8.496445e-10, 8.594247e-10, 8.6951274e-10,
+	8.799301e-10, 8.9070046e-10, 9.018503e-10, 9.134092e-10,
+	9.254101e-10, 9.378904e-10, 9.508923e-10, 9.644638e-10,
+	9.786603e-10, 9.935448e-10, 1.0091913e-09, 1.025686e-09,
+	1.0431306e-09, 1.0616465e-09, 1.08138e-09, 1.1025096e-09,
+	1.1252564e-09, 1.1498986e-09, 1.1767932e-09, 1.206409e-09,
+	1.2393786e-09, 1.276585e-09, 1.3193139e-09, 1.3695435e-09,
+	1.4305498e-09, 1.508365e-09, 1.6160854e-09, 1.7921248e-09,
+}
+var fe = [256]float32{
+	1, 0.9381437, 0.90046996, 0.87170434, 0.8477855, 0.8269933,
+	0.8084217, 0.7915276, 0.77595687, 0.7614634, 0.7478686,
+	0.7350381, 0.72286767, 0.71127474, 0.70019263, 0.6895665,
+	0.67935055, 0.6695063, 0.66000086, 0.65080583, 0.6418967,
+	0.63325197, 0.6248527, 0.6166822, 0.60872537, 0.60096896,
+	0.5934009, 0.58601034, 0.5787874, 0.57172304, 0.5648092,
+	0.5580383, 0.5514034, 0.5448982, 0.5385169, 0.53225386,
+	0.5261042, 0.52006316, 0.5141264, 0.50828975, 0.5025495,
+	0.496902, 0.49134386, 0.485872, 0.48048335, 0.4751752,
+	0.46994483, 0.46478975, 0.45970762, 0.45469615, 0.44975325,
+	0.44487688, 0.44006512, 0.43531612, 0.43062815, 0.42599955,
+	0.42142874, 0.4169142, 0.41245446, 0.40804818, 0.403694,
+	0.3993907, 0.39513698, 0.39093173, 0.38677382, 0.38266218,
+	0.37859577, 0.37457356, 0.37059465, 0.3666581, 0.362763,
+	0.35890847, 0.35509375, 0.351318, 0.3475805, 0.34388044,
+	0.34021714, 0.3365899, 0.33299807, 0.32944095, 0.32591796,
+	0.3224285, 0.3189719, 0.31554767, 0.31215525, 0.30879408,
+	0.3054636, 0.3021634, 0.29889292, 0.2956517, 0.29243928,
+	0.28925523, 0.28609908, 0.28297043, 0.27986884, 0.27679393,
+	0.2737453, 0.2707226, 0.2677254, 0.26475343, 0.26180625,
+	0.25888354, 0.25598502, 0.2531103, 0.25025907, 0.24743107,
+	0.24462597, 0.24184346, 0.23908329, 0.23634516, 0.23362878,
+	0.23093392, 0.2282603, 0.22560766, 0.22297576, 0.22036438,
+	0.21777324, 0.21520215, 0.21265087, 0.21011916, 0.20760682,
+	0.20511365, 0.20263945, 0.20018397, 0.19774707, 0.19532852,
+	0.19292815, 0.19054577, 0.1881812, 0.18583426, 0.18350479,
+	0.1811926, 0.17889754, 0.17661946, 0.17435817, 0.17211354,
+	0.1698854, 0.16767362, 0.16547804, 0.16329853, 0.16113494,
+	0.15898713, 0.15685499, 0.15473837, 0.15263714, 0.15055119,
+	0.14848037, 0.14642459, 0.14438373, 0.14235765, 0.14034624,
+	0.13834943, 0.13636707, 0.13439907, 0.13244532, 0.13050574,
+	0.1285802, 0.12666863, 0.12477092, 0.12288698, 0.12101672,
+	0.119160056, 0.1173169, 0.115487166, 0.11367077, 0.11186763,
+	0.11007768, 0.10830083, 0.10653701, 0.10478614, 0.10304816,
+	0.101323, 0.09961058, 0.09791085, 0.09622374, 0.09454919,
+	0.09288713, 0.091237515, 0.08960028, 0.087975375, 0.08636274,
+	0.08476233, 0.083174095, 0.081597984, 0.08003395, 0.07848195,
+	0.076941945, 0.07541389, 0.07389775, 0.072393484, 0.07090106,
+	0.069420435, 0.06795159, 0.066494495, 0.06504912, 0.063615434,
+	0.062193416, 0.060783047, 0.059384305, 0.057997175,
+	0.05662164, 0.05525769, 0.053905312, 0.052564494, 0.051235236,
+	0.049917534, 0.048611384, 0.047316793, 0.046033762, 0.0447623,
+	0.043502413, 0.042254124, 0.041017443, 0.039792392,
+	0.038578995, 0.037377283, 0.036187284, 0.035009038,
+	0.033842582, 0.032687962, 0.031545233, 0.030414443, 0.02929566,
+	0.02818895, 0.027094385, 0.026012046, 0.024942026, 0.023884421,
+	0.022839336, 0.021806888, 0.020787204, 0.019780423, 0.0187867,
+	0.0178062, 0.016839107, 0.015885621, 0.014945968, 0.014020392,
+	0.013109165, 0.012212592, 0.011331013, 0.01046481, 0.009614414,
+	0.008780315, 0.007963077, 0.0071633533, 0.006381906,
+	0.0056196423, 0.0048776558, 0.004157295, 0.0034602648,
+	0.0027887989, 0.0021459677, 0.0015362998, 0.0009672693,
+	0.00045413437,
+}
diff --git a/rand/modulo_test.go b/rand/modulo_test.go
new file mode 100644
index 0000000..da963c7
--- /dev/null
+++ b/rand/modulo_test.go
@@ -0,0 +1,50 @@
+// 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.
+
+// This file validates that the calculation in Uint64n corrects for
+// possible bias.
+
+package rand
+
+import (
+	"testing"
+)
+
+// modSource is used to probe the upper region of uint64 space. It
+// generates values sequentially in [maxUint64-15,maxUint64]. With
+// modEdge == 15 and maxUint64 == 1<<64-1 == 18446744073709551615,
+// this means that Uint64n(10) will repeatedly probe the top range.
+// We thus expect a bias to result unless the calculation in Uint64n
+// gets the edge condition right. We test this by calling Uint64n 100
+// times; the results should be perfectly evenly distributed across
+// [0,10).
+type modSource uint64
+
+const modEdge = 15
+
+func (m *modSource) Seed(uint64) {}
+
+// Uint64 returns a non-pseudo-random 64-bit unsigned integer as a uint64.
+func (m *modSource) Uint64() uint64 {
+	if *m > modEdge {
+		*m = 0
+	}
+	r := maxUint64 - *m
+	*m++
+	return uint64(r)
+}
+
+func TestUint64Modulo(t *testing.T) {
+	var src modSource
+	rng := New(&src)
+	var result [10]uint64
+	for i := 0; i < 100; i++ {
+		result[rng.Uint64n(10)]++
+	}
+	for _, r := range result {
+		if r != 10 {
+			t.Fatal(result)
+		}
+	}
+}
diff --git a/rand/normal.go b/rand/normal.go
new file mode 100644
index 0000000..ba4ea54
--- /dev/null
+++ b/rand/normal.go
@@ -0,0 +1,157 @@
+// 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 (
+	"math"
+)
+
+/*
+ * Normal distribution
+ *
+ * See "The Ziggurat Method for Generating Random Variables"
+ * (Marsaglia & Tsang, 2000)
+ * http://www.jstatsoft.org/v05/i08/paper [pdf]
+ */
+
+const (
+	rn = 3.442619855899
+)
+
+func absInt32(i int32) uint32 {
+	if i < 0 {
+		return uint32(-i)
+	}
+	return uint32(i)
+}
+
+// NormFloat64 returns a normally distributed float64 in the range
+// [-math.MaxFloat64, +math.MaxFloat64] with
+// standard normal distribution (mean = 0, stddev = 1).
+// To produce a different normal distribution, callers can
+// adjust the output using:
+//
+//  sample = NormFloat64() * desiredStdDev + desiredMean
+//
+func (r *Rand) NormFloat64() float64 {
+	for {
+		j := int32(r.Uint32()) // Possibly negative
+		i := j & 0x7F
+		x := float64(j) * float64(wn[i])
+		if absInt32(j) < kn[i] {
+			// This case should be hit better than 99% of the time.
+			return x
+		}
+
+		if i == 0 {
+			// This extra work is only required for the base strip.
+			for {
+				x = -math.Log(r.Float64()) * (1.0 / rn)
+				y := -math.Log(r.Float64())
+				if y+y >= x*x {
+					break
+				}
+			}
+			if j > 0 {
+				return rn + x
+			}
+			return -rn - x
+		}
+		if fn[i]+float32(r.Float64())*(fn[i-1]-fn[i]) < float32(math.Exp(-.5*x*x)) {
+			return x
+		}
+	}
+}
+
+var kn = [128]uint32{
+	0x76ad2212, 0x0, 0x600f1b53, 0x6ce447a6, 0x725b46a2,
+	0x7560051d, 0x774921eb, 0x789a25bd, 0x799045c3, 0x7a4bce5d,
+	0x7adf629f, 0x7b5682a6, 0x7bb8a8c6, 0x7c0ae722, 0x7c50cce7,
+	0x7c8cec5b, 0x7cc12cd6, 0x7ceefed2, 0x7d177e0b, 0x7d3b8883,
+	0x7d5bce6c, 0x7d78dd64, 0x7d932886, 0x7dab0e57, 0x7dc0dd30,
+	0x7dd4d688, 0x7de73185, 0x7df81cea, 0x7e07c0a3, 0x7e163efa,
+	0x7e23b587, 0x7e303dfd, 0x7e3beec2, 0x7e46db77, 0x7e51155d,
+	0x7e5aabb3, 0x7e63abf7, 0x7e6c222c, 0x7e741906, 0x7e7b9a18,
+	0x7e82adfa, 0x7e895c63, 0x7e8fac4b, 0x7e95a3fb, 0x7e9b4924,
+	0x7ea0a0ef, 0x7ea5b00d, 0x7eaa7ac3, 0x7eaf04f3, 0x7eb3522a,
+	0x7eb765a5, 0x7ebb4259, 0x7ebeeafd, 0x7ec2620a, 0x7ec5a9c4,
+	0x7ec8c441, 0x7ecbb365, 0x7ece78ed, 0x7ed11671, 0x7ed38d62,
+	0x7ed5df12, 0x7ed80cb4, 0x7eda175c, 0x7edc0005, 0x7eddc78e,
+	0x7edf6ebf, 0x7ee0f647, 0x7ee25ebe, 0x7ee3a8a9, 0x7ee4d473,
+	0x7ee5e276, 0x7ee6d2f5, 0x7ee7a620, 0x7ee85c10, 0x7ee8f4cd,
+	0x7ee97047, 0x7ee9ce59, 0x7eea0eca, 0x7eea3147, 0x7eea3568,
+	0x7eea1aab, 0x7ee9e071, 0x7ee98602, 0x7ee90a88, 0x7ee86d08,
+	0x7ee7ac6a, 0x7ee6c769, 0x7ee5bc9c, 0x7ee48a67, 0x7ee32efc,
+	0x7ee1a857, 0x7edff42f, 0x7ede0ffa, 0x7edbf8d9, 0x7ed9ab94,
+	0x7ed7248d, 0x7ed45fae, 0x7ed1585c, 0x7ece095f, 0x7eca6ccb,
+	0x7ec67be2, 0x7ec22eee, 0x7ebd7d1a, 0x7eb85c35, 0x7eb2c075,
+	0x7eac9c20, 0x7ea5df27, 0x7e9e769f, 0x7e964c16, 0x7e8d44ba,
+	0x7e834033, 0x7e781728, 0x7e6b9933, 0x7e5d8a1a, 0x7e4d9ded,
+	0x7e3b737a, 0x7e268c2f, 0x7e0e3ff5, 0x7df1aa5d, 0x7dcf8c72,
+	0x7da61a1e, 0x7d72a0fb, 0x7d30e097, 0x7cd9b4ab, 0x7c600f1a,
+	0x7ba90bdc, 0x7a722176, 0x77d664e5,
+}
+var wn = [128]float32{
+	1.7290405e-09, 1.2680929e-10, 1.6897518e-10, 1.9862688e-10,
+	2.2232431e-10, 2.4244937e-10, 2.601613e-10, 2.7611988e-10,
+	2.9073963e-10, 3.042997e-10, 3.1699796e-10, 3.289802e-10,
+	3.4035738e-10, 3.5121603e-10, 3.616251e-10, 3.7164058e-10,
+	3.8130857e-10, 3.9066758e-10, 3.9975012e-10, 4.08584e-10,
+	4.1719309e-10, 4.2559822e-10, 4.338176e-10, 4.418672e-10,
+	4.497613e-10, 4.5751258e-10, 4.651324e-10, 4.7263105e-10,
+	4.8001775e-10, 4.87301e-10, 4.944885e-10, 5.015873e-10,
+	5.0860405e-10, 5.155446e-10, 5.2241467e-10, 5.2921934e-10,
+	5.359635e-10, 5.426517e-10, 5.4928817e-10, 5.5587696e-10,
+	5.624219e-10, 5.6892646e-10, 5.753941e-10, 5.818282e-10,
+	5.882317e-10, 5.946077e-10, 6.00959e-10, 6.072884e-10,
+	6.135985e-10, 6.19892e-10, 6.2617134e-10, 6.3243905e-10,
+	6.386974e-10, 6.449488e-10, 6.511956e-10, 6.5744005e-10,
+	6.6368433e-10, 6.699307e-10, 6.7618144e-10, 6.824387e-10,
+	6.8870465e-10, 6.949815e-10, 7.012715e-10, 7.075768e-10,
+	7.1389966e-10, 7.202424e-10, 7.266073e-10, 7.329966e-10,
+	7.394128e-10, 7.4585826e-10, 7.5233547e-10, 7.58847e-10,
+	7.653954e-10, 7.719835e-10, 7.7861395e-10, 7.852897e-10,
+	7.920138e-10, 7.987892e-10, 8.0561924e-10, 8.125073e-10,
+	8.194569e-10, 8.2647167e-10, 8.3355556e-10, 8.407127e-10,
+	8.479473e-10, 8.55264e-10, 8.6266755e-10, 8.7016316e-10,
+	8.777562e-10, 8.8545243e-10, 8.932582e-10, 9.0117996e-10,
+	9.09225e-10, 9.174008e-10, 9.2571584e-10, 9.341788e-10,
+	9.427997e-10, 9.515889e-10, 9.605579e-10, 9.697193e-10,
+	9.790869e-10, 9.88676e-10, 9.985036e-10, 1.0085882e-09,
+	1.0189509e-09, 1.0296151e-09, 1.0406069e-09, 1.0519566e-09,
+	1.063698e-09, 1.0758702e-09, 1.0885183e-09, 1.1016947e-09,
+	1.1154611e-09, 1.1298902e-09, 1.1450696e-09, 1.1611052e-09,
+	1.1781276e-09, 1.1962995e-09, 1.2158287e-09, 1.2369856e-09,
+	1.2601323e-09, 1.2857697e-09, 1.3146202e-09, 1.347784e-09,
+	1.3870636e-09, 1.4357403e-09, 1.5008659e-09, 1.6030948e-09,
+}
+var fn = [128]float32{
+	1, 0.9635997, 0.9362827, 0.9130436, 0.89228165, 0.87324303,
+	0.8555006, 0.8387836, 0.8229072, 0.8077383, 0.793177,
+	0.7791461, 0.7655842, 0.7524416, 0.73967725, 0.7272569,
+	0.7151515, 0.7033361, 0.69178915, 0.68049186, 0.6694277,
+	0.658582, 0.6479418, 0.63749546, 0.6272325, 0.6171434,
+	0.6072195, 0.5974532, 0.58783704, 0.5783647, 0.56903,
+	0.5598274, 0.5507518, 0.54179835, 0.5329627, 0.52424055,
+	0.5156282, 0.50712204, 0.49871865, 0.49041483, 0.48220766,
+	0.4740943, 0.46607214, 0.4581387, 0.45029163, 0.44252872,
+	0.43484783, 0.427247, 0.41972435, 0.41227803, 0.40490642,
+	0.39760786, 0.3903808, 0.3832238, 0.37613547, 0.36911446,
+	0.3621595, 0.35526937, 0.34844297, 0.34167916, 0.33497685,
+	0.3283351, 0.3217529, 0.3152294, 0.30876362, 0.30235484,
+	0.29600215, 0.28970486, 0.2834622, 0.2772735, 0.27113807,
+	0.2650553, 0.25902456, 0.2530453, 0.24711695, 0.241239,
+	0.23541094, 0.22963232, 0.2239027, 0.21822165, 0.21258877,
+	0.20700371, 0.20146611, 0.19597565, 0.19053204, 0.18513499,
+	0.17978427, 0.17447963, 0.1692209, 0.16400786, 0.15884037,
+	0.15371831, 0.14864157, 0.14361008, 0.13862377, 0.13368265,
+	0.12878671, 0.12393598, 0.119130544, 0.11437051, 0.10965602,
+	0.104987256, 0.10036444, 0.095787846, 0.0912578, 0.08677467,
+	0.0823389, 0.077950984, 0.073611505, 0.06932112, 0.06508058,
+	0.06089077, 0.056752663, 0.0526674, 0.048636295, 0.044660863,
+	0.040742867, 0.03688439, 0.033087887, 0.029356318,
+	0.025693292, 0.022103304, 0.018592102, 0.015167298,
+	0.011839478, 0.008624485, 0.005548995, 0.0026696292,
+}
diff --git a/rand/race_test.go b/rand/race_test.go
new file mode 100644
index 0000000..376224f
--- /dev/null
+++ b/rand/race_test.go
@@ -0,0 +1,48 @@
+// Copyright 2016 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 (
+	"sync"
+	"testing"
+)
+
+// TestConcurrent exercises the rand API concurrently, triggering situations
+// where the race detector is likely to detect issues.
+func TestConcurrent(t *testing.T) {
+	const (
+		numRoutines = 10
+		numCycles   = 10
+	)
+	var wg sync.WaitGroup
+	defer wg.Wait()
+	wg.Add(numRoutines)
+	for i := 0; i < numRoutines; i++ {
+		go func(i int) {
+			defer wg.Done()
+			buf := make([]byte, 997)
+			for j := 0; j < numCycles; j++ {
+				var seed uint64
+				seed += uint64(ExpFloat64())
+				seed += uint64(Float32())
+				seed += uint64(Float64())
+				seed += uint64(Intn(Int()))
+				seed += uint64(Int31n(Int31()))
+				seed += uint64(Int63n(Int63()))
+				seed += uint64(NormFloat64())
+				seed += uint64(Uint32())
+				seed += uint64(Uint64())
+				for _, p := range Perm(10) {
+					seed += uint64(p)
+				}
+				Read(buf)
+				for _, b := range buf {
+					seed += uint64(b)
+				}
+				Seed(uint64(i*j) * seed)
+			}
+		}(i)
+	}
+}
diff --git a/rand/rand.go b/rand/rand.go
new file mode 100644
index 0000000..47ccd1f
--- /dev/null
+++ b/rand/rand.go
@@ -0,0 +1,334 @@
+// 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 "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 {
+	var rng PCGSource
+	rng.Seed(seed)
+	return &rng
+}
+
+// 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}
+}
+
+// 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 {
+	return float64(r.Uint64n(1<<53)) / (1 << 53)
+	// 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
+}
+
+// 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.
+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.Uint64, &r.readVal, &r.readPos)
+}
+
+func read(p []byte, uint64 func() uint64, readVal *uint64, readPos *int8) (n int, err error) {
+	pos := *readPos
+	val := *readVal
+	for n = 0; n < len(p); n++ {
+		if pos == 0 {
+			val = 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)})
+
+// 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) }
+
+// 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.
+// It is just a standard Source with its operations protected by a sync.Mutex.
+type LockedSource struct {
+	lk  sync.Mutex
+	src Source
+}
+
+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.Uint64, readVal, readPos)
+	s.lk.Unlock()
+	return
+}
diff --git a/rand/rand_test.go b/rand/rand_test.go
new file mode 100644
index 0000000..c81293f
--- /dev/null
+++ b/rand/rand_test.go
@@ -0,0 +1,548 @@
+// 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"
+	"testing"
+	"testing/iotest"
+)
+
+const (
+	numTestSamples = 10000
+)
+
+type statsResults struct {
+	mean        float64
+	stddev      float64
+	closeEnough float64
+	maxError    float64
+}
+
+func max(a, b float64) float64 {
+	if a > b {
+		return a
+	}
+	return b
+}
+
+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 (this *statsResults) checkSimilarDistribution(expected *statsResults) error {
+	if !nearEqual(this.mean, expected.mean, expected.closeEnough, expected.maxError) {
+		s := fmt.Sprintf("mean %v != %v (allowed error %v, %v)", this.mean, expected.mean, expected.closeEnough, expected.maxError)
+		fmt.Println(s)
+		return errors.New(s)
+	}
+	if !nearEqual(this.stddev, expected.stddev, 0, expected.maxError) {
+		s := fmt.Sprintf("stddev %v != %v (allowed error %v, %v)", this.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.Errorf(err.Error())
+	}
+}
+
+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
+			}
+		}
+	}
+}
+
+//
+// Table generation tests
+//
+
+func initNorm() (testKn []uint32, testWn, testFn []float32) {
+	const m1 = 1 << 31
+	var (
+		dn float64 = rn
+		tn         = dn
+		vn float64 = 9.91256303526217e-3
+	)
+
+	testKn = make([]uint32, 128)
+	testWn = make([]float32, 128)
+	testFn = make([]float32, 128)
+
+	q := vn / math.Exp(-0.5*dn*dn)
+	testKn[0] = uint32((dn / q) * m1)
+	testKn[1] = 0
+	testWn[0] = float32(q / m1)
+	testWn[127] = float32(dn / m1)
+	testFn[0] = 1.0
+	testFn[127] = float32(math.Exp(-0.5 * dn * dn))
+	for i := 126; i >= 1; i-- {
+		dn = math.Sqrt(-2.0 * math.Log(vn/dn+math.Exp(-0.5*dn*dn)))
+		testKn[i+1] = uint32((dn / tn) * m1)
+		tn = dn
+		testFn[i] = float32(math.Exp(-0.5 * dn * dn))
+		testWn[i] = float32(dn / m1)
+	}
+	return
+}
+
+func initExp() (testKe []uint32, testWe, testFe []float32) {
+	const m2 = 1 << 32
+	var (
+		de float64 = re
+		te         = de
+		ve float64 = 3.9496598225815571993e-3
+	)
+
+	testKe = make([]uint32, 256)
+	testWe = make([]float32, 256)
+	testFe = make([]float32, 256)
+
+	q := ve / math.Exp(-de)
+	testKe[0] = uint32((de / q) * m2)
+	testKe[1] = 0
+	testWe[0] = float32(q / m2)
+	testWe[255] = float32(de / m2)
+	testFe[0] = 1.0
+	testFe[255] = float32(math.Exp(-de))
+	for i := 254; i >= 1; i-- {
+		de = -math.Log(ve/de + math.Exp(-de))
+		testKe[i+1] = uint32((de / te) * m2)
+		te = de
+		testFe[i] = float32(math.Exp(-de))
+		testWe[i] = float32(de / m2)
+	}
+	return
+}
+
+// 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 TestNormTables(t *testing.T) {
+	testKn, testWn, testFn := initNorm()
+	if i := compareUint32Slices(kn[0:], testKn); i >= 0 {
+		t.Errorf("kn disagrees at index %v; %v != %v", i, kn[i], testKn[i])
+	}
+	if i := compareFloat32Slices(wn[0:], testWn); i >= 0 {
+		t.Errorf("wn disagrees at index %v; %v != %v", i, wn[i], testWn[i])
+	}
+	if i := compareFloat32Slices(fn[0:], testFn); i >= 0 {
+		t.Errorf("fn disagrees at index %v; %v != %v", i, fn[i], testFn[i])
+	}
+}
+
+func TestExpTables(t *testing.T) {
+	testKe, testWe, testFe := initExp()
+	if i := compareUint32Slices(ke[0:], testKe); i >= 0 {
+		t.Errorf("ke disagrees at index %v; %v != %v", i, ke[i], testKe[i])
+	}
+	if i := compareFloat32Slices(we[0:], testWe); i >= 0 {
+		t.Errorf("we disagrees at index %v; %v != %v", i, we[i], testWe[i])
+	}
+	if i := compareFloat32Slices(fe[0:], testFe); i >= 0 {
+		t.Errorf("fe disagrees at index %v; %v != %v", i, fe[i], testFe[i])
+	}
+}
+
+func hasSlowFloatingPoint() bool {
+	switch runtime.GOARCH {
+	case "arm":
+		return os.Getenv("GOARM") == "5"
+	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)
+	}
+}
+
+// Benchmarks
+
+func BenchmarkSource(b *testing.B) {
+	rng := NewSource(0)
+	for n := b.N; n > 0; n-- {
+		rng.Uint64()
+	}
+}
+
+func BenchmarkInt63Threadsafe(b *testing.B) {
+	for n := b.N; n > 0; n-- {
+		Int63()
+	}
+}
+
+func BenchmarkInt63Unthreadsafe(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Int63()
+	}
+}
+
+func BenchmarkIntn1000(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Intn(1000)
+	}
+}
+
+func BenchmarkInt63n1000(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Int63n(1000)
+	}
+}
+
+func BenchmarkInt31n1000(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Int31n(1000)
+	}
+}
+
+func BenchmarkFloat32(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Float32()
+	}
+}
+
+func BenchmarkFloat64(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Float64()
+	}
+}
+
+func BenchmarkPerm3(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Perm(3)
+	}
+}
+
+func BenchmarkPerm30(b *testing.B) {
+	r := New(NewSource(1))
+	for n := b.N; n > 0; n-- {
+		r.Perm(30)
+	}
+}
+
+func BenchmarkRead3(b *testing.B) {
+	r := New(NewSource(1))
+	buf := make([]byte, 3)
+	b.ResetTimer()
+	for n := b.N; n > 0; n-- {
+		r.Read(buf)
+	}
+}
+
+func BenchmarkRead64(b *testing.B) {
+	r := New(NewSource(1))
+	buf := make([]byte, 64)
+	b.ResetTimer()
+	for n := b.N; n > 0; n-- {
+		r.Read(buf)
+	}
+}
+
+func BenchmarkRead1000(b *testing.B) {
+	r := New(NewSource(1))
+	buf := make([]byte, 1000)
+	b.ResetTimer()
+	for n := b.N; n > 0; n-- {
+		r.Read(buf)
+	}
+}
diff --git a/rand/regress_test.go b/rand/regress_test.go
new file mode 100644
index 0000000..ace8107
--- /dev/null
+++ b/rand/regress_test.go
@@ -0,0 +1,490 @@
+// Copyright 2014 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.
+
+// Test that random number sequences generated by a specific seed
+// do not change from version to version.
+//
+// If the generator changes, the golden outputs need updating, and
+// client programs may break. Although the desire for compatibility
+// is not as stringent as in the original math/rand package,
+// when possible avoid changing the generator.
+
+package rand_test
+
+import (
+	"flag"
+	"fmt"
+	"reflect"
+	"testing"
+
+	. "golang.org/x/exp/rand"
+)
+
+var printgolden = flag.Bool("printgolden", false, "print golden results for regression test")
+
+// TestSource verifies that the output of the default Source is locked down.
+func TestSourceRegress(t *testing.T) {
+	src := NewSource(1)
+	var got [20]uint64
+	for i := range got {
+		got[i] = src.Uint64()
+	}
+	want := [20]uint64{
+		0x34e936394905d167,
+		0x817c0ef62fe4c731,
+		0x937987e6e24f5a40,
+		0x0c0a8307fe226199,
+		0xf96363568d8bab56,
+		0xbaef3af36bd02620,
+		0x8f18e416eb6b936b,
+		0x05a43fc149f3a67a,
+		0xdab012eb3ce01697,
+		0xf76c495a133c6aa9,
+		0x304b24c5040ce457,
+		0x47d77e0abb413159,
+		0x52a810fa9e452f04,
+		0x2d24b66380cf4780,
+		0x5ec7691b92018ef5,
+		0x5076dfa749261ea0,
+		0xac8f11ad3941d213,
+		0x13fa8d67de91db25,
+		0xb50883a9893274eb,
+		0xeb8f59263f9109ac,
+	}
+	if got != want {
+		t.Errorf("got:\n\t%#016x\nwant:\n\t%#016x", got, want)
+		if *printgolden {
+			for _, x := range got {
+				fmt.Printf("\t\t%#016x,\n", x)
+			}
+		}
+	}
+}
+
+// TestRegress validates that the output stream is locked down, for instance so
+// optimizations do not change the output. It iterates over methods of the
+// Rand type to find functions to evaluate and checks the first 20 results
+// against the golden results.
+func TestRegress(t *testing.T) {
+	var int32s = []int32{1, 10, 32, 1 << 20, 1<<20 + 1, 1000000000, 1 << 30, 1<<31 - 2, 1<<31 - 1}
+	var int64s = []int64{1, 10, 32, 1 << 20, 1<<20 + 1, 1000000000, 1 << 30, 1<<31 - 2, 1<<31 - 1, 1000000000000000000, 1 << 60, 1<<63 - 2, 1<<63 - 1}
+	var uint64s = []uint64{1, 10, 32, 1 << 20, 1<<20 + 1, 1000000000, 1 << 30, 1<<31 - 2, 1<<31 - 1, 1000000000000000000, 1 << 60, 1<<64 - 2, 1<<64 - 1}
+	var permSizes = []int{0, 1, 5, 8, 9, 10, 16}
+	var readBufferSizes = []int{1, 7, 8, 9, 10}
+	r := New(NewSource(0))
+
+	rv := reflect.ValueOf(r)
+	n := rv.NumMethod()
+	p := 0
+	if *printgolden {
+		fmt.Printf("var regressGolden = []interface{}{\n")
+	}
+	for i := 0; i < n; i++ {
+		m := rv.Type().Method(i)
+		mv := rv.Method(i)
+		mt := mv.Type()
+		if mt.NumOut() == 0 {
+			continue
+		}
+		r.Seed(0)
+		if *printgolden && i > 0 {
+			fmt.Println()
+		}
+		for repeat := 0; repeat < 20; repeat++ {
+			var args []reflect.Value
+			var argstr string
+			if mt.NumIn() == 1 {
+				var x interface{}
+				switch mt.In(0).Kind() {
+				default:
+					t.Fatalf("unexpected argument type for r.%s", m.Name)
+
+				case reflect.Int:
+					if m.Name == "Perm" {
+						x = permSizes[repeat%len(permSizes)]
+						break
+					}
+					big := int64s[repeat%len(int64s)]
+					if int64(int(big)) != big {
+						r.Int63n(big) // what would happen on 64-bit machine, to keep stream in sync
+						if *printgolden {
+							fmt.Printf("\tskipped, // must run printgolden on 64-bit machine\n")
+						}
+						p++
+						continue
+					}
+					x = int(big)
+
+				case reflect.Int32:
+					x = int32s[repeat%len(int32s)]
+
+				case reflect.Int64:
+					x = int64s[repeat%len(int64s)]
+
+				case reflect.Uint64:
+					x = uint64s[repeat%len(uint64s)]
+
+				case reflect.Slice:
+					if m.Name == "Read" {
+						n := readBufferSizes[repeat%len(readBufferSizes)]
+						x = make([]byte, n)
+					}
+				}
+				argstr = fmt.Sprint(x)
+				args = append(args, reflect.ValueOf(x))
+			}
+
+			var out interface{}
+			out = mv.Call(args)[0].Interface()
+			if m.Name == "Int" || m.Name == "Intn" {
+				out = int64(out.(int))
+			}
+			if m.Name == "Read" {
+				out = args[0].Interface().([]byte)
+			}
+			if *printgolden {
+				var val string
+				big := int64(1 << 60)
+				if int64(int(big)) != big && (m.Name == "Int" || m.Name == "Intn") {
+					// 32-bit machine cannot print 64-bit results
+					val = "truncated"
+				} else if reflect.TypeOf(out).Kind() == reflect.Slice {
+					val = fmt.Sprintf("%#v", out)
+				} else {
+					val = fmt.Sprintf("%T(%v)", out, out)
+				}
+				fmt.Printf("\t%s, // %s(%s)\n", val, m.Name, argstr)
+			} else {
+				want := regressGolden[p]
+				if m.Name == "Int" {
+					want = int64(int(uint(want.(int64)) << 1 >> 1))
+				}
+				if !reflect.DeepEqual(out, want) {
+					t.Errorf("r.%s(%s) = %v, want %v", m.Name, argstr, out, want)
+				}
+			}
+			p++
+		}
+	}
+	if *printgolden {
+		fmt.Printf("}\n")
+	}
+}
+
+var regressGolden = []interface{}{
+	float64(0.6279600685109523),   // ExpFloat64()
+	float64(0.16198826513357806),  // ExpFloat64()
+	float64(0.007880404652650552), // ExpFloat64()
+	float64(0.41649788761745654),  // ExpFloat64()
+	float64(1.6958707787276301),   // ExpFloat64()
+	float64(2.7227327706138036),   // ExpFloat64()
+	float64(2.4235600263079657),   // ExpFloat64()
+	float64(1.277967771105418),    // ExpFloat64()
+	float64(0.7111660437031769),   // ExpFloat64()
+	float64(0.23090401427981888),  // ExpFloat64()
+	float64(1.4746763588379928),   // ExpFloat64()
+	float64(1.4868726779832278),   // ExpFloat64()
+	float64(0.1686257242078103),   // ExpFloat64()
+	float64(0.2732721816228957),   // ExpFloat64()
+	float64(0.4644536065869748),   // ExpFloat64()
+	float64(0.01319850986379164),  // ExpFloat64()
+	float64(0.7184492551742854),   // ExpFloat64()
+	float64(0.1913536422195827),   // ExpFloat64()
+	float64(0.16034475958495667),  // ExpFloat64()
+	float64(0.40599859014785644),  // ExpFloat64()
+
+	float32(0.7979972),   // Float32()
+	float32(0.7725961),   // Float32()
+	float32(0.21894403),  // Float32()
+	float32(0.96194494),  // Float32()
+	float32(0.2915732),   // Float32()
+	float32(0.59569645),  // Float32()
+	float32(0.99596655),  // Float32()
+	float32(0.4979039),   // Float32()
+	float32(0.98148686),  // Float32()
+	float32(0.01380035),  // Float32()
+	float32(0.086487144), // Float32()
+	float32(0.6114401),   // Float32()
+	float32(0.71081316),  // Float32()
+	float32(0.6342346),   // Float32()
+	float32(0.008082573), // Float32()
+	float32(0.33020085),  // Float32()
+	float32(0.032625034), // Float32()
+	float32(0.9278005),   // Float32()
+	float32(0.34497985),  // Float32()
+	float32(0.66506875),  // Float32()
+
+	float64(0.797997151016231),    // Float64()
+	float64(0.7725961454373316),   // Float64()
+	float64(0.21894402538580782),  // Float64()
+	float64(0.9619449481780457),   // Float64()
+	float64(0.2915731877602916),   // Float64()
+	float64(0.5956964580775652),   // Float64()
+	float64(0.9959665347028619),   // Float64()
+	float64(0.49790390966591147),  // Float64()
+	float64(0.9814868602566785),   // Float64()
+	float64(0.013800350332924483), // Float64()
+	float64(0.08648714463652596),  // Float64()
+	float64(0.6114401479210267),   // Float64()
+	float64(0.7108131531183706),   // Float64()
+	float64(0.6342346133706837),   // Float64()
+	float64(0.008082572853887138), // Float64()
+	float64(0.3302008651926287),   // Float64()
+	float64(0.03262503454637655),  // Float64()
+	float64(0.9278004634858956),   // Float64()
+	float64(0.3449798628384906),   // Float64()
+	float64(0.665068719316529),    // Float64()
+
+	int64(5474557666971700975), // Int()
+	int64(5591422465364813936), // Int()
+	int64(74029666500212977),   // Int()
+	int64(8088122161323000979), // Int()
+	int64(7298457654139700474), // Int()
+	int64(1590632625527662686), // Int()
+	int64(9052198920789078554), // Int()
+	int64(7381380909356947872), // Int()
+	int64(1738222704626512495), // Int()
+	int64(3278744831230954970), // Int()
+	int64(7062423222661652521), // Int()
+	int64(6715870808026712034), // Int()
+	int64(528819992478005418),  // Int()
+	int64(2284534088986354339), // Int()
+	int64(945828723091990082),  // Int()
+	int64(3813019469742317492), // Int()
+	int64(1369388146907482806), // Int()
+	int64(7367238674766648970), // Int()
+	int64(8217673022687244206), // Int()
+	int64(3185531743396549562), // Int()
+
+	int32(1711064216), // Int31()
+	int32(650927245),  // Int31()
+	int32(8618187),    // Int31()
+	int32(941581344),  // Int31()
+	int32(1923394120), // Int31()
+	int32(1258915833), // Int31()
+	int32(1053814650), // Int31()
+	int32(859305834),  // Int31()
+	int32(1276097579), // Int31()
+	int32(1455437958), // Int31()
+	int32(1895916096), // Int31()
+	int32(781830261),  // Int31()
+	int32(61562749),   // Int31()
+	int32(265954771),  // Int31()
+	int32(1183850779), // Int31()
+	int32(443893888),  // Int31()
+	int32(1233159585), // Int31()
+	int32(857659461),  // Int31()
+	int32(956663049),  // Int31()
+	int32(370844703),  // Int31()
+
+	int32(0),          // Int31n(1)
+	int32(6),          // Int31n(10)
+	int32(17),         // Int31n(32)
+	int32(1000595),    // Int31n(1048576)
+	int32(424333),     // Int31n(1048577)
+	int32(382438494),  // Int31n(1000000000)
+	int32(902738458),  // Int31n(1073741824)
+	int32(1204933878), // Int31n(2147483646)
+	int32(1376191263), // Int31n(2147483647)
+	int32(0),          // Int31n(1)
+	int32(9),          // Int31n(10)
+	int32(2),          // Int31n(32)
+	int32(440490),     // Int31n(1048576)
+	int32(176312),     // Int31n(1048577)
+	int32(946765890),  // Int31n(1000000000)
+	int32(665034676),  // Int31n(1073741824)
+	int32(1947285452), // Int31n(2147483646)
+	int32(1702344608), // Int31n(2147483647)
+	int32(0),          // Int31n(1)
+	int32(2),          // Int31n(10)
+
+	int64(5474557666971700975), // Int63()
+	int64(5591422465364813936), // Int63()
+	int64(74029666500212977),   // Int63()
+	int64(8088122161323000979), // Int63()
+	int64(7298457654139700474), // Int63()
+	int64(1590632625527662686), // Int63()
+	int64(9052198920789078554), // Int63()
+	int64(7381380909356947872), // Int63()
+	int64(1738222704626512495), // Int63()
+	int64(3278744831230954970), // Int63()
+	int64(7062423222661652521), // Int63()
+	int64(6715870808026712034), // Int63()
+	int64(528819992478005418),  // Int63()
+	int64(2284534088986354339), // Int63()
+	int64(945828723091990082),  // Int63()
+	int64(3813019469742317492), // Int63()
+	int64(1369388146907482806), // Int63()
+	int64(7367238674766648970), // Int63()
+	int64(8217673022687244206), // Int63()
+	int64(3185531743396549562), // Int63()
+
+	int64(0),                   // Int63n(1)
+	int64(6),                   // Int63n(10)
+	int64(17),                  // Int63n(32)
+	int64(1000595),             // Int63n(1048576)
+	int64(424333),              // Int63n(1048577)
+	int64(382438494),           // Int63n(1000000000)
+	int64(902738458),           // Int63n(1073741824)
+	int64(1204933878),          // Int63n(2147483646)
+	int64(1376191263),          // Int63n(2147483647)
+	int64(502116868085730778),  // Int63n(1000000000000000000)
+	int64(144894195020570665),  // Int63n(1152921504606846976)
+	int64(6715870808026712034), // Int63n(9223372036854775806)
+	int64(528819992478005418),  // Int63n(9223372036854775807)
+	int64(0),                   // Int63n(1)
+	int64(0),                   // Int63n(10)
+	int64(20),                  // Int63n(32)
+	int64(854710),              // Int63n(1048576)
+	int64(649893),              // Int63n(1048577)
+	int64(687244206),           // Int63n(1000000000)
+	int64(836883386),           // Int63n(1073741824)
+
+	int64(0),                   // Intn(1)
+	int64(6),                   // Intn(10)
+	int64(17),                  // Intn(32)
+	int64(1000595),             // Intn(1048576)
+	int64(424333),              // Intn(1048577)
+	int64(382438494),           // Intn(1000000000)
+	int64(902738458),           // Intn(1073741824)
+	int64(1204933878),          // Intn(2147483646)
+	int64(1376191263),          // Intn(2147483647)
+	int64(502116868085730778),  // Intn(1000000000000000000)
+	int64(144894195020570665),  // Intn(1152921504606846976)
+	int64(6715870808026712034), // Intn(9223372036854775806)
+	int64(528819992478005418),  // Intn(9223372036854775807)
+	int64(0),                   // Intn(1)
+	int64(0),                   // Intn(10)
+	int64(20),                  // Intn(32)
+	int64(854710),              // Intn(1048576)
+	int64(649893),              // Intn(1048577)
+	int64(687244206),           // Intn(1000000000)
+	int64(836883386),           // Intn(1073741824)
+
+	float64(-0.5410658516792047),  // NormFloat64()
+	float64(0.615296849055287),    // NormFloat64()
+	float64(0.007477442280032887), // NormFloat64()
+	float64(1.3443892057169684),   // NormFloat64()
+	float64(-0.17508902754863512), // NormFloat64()
+	float64(-2.03494397556937),    // NormFloat64()
+	float64(2.5213558871972306),   // NormFloat64()
+	float64(1.4572921639613627),   // NormFloat64()
+	float64(-1.5164961164210644),  // NormFloat64()
+	float64(-0.4861150771891445),  // NormFloat64()
+	float64(-0.8699409548614199),  // NormFloat64()
+	float64(1.6271559815452794),   // NormFloat64()
+	float64(0.1659465769926195),   // NormFloat64()
+	float64(0.2921716191987018),   // NormFloat64()
+	float64(-1.2550269636927838),  // NormFloat64()
+	float64(0.11257973349467548),  // NormFloat64()
+	float64(0.5437525915836436),   // NormFloat64()
+	float64(0.781754430770282),    // NormFloat64()
+	float64(0.5201256313962235),   // NormFloat64()
+	float64(1.3826174159276245),   // NormFloat64()
+
+	[]int{},                                                     // Perm(0)
+	[]int{0},                                                    // Perm(1)
+	[]int{0, 2, 3, 1, 4},                                        // Perm(5)
+	[]int{5, 6, 3, 7, 4, 2, 0, 1},                               // Perm(8)
+	[]int{8, 4, 5, 2, 7, 3, 0, 6, 1},                            // Perm(9)
+	[]int{6, 1, 5, 3, 2, 9, 7, 0, 8, 4},                         // Perm(10)
+	[]int{12, 5, 1, 9, 15, 7, 13, 6, 10, 11, 8, 0, 4, 2, 14, 3}, // Perm(16)
+	[]int{},                                                     // Perm(0)
+	[]int{0},                                                    // Perm(1)
+	[]int{0, 2, 3, 4, 1},                                        // Perm(5)
+	[]int{3, 2, 7, 4, 0, 6, 5, 1},                               // Perm(8)
+	[]int{0, 6, 2, 1, 3, 7, 5, 8, 4},                            // Perm(9)
+	[]int{2, 5, 6, 4, 7, 3, 0, 8, 1, 9},                         // Perm(10)
+	[]int{3, 6, 5, 4, 9, 15, 13, 7, 1, 11, 10, 8, 12, 0, 2, 14}, // Perm(16)
+	[]int{},                             // Perm(0)
+	[]int{0},                            // Perm(1)
+	[]int{2, 4, 3, 1, 0},                // Perm(5)
+	[]int{1, 6, 7, 5, 4, 3, 2, 0},       // Perm(8)
+	[]int{7, 6, 8, 2, 0, 1, 3, 4, 5},    // Perm(9)
+	[]int{2, 9, 7, 1, 5, 4, 0, 6, 8, 3}, // Perm(10)
+
+	[]byte{0xef}, // Read([0])
+	[]byte{0x4e, 0x3d, 0x52, 0x31, 0x89, 0xf9, 0xcb},                   // Read([0 0 0 0 0 0 0])
+	[]byte{0x70, 0x68, 0x35, 0x8d, 0x1b, 0xb9, 0x98, 0x4d},             // Read([0 0 0 0 0 0 0 0])
+	[]byte{0xf1, 0xf8, 0x95, 0xe6, 0x96, 0x1, 0x7, 0x1, 0x93},          // Read([0 0 0 0 0 0 0 0 0])
+	[]byte{0x44, 0x9f, 0xc5, 0x40, 0xc8, 0x3e, 0x70, 0xfa, 0x44, 0x3a}, // Read([0 0 0 0 0 0 0 0 0 0])
+	[]byte{0x4b}, // Read([0])
+	[]byte{0x91, 0x54, 0x49, 0xe5, 0x5e, 0x28, 0xb9},                   // Read([0 0 0 0 0 0 0])
+	[]byte{0x4, 0xf2, 0xf, 0x13, 0x96, 0x1a, 0xb2, 0xce},               // Read([0 0 0 0 0 0 0 0])
+	[]byte{0x35, 0xf5, 0xde, 0x9f, 0x7d, 0xa0, 0x19, 0x12, 0x2e},       // Read([0 0 0 0 0 0 0 0 0])
+	[]byte{0xd4, 0xee, 0x6f, 0x66, 0x6f, 0x32, 0xc8, 0x21, 0x57, 0x68}, // Read([0 0 0 0 0 0 0 0 0 0])
+	[]byte{0x1f},                                                      // Read([0])
+	[]byte{0x98, 0xda, 0x4d, 0xab, 0x6e, 0xd, 0x71},                   // Read([0 0 0 0 0 0 0])
+	[]byte{0x80, 0xad, 0x29, 0xa0, 0x37, 0xb0, 0x80, 0xc4},            // Read([0 0 0 0 0 0 0 0])
+	[]byte{0x2, 0xe2, 0xe2, 0x7, 0xd9, 0xed, 0xea, 0x90, 0x33},        // Read([0 0 0 0 0 0 0 0 0])
+	[]byte{0x5d, 0xaa, 0xb8, 0xc6, 0x39, 0xfb, 0xbe, 0x56, 0x7, 0xa3}, // Read([0 0 0 0 0 0 0 0 0 0])
+	[]byte{0x62}, // Read([0])
+	[]byte{0x4d, 0x63, 0xa6, 0x4b, 0xb4, 0x1f, 0x42},                // Read([0 0 0 0 0 0 0])
+	[]byte{0x66, 0x42, 0x62, 0x36, 0x42, 0x20, 0x8d, 0xb4},          // Read([0 0 0 0 0 0 0 0])
+	[]byte{0x9f, 0xa3, 0x67, 0x1, 0x91, 0xea, 0x34, 0xb6, 0xa},      // Read([0 0 0 0 0 0 0 0 0])
+	[]byte{0xd, 0xa8, 0x43, 0xb, 0x1, 0x93, 0x8a, 0x56, 0xfc, 0x98}, // Read([0 0 0 0 0 0 0 0 0 0])
+
+	uint32(3422128433), // Uint32()
+	uint32(1301854491), // Uint32()
+	uint32(17236374),   // Uint32()
+	uint32(1883162688), // Uint32()
+	uint32(3846788241), // Uint32()
+	uint32(2517831666), // Uint32()
+	uint32(2107629301), // Uint32()
+	uint32(1718611668), // Uint32()
+	uint32(2552195159), // Uint32()
+	uint32(2910875917), // Uint32()
+	uint32(3791832192), // Uint32()
+	uint32(1563660522), // Uint32()
+	uint32(123125499),  // Uint32()
+	uint32(531909542),  // Uint32()
+	uint32(2367701558), // Uint32()
+	uint32(887787777),  // Uint32()
+	uint32(2466319171), // Uint32()
+	uint32(1715318922), // Uint32()
+	uint32(1913326099), // Uint32()
+	uint32(741689406),  // Uint32()
+
+	uint64(14697929703826476783), // Uint64()
+	uint64(5591422465364813936),  // Uint64()
+	uint64(74029666500212977),    // Uint64()
+	uint64(8088122161323000979),  // Uint64()
+	uint64(16521829690994476282), // Uint64()
+	uint64(10814004662382438494), // Uint64()
+	uint64(9052198920789078554),  // Uint64()
+	uint64(7381380909356947872),  // Uint64()
+	uint64(10961594741481288303), // Uint64()
+	uint64(12502116868085730778), // Uint64()
+	uint64(16285795259516428329), // Uint64()
+	uint64(6715870808026712034),  // Uint64()
+	uint64(528819992478005418),   // Uint64()
+	uint64(2284534088986354339),  // Uint64()
+	uint64(10169200759946765890), // Uint64()
+	uint64(3813019469742317492),  // Uint64()
+	uint64(10592760183762258614), // Uint64()
+	uint64(7367238674766648970),  // Uint64()
+	uint64(8217673022687244206),  // Uint64()
+	uint64(3185531743396549562),  // Uint64()
+
+	uint64(0),                   // Uint64n(1)
+	uint64(6),                   // Uint64n(10)
+	uint64(17),                  // Uint64n(32)
+	uint64(1000595),             // Uint64n(1048576)
+	uint64(424333),              // Uint64n(1048577)
+	uint64(382438494),           // Uint64n(1000000000)
+	uint64(902738458),           // Uint64n(1073741824)
+	uint64(1204933878),          // Uint64n(2147483646)
+	uint64(1376191263),          // Uint64n(2147483647)
+	uint64(502116868085730778),  // Uint64n(1000000000000000000)
+	uint64(144894195020570665),  // Uint64n(1152921504606846976)
+	uint64(6715870808026712034), // Uint64n(18446744073709551614)
+	uint64(528819992478005418),  // Uint64n(18446744073709551615)
+	uint64(0),                   // Uint64n(1)
+	uint64(0),                   // Uint64n(10)
+	uint64(20),                  // Uint64n(32)
+	uint64(854710),              // Uint64n(1048576)
+	uint64(649893),              // Uint64n(1048577)
+	uint64(687244206),           // Uint64n(1000000000)
+	uint64(836883386),           // Uint64n(1073741824)
+}
diff --git a/rand/rng.go b/rand/rng.go
new file mode 100644
index 0000000..d442985
--- /dev/null
+++ b/rand/rng.go
@@ -0,0 +1,91 @@
+// 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/bits"
+
+// PCGSource is an implementation of a 64-bit permuted congruential
+// generator as defined in
+//
+// 	PCG: A Family of Simple Fast Space-Efficient Statistically Good
+// 	Algorithms for Random Number Generation
+// 	Melissa E. O’Neill, Harvey Mudd College
+// 	http://www.pcg-random.org/pdf/toms-oneill-pcg-family-v1.02.pdf
+//
+// The generator here is the congruential generator PCG XSL RR 128/64 (LCG)
+// as found in the software available at http://www.pcg-random.org/.
+// It has period 2^128 with 128 bits of state, producing 64-bit values.
+// Is state is represented by two uint64 words.
+
+type PCGSource struct {
+	low  uint64
+	high uint64
+}
+
+const (
+	maxUint32 = (1 << 32) - 1
+
+	multiplier = 47026247687942121848144207491837523525
+
+	increment = 117397592171526113268558934119004209487
+	incHigh   = increment >> 64
+	incLow    = increment & maxUint64
+
+	// TODO: Use these?
+	initializer = 245720598905631564143578724636268694099
+	initHigh    = initializer >> 64
+	initLow     = initializer & maxUint64
+)
+
+// Seed uses the provided seed value to initialize the generator to a deterministic state.
+func (pcg *PCGSource) Seed(seed uint64) {
+	pcg.low = seed
+	pcg.high = seed // TODO: What is right?
+}
+
+// Uint64 returns a pseudo-random 64-bit unsigned integer as a uint64.
+func (pcg *PCGSource) Uint64() uint64 {
+	pcg.multiply()
+	pcg.add()
+	// XOR high and low 64 bits together and rotate right by high 6 bits of state.
+	return bits.RotateLeft64(pcg.high^pcg.low, -int(pcg.high>>58))
+}
+
+func (pcg *PCGSource) add() {
+	old := pcg.low
+	pcg.low += incLow
+	if pcg.low < old {
+		// Carry occurred.
+		pcg.high++
+	}
+	pcg.high += incHigh
+}
+
+func (pcg *PCGSource) multiply() {
+	// Break each lower word into two separate 32-bit 'digits' each stored
+	// in a 64-bit word with 32 high zero bits.  This allows the overflow
+	// into the high word to be computed.
+	s0 := (pcg.low >> 00) & maxUint32
+	s1 := (pcg.low >> 32) & maxUint32
+
+	const (
+		m0    = (multiplier >> 00) & maxUint32
+		m1    = (multiplier >> 32) & maxUint32
+		mLow  = multiplier & (1<<64 - 1)
+		mHigh = multiplier >> 64 & (1<<64 - 1)
+	)
+
+	high := pcg.low*mHigh + pcg.high*mLow
+	s0m0 := s0 * m0
+	s0m1 := s0 * m1
+	s1m0 := s1 * m0
+	s1m1 := s1 * m1
+	high += (s0m1 >> 32) + (s1m0 >> 32)
+	carry := (s0m1 & maxUint32) + (s1m0 & maxUint32) + s0m0>>32
+	high += (carry >> 32)
+
+	pcg.low *= mLow
+	pcg.high = high + s1m1
+}
diff --git a/rand/zipf.go b/rand/zipf.go
new file mode 100644
index 0000000..f04c814
--- /dev/null
+++ b/rand/zipf.go
@@ -0,0 +1,77 @@
+// 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.
+
+// W.Hormann, G.Derflinger:
+// "Rejection-Inversion to Generate Variates
+// from Monotone Discrete Distributions"
+// http://eeyore.wu-wien.ac.at/papers/96-04-04.wh-der.ps.gz
+
+package rand
+
+import "math"
+
+// A Zipf generates Zipf distributed variates.
+type Zipf struct {
+	r            *Rand
+	imax         float64
+	v            float64
+	q            float64
+	s            float64
+	oneminusQ    float64
+	oneminusQinv float64
+	hxm          float64
+	hx0minusHxm  float64
+}
+
+func (z *Zipf) h(x float64) float64 {
+	return math.Exp(z.oneminusQ*math.Log(z.v+x)) * z.oneminusQinv
+}
+
+func (z *Zipf) hinv(x float64) float64 {
+	return math.Exp(z.oneminusQinv*math.Log(z.oneminusQ*x)) - z.v
+}
+
+// NewZipf returns a Zipf variate generator.
+// The generator generates values k ∈ [0, imax]
+// such that P(k) is proportional to (v + k) ** (-s).
+// Requirements: s > 1 and v >= 1.
+func NewZipf(r *Rand, s float64, v float64, imax uint64) *Zipf {
+	z := new(Zipf)
+	if s <= 1.0 || v < 1 {
+		return nil
+	}
+	z.r = r
+	z.imax = float64(imax)
+	z.v = v
+	z.q = s
+	z.oneminusQ = 1.0 - z.q
+	z.oneminusQinv = 1.0 / z.oneminusQ
+	z.hxm = z.h(z.imax + 0.5)
+	z.hx0minusHxm = z.h(0.5) - math.Exp(math.Log(z.v)*(-z.q)) - z.hxm
+	z.s = 1 - z.hinv(z.h(1.5)-math.Exp(-z.q*math.Log(z.v+1.0)))
+	return z
+}
+
+// Uint64 returns a value drawn from the Zipf distribution described
+// by the Zipf object.
+func (z *Zipf) Uint64() uint64 {
+	if z == nil {
+		panic("rand: nil Zipf")
+	}
+	k := 0.0
+
+	for {
+		r := z.r.Float64() // r on [0,1]
+		ur := z.hxm + r*z.hx0minusHxm
+		x := z.hinv(ur)
+		k = math.Floor(x + 0.5)
+		if k-x <= z.s {
+			break
+		}
+		if ur >= z.h(k+0.5)-math.Exp(-math.Log(k+z.v)*z.q) {
+			break
+		}
+	}
+	return uint64(k)
+}