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 // MarshalBinary implements the encoding.BinaryMarshaler interface. 35 func (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 43 var errUnmarshalPCG = errors.New("invalid PCG encoding") 44 45 // UnmarshalBinary implements the encoding.BinaryUnmarshaler interface. 46 func (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 55 func (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. 82 func (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