function [probM, expected,P] = solve_hhh_tthh(p_coin)
% Solve the HHH vs TTHH coin-tossing problem,
% using a coin that lands heads with probability p_coin
% (with prob = 1/2 as the default).  
% 
% The matrix P contains the  probabilities for
% transitions between the states, with rows and columns
% in the order of the states 
% 	S, H, T, HH, TT, TTH, Mwins, Rwins.
% 
% The return values give the probability that M wins, 
% starting from each particular state, and the expected 
% number of tosses needed to complete the game.

if(nargin <1 )
	p_coin = 1/2
end

states = {'S' 'H' 'T' 'HH' 'TT' 'TTH' 'Mwins' 'Rwins'};

P= [0 p_coin (1-p_coin) 0 0 0 0 0;   % S
	 0 0 (1-p_coin) p_coin 0 0 0 0;   % H
	 0 p_coin 0 0 (1-p_coin) 0 0 0;   % T
	 0 0 (1-p_coin) 0 0 0 p_coin 0;	% HH
	 0 0 0 0 (1-p_coin) p_coin 0 0;	% TT
	 0 0 (1-p_coin) 0 0 0 0 p_coin;	% TTH
	 0 0  0 0 0 0 1 0;	 	% Mwins
	 0 0  0 0 0 0 0 1;	 	% Rwins
	 ];
	 
% Solve the equations P * probM = probM 
% with boundary conditions probM(7) = 1, probM(8) = 0

probM = [eye(8) - P; P(7:8,:)]\[0 0 0 0 0 0 0 0 1 0]';

% Equations 7 and 8 are trivial.  Maybe better to discard them:
% 
% probM = [eye(6,8) - P(1:6,1:8); P(7:8,:)]\[0 0 0 0 0 0 1 0]';
% 
% Alternatively, solve only for the nontrivial probM values:
% 
% probM = [eye(6) - P(1:6,1:6)]\[0 0 0 1/2 0 0 ]';
% 


% The expected values solve the equations P*(1+expected) = expected
% with the boundary conditions expected(7) = expected(8) = 0.

expected = [eye(6,8) - P(1:6,1:8); P(7:8,:)]\[1 1 1 1 1 1 0 0]';



% print a little table of the results

fprintf (1, '%10s %10s %10s\n','state','probM','expected');
for i =1:8
fprintf (1, '%10s %10.3f %10.3f\n', states{i}, probM(i),expected(i));
end
 
