From 1a57267a17c2fc17fb6e104846fabc3e363c326c Mon Sep 17 00:00:00 2001 From: Emile Date: Fri, 16 Aug 2024 19:50:26 +0200 Subject: initial commit --- vendor/modernc.org/mathutil/rnd.go | 383 +++++++++++++++++++++++++++++++++++++ 1 file changed, 383 insertions(+) create mode 100644 vendor/modernc.org/mathutil/rnd.go (limited to 'vendor/modernc.org/mathutil/rnd.go') diff --git a/vendor/modernc.org/mathutil/rnd.go b/vendor/modernc.org/mathutil/rnd.go new file mode 100644 index 0000000..598bad4 --- /dev/null +++ b/vendor/modernc.org/mathutil/rnd.go @@ -0,0 +1,383 @@ +// Copyright (c) 2014 The mathutil 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 mathutil // import "modernc.org/mathutil" + +import ( + "fmt" + "math" + "math/big" +) + +// FC32 is a full cycle PRNG covering the 32 bit signed integer range. +// In contrast to full cycle generators shown at e.g. http://en.wikipedia.org/wiki/Full_cycle, +// this code doesn't produce values at constant delta (mod cycle length). +// The 32 bit limit is per this implementation, the algorithm used has no intrinsic limit on the cycle size. +// Properties include: +// - Adjustable limits on creation (hi, lo). +// - Positionable/randomly accessible (Pos, Seek). +// - Repeatable (deterministic). +// - Can run forward or backward (Next, Prev). +// - For a billion numbers cycle the Next/Prev PRN can be produced in cca 100-150ns. +// That's like 5-10 times slower compared to PRNs generated using the (non FC) rand package. +type FC32 struct { + cycle int64 // On average: 3 * delta / 2, (HQ: 2 * delta) + delta int64 // hi - lo + factors [][]int64 // This trades some space for hopefully a bit of speed (multiple adding vs multiplying). + lo int + mods []int // pos % set + pos int64 // Within cycle. + primes []int64 // Ordered. ∏ primes == cycle. + set []int64 // Reordered primes (magnitude order bases) according to seed. +} + +// NewFC32 returns a newly created FC32 adjusted for the closed interval [lo, hi] or an Error if any. +// If hq == true then trade some generation time for improved (pseudo)randomness. +func NewFC32(lo, hi int, hq bool) (r *FC32, err error) { + if lo > hi { + return nil, fmt.Errorf("invalid range %d > %d", lo, hi) + } + + if uint64(hi)-uint64(lo) > math.MaxUint32 { + return nil, fmt.Errorf("range out of int32 limits %d, %d", lo, hi) + } + + delta := int64(hi) - int64(lo) + // Find the primorial covering whole delta + n, set, p := int64(1), []int64{}, uint32(2) + if hq { + p++ + } + for { + set = append(set, int64(p)) + n *= int64(p) + if n > delta { + break + } + p, _ = NextPrime(p) + } + + // Adjust the set so n ∊ [delta, 2 * delta] (HQ: [delta, 3 * delta]) + // while keeping the cardinality of the set (correlates with the statistic "randomness quality") + // at max, i.e. discard atmost one member. + i := -1 // no candidate prime + if n > 2*(delta+1) { + for j, p := range set { + q := n / p + if q < delta+1 { + break + } + + i = j // mark the highest candidate prime set index + } + } + if i >= 0 { // shrink the inner cycle + n = n / set[i] + set = delete(set, i) + } + r = &FC32{ + cycle: n, + delta: delta, + factors: make([][]int64, len(set)), + lo: lo, + mods: make([]int, len(set)), + primes: set, + } + r.Seed(1) // the default seed should be always non zero + return +} + +// Cycle reports the length of the inner FCPRNG cycle. +// Cycle is atmost the double (HQ: triple) of the generator period (hi - lo + 1). +func (r *FC32) Cycle() int64 { + return r.cycle +} + +// Next returns the first PRN after Pos. +func (r *FC32) Next() int { + return r.step(1) +} + +// Pos reports the current position within the inner cycle. +func (r *FC32) Pos() int64 { + return r.pos +} + +// Prev return the first PRN before Pos. +func (r *FC32) Prev() int { + return r.step(-1) +} + +// Seed uses the provided seed value to initialize the generator to a deterministic state. +// A zero seed produces a "canonical" generator with worse randomness than for most non zero seeds. +// Still, the FC property holds for any seed value. +func (r *FC32) Seed(seed int64) { + u := uint64(seed) + r.set = mix(r.primes, &u) + n := int64(1) + for i, p := range r.set { + k := make([]int64, p) + v := int64(0) + for j := range k { + k[j] = v + v += n + } + n *= p + r.factors[i] = mix(k, &u) + } +} + +// Seek sets Pos to |pos| % Cycle. +func (r *FC32) Seek(pos int64) { //vet:ignore + if pos < 0 { + pos = -pos + } + pos %= r.cycle + r.pos = pos + for i, p := range r.set { + r.mods[i] = int(pos % p) + } +} + +func (r *FC32) step(dir int) int { + for { // avg loops per step: 3/2 (HQ: 2) + y := int64(0) + pos := r.pos + pos += int64(dir) + switch { + case pos < 0: + pos = r.cycle - 1 + case pos >= r.cycle: + pos = 0 + } + r.pos = pos + for i, mod := range r.mods { + mod += dir + p := int(r.set[i]) + switch { + case mod < 0: + mod = p - 1 + case mod >= p: + mod = 0 + } + r.mods[i] = mod + y += r.factors[i][mod] + } + if y <= r.delta { + return int(y) + r.lo + } + } +} + +func delete(set []int64, i int) (y []int64) { + for j, v := range set { + if j != i { + y = append(y, v) + } + } + return +} + +func mix(set []int64, seed *uint64) (y []int64) { + for len(set) != 0 { + *seed = rol(*seed) + i := int(*seed % uint64(len(set))) + y = append(y, set[i]) + set = delete(set, i) + } + return +} + +func rol(u uint64) (y uint64) { + y = u << 1 + if int64(u) < 0 { + y |= 1 + } + return +} + +// FCBig is a full cycle PRNG covering ranges outside of the int32 limits. +// For more info see the FC32 docs. +// Next/Prev PRN on a 1e15 cycle can be produced in about 2 µsec. +type FCBig struct { + cycle *big.Int // On average: 3 * delta / 2, (HQ: 2 * delta) + delta *big.Int // hi - lo + factors [][]*big.Int // This trades some space for hopefully a bit of speed (multiple adding vs multiplying). + lo *big.Int + mods []int // pos % set + pos *big.Int // Within cycle. + primes []int64 // Ordered. ∏ primes == cycle. + set []int64 // Reordered primes (magnitude order bases) according to seed. +} + +// NewFCBig returns a newly created FCBig adjusted for the closed interval [lo, hi] or an Error if any. +// If hq == true then trade some generation time for improved (pseudo)randomness. +func NewFCBig(lo, hi *big.Int, hq bool) (r *FCBig, err error) { + if lo.Cmp(hi) > 0 { + return nil, fmt.Errorf("invalid range %d > %d", lo, hi) + } + + delta := big.NewInt(0) + delta.Add(delta, hi).Sub(delta, lo) + + // Find the primorial covering whole delta + n, set, pp, p := big.NewInt(1), []int64{}, big.NewInt(0), uint32(2) + if hq { + p++ + } + for { + set = append(set, int64(p)) + pp.SetInt64(int64(p)) + n.Mul(n, pp) + if n.Cmp(delta) > 0 { + break + } + p, _ = NextPrime(p) + } + + // Adjust the set so n ∊ [delta, 2 * delta] (HQ: [delta, 3 * delta]) + // while keeping the cardinality of the set (correlates with the statistic "randomness quality") + // at max, i.e. discard atmost one member. + dd1 := big.NewInt(1) + dd1.Add(dd1, delta) + dd2 := big.NewInt(0) + dd2.Lsh(dd1, 1) + i := -1 // no candidate prime + if n.Cmp(dd2) > 0 { + q := big.NewInt(0) + for j, p := range set { + pp.SetInt64(p) + q.Set(n) + q.Div(q, pp) + if q.Cmp(dd1) < 0 { + break + } + + i = j // mark the highest candidate prime set index + } + } + if i >= 0 { // shrink the inner cycle + pp.SetInt64(set[i]) + n.Div(n, pp) + set = delete(set, i) + } + r = &FCBig{ + cycle: n, + delta: delta, + factors: make([][]*big.Int, len(set)), + lo: lo, + mods: make([]int, len(set)), + pos: big.NewInt(0), + primes: set, + } + r.Seed(1) // the default seed should be always non zero + return +} + +// Cycle reports the length of the inner FCPRNG cycle. +// Cycle is atmost the double (HQ: triple) of the generator period (hi - lo + 1). +func (r *FCBig) Cycle() *big.Int { + return r.cycle +} + +// Next returns the first PRN after Pos. +func (r *FCBig) Next() *big.Int { + return r.step(1) +} + +// Pos reports the current position within the inner cycle. +func (r *FCBig) Pos() *big.Int { + return r.pos +} + +// Prev return the first PRN before Pos. +func (r *FCBig) Prev() *big.Int { + return r.step(-1) +} + +// Seed uses the provided seed value to initialize the generator to a deterministic state. +// A zero seed produces a "canonical" generator with worse randomness than for most non zero seeds. +// Still, the FC property holds for any seed value. +func (r *FCBig) Seed(seed int64) { + u := uint64(seed) + r.set = mix(r.primes, &u) + n := big.NewInt(1) + v := big.NewInt(0) + pp := big.NewInt(0) + for i, p := range r.set { + k := make([]*big.Int, p) + v.SetInt64(0) + for j := range k { + k[j] = big.NewInt(0) + k[j].Set(v) + v.Add(v, n) + } + pp.SetInt64(p) + n.Mul(n, pp) + r.factors[i] = mixBig(k, &u) + } +} + +// Seek sets Pos to |pos| % Cycle. +func (r *FCBig) Seek(pos *big.Int) { + r.pos.Set(pos) + r.pos.Abs(r.pos) + r.pos.Mod(r.pos, r.cycle) + mod := big.NewInt(0) + pp := big.NewInt(0) + for i, p := range r.set { + pp.SetInt64(p) + r.mods[i] = int(mod.Mod(r.pos, pp).Int64()) + } +} + +func (r *FCBig) step(dir int) (y *big.Int) { + y = big.NewInt(0) + d := big.NewInt(int64(dir)) + for { // avg loops per step: 3/2 (HQ: 2) + r.pos.Add(r.pos, d) + switch { + case r.pos.Sign() < 0: + r.pos.Add(r.pos, r.cycle) + case r.pos.Cmp(r.cycle) >= 0: + r.pos.SetInt64(0) + } + for i, mod := range r.mods { + mod += dir + p := int(r.set[i]) + switch { + case mod < 0: + mod = p - 1 + case mod >= p: + mod = 0 + } + r.mods[i] = mod + y.Add(y, r.factors[i][mod]) + } + if y.Cmp(r.delta) <= 0 { + y.Add(y, r.lo) + return + } + y.SetInt64(0) + } +} + +func deleteBig(set []*big.Int, i int) (y []*big.Int) { + for j, v := range set { + if j != i { + y = append(y, v) + } + } + return +} + +func mixBig(set []*big.Int, seed *uint64) (y []*big.Int) { + for len(set) != 0 { + *seed = rol(*seed) + i := int(*seed % uint64(len(set))) + y = append(y, set[i]) + set = deleteBig(set, i) + } + return +} -- cgit 1.4.1