exp-spec.txt

/*
 *      $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);
}