3 files changed, 206 insertions, 0 deletions
diff --git a/api/next/61716.txt b/api/next/61716.txt
index b84e7e1147..ac974a6117 100644
--- a/api/next/61716.txt
+++ b/api/next/61716.txt
@@ -9,6 +9,7 @@ pkg math/rand/v2, func Int64N(int64) int64 #61716
 pkg math/rand/v2, func IntN(int) int #61716
 pkg math/rand/v2, func N[$0 intType]($0) $0 #61716
 pkg math/rand/v2, func New(Source) *Rand #61716
+pkg math/rand/v2, func NewPCG(uint64, uint64) *PCG #61716
 pkg math/rand/v2, func NewSource(int64) Source #61716
 pkg math/rand/v2, func NewZipf(*Rand, float64, float64, uint64) *Zipf #61716
 pkg math/rand/v2, func NormFloat64() float64 #61716
@@ -19,6 +20,10 @@ pkg math/rand/v2, func Uint32N(uint32) uint32 #61716
 pkg math/rand/v2, func Uint64() uint64 #61716
 pkg math/rand/v2, func Uint64N(uint64) uint64 #61716
 pkg math/rand/v2, func UintN(uint) uint #61716
+pkg math/rand/v2, method (*PCG) MarshalBinary() ([]uint8, error) #61716
+pkg math/rand/v2, method (*PCG) Seed(uint64, uint64) #61716
+pkg math/rand/v2, method (*PCG) Uint64() uint64 #61716
+pkg math/rand/v2, method (*PCG) UnmarshalBinary([]uint8) error #61716
 pkg math/rand/v2, method (*Rand) ExpFloat64() float64 #61716
 pkg math/rand/v2, method (*Rand) Float32() float32 #61716
 pkg math/rand/v2, method (*Rand) Float64() float64 #61716
@@ -37,6 +42,7 @@ pkg math/rand/v2, method (*Rand) Uint64() uint64 #61716
 pkg math/rand/v2, method (*Rand) Uint64N(uint64) uint64 #61716
 pkg math/rand/v2, method (*Rand) UintN(uint) uint #61716
 pkg math/rand/v2, method (*Zipf) Uint64() uint64 #61716
+pkg math/rand/v2, type PCG struct #61716
 pkg math/rand/v2, type Rand struct #61716
 pkg math/rand/v2, type Source interface { Uint64 } #61716
 pkg math/rand/v2, type Source interface, Uint64() uint64 #61716
diff --git a/src/math/rand/v2/pcg.go b/src/math/rand/v2/pcg.go
new file mode 100644
index 0000000000..77708d799e
--- /dev/null
+++ b/src/math/rand/v2/pcg.go
@@ -0,0 +1,121 @@
+// Copyright 2023 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 (
+	"errors"
+	"math/bits"
+)
+
+// https://numpy.org/devdocs/reference/random/upgrading-pcg64.html
+// https://github.com/imneme/pcg-cpp/commit/871d0494ee9c9a7b7c43f753e3d8ca47c26f8005
+
+// A PCG is a PCG generator with 128 bits of internal state.
+// A zero PCG is equivalent to NewPCG(0, 0).
+type PCG struct {
+	hi uint64
+	lo uint64
+}
+
+// NewPCG returns a new PCG seeded with the given values.
+func NewPCG(seed1, seed2 uint64) *PCG {
+	return &PCG{seed1, seed2}
+}
+
+// Seed resets the PCG to behave the same way as NewPCG(seed1, seed2).
+func (p *PCG) Seed(seed1, seed2 uint64) {
+	p.hi = seed1
+	p.lo = seed2
+}
+
+// binary.bigEndian.Uint64, copied to avoid dependency
+func beUint64(b []byte) uint64 {
+	_ = b[7] // bounds check hint to compiler; see golang.org/issue/14808
+	return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 |
+		uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56
+}
+
+// binary.bigEndian.PutUint64, copied to avoid dependency
+func bePutUint64(b []byte, v uint64) {
+	_ = b[7] // early bounds check to guarantee safety of writes below
+	b[0] = byte(v >> 56)
+	b[1] = byte(v >> 48)
+	b[2] = byte(v >> 40)
+	b[3] = byte(v >> 32)
+	b[4] = byte(v >> 24)
+	b[5] = byte(v >> 16)
+	b[6] = byte(v >> 8)
+	b[7] = byte(v)
+}
+
+// MarshalBinary implements the encoding.BinaryMarshaler interface.
+func (p *PCG) MarshalBinary() ([]byte, error) {
+	b := make([]byte, 20)
+	copy(b, "pcg:")
+	bePutUint64(b[4:], p.hi)
+	bePutUint64(b[4+8:], p.lo)
+	return b, nil
+}
+
+var errUnmarshalPCG = errors.New("invalid PCG encoding")
+
+// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface.
+func (p *PCG) UnmarshalBinary(data []byte) error {
+	if len(data) != 20 || string(data[:4]) != "pcg:" {
+		return errUnmarshalPCG
+	}
+	p.hi = beUint64(data[4:])
+	p.lo = beUint64(data[4+8:])
+	return nil
+}
+
+func (p *PCG) next() (hi, lo uint64) {
+	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L161
+	//
+	// Numpy's PCG multiplies by the 64-bit value cheapMul
+	// instead of the 128-bit value used here and in the official PCG code.
+	// This does not seem worthwhile, at least for Go: not having any high
+	// bits in the multiplier reduces the effect of low bits on the highest bits,
+	// and it only saves 1 multiply out of 3.
+	// (On 32-bit systems, it saves 1 out of 6, since Mul64 is doing 4.)
+	const (
+		mulHi = 2549297995355413924
+		mulLo = 4865540595714422341
+		incHi = 6364136223846793005
+		incLo = 1442695040888963407
+	)
+
+	// state = state * mul + inc
+	hi, lo = bits.Mul64(p.lo, mulLo)
+	hi += p.hi*mulLo + p.lo*mulHi
+	lo, c := bits.Add64(lo, incLo, 0)
+	hi, _ = bits.Add64(hi, incHi, c)
+	p.lo = lo
+	p.hi = hi
+	return hi, lo
+}
+
+// Uint64 return a uniformly-distributed random uint64 value.
+func (p *PCG) Uint64() uint64 {
+	hi, lo := p.next()
+
+	// XSL-RR would be
+	//	hi, lo := p.next()
+	//	return bits.RotateLeft64(lo^hi, -int(hi>>58))
+	// but Numpy uses DXSM and O'Neill suggests doing the same.
+	// See https://github.com/golang/go/issues/21835#issuecomment-739065688
+	// and following comments.
+
+	// DXSM "double xorshift multiply"
+	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L1015
+
+	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L176
+	const cheapMul = 0xda942042e4dd58b5
+	hi ^= hi >> 32
+	hi *= cheapMul
+	hi ^= hi >> 48
+	hi *= (lo | 1)
+	return hi
+}
diff --git a/src/math/rand/v2/pcg_test.go b/src/math/rand/v2/pcg_test.go
new file mode 100644
index 0000000000..db866c8c85
--- /dev/null
+++ b/src/math/rand/v2/pcg_test.go
@@ -0,0 +1,79 @@
+// Copyright 2023 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 (
+	. "math/rand/v2"
+	"testing"
+)
+
+func BenchmarkPCG_DXSM(b *testing.B) {
+	var p PCG
+	var t uint64
+	for n := b.N; n > 0; n-- {
+		t += p.Uint64()
+	}
+	Sink = t
+}
+
+func TestPCGMarshal(t *testing.T) {
+	var p PCG
+	const (
+		seed1 = 0x123456789abcdef0
+		seed2 = 0xfedcba9876543210
+		want  = "pcg:\x12\x34\x56\x78\x9a\xbc\xde\xf0\xfe\xdc\xba\x98\x76\x54\x32\x10"
+	)
+	p.Seed(seed1, seed2)
+	data, err := p.MarshalBinary()
+	if string(data) != want || err != nil {
+		t.Errorf("MarshalBinary() = %q, %v, want %q, nil", data, err, want)
+	}
+
+	q := PCG{}
+	if err := q.UnmarshalBinary([]byte(want)); err != nil {
+		t.Fatalf("UnmarshalBinary(): %v", err)
+	}
+	if q != p {
+		t.Fatalf("after round trip, q = %#x, but p = %#x", q, p)
+	}
+
+	qu := q.Uint64()
+	pu := p.Uint64()
+	if qu != pu {
+		t.Errorf("after round trip, q.Uint64() = %#x, but p.Uint64() = %#x", qu, pu)
+	}
+}
+
+func TestPCG(t *testing.T) {
+	p := NewPCG(1, 2)
+	want := []uint64{
+		0xc4f5a58656eef510,
+		0x9dcec3ad077dec6c,
+		0xc8d04605312f8088,
+		0xcbedc0dcb63ac19a,
+		0x3bf98798cae97950,
+		0xa8c6d7f8d485abc,
+		0x7ffa3780429cd279,
+		0x730ad2626b1c2f8e,
+		0x21ff2330f4a0ad99,
+		0x2f0901a1947094b0,
+		0xa9735a3cfbe36cef,
+		0x71ddb0a01a12c84a,
+		0xf0e53e77a78453bb,
+		0x1f173e9663be1e9d,
+		0x657651da3ac4115e,
+		0xc8987376b65a157b,
+		0xbb17008f5fca28e7,
+		0x8232bd645f29ed22,
+		0x12be8f07ad14c539,
+		0x54908a48e8e4736e,
+	}
+
+	for i, x := range want {
+		if u := p.Uint64(); u != x {
+			t.Errorf("PCG #%d = %#x, want %#x", i, u, x)
+		}
+	}
+}