half pills

You have a prescription to take one half of a pill per day for 20 days but the pharmacist (who is too busy to divide pills for you) gives you 10 whole pills in a bottle. On day 1, you remove a pill from the bottle, break it into two half pills, take one, and return the other half pill to the bottle.On all subsequent days you shake the bottle thoroughly and pour something out  whatever comes out first - either a half pill or a whole pill; if it' s a half pill you take it and you' re done for that day; if it' s a whole pill, you split it into two half - pills, take one, and put the other back in the bottle, exactly like you did on day 1. On day 20 there can be only one half pill left in the bottle, but on day 19 there are two possibilities : either there is one whole pill or there are two half - pills left in the bottle.What is the probability that there are two half - pills in the bottle on day 19?

Solution by Michael A. Gottlieb
One way to solve this problem approximately is to use difference equations: If there are w whole pills and h half pills in the bottle on a given day, then the probability of drawing a whole pill equals w/(w+h). One can think of this probability as the “mean whole pills” leaving the bottle by considering w and h not as integer numbers of pills but as real-valued (mean) quantities of pills. Then one can approximate that on the following day there will be, on average, w + dw whole pills in the bottle, with dw = - w/(w+h). The case for half pills is similar, but a bit more complicated because if we draw a whole pill it adds a new half pill to the bottle, while if we draw a half pill it subtracts one, so dh = w/(w+h) - h/(w+h) = (w-h)/(w+h). Thus an approximate solution can be found by the following procedure:

w = W (initial number of whole pills)
h = H   (initial number of half pills)

Repeat (2W +H - 2) times:
      w = w - w/(w+h)
      h = h + (w-h)/(w+h)	
	
The probability of two half pills on the second to last day is then h/(w+h).

w = W = 10.0;  h = H = 0.0;
Do[
	dw = -w/(w+h);
	dh =  (w-h)/(w+h);
	w +=dw;
	h  +=dh;
,{2W + H-2}
];
h/(w+h)

The true probabilty is (21937801980489931824271/28810829817907200000000) ≈ 0.761443, which differs from the value calculated here by about 0.033.  This method is pretty fast, but it’s running time is still proportional to 2W+H-2 (the number of times it iterates). It always overestimates the true probability and works better for larger numbers of pills. For example, when W=30 (and H=0), it yields a probability of 0.836144, wheras the true probability is 0.820834, a difference of 0.015.