function [Pexact,Pa6exact,Pnumber,Pa6number]=crapsmof65_ace6(number, H)
% format [Pexact,Pa6exact,Pnumber,Pa6number]=crapsmof65_ace6(number, H)
% input: number=value of 1st roll
%        H    ='H' for Pnumber & odds of making 4, 6, 8, or 10 the hard way
% output: Pexact= Exact probability of winning, given 1st roll = number
%         Pnumber=Probability of rolling that number by rolling 2 fair dice.
%         Pa6=Probability of rolling that number with loaded ace 6 flats.
% Using the 'method of reduced sample space' from Mosteller, F. 1965. Fifty
%   challenging problems in probability with solutions, p. 9
% see craps.m for a more laborious solution using absorbing Markov chains
% written 6/2/03 by E. Gallagher for EEOS601, revised 1/12/11;
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% See Larsen & Marx Example 3.3.6 Ace-six flats are foreshortened so that
% the p(1)=p(6)=1/4 & p(2)=p(3)=p(4)=p(5)=1/8;
if nargin > 1 
    if H=='H'
        if mod(number,2) | number ==2 | number==12 
            error(...
              'Can only make even numbers between 4 & 10 the ''hard'' way')
            return
        end
    end
end
dief=1:6; % create a sample space of two dice rolls; while simultaneously
         % creating a matrix of probabilities corresponding to the faces
diep=[1/4 repmat(1/8,1,4) 1/4];
die1f=repmat(dief,6,1); % create a matrix, with each row equal to 1 to 6;
die1p=repmat(diep,6,1); % create a matrix, with each row equal to 1 to 6;
die2f=repmat(dief',1,6); % create a matrix, with each column equal to 1 to 6;
die2p=repmat(diep',1,6); % create a matrix, with each column equal to 1 to 6;
sumdicef=die1f+die2f;  % This gives the sample space, the 36 sums of two die.
sumdicep=die1p.*die2p;  % This gives the sample space,the 36 elementwise 
                        % products two die.
i=find(sumdicef==number);  % find the indices for each value equal to number;
lengthi=length(i);     % find how many sums of 2 die equal the number in input;
if lengthi==0          % only true if the number isn't in the sample space.
    disp('Your number isn''t possible, not in the sample space')
    number
    % The following odd syntax will print the number
    fprintf('Your number, %6.0f, isn''t possible, not in the sample space\n',number)
    return % This stops the program after displaying the message;
end
[r,c]=size(sumdicef); % find size of sample space
Pnumber=lengthi/(r*c); % What is the probability of rolling your number with 1 roll
                       % of two dice
Pa6number=sum(sumdicep(i)); % sum the probabilities corresponding to the summed faces
% Is the craps shooter an immediate winner?
if (number==7 | number==11) % 2 logical statements using the | OR command.
    Pexact=1;              % Matlab performs statements in this if section only
                            % if the logical statement number equals 7 or equals 11
                            % is true.
    Pa6exact=1;
    return
end
% Is the craps shooter an immediate loser?
if (number==2 | number==3 | number==12)
    Pexact=0;
    Pa6exact=0;
    return  % exit the function and return to the command window.
end
% Exact probability of making a point - no need for sum of geometric series
% or Markov chains
% Use Mosteller's 'Method of reduced sample space: 
% all we need is Probability of rolling a given
% number - which wins - and the probability of rolling a 7, which loses.
%  We don't need any of the other probabilities 
% (e.g., probability of rolling 2, 3 ...)
%  First find probability of rolling a 7, which results in a loss.
j=find(sumdicef==7);P7=length(j)/(r*c); % P7 is the probability of rolling 7
Pa6_7=sum(sumdicep(j));
if nargin>1 & H=='H'
    ih=find(sumdicef==number & die1f==number/2);
    lengthih=length(ih);
    Pnumberh=lengthih/(r*c); % Another way to lose is to roll a soft number
    Pnumbera6h=sum(sumdicep(ih))
    Pexact=Pnumberh/(P7+Pnumber); % You win or lose if you roll your hard 
                                  % number (4,6,8)
    Pa6exact=Pnumbera6h/(Pa6_7+Pa6number);
else
    Pexact=Pnumber/(Pnumber+P7);
    Pa6exact=Pa6number/(Pa6_7+Pa6number);
end
