1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// // 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 variate generator.
// // The generator generates values k ∈ [0, imax]
// // such that P(k) is proportional to (v + k) ** (-s).
// // Requirements: s > 1 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 distribution described
// // by the Zipf object.
// func (z *Zipf) Uint64() uint64 {
// if z == nil {
// panic("rand: nil Zipf")
// }
// 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)
// }