## 2008-03-13

### Bitwise Dice

I'm rereading this paper by Knuth and Yao [1], which could be called "How to throw a die when all you have is a penny", or something. Its premise is, how can you generate a random variable with arbitrary distribution when all you have is a perfectly fair binary random variate? It begins with this blindingly eye-catching diagram:

Every round node is a coin-toss (a bit-flip), and transitions are taken depending on the value obtained, until one of the dice nodes are reached. What's the expected number of flips needed to get a die? As explained in the paper, it's T = 3 + 2/8⋅T, which solves for T = 4. After this preliminaries, they show us another diagram, but not before observing that [t]he reader may have noticed a source of inefficiency in [the figure]: When three equal values are obtained, we could have used this common result as the value of the next flip, since 000 and 111 are equally probable.

I certainly hadn't; the paper states that the average number of bit-flips needed is 11/3, and that that is minimal. The rest of the paper is devoted to prove constructively that this is so, by giving an algorithm that minimizes the number of bit-flips required to compute an arbitrary random variate. What it doesn't show is how to arrive at this number of flips. Applying the method cited above, we have T = 1 + U and U = 2 + 1/4⋅U, for U = 8/3, T = 11/3.

All fine and dandy but, why am I writing this? For one, to show you pretty finite state machines (actually, Markov chains) that generate dice throws! You might think, what a roundabout way to do it, and counter with:

```let die () = 1 + Random.int 6
```

One problem I see is that, as usually implemented, it is easier to find in the wild fairer bits than dice; or rather, usual implementations of integer random variates tend to take the lower bits from a linear congruential generator, which is a rather lousy way to do it. It is better to do:

```let bit () = Random.float 1. < 0.5
```

which, as it takes into account the more significant bits first, is more apt to be fairer. Another good source of random bits is a linear feedback shift register; the table in [2, 376] shows that x31 + x3 + 1 is a primitive polynomial, modulo 2, which gives the following generator with period 231 - 1:

```let seed = ref 0x33c6ef37

let set_seed s =
seed := if s == 0 then 0x33c6ef37 else s

let rand_bit () =
let b = !seed land 1 in
seed := (!seed lsr 1) lxor (-b land 0x40000004);
b
```

Of course, if you can control your source of bits this disquisition is probably moot; if not, it might be something you'd want to account for. Anyway, back to the bones. There are many ways to write a FSM, but most involve using a table. Instead of falling so low, I'll show you a little functional magic. A state is a finite map from an input symbol to an edge or transition; in turn, an edge can be a way to get to the next state or a final, result symbol:

```type ('a, 'b) state = 'a -> ('a, 'b) edge
and ('a, 'b) edge  = Halt of 'b | Next of ('a, 'b) state
```

I'll represent the second machine as a value of type `(bool, int) state`, with maximally-shared states via value recursion. I number the states from left to right, and from top to bottom:

```let die_fsm : (bool, int) state =
let
rec st0 = function false -> Next st1 | true -> Next st2
and st1 = function false -> Next st3 | true -> Next st4
and st2 = function false -> Next st5 | true -> Next st6
and st3 = function false -> Next st1 | true -> Halt 1
and st4 = function false -> Halt 2   | true -> Halt 3
and st5 = function false -> Halt 4   | true -> Halt 5
and st6 = function false -> Halt 6   | true -> Next st2
in st0
```

The state machine is run simply via the following function:

```let rec run_fsm st input =
match st (input ()) with
| Next st -> run_fsm st input
| Halt x  -> x
```

That's it, really! I can throw my die like so:

```let die () = run_fsm die_fsm (fun () -> rand_bit () == 1)
```

Let's see how fair our die is. For that, I'll compute a histogram of recorded values in an experiment involving n throws:

```let test_die n =
let hist = Array.make 6 0 in
for i = 1 to n do
let d = die () - 1 in
hist.(d) <- hist.(d) + 1
done;
Array.map (fun s -> float s /. float n) hist
```

You can check that the die is fair:

```# test_die 100_000 ;;
- : float array = [|0.16379; 0.16713; 0.16702; 0.16768; 0.1661; 0.16828|]
```

And of course, let's make sure I haven't made a silly mistake:

```# Array.fold_left (+.) 0. (test_die 100_000) ;;
- : float = 1.
```

Now, shoot!

#### References

1. Knuth, D.E.K. and Yao, A. C. The Complexity of Nonuniform Random Number Generation. Reprinted in Donald E. Knuth, Selected Papers on Analysis of Algorithms (Stanford: CSLI lecture notes 102, 2000).
2. Bruce Schneier, Applied Cryptography Second Edition: protocols, algorithms, and source code in C (New York: John Wiley & Sons, Inc. 1996)