initial commit

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