gtsocial-umbx

rnd.go (9303B)

---

 // 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
 }