source: code/trunk/vendor/modernc.org/mathutil/rnd.go@ 822

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