forked from perfsonar/owamp
-
Notifications
You must be signed in to change notification settings - Fork 0
/
exp-spec.txt
415 lines (324 loc) · 10.6 KB
/
exp-spec.txt
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
/*
* $Id$
*/
/************************************************************************
* *
* Copyright (C) 2002 *
* Internet2 *
* All Rights Reserved *
* *
************************************************************************/
/*
* File: exp-spec.txt
*
* Author: Anatoly Karp
* Internet2
*
* Date: Sun Jun 20 12:29:20 MDT 2002
*
* Description: Part of OWAMP specification describing
* generation of uniform(0,1) and
* exponential (mean 1) random variates.
*/
1. Introduction.
This part of the document describes in detail the way exponential
random quantities used in the protocol are generated. While there is
a fair number of algorithms for generating exponential random
variables, most of them rely on having logarithmic function as a
primitive, resulting in potentially different values, depending on the
particular implementation of the math library.
This RFC adapts algorithm 3.4.1.S in [1] which is free of the above
mentioned problem, and guarantees the same output on any
implementation. The algorithm belongs to the 'ziggurat' family
developed in the 1970s by G.Marsaglia, M.Sibuya and J.H.Ahrens. It
replaces the use of logarithmic function by clever bit manipulation,
still producing the exponential variates on output.
2. High-level description of the algorithm.
For ease of exposition, the algorithm is first described with all
arithmetic operations being interpreted in their natural sense. Later,
exact details on data types, arithmetic, and generation of the uniform
random variates used by the algorithm are given.
Algorithm S: Given a real positive number 'mu', produce an exponential
random variate with mean 'mu'.
First, the constants
Q[k] = (ln2)/(1!) + (ln2)^2/(2!) + ... + (ln2)^k/(k!), 1 <= k <= 11
are computed in advance.
S1. [Get U and shift.] Generate a 32-bit uniform random binary fraction
U = (.b0 b1 b2 ... b31) [note the decimal point]
Locate the first zero bit b_j, and shift off the leading (j+1) bits,
setting U <- (.b_{j+1} ... b31)
NOTE: in the rare case that the zero has not been found it is prescribed
that the algorithm return (mu*32*ln2).
S2. [Immediate acceptance?] If U < ln2, set X <- mu*(j*ln2 + U) and terminate
the algorithm. (Note that Q[1] = ln2.)
S3. [Minimize.] Find the least k >= 2 such that U < Q[k]. Generate
k new uniform random binary fractions U1,...,Uk and set
V <- min(U1,...,Uk).
S4. [Deliver the answer.] Set X <- mu*(j + V)*ln2.
3. Data types, representation and arithmetic.
The high-level algorithm operates on real numbers, that is objects typically
represented as floating point numbers. This specification prescribes
that unsigned 64-bit integers be used instead.
typedef u_int64_t OWPnum64;
OWPnum64 integers are interpreted as real numbers by placing the
decimal point after the first 32 bits. In other words, coneceptually
the interpretation is given by the map:
OWPnum64 u;
u ---> (double)u / (2**32) (1)
The algorithm produces a sequence of such OWPnum64 integers which is
guaranteed to be the same on any implementation. Any further interpretation
(such as given by (1)) is done by the application, and is not part of
this specification.
a) This RFC specifies that the OWPNum64 representations of the
first 11 values of the Q array in the high-level algorithm be as
follows:
#1 0xB17217F8,
#2 0xEEF193F7,
#3 0xFD271862,
#4 0xFF9D6DD0,
#5 0xFFF4CFD0,
#6 0xFFFEE819,
#7 0xFFFFE7FF,
#8 0xFFFFFE2B,
#9 0xFFFFFFE0,
#10 0xFFFFFFFE,
#11 0xFFFFFFFF
(For example, Q[1] = ln2 is indeed approximated by 0xB17217F8/(2**32) =
0.693147180601954);
b) small integer 'j' in the high-level algorithm is represented
as OWPnum64 value j * (2**32);
c) operation of addition is done as usual on OWPNum64 numbers;
however, the operation of multiplication in the high-level algorithm
should be replaced by
(u, v) ---> (u * v) >> 32
NOTE: correct implementation should compute (u * v) exactly. For
example, a fragment of unsigned 128-bit arithmetic can be
implemented for this purpose (see sample implementation below)
4. Uniform random quantities.
The procedure for obtaining a sequence of 32-bit random numbers (such
as 'U' in algorithm S) relies on using AES encryption in counter
mode. To describe the exact working of the algorithm we introduce two
primitives from rijndael. Their prototypes and specification are given
below, and they are provided by every rijndael implementation, such as
the reference implementation available at
http://www.esat.kuleuven.ac.be/~rijmen/rijndael/rijndaelref.zip:
1. this function initializes a rijndael key with bytes from 'seed'
KeyInit(BYTE seed[16]);
2. this function encrypts the 16-octet block 'inblock' with
the 'key' returning a 16-octet encrypted block. Here 'keyInstance'
is an opaque type used to represent rijndael keys.
BlockEncrypt(keyInstance key, BYTE inblock[16]);
Algorithm Unif: given a 16-octet quantity 'seed', produce a sequence
of unsigned 32-bit pseudo-random uniformly distributed integers. In OWAMP,
the shared secret from Control protocol plays the role of 'seed'.
U1. [Initialize rijndael key]
key <- KeyInit(seed)
[Initialize an unsigned 16-octet (network byte order) counter]
c <- 0
U2. [Need more random bytes?]
Set i <- c mod 4. If (i == 0) set s <- BlockEncrypt(key, c)
U3. [Increment the counter as unsigned 16-octet quantity]
c <- c + 1
U4. [Do output] Output the i_th quartet of octets from s starting
from high-order octets, converted to native byte order and
represented as OWPNum64 value (as in 3.b).
U5. [Loop] Go to step U2.
5. Appendix. Sample implementation.
/*
** Example usage: generate a stream of exponential (mean 1)
** random quantities (ignoring error checking during initialization).
** Assume that a 16-octet 'seed' has been initialized
** (as the shared secret in OWAMP, for example)
** unsigned char seed[16];
** OWPrand_context64 next;
** (initialize state)
** OWPrand_context64_init(&next, seed);
** (generate a sequence of exponential variates)
** while (1) {
** OWPnum64 num = OWPexp_rand64(&next);
<do something with num here>
...
** }
*/
#include <stdlib.h>
#include <stdio.h>
#include "rijndael-alg-fst.h"
typedef u_int64_t OWPnum64;
/* (K - 1) is the first k such that Q[k] > 1 - 1/(2^32). */
#define K 12
#define MASK32(x) ((x) & 0xFFFFFFFF)
typedef struct OWPrand_context64 {
unsigned char counter[16]; /* 16-octet counter (network byte order) */
keyInstance key; /* key used to encrypt the counter. */
BYTE out[16]; /* the encrypted block is kept there. */
} OWPrand_context64;
/*
** The array has been computed according to the formula:
**
** Q[k] = (ln2)/(1!) + (ln2)^2/(2!) + ... + (ln2)^k/(k!)
**
** as described in algorithm S. (The values below have been
** multiplied by 2^32 and rounded to the nearest integer.)
*/
static OWPnum64 Q[K] = {
0, /* Placeholder - so array indices start from 1. */
0xB17217F8,
0xEEF193F7,
0xFD271862,
0xFF9D6DD0,
0xFFF4CFD0,
0xFFFEE819,
0xFFFFE7FF,
0xFFFFFE2B,
0xFFFFFFE0,
0xFFFFFFFE,
0xFFFFFFFF
};
/* this element represents ln2 */
#define LN2 Q[1]
/*
** Convert an unsigned 32-bit integer into a OWPnum64 number..
*/
OWPnum64
OWPulong2num64(u_int32_t a)
{
return ((u_int64_t)1 << 32) * a;
}
/*
** Arithmetic functions on OWPnum64 numbers.
*/
/*
** Addition.
*/
OWPnum64
OWPnum64_add(OWPnum64 x, OWPnum64 y)
{
return x + y;
}
/*
** Multiplication. Allows overflow. Straightforward implementation
** of Knuth vol.2 Algorithm 4.3.1.M (p.268)
*/
OWPnum64
OWPnum64_mul(OWPnum64 x, OWPnum64 y)
{
unsigned long w[4];
u_int64_t xdec[2];
u_int64_t ydec[2];
int i, j;
u_int64_t k, t;
OWPnum64 ret;
xdec[0] = x & 0xFFFFFFFF;
xdec[1] = x >> 32;
ydec[0] = y & 0xFFFFFFFF;
ydec[1] = y >> 32;
for (j = 0; j < 4; j++)
w[j] = 0;
for (j = 0; j < 2; j++) {
k = 0;
for (i = 0; ; ) {
t = k + (xdec[i]*ydec[j]) + w[i + j];
w[i + j] = t%((u_int64_t)1 << 32);
k = t/((u_int64_t)1 << 32);
if (++i < 2)
continue;
else {
w[j + 2] = k;
break;
}
}
}
ret = w[2];
ret <<= 32;
return w[1] + ret;
}
/*
** Seed the random number generator using a 16-byte quantity 'seed'
** (== the shared secret in OWAMP). This function implements step U1
** of algorithm Unif.
*/
void
OWPrand_context64_init(OWPrand_context64 *next, BYTE *seed)
{
int i;
/* Initialize the key */
rijndaelKeyInit(next->key, seed);
/* Initialize the counter with zeros */
memset(next->out, 0, 16);
for (i = 0; i < 16; i++)
next->counter[i] = 0UL;
}
/*
** Random number generating functions.
*/
/*
** Generate and return a 32-bit uniform random string (saved in the lower
** half of the OWPnum64. This function implements steps U2-U4 of the
** alsgorithm Unif.
*/
OWPnum64
OWPunif_rand64(OWPrand_context64 *next)
{
int j;
/* step U2 */
u_int8_t i = next->counter[15] & (u_int8_t)3;
if (!i)
rijndaelEncrypt(next->key, next->counter, next->out);
/* Step U3. Increment next.counter as a 16-octet single quantity
in network byte order for AES counter mode. */
for (j = 15; j >= 0; j--)
if (++next->counter[j])
break;
/* Step U4. Do output. */
return OWPraw2num64((next->out) + 4*i);
}
/*
** This function converts a 32-bit binary string (network byte order)
** into a OWPnum64 number (32 least significant bits). It is used in
** step U4 of algorithm Unif.
*/
static OWPnum64
OWPraw2num64(const unsigned char *raw)
{
return ((u_int32_t)(raw[0]) << 24)
+ ((u_int32_t)(raw[1]) << 16)
+ ((u_int32_t)(raw[2]) << 8)
+ (u_int32_t)raw[3];
}
/*
** Generate a mean 1 exponential deviate.
*/
OWPnum64
OWPexp_rand64(OWPrand_context64 *next)
{
unsigned long i, k;
u_int32_t j = 0;
OWPnum64 U, V, J, tmp;
u_int32_t mask = 0x80000000; /* see if first bit in the lower
32 bits is zero */
/* Step S1. Get U and shift */
U = OWPunif_rand64(next);
while (U & mask && j < 32){ /* shift until find first '0' */
U <<= 1;
j++;
}
/* remove the '0' itself */
U <<= 1;
U &= 0xFFFFFFFF; /* Keep only the fractional part. */
J = OWPulong2num64(j);
/* Step S2. Immediate acceptance? */
if (U < LN2) /* return (j*ln2 + U) */
return OWPnum64_add(OWPnum64_mul(J, LN2), U);
/* Step S3. Minimize. */
for (k = 2; k < K; k++)
if (U < Q[k])
break;
V = OWPunif_rand64(next);
for (i = 2; i <= k; i++){
tmp = OWPunif_rand64(next);
if (tmp < V)
V = tmp;
}
/* Step S4. Return (j+V)*ln2 */
return OWPnum64_mul(OWPnum64_add(J, V), LN2);
}