1#include <math.h>
2#include <string.h>
3#include <inttypes.h>
4#include <stdio.h>
5#include <unistd.h>
6#include <sys/types.h>
7#include <fcntl.h>
8#include "ieee754.h"
9#include "../log.h"
10#include "zipf.h"
11#include "../minmax.h"
12#include "../hash.h"
13
14#define ZIPF_MAX_GEN	10000000
15
16static void zipf_update(struct zipf_state *zs)
17{
18	unsigned long to_gen;
19	unsigned int i;
20
21	/*
22	 * It can become very costly to generate long sequences. Just cap it at
23	 * 10M max, that should be doable in 1-2s on even slow machines.
24	 * Precision will take a slight hit, but nothing major.
25	 */
26	to_gen = min(zs->nranges, ZIPF_MAX_GEN);
27
28	for (i = 0; i < to_gen; i++)
29		zs->zetan += pow(1.0 / (double) (i + 1), zs->theta);
30}
31
32static void shared_rand_init(struct zipf_state *zs, unsigned long nranges,
33			     unsigned int seed)
34{
35	memset(zs, 0, sizeof(*zs));
36	zs->nranges = nranges;
37
38	init_rand_seed(&zs->rand, seed);
39	zs->rand_off = __rand(&zs->rand);
40}
41
42void zipf_init(struct zipf_state *zs, unsigned long nranges, double theta,
43	       unsigned int seed)
44{
45	shared_rand_init(zs, nranges, seed);
46
47	zs->theta = theta;
48	zs->zeta2 = pow(1.0, zs->theta) + pow(0.5, zs->theta);
49
50	zipf_update(zs);
51}
52
53unsigned long long zipf_next(struct zipf_state *zs)
54{
55	double alpha, eta, rand_uni, rand_z;
56	unsigned long long n = zs->nranges;
57	unsigned long long val;
58
59	alpha = 1.0 / (1.0 - zs->theta);
60	eta = (1.0 - pow(2.0 / n, 1.0 - zs->theta)) / (1.0 - zs->zeta2 / zs->zetan);
61
62	rand_uni = (double) __rand(&zs->rand) / (double) FRAND_MAX;
63	rand_z = rand_uni * zs->zetan;
64
65	if (rand_z < 1.0)
66		val = 1;
67	else if (rand_z < (1.0 + pow(0.5, zs->theta)))
68		val = 2;
69	else
70		val = 1 + (unsigned long long)(n * pow(eta*rand_uni - eta + 1.0, alpha));
71
72	return (__hash_u64(val - 1) + zs->rand_off) % zs->nranges;
73}
74
75void pareto_init(struct zipf_state *zs, unsigned long nranges, double h,
76		 unsigned int seed)
77{
78	shared_rand_init(zs, nranges, seed);
79	zs->pareto_pow = log(h) / log(1.0 - h);
80}
81
82unsigned long long pareto_next(struct zipf_state *zs)
83{
84	double rand = (double) __rand(&zs->rand) / (double) FRAND_MAX;
85	unsigned long long n = zs->nranges - 1;
86
87	return (__hash_u64(n * pow(rand, zs->pareto_pow)) + zs->rand_off) % zs->nranges;
88}
89