1 /*
2 * rand-gamma.c
3 *
4 * Custom variable implementation that returns gamma distributed random
5 * numbers.
6 *
7 */
8
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <math.h>
12
13 #include "cvar.h"
14 #include "cvar_trace.h"
15 #include "cvar_tokens.h"
16
17 #define VERSION "0.1.1 (alpha)"
18
19 #define USAGE_LEN 2048
20 static char usage[USAGE_LEN + 1];
21
22 #define MEAN "mean"
23 #define GAMMA "gamma"
24
25 #define MEAN_DEFAULT 4096.0
26 #define GAMMA_DEFAULT 1.5
27
28 typedef struct handle {
29 double mean;
30 double scaledmean;
31 double gamma;
32 } handle_t;
33
34 /******* Part from the original FB distribution ****************/
35
36 /*
37 * This is valid for 0 < a <= 1
38 *
39 * From Knuth Volume 2, 3rd edition, pages 586 - 587.
40 */
41 static double
gamma_dist_knuth_algG(double a,double (* src)(unsigned short *),unsigned short * xi)42 gamma_dist_knuth_algG(double a, double (*src)(unsigned short *),
43 unsigned short *xi)
44 {
45 double p, U, V, X, q;
46
47 p = M_E/(a + M_E);
48 G2:
49 /* get a random number U */
50 U = (*src)(xi);
51
52 do {
53 /* get a random number V */
54 V = (*src)(xi);
55
56 } while (V == 0);
57
58 if (U < p) {
59 X = pow(V, 1/a);
60 /* q = e^(-X) */
61 q = exp(-X);
62 } else {
63 X = 1 - log(V);
64 q = pow(X, a-1);
65 }
66
67 /*
68 * X now has density g, and q = f(X)/cg(X)
69 */
70
71 /* get a random number U */
72 U = (*src)(xi);
73
74 if (U >= q)
75 goto G2;
76 return (X);
77 }
78
79 /*
80 * This is valid for a > 1
81 *
82 * From Knuth Volume 2, 3rd edition, page 134.
83 */
84 static double
gamma_dist_knuth_algA(double a,double (* src)(unsigned short *),unsigned short * xi)85 gamma_dist_knuth_algA(double a, double (*src)(unsigned short *),
86 unsigned short *xi)
87 {
88 double U, Y, X, V;
89
90 A1:
91 /* get a random number U */
92 U = (*src)(xi);
93
94 Y = tan(M_PI*U);
95 X = (sqrt((2*a) - 1) * Y) + a - 1;
96
97 if (X <= 0)
98 goto A1;
99
100 /* get a random number V */
101 V = (*src)(xi);
102
103 if (V > ((1 + (Y*Y)) * exp((a-1) * log(X/(a-1)) - sqrt(2*a -1) * Y)))
104 goto A1;
105
106 return (X);
107 }
108
109 /*
110 * fetch a uniformly distributed random number using the drand48 generator
111 */
112 /* ARGSUSED */
113 static double
default_src(unsigned short * xi)114 default_src(unsigned short *xi)
115 {
116 return (drand48());
117 }
118
119 /*
120 * Sample the gamma distributed random variable with gamma 'a' and
121 * result mulitplier 'b', which is usually mean/gamma. Uses the default
122 * drand48 random number generator as input
123 */
124 static double
gamma_dist_knuth(double a,double b)125 gamma_dist_knuth(double a, double b)
126 {
127 if (a <= 1.0)
128 return (b * gamma_dist_knuth_algG(a, default_src, NULL));
129 else
130 return (b * gamma_dist_knuth_algA(a, default_src, NULL));
131 }
132
133 /*
134 * Sample the gamma distributed random variable with gamma 'a' and
135 * multiplier 'b', which is mean / gamma adjusted for the specified
136 * minimum value. The suppled random number source function is
137 * used to optain the uniformly distributed random numbers.
138 */
139 static double
gamma_dist_knuth_src(double a,double b,double (* src)(unsigned short *),unsigned short * xi)140 gamma_dist_knuth_src(double a, double b,
141 double (*src)(unsigned short *), unsigned short *xi)
142 {
143 if (a <= 1.0)
144 return (b * gamma_dist_knuth_algG(a, src, xi));
145 else
146 return (b * gamma_dist_knuth_algA(a, src, xi));
147 }
148 /*************************************************************/
149
cvar_alloc_handle(const char * cvar_parameters,void * (* cvar_malloc)(size_t size),void (* cvar_free)(void * cvar_ptr))150 void *cvar_alloc_handle(const char *cvar_parameters,
151 void *(*cvar_malloc)(size_t size), void (*cvar_free)(void *cvar_ptr))
152 {
153 cvar_token_t *list_head;;
154 cvar_token_t *t;
155 handle_t handle;
156 handle_t *state = NULL;
157 int ret = 0;
158
159 cvar_trace("entry");
160
161 /* Tokenize parameters supplied by filebench. */
162 list_head = NULL;
163 ret = tokenize(cvar_parameters, DEFAULT_PARAMETER_DELIMITER,
164 DEFAULT_KEY_VALUE_DELIMITER, &list_head);
165 if (ret)
166 goto out;
167
168 /* Get the values for mean and gamma. */
169 t = find_token(list_head, MEAN);
170 if (t && t->value) {
171 t->used = 1;
172 handle.mean = atof(t->value);
173 } else
174 handle.mean = MEAN_DEFAULT;
175
176 t = find_token(list_head, GAMMA);
177 if (t && t->value) {
178 t->used = 1;
179 handle.gamma = atof(t->value);
180 } else
181 handle.gamma= GAMMA_DEFAULT;
182
183 if (!handle.gamma) {
184 cvar_log_error("Invalid parameter values: mean = %lf "
185 "and gamma = %lf. gamma must be greater "
186 "than 0", handle.mean, handle.gamma);
187 goto out;
188 }
189
190 handle.scaledmean = handle.mean / handle.gamma;
191
192 cvar_trace("mean = %lf, gamma = %lf", handle.mean, handle.gamma);
193
194 t = unused_tokens(list_head);
195 if (t) {
196 cvar_log_error("Unsupported parameter %s", t->key);
197 goto out;
198 }
199
200 /* All set. Now allocate space for the handle in the shared segment and
201 * copy the state over. */
202 state = (handle_t *)cvar_malloc(sizeof(handle_t));
203 if (!state) {
204 cvar_log_error("Out of memory");
205 goto out;
206 }
207
208 *state = handle;
209
210 out:
211 free_tokens(list_head);
212
213 cvar_trace("exit");
214 return state;
215 }
216
cvar_revalidate_handle(void * cvar_handle)217 int cvar_revalidate_handle(void *cvar_handle)
218 {
219 handle_t *h = (handle_t *) cvar_handle;
220
221 return 0;
222 }
223
cvar_next_value(void * cvar_handle,double * value)224 int cvar_next_value(void *cvar_handle, double *value)
225 {
226 handle_t *h = (handle_t *) cvar_handle;
227
228 if (!h) {
229 cvar_log_error("NULL cvar_handle");
230 return -1;
231 }
232
233 if (!value) {
234 cvar_log_error("NULL value");
235 return -1;
236 }
237
238 *value = gamma_dist_knuth(h->gamma, h->scaledmean);
239
240 return 0;
241 }
242
cvar_free_handle(void * handle,void (* cvar_free)(void * ptr))243 void cvar_free_handle(void *handle, void (*cvar_free)(void *ptr))
244 {
245 cvar_free(handle);
246 }
247
cvar_usage()248 const char *cvar_usage()
249 {
250 int offset;
251
252 if (usage[0])
253 return usage;
254
255 offset = 0;
256
257 offset += snprintf(usage + offset, USAGE_LEN - offset,
258 "\tparameter\tdefault\n");
259 offset += snprintf(usage + offset, USAGE_LEN - offset,
260 "\t---------\t-------\n");
261
262 offset += snprintf(usage + offset, USAGE_LEN - offset,
263 "\t%s\t\t%.1f\n", MEAN, MEAN_DEFAULT);
264 offset += snprintf(usage + offset, USAGE_LEN - offset,
265 "\t%s\t\t%.1f\n", GAMMA, GAMMA_DEFAULT);
266
267 offset += snprintf(usage + offset, USAGE_LEN - offset,
268 "Use '%c' to delimit parameters and '%c' to delimit key-value "
269 "pairs.\n", DEFAULT_PARAMETER_DELIMITER,
270 DEFAULT_KEY_VALUE_DELIMITER);
271 offset += snprintf(usage + offset, USAGE_LEN - offset,
272 "Example: '%s%c%.1f%c%s%c%.1f'",
273 MEAN, DEFAULT_KEY_VALUE_DELIMITER, MEAN_DEFAULT,
274 DEFAULT_PARAMETER_DELIMITER,
275 GAMMA, DEFAULT_KEY_VALUE_DELIMITER, GAMMA_DEFAULT);
276
277 return usage;
278 }
279
cvar_version()280 const char *cvar_version()
281 {
282 return VERSION;
283 }
284