A Quantum Diaries Survivor

Many thanks to Riqie Arneberg, who  criticized me for blindly trusting some information on a web site about the odds of picking up 13 cards in a 52-card deck and ending up with 28 or more high-card points (HCP). She was right: if you do not know something and your life is not at stake, go search the web; but if you can find out the solution of your problem by yourself, trust nobody…

For those of you that do not know what is a HCP count: Aces count 4 points, Kings 3, Queens 2 points, and Jacks 1. All other cards do not count. The hand I was dealt yesterday was a whooping 28 HCP, since I had 3 aces, 4 kings, and 2 queens.

So it took me ten minutes of fun this morning to write the following piece of code:

      program prob_hcp
c
c     Computes odds of getting a given HCP in a 13-card hand
c     ——————————————————
      real rndm, dummy
      integer i,j
      integer hcp
      integer x
      real sum_hcp
      real total
      real seen(40)
      integer dealt(13)
c    
c     Start of code - some initializations
c     ————————————
      total = 0
      do hcp = 1, 40
        seen(hcp) = 0
      enddo
c
c     Loop here on hands dealt
c     ————————
 10   sum_hcp = 0
c
c     Refresh deck
c     ————
      do j=1, 13
        dealt(j)=0
      enddo
c
c     Deal hand
c     ———
      do i = 1, 13
c
c       Deal a new card
c       —————
 20     x = 1+ifix(rndm(dummy)*52)
c
c       Check the card is still in the deck
c       ———————————–
        if ( i.gt.1 ) then
          do j = 1, i-1
            if ( x.eq.dealt(j) ) goto 20
          enddo
        endif
c
c       Update deck
c       ———–
        dealt(i)=x
c
c       Add HCP of cards to sum
c       ———————–
        if ( x.le.4 ) then         ! an Ace
          sum_hcp=sum_hcp+4
        else if ( x.le.8 ) then    ! a King
          sum_hcp=sum_hcp+3
        else if ( x.le.12 ) then   ! a Queen
          sum_hcp=sum_hcp+2
        else if ( x.le.16 ) then   ! a Jack
          sum_hcp=sum_hcp+1
        endif
      enddo
c
c     Update total hands dealt
c     ————————
      total=total+1
c
c     Print out results
c     —————–
      if ( mod(total,100000).eq.0 ) then
        print *, ‘****** Total hands dealt = ‘, total, ‘ ******’
      endif
c
c     Print probability that hand has at least hcp points
c     —————————————————
      do hcp = 1, 40
        if ( sum_hcp.ge.hcp ) seen(hcp)=seen(hcp)+1
        if (mod(total,100000).eq.0 ) print *, hcp, seen(hcp)/total
      enddo
      goto 10
c
c     That’s all folks!
c     —————–
      end

Here is the output after 16.7 million hands dealt: 

 ——— k =   16700000.——–
 1  0.99636215
 2  0.988477051
 3  0.974916101
 4  0.950280249
 5  0.911806643
 6  0.86003077
 7  0.79450506
 8  0.714180052
 9  0.625258327
 10  0.531646669
 11  0.437651783
 12  0.348090291
 13  0.267864734
 14  0.198694006
 15  0.141775325
 16  0.0975567698
 17  0.064420417
 18  0.0407968275
 19  0.0247638319
 20  0.0144373057
 21  0.00801041909
 22  0.00424473034
 23  0.00214544917
 24  0.00102688628
 25  0.000462155702
 26  0.000193892221
 27  7.59880204E-05
 28  2.86826344E-05
 29  9.70059864E-06
 30  2.75449111E-06
 31  5.38922166E-07
 32  5.9880243E-08

It transpires that P(>=28HCP) = 28.68E-6, or a bit less than three times in a hundred thousand.