% octave/matlab function to compute probability that M blackjack trials will produce
% an estimate which rounds to 8, where 8/169 is the chance of getting blackjack
% on two cards drawn from an infinite deck:
%  169 = 13*13 possibilities, A...10,J,Q,K each card
%  8 ways to get blackjack: A and 10,J,Q,K or 10,J,Q,K and A
% est = 169*S/M for S blackjacks in M Monte-Carlo trials
% P( 7.5 <= est <= 8.5) = P( 7.5*M/169 <= S <= 8.5*M/169) = P( low <= S <= high)
%
function p = bj(M)
 low = ceil(7.5*M/169);
 high = floor(8.5*M/169);
 p = sum( binopdf([low:high],M,8/169) );
end