==> 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