summaryrefslogtreecommitdiffstats
path: root/vendor/modernc.org/mathutil/rnd.go
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/modernc.org/mathutil/rnd.go')
-rw-r--r--vendor/modernc.org/mathutil/rnd.go383
1 files changed, 383 insertions, 0 deletions
diff --git a/vendor/modernc.org/mathutil/rnd.go b/vendor/modernc.org/mathutil/rnd.go
new file mode 100644
index 00000000..598bad4f
--- /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
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-29 11:08:54 +0000