Subversion Repositories openrisc

// Copyright 2009 The Go Authors. All rights reserved.

// Use of this source code is governed by a BSD-style

// license that can be found in the LICENSE file.

// W.Hormann, G.Derflinger:

// "Rejection-Inversion to Generate Variates

// from Monotone Discrete Distributions"

// http://eeyore.wu-wien.ac.at/papers/96-04-04.wh-der.ps.gz

package rand

import "math"

// A Zipf generates Zipf distributed variates.

type Zipf struct {

        r            *Rand

        imax         float64

        v            float64

        q            float64

        s            float64

        oneminusQ    float64

        oneminusQinv float64

        hxm          float64

        hx0minusHxm  float64

}

func (z *Zipf) h(x float64) float64 {

        return math.Exp(z.oneminusQ*math.Log(z.v+x)) * z.oneminusQinv

}

func (z *Zipf) hinv(x float64) float64 {

        return math.Exp(z.oneminusQinv*math.Log(z.oneminusQ*x)) - z.v

}

// NewZipf returns a Zipf generating variates p(k) on [0, imax]

// proportional to (v+k)**(-s) where s>1 and k>=0, and v>=1.

//

func NewZipf(r *Rand, s float64, v float64, imax uint64) *Zipf {

        z := new(Zipf)

        if s <= 1.0 || v < 1 {

                return nil

        }

        z.r = r

        z.imax = float64(imax)

        z.v = v

        z.q = s

        z.oneminusQ = 1.0 - z.q

        z.oneminusQinv = 1.0 / z.oneminusQ

        z.hxm = z.h(z.imax + 0.5)

        z.hx0minusHxm = z.h(0.5) - math.Exp(math.Log(z.v)*(-z.q)) - z.hxm

        z.s = 1 - z.hinv(z.h(1.5)-math.Exp(-z.q*math.Log(z.v+1.0)))

        return z

}

// Uint64 returns a value drawn from the Zipf distributed described

// by the Zipf object.

func (z *Zipf) Uint64() uint64 {

        k := 0.0

        for {

                r := z.r.Float64() // r on [0,1]

                ur := z.hxm + r*z.hx0minusHxm

                x := z.hinv(ur)

                k = math.Floor(x + 0.5)

                if k-x <= z.s {

                        break

                }

                if ur >= z.h(k+0.5)-math.Exp(-math.Log(k+z.v)*z.q) {

                        break

                }

        }

        return uint64(k)

}