Dice Programs

(Knuth & Yao)

This case study considers two probabilistic programs, due to Knuth and Yao [KY76], which model fair dice using only fair coins. We have reimplemented work done by Joe Hurd, the key difference being that the latter uses a theorem prover (HOL) whereas we use a model checker (PRISM) to show the correctness of such probabilistic programs.

Value of a die

The following program models a die using only fair coins. Starting at the root vertex (state 0), one repeatedly tosses a coin. Every time heads appears, one takes the upper branch and when tails appears, the lower branch. This continues until the value of the die is decided.

The PRISM code for this program is as follows.

// Knuth's model of a fair die using only fair coins
dtmc

module die

	// local state
	s : [0..7] init 0;
	// value of the dice
	d : [0..6] init 0;
	
	[] s=0 -> 0.5 : (s'=1) + 0.5 : (s'=2);
	[] s=1 -> 0.5 : (s'=3) + 0.5 : (s'=4);
	[] s=2 -> 0.5 : (s'=5) + 0.5 : (s'=6);
	[] s=3 -> 0.5 : (s'=1) + 0.5 : (s'=7) & (d'=1);
	[] s=4 -> 0.5 : (s'=7) & (d'=2) + 0.5 : (s'=7) & (d'=3);
	[] s=5 -> 0.5 : (s'=7) & (d'=4) + 0.5 : (s'=7) & (d'=5);
	[] s=6 -> 0.5 : (s'=2) + 0.5 : (s'=7) & (d'=6);
	[] s=7 -> (s'=7);
	
endmodule

rewards "coin_flips"
	[] s<7 : 1;
endrewards

This model has 13 states and in PRISM it takes 4 iterations to find these reachable states.

To prove the above program is correct, we show that the probability of reaching a state where d=k for k=1,...,6 is 1/6. In PRISM this corresponds to calculating the probability of satisfying the formula

P=? [ F s=7 & d=k ]

for k=1..6.

Performing this verification in PRISM using iterative numerical methods, we find that the probability for each k is indeed 1/6 (up to an accuracy of six decimal places), with each case requiring 22 iterations.

As shown in [KY76], this program takes on average 11/3 coin tosses to output a dice throw and this value is optimal.

The expected time can be calculated by in PRISM through the formula

R=? [ F s=7 ]

Performing this verification in PRISM using iterative numerical methods, we find that the expected number of coin tosses is indeed 11/3 (up to an accuracy of six decimal places), requiring 21 iterations.

Sum of two dice

To generate the sum of two dice throws using the above program using PRISM, we combine two such processes in asynchronous parallel composition as follows:

// sum of two dice as the asynchronous parallel composition of
// two copies of Knuth's model of a fair die using only fair coins

mdp

module die1

	// local state
	s1 : [0..7] init 0;
	// value of the dice
	d1 : [0..6] init 0;
	
	[] s1=0 -> 0.5 : (s1'=1) + 0.5 : (s1'=2);
	[] s1=1 -> 0.5 : (s1'=3) + 0.5 : (s1'=4);
	[] s1=2 -> 0.5 : (s1'=5) + 0.5 : (s1'=6);
	[] s1=3 -> 0.5 : (s1'=1) + 0.5 : (s1'=7) & (d1'=1);
	[] s1=4 -> 0.5 : (s1'=7) & (d1'=2) + 0.5 : (s1'=7) & (d1'=3);
	[] s1=5 -> 0.5 : (s1'=7) & (d1'=4) + 0.5 : (s1'=7) & (d1'=5);
	[] s1=6 -> 0.5 : (s1'=2) + 0.5 : (s1'=7) & (d1'=6);
	[] s1=7 & s2=7 -> (s1'=7);
	
endmodule

module die2 = die1 [ s1=s2, s2=s1, d1=d2 ] endmodule

rewards "coin flips"
	[] s1<7 | s2<7 : 1;
endrewards

To prove the above program is correct for calculating the sum of two dice, we show that, both the minimum and maximum probability of reaching a state where s1=7, s2=7 and d1+d2=k for k=2,...,12 is:

k:		probability:
2		1/36
3		1/18
4		3/36
5		1/9
6		5/36
7		1/6
8		5/36
9		1/9
10		3/36
11		1/18
12		1/36

In PRISM this corresponds to verifying the formulae:

Pmin=?[ F s1=7 & s2=7 & d1+d2=k] and Pmax=?[ F s1=7 & s2=7 & d1+d2=k]

for k=2,...,12.

Performing this verification in PRISM we find that the probabilities for each k is as shown in the table and each calculation requires 29 iterations.

Again we can calculate the expected number of coin flips, in this case the minimum and maximum number of coin flips can be calculated through the following formulae:

Rmin=?[ F s1=7 & s2=7 ] and Rmax=?[ F s1=7 & s2=7 ]

Verifying these properties we find that the both the minimum and maximum number of coin flips to find the sum of two dice is 22/3 (two times the expected time to find the value of one dice).

We note that this approach is not optimal and [KY76] shows that the optimal program has an expected 79/18 coin flips and is realized by the following program:

The PRISM code for this program is given below. Note that verifying the formula R=?[ F s=34 ] with PRISM does indeed show that the expected number of coin flips in this case is 79/18 (up to an accuracy of six decimal places), requiring 21 iterations.

// optimal program for the sum of two dice

dtmc

module sum_of_two_dice
	
	// local state
	s : [0..34] init 0;
	// total of two dice
	d : [0..12] init 0;
	
	[] s=0  -> 0.5 : (s'=1)  + 0.5 : (s'=2);
	[] s=1  -> 0.5 : (s'=3)  + 0.5 : (s'=4);
	[] s=2  -> 0.5 : (s'=5)  + 0.5 : (s'=6);
	[] s=3  -> 0.5 : (s'=7)  + 0.5 : (s'=34) & (d'=6);
	[] s=4  -> 0.5 : (s'=8)  + 0.5 : (s'=34) & (d'=7);
	[] s=5  -> 0.5 : (s'=9)  + 0.5 : (s'=34) & (d'=8);
	[] s=6  -> 0.5 : (s'=10) + 0.5 : (s'=11);
	[] s=7  -> 0.5 : (s'=12) + 0.5 : (s'=34) & (d'=4);
	[] s=8  -> 0.5 : (s'=13) + 0.5 : (s'=34) & (d'=5);
	[] s=9  -> 0.5 : (s'=14) + 0.5 : (s'=34) & (d'=9);
	[] s=10 -> 0.5 : (s'=15) + 0.5 : (s'=34) & (d'=10);
	[] s=11 -> 0.5 : (s'=16) + 0.5 : (s'=17);
	[] s=12 -> 0.5 : (s'=18) + 0.5 : (s'=34) & (d'=3);
	[] s=13 -> 0.5 : (s'=19) + 0.5 : (s'=34) & (d'=5);
	[] s=14 -> 0.5 : (s'=20) + 0.5 : (s'=34) & (d'=7);
	[] s=15 -> 0.5 : (s'=21) + 0.5 : (s'=34) & (d'=9);
	[] s=16 -> 0.5 : (s'=22) + 0.5 : (s'=34) & (d'=11);
	[] s=17 -> 0.5 : (s'=23) + 0.5 : (s'=24);
	[] s=18 -> 0.5 : (s'=25) + 0.5 : (s'=34) & (d'=2);
	[] s=19 -> 0.5 : (s'=26) + 0.5 : (s'=34) & (d'=3);
	[] s=20 -> 0.5 : (s'=27) + 0.5 : (s'=34) & (d'=4);
	[] s=21 -> 0.5 : (s'=34) & (d'=5) + 0.5 : (s'=34) & (d'=9);
	[] s=22 -> 0.5 : (s'=28) + 0.5 : (s'=34) & (d'=10);
	[] s=23 -> 0.5 : (s'=29) + 0.5 : (s'=34) & (d'=11);
	[] s=24 -> 0.5 : (s'=30) + 0.5 : (s'=34) & (d'=12);
	[] s=25 -> 0.5 : (s'=1)  + 0.5 : (s'=34) & (d'=2);
	[] s=26 -> 0.5 : (s'=31) + 0.5 : (s'=34) & (d'=3);
	[] s=27 -> 0.5 : (s'=32) + 0.5 : (s'=34) & (d'=6);
	[] s=28 -> 0.5 : (s'=34) & (d'=7) + 0.5 : (s'=34) & (d'=8);
	[] s=29 -> 0.5 : (s'=33) + 0.5 : (s'=34) & (d'=11);
	[] s=30 -> 0.5 : (s'=2)  + 0.5 : (s'=34) & (d'=12);
	[] s=31 -> 0.5 : (s'=34) & (d'=2) + 0.5 : (s'=34) & (d'=4);
	[] s=32 -> 0.5 : (s'=34) & (d'=6) + 0.5 : (s'=34) & (d'=8);
	[] s=33 -> 0.5 : (s'=34) & (d'=10) + 0.5 : (s'=34) & (d'=12);
	[] s=34 -> (s'=34);
	
endmodule

rewards "coin_flips"
	[] s<34 : 1;
endrewards

Site hosted at the Department of Computer Science, University of Oxford