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

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

[ Editor: Conway's algorithm for the odds that one sequence appears
before another, from a 1998 thread "coin flipper & independent sets". ]

Nick Wedd <Nick@maproom.co.uk> said:
> I propose to Jamie that we play the following game.  We toss a fair coin
> until either
> THHHT comes up, when Jamie wins and I pay him a dollar;  or
> HTHHH comes up, when I win and Jamie pays me a dollar.
> These two are equally probable, so the game is, I allege, fair.
>
> ObPuzzle:  what is the value to me of playing this game once?

In "Time Travel and Other Mathematical Bewilderments", Martin Gardner
describes a procedure for solving this problem, attributed to J. H.
Conway.  He says "I have no idea why it works.  It just cranks out the
answer as if by magic, like so many of Conway's other algorithms."

Let A = THHHT and B = HTHHH.  Write the four sets of pairs of A and B, one
line above the other.  Write a binary digit on top of each of the pairs,
where putting a 1 in the n'th place means "The last n digits of the top
row are the same as the first n digits of the bottom row."  So, it looks
like this:

   10001
A: THHHT
A: THHHT

   00010
A: THHHT
B: HTHHH

   01000
B: HTHHH
A: THHHT

   10001
B: HTHHH
B: HTHHH

Take each number as an n-digit binary number.  So AA=17, AB=2, BA=8, and
BB=17.  The odds that B will come up before A is AA-AB:BB-BA = 15:9 = 5:3.
(This corresponds to the 5/8 probability that Glenn came up with
independently.)

ObPuzzle: How would the mathematics change if you had more than two people
with winning combinations?  How would they change if the two strings were
of different lengths?

-- Matthew Daly <mwdaly@pobox.com>

[ Editor: the n'th place is counted from the right, as the four
examples show.  Glenn's posting (probably Glenn Rhoads's) is not in the
collected material.  With B (Wedd's sequence) first with probability
5/8, the game is worth 5/8 - 3/8 = $0.25 to Wedd per play. ]