1// Copyright 2023 The Go Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package rand
6
7import (
8	"errors"
9	"internal/byteorder"
10	"math/bits"
11)
12
13// https://numpy.org/devdocs/reference/random/upgrading-pcg64.html
14// https://github.com/imneme/pcg-cpp/commit/871d0494ee9c9a7b7c43f753e3d8ca47c26f8005
15
16// A PCG is a PCG generator with 128 bits of internal state.
17// A zero PCG is equivalent to NewPCG(0, 0).
18type PCG struct {
19	hi uint64
20	lo uint64
21}
22
23// NewPCG returns a new PCG seeded with the given values.
24func NewPCG(seed1, seed2 uint64) *PCG {
25	return &PCG{seed1, seed2}
26}
27
28// Seed resets the PCG to behave the same way as NewPCG(seed1, seed2).
29func (p *PCG) Seed(seed1, seed2 uint64) {
30	p.hi = seed1
31	p.lo = seed2
32}
33
34// MarshalBinary implements the encoding.BinaryMarshaler interface.
35func (p *PCG) MarshalBinary() ([]byte, error) {
36	b := make([]byte, 20)
37	copy(b, "pcg:")
38	byteorder.BePutUint64(b[4:], p.hi)
39	byteorder.BePutUint64(b[4+8:], p.lo)
40	return b, nil
41}
42
43var errUnmarshalPCG = errors.New("invalid PCG encoding")
44
45// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface.
46func (p *PCG) UnmarshalBinary(data []byte) error {
47	if len(data) != 20 || string(data[:4]) != "pcg:" {
48		return errUnmarshalPCG
49	}
50	p.hi = byteorder.BeUint64(data[4:])
51	p.lo = byteorder.BeUint64(data[4+8:])
52	return nil
53}
54
55func (p *PCG) next() (hi, lo uint64) {
56	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L161
57	//
58	// Numpy's PCG multiplies by the 64-bit value cheapMul
59	// instead of the 128-bit value used here and in the official PCG code.
60	// This does not seem worthwhile, at least for Go: not having any high
61	// bits in the multiplier reduces the effect of low bits on the highest bits,
62	// and it only saves 1 multiply out of 3.
63	// (On 32-bit systems, it saves 1 out of 6, since Mul64 is doing 4.)
64	const (
65		mulHi = 2549297995355413924
66		mulLo = 4865540595714422341
67		incHi = 6364136223846793005
68		incLo = 1442695040888963407
69	)
70
71	// state = state * mul + inc
72	hi, lo = bits.Mul64(p.lo, mulLo)
73	hi += p.hi*mulLo + p.lo*mulHi
74	lo, c := bits.Add64(lo, incLo, 0)
75	hi, _ = bits.Add64(hi, incHi, c)
76	p.lo = lo
77	p.hi = hi
78	return hi, lo
79}
80
81// Uint64 return a uniformly-distributed random uint64 value.
82func (p *PCG) Uint64() uint64 {
83	hi, lo := p.next()
84
85	// XSL-RR would be
86	//	hi, lo := p.next()
87	//	return bits.RotateLeft64(lo^hi, -int(hi>>58))
88	// but Numpy uses DXSM and O'Neill suggests doing the same.
89	// See https://github.com/golang/go/issues/21835#issuecomment-739065688
90	// and following comments.
91
92	// DXSM "double xorshift multiply"
93	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L1015
94
95	// https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L176
96	const cheapMul = 0xda942042e4dd58b5
97	hi ^= hi >> 32
98	hi *= cheapMul
99	hi ^= hi >> 48
100	hi *= (lo | 1)
101	return hi
102}
103