blob: c4e7b7d93e9ab20012aa08106636d6f449e7f166 [file] [log] [blame]
 // 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 := float64(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) }