# 406 probability/flips/waiting.time.p

## Description

This article is from the Puzzles FAQ, by Chris Cole chris@questrel.questrel.com and Matthew Daly mwdaly@pobox.com with numerous contributions by others.

# 406 probability/flips/waiting.time.p

Compute the expected waiting time for a sequence of coin flips, or the
probabilty that one sequence of coin flips will occur before another.

probability/flips/waiting.time.s

Here's a C program I had lying around that is relevant to the
current discussion of coin-flipping. The algorithm is N^3 (for N flips)
but it could certainly be improved. Compile with

cc -o flip flip.c -lm

-- Guy

```_________________ Cut here ___________________

#include <stdio.h>
#include <math.h>

char *progname;              /* Program name */

#define NOT(c) (('H' + 'T') - (c))

/* flip.c -- a program to compute the expected waiting time for a sequence
of coin flips, or the probabilty that one sequence
of coin flips will occur before another.

Guy Jacobson, 11/1/90
*/

main (ac, av) int ac; char **av;
{
char *f1, *f2, *parseflips ();
double compute ();

progname = av[0];

if (ac == 2) {
f1 = parseflips (av[1]);
printf ("Expected number of flips until %s = %.1f\n",
f1, compute (f1, NULL));
}
else if (ac == 3) {

f1 = parseflips (av[1]);
f2 = parseflips (av[2]);

if (strcmp (f1, f2) == 0) {
printf ("Can't use the same flip sequence.\n");
exit (1);
}
printf ("Probability of flipping %s before %s = %.1f%%\n",
av[1], av[2], compute (f1, f2) * 100.0);
}
else
usage ();
}

char *parseflips (s) char *s;
{
char *f = s;

while (*s)
if (*s == 'H' || *s == 'h')
*s++ = 'H';
else if (*s == 'T' || *s == 't')
*s++ = 'T';
else
usage ();

return f;
}

usage ()
{
printf ("usage: %s {HT}^n\n", progname);
printf ("\tto get the expected waiting time, or\n");
printf ("usage: %s s1 s2\n\t(where s1, s2 in {HT}^n for some fixed n)\n",
progname);
printf ("\tto get the probability that s1 will occur before s2\n");
exit (1);
}

/*
compute --  if f2 is non-null, compute the probability that flip
sequence f1 will occur before f2.  With null f2, compute
the expected waiting time until f1 is flipped

technique:

Build a DFA to recognize (H+T)*f1 [or (H+T)*(f1+f2) when f2
is non-null].  Randomly flipping coins is a Markov process on the
graph of this DFA.  We can solve for the probability that f1 precedes
f2 or the expected waiting time for f1 by setting up a linear system
of equations relating the values of these unknowns starting from each
state of the DFA.  Solve this linear system by Gaussian Elimination.
*/

typedef struct state {
char *s;                /* pointer to substring string matched */
int len;                /* length of substring matched */
int backup;             /* number of one of the two next states */
} state;

double compute (f1, f2) char *f1, *f2;
{
double solvex0 ();
int i, j, n1, n;

state *dfa;
int nstates;

char *malloc ();

n = n1 = strlen (f1);
if (f2)
n += strlen (f2); /* n + 1 states in the DFA */

dfa = (state *) malloc ((unsigned) ((n + 1) * sizeof (state)));

if (!dfa) {
printf ("Ouch, out of memory!\n");
exit (1);
}

/* set up the backbone of the DFA */

for (i = 0; i <= n; i++) {
dfa[i].s = (i <= n1) ? f1 : f2;
dfa[i].len = (i <= n1) ? i : i - n1;
}

/* for i not a final state, one next state of i is simply
i+1 (this corresponds to another matching character of dfs[i].s
The other next state (the backup state) is now computed.
It is the state whose substring matches the longest suffix
with the last character changed */

for (i = 0; i <= n; i++) {
dfa[i].backup = 0;
for (j = 1; j <= n; j++)
if ((dfa[j].len > dfa[dfa[i].backup].len)
&& dfa[i].s[dfa[i].len] == NOT (dfa[j].s[dfa[j].len - 1])
&& strncmp (dfa[j].s, dfa[i].s + dfa[i].len - dfa[j].len + 1,
dfa[j].len - 1) == 0)
dfa[i].backup = j;
}

/* our dfa has n + 1 states, so build a system n + 1 equations
in n + 1 unknowns */

eqsystem (n + 1);

for (i = 0; i < n; i++)
if (i == n1)
equation (1.0, n1, 0.0, 0, 0.0, 0, -1.0);
else
equation (1.0, i, -0.5, i + 1, -0.5, dfa[i].backup, f2 ? 0.0 : -1.0);
equation (1.0, n, 0.0, 0, 0.0, 0, 0.0);

free (dfa);

return solvex0 ();
}

/* a simple gaussian elimination equation solver */

double *m, **M;
int rank;
int neq = 0;

/* create an n by n system of linear equations.  allocate space
for the matrix m, filled with zeroes and the dope vector M */

eqsystem (n) int n;
{
char *calloc ();
int i;

m = (double *) calloc (n * (n + 1), sizeof (double));
M = (double **) calloc (n, sizeof (double *));

if (!m || !M) {
printf ("Ouch, out of memory!\n");
exit (1);
}

for (i = 0; i < n; i++)
M[i] = &m[i * (n + 1)];
rank = n;
neq = 0;
}

/* add a new equation a * x_na + b * x_nb + c * x_nc + d = 0.0
(note that na, nb, and nc are not necessarily all distinct.) */

equation (a, na, b, nb, c, nc, d) double a, b, c, d; int na, nb, nc;
{
double *eq = M[neq++]; /* each row is an equation */
eq[na + 1] += a;
eq[nb + 1] += b;
eq[nc + 1] += c;
eq[0] = d;             /* column zero holds the constant term */
}

/* solve for the value of variable x_0.  This will go nuts if
therer are errors (for example, if m is singular) */

double solvex0 ()
{
register i, j, jmax, k;
register double  max, val;
register double *maxrow, *row;

for (i = rank; i > 0; --i) {      /* for each variable */

/* find pivot element--largest value in ith column*/
max = 0.0;
for (j = 0; j < i; j++)
if (fabs (M[j][i]) > fabs (max)) {
max = M[j][i];
jmax = j;
}
/* swap pivot row with last row using dope vectors */
maxrow = M[jmax];
M[jmax] = M[i - 1];
M[i - 1] = maxrow;

/* normalize pivot row */
max = 1.0 / max;
for (k = 0; k <= i; k++)
maxrow[k] *= max;

/* now eliminate variable i by subtracting multiples of pivot row */
for (j = 0; j < i - 1; j++) {
row = M[j];
if (val = row[i])              /* if variable i is in this eq */
for (k = 0; k <= i; k++)
row[k] -= maxrow[k] * val;
}
}

/* the value of x0 is now in constant column of first row
we only need x0, so no need to back-substitute */

val = -M[0][0];

free (M);
free (m);

return val;
}

_________________________________________________________________
Guy Jacobson   (201) 582-6558              AT&T Bell Laboratories
uucp:  {att,ucbvax}!ulysses!guy	   600 Mountain Avenue
internet:  guy@ulysses.att.com         Murray Hill NJ, 07974
```

Continue to: