pcg.go

     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(make([]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 || string(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
View as plain text