# program that compute the approximate probability of
# drawing a k-card hand having at least one ace
#
# input:  k --- number of cards in a hand
#         n --- number of trials for the experiment
#
# This program require draw_k_card package as subprogram
#
get_1_ace:=proc(k::integer, # number of cards drawn
                     n::integer  # number of trials
                    )
   local count, i, j, hand, flag, cnt;

   count:=0;   # initialize counting
   for i from 1 to n do

      #
      # simulate the experiment
      #
      draw_k_card(k,hand);  # draw a k card hand

      #
      # checking
      #
      cnt:=0;              # initialize the flag=0, meaning no ace
      for j from 1 to k do
         if hand[j][1] = "ace" then
            cnt:=cnt+1;        # raise the flag, an ace is in the hand
         fi;
      od;
      
      if cnt=1 then
         count:=count+1;    # count
      fi;
   od;

   evalf(count/n);          # calculate the probability
end;

#
################## subprogram ###################################
#
# program that drawing k cards, without replacement
# out a standard deck
#
#  input:  k    --- the number of balls to be drawn
#  output: hand --- array that contains the names of cards drawn
#
draw_k_card:=proc(
                  k::integer, # number of cards to draw
               hand::evaln    # output: the drawn hand
                  )
   local deck, i, m, n, temp;

   if k < 1 or k > 52 then   # quick exit if k < 1 or k > 52
      print(`1st argument must be between 1 and 52`);
      RETURN();
   fi;

   hand:=array(1..k);       # open the space for output
   deck:=array(1..52,       # making a pool of numbers to draw
         [seq(i,i=1..52)]);

   n:=52;                   # the number of cards remaining

   for i from 1 to k do     # loop to draw card 1 by 1

      m:=int_ran(1,n);      # draw a number from 1 to n
      hand[i]:=deck[m];

      if m<>n then          # move the number to the end
         temp:=deck[m];     # of the pool to avoid redrawn
         deck[m]:=deck[n];
         deck[n]:=temp;
      fi;

      n:=n-1;               # update n

   od;

   #
   # get the names of the cards
   #
   for i from 1 to k do
      m:=hand[i];         # get the number
      hand[i]:=card(m);   # call subprogram that identifies the card
   od;

end;
################# end of the main program #############################

###################### subprogram #####################################
#
# function that identify the card name given its number
# input: m --- a number from 1 to 52
# output: an array of two entries that identifies the card
#
card:=proc(m::integer)
   local rank, suit;
   if m>0 and m<= 52 then  # quick exit if m is beyond 1-52

      rank:= m mod 13;     # get numerical rank (1-13)
      suit:=ceil(m/13);    # get numerical suit (1-4)
      #
      # identifies the actual suit as a string
      #
      if suit=1 then
         suit:="heart";
      elif suit=2 then
         suit:="spade";
      elif suit=3 then
         suit:="diamond";
      else
         suit:="club";
      fi;
      #
      # identifies the actual rank
      #
      if rank=1 then
         rank:="ace";
      elif rank=11 then
         rank:="jack"
      elif rank=12 then
         rank:="queen"
      elif rank=0 then
         rank:="king"
      fi;                # remaining cases are rank=2,3,4,...,10,
                         # already defined
      #
      # output the card name as an array
      #
      array(1..2,[rank,suit])
   fi;
end;
############# end of subprogram card #################################


###################### subprogram #####################################
#
#
# function that generates a random integer between m and n, inclusive
#
int_ran:=proc(m::integer,n::integer)
   round(evalf(m-0.5+(n-m+1)*rand()/999999999999))
end;
############# end of subprogram int_ran ###############################