# A Big Problem That Doesn't Need a Bignum

When I finally wrapped up the mathematical modeling of the k-consecutive heads coin toss problem I had a concise formula:

To get the exact result I needed, I used the BigInteger and BigDecimal Java classes to handle those big numbers, and had to wait about an hour to get the result. As it turns out, there is a much faster and easier way to do things

#### Simplifying the Computation

To solve this problem in a more reasonable matter, take a look at the part of the equation where all the real work is done:*p*we are solving for is the probability that

*no*runs of

*k*consecutive heads will occur in

*n*tosses of a fair coin. Although we might be able to find a closed-form equation that defines the k-step Fibonacci number, it is going to be a complicated polynomial, and the chances of us simplifying this equation are effectively nil.

So at first blush I seem to be stuck with calculating the top and bottom parts of the fraction independently. But a little inspection of the equation shows that I can actually calculate *p* in an iterative fashion, starting with *n*=1, and working my way up to the desired result.

#### Calculating the Fibonacci Number

The first step in the calculation is writing the code to calculate k-step Fibonacci numbers. Recall that the k-step Fibonacci number*n*is found by summing the values of the previous

*k*numbers - a recursive definition. The base cases for all values of

*k*are: Fib

_{k}(n) is 1 when n is 1, and 0 when n is less than 1.

To determine the value of Fib_{k} in an iterative fashion, I make use of the the handy C++ container `deque<T>`

. To calculate a *k*-step Fibonacci number, I make use of a `deque<double>`

of size *k* that contains the previous k Fibonacci numbers. With the data structure set up like that, calculating the next Fibonacci number is simple:

q.push_front( accumulate( q.begin(), q.end(), 0 ) ); q.pop_back();Note that after this calculation, q[ 0 ] contains the next new Fibonacci number, and the previous

*k-1*numbers are still stored in sequence in the

`deque`

. The call to `pop_back()`

discards Fib_{k}(n-k), which is not going to be needed for the next calculation.

Add the code to initialize the `deque`

, and a loop, and we can print out Fib_{k}(n) for values 1 through *n*:

deque<double> q(k, 0); q[0] = 1.0; for ( int i = 1 ; i < n ; i++ ) { cout << q[ 0 ] << " "; q.push_front( accumulate( q.begin(), q.end(), 0.0 ) ); q.pop_back(); } cout << q[ 0 ] << endl;Using the familiar value of

*k*=2, and

*n*=20, I get the following output:

1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765The

`deque`

is perfect for this class - it lets us simulate a rolling window through the Fibonacci sequence by pushing new values on one end and popping aged ones off the other end. All this while providing constant time access to indexed elements. (Admittedly this code could be written just as efficiently with a `list<double>`

.)
#### Divide By Two

This is a nice simple calculation, but it only gets me halfway to the goal. I'm calculating the top part of the fraction that computes*p*, but I need to divide by the bottom value as well. As it turns out, I can accomplish this in a way that allows me to iterate to the next value: before I add the values of the previous Fibonacci numbers to

`q[0]`

, I just divide all of them by two. The resulting loop ends up looking like this:
deque<double> q(k, 0); q[0] = 1.0; q[1] = 1.0; for ( int i = 1 ; i < n ; i++ ) { q.push_front(0.0); for ( int j = 1 ; j <= k ; j++ ) { q[ j ] /= 2.0; q[ 0 ] += q[ j ]; } q.pop_back(); } cout << q[ 0 ] << endl;Note that I've slightly changed the initialization - the

`deque<double>`

is initialized with Fib_{k}(2) - this is because the top half of the fraction uses Fib

_{k}(n+2), and the bottom uses 2

^{n}- the value of n for the Fibonacci calculation always has to be two ahead of the exponent of 2.

This calculation allows me to calculate the probabilities for large numbers of tosses in a matter of seconds, instead of the hour-plus it took using BigInteger. Using this code I was able to get the final word on just when we could be certain that a certain number of coin tosses would produce a run of 20 heads. If we want to be modestly certain, with 90% probability I can say that 4,828,842 coin tosses will produce a run of 20 heads.

If my life depended on it, and I wanted to get out to five nines, I can say that 24,144,137 coin tosses will produce a run of 20 heads with 99.999% probability. The output of kfib.cpp is shown here:

There is a 90% probability of seeing 20 consecutive heads in 4828842 tosses of a fair coin There is a 99% probability of seeing 20 consecutive heads in 9657666 tosses of a fair coin There is a 99.9% probability of seeing 20 consecutive heads in 14486490 tosses of a fair coin There is a 99.99% probability of seeing 20 consecutive heads in 19315313 tosses of a fair coin There is a 99.999% probability of seeing 20 consecutive heads in 24144137 tosses of a fair coinAnd that's the last thing I am going to say about it.

#### Note for DDJ Readers

I read all the comments posted on my DDJ blog articles, but you should know that unfortunately, this is done without the support of a notification system. So if you don't see a response right away, be patient - I'll be making a pass through soon.