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 5 package rand
6 7 import (
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).
18 type PCG struct {
19 hi uint64
20 lo uint64
21 }
22 23 // NewPCG returns a new PCG seeded with the given values.
24 func 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).
29 func (p *PCG) Seed(seed1, seed2 uint64) {
30 p.hi = seed1
31 p.lo = seed2
32 }
33 34 // AppendBinary implements the [encoding.BinaryAppender] interface.
35 func (p *PCG) AppendBinary(b []byte) ([]byte, error) {
36 b = append(b, "pcg:"...)
37 b = byteorder.BEAppendUint64(b, p.hi)
38 b = byteorder.BEAppendUint64(b, p.lo)
39 return b, nil
40 }
41 42 // MarshalBinary implements the [encoding.BinaryMarshaler] interface.
43 func (p *PCG) MarshalBinary() ([]byte, error) {
44 return p.AppendBinary([]byte{:0:20})
45 }
46 47 var errUnmarshalPCG = errors.New("invalid PCG encoding")
48 49 // UnmarshalBinary implements the [encoding.BinaryUnmarshaler] interface.
50 func (p *PCG) UnmarshalBinary(data []byte) error {
51 if len(data) != 20 || data[:4] != "pcg:" {
52 return errUnmarshalPCG
53 }
54 p.hi = byteorder.BEUint64(data[4:])
55 p.lo = byteorder.BEUint64(data[4+8:])
56 return nil
57 }
58 59 func (p *PCG) next() (hi, lo uint64) {
60 // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L161
61 //
62 // Numpy's PCG multiplies by the 64-bit value cheapMul
63 // instead of the 128-bit value used here and in the official PCG code.
64 // This does not seem worthwhile, at least for Go: not having any high
65 // bits in the multiplier reduces the effect of low bits on the highest bits,
66 // and it only saves 1 multiply out of 3.
67 // (On 32-bit systems, it saves 1 out of 6, since Mul64 is doing 4.)
68 const (
69 mulHi = 2549297995355413924
70 mulLo = 4865540595714422341
71 incHi = 6364136223846793005
72 incLo = 1442695040888963407
73 )
74 75 // state = state * mul + inc
76 hi, lo = bits.Mul64(p.lo, mulLo)
77 hi += p.hi*mulLo + p.lo*mulHi
78 lo, c := bits.Add64(lo, incLo, 0)
79 hi, _ = bits.Add64(hi, incHi, c)
80 p.lo = lo
81 p.hi = hi
82 return hi, lo
83 }
84 85 // Uint64 return a uniformly-distributed random uint64 value.
86 func (p *PCG) Uint64() uint64 {
87 hi, lo := p.next()
88 89 // XSL-RR would be
90 // hi, lo := p.next()
91 // return bits.RotateLeft64(lo^hi, -int(hi>>58))
92 // but Numpy uses DXSM and O'Neill suggests doing the same.
93 // See https://github.com/golang/go/issues/21835#issuecomment-739065688
94 // and following comments.
95 96 // DXSM "double xorshift multiply"
97 // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L1015
98 99 // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L176
100 const cheapMul = 0xda942042e4dd58b5
101 hi ^= hi >> 32
102 hi *= cheapMul
103 hi ^= hi >> 48
104 hi *= (lo | 1)
105 return hi
106 }
107