C compute the odds of getting all four suits after a deal of n cards, where n is from 4 to 39.

	real rfact
	real*8 rfact1, prob(1000), tprob(40), cumprob
	dimension icomb(1000)

	print *, ' Begin computation of probability of reaching all 4 '
	print *, '   suits after receiving n cards.'

c first find the number of combinations to reach n-1 cards with only 3 suits
	do n = 4, 40
c	do n = 4, 6
	  nm1 = n - 1
	  ic = 0
	  iend = n-3
	  if(iend .gt. 13) iend = 13
	  jend = n-3 
	  if(jend .gt. 13) jend = 13
	  kend = n-3
	  if(kend .gt. 13) kend = 13	
	  do i = 1, iend
	    do j = 1, jend
	      do k = 1, kend
		if (i+j+k .eq. n-1) then
c sort small to large, i, j, k
		  ni = i
		  nj = j
		  nk = k
		  if (ni .gt. nj) then
		    nit = ni
		    ni = nj
		    nj = nit
		  endif
		  if (ni .gt. nk) then
		    nit = ni
		    ni = nk
		    nk = nit
		  endif
		  if (nj .gt. nk) then
		    njt = nj
		    nj = nk
		    nk = njt
		  endif
		  ic = ic + 1
		  icomb(ic) = 10000*ni+100*nj + nk
		endif
		
	      enddo
	    enddo
	  enddo
	  icf = 1
	  do i = 1, ic
	    do j = i+1, ic
  	      if (icomb(j) .eq. icomb(i)) icomb(j) = 0
	    enddo
	  enddo
	  nprob = 0
	  tprob(n) = 0.
	  do i = 1, ic
	    if(icomb(i) .ne. 0) then
		nprob = nprob + 1
c compute probability
		n1 = icomb(i)/10000
		n2 = (icomb(i) - n1*10000)/100
		n3 = icomb(i) - n1*10000 - n2*100
c		print *, ' *** n, n1, n2, n3', n, n1, n2, n3, '***'
		prob(nprob) = rfact1(nm1,n1)*rfact1(nm1-n1,n2)*
     *            rfact1(nm1-n1-n2,n3)*4*3*2/
     *            (rfact(n1)*rfact(n2)*rfact(n3))
		if ((n1 .eq. n2) .and. (n1 .eq. n3)) then
		  prob(nprob) = prob(nprob)/6.
		else if( (n1 .eq. n2) .or. (n1 .eq. n3) .or. 
     *	          (n2 .eq. n3)) then
	          prob(nprob) = prob(nprob)/2.
		endif
c234567890123456789012345678901234567890c234567890123456789012345678901234567890
c        1         2         3         4         5         6         7

c		print *, 'n, icomb, # of ways to fill ', n, icomb(i), 
c     *            prob(nprob)	
		prob(nprob) = prob(nprob)*rfact1(13,n1)*rfact1(13,n2)
     *            *rfact1(13,n3)*13/rfact1(52,n)
		tprob(n) = tprob(n) + prob(nprob)
c	      print *, 'n, icomb, prob ', n, icomb(i), prob(nprob)
	    endif
	  enddo
c	  print *, 'n, total probability ', n, tprob(n)
	enddo
	cumprob = 0.
	do n = 4, 40
	  cumprob = cumprob + tprob(n)
	  print *, n, tprob(n), cumprob
	enddo
	stop
	end

	real function rfact(n)
	rfact = 1
	if (n .lt. 0) then
	  print *, ' can not find factorial of a negative number'
	  stop
	else if(n .eq. 0) then
	  rfact = 1
	else if(n .eq. 1) then
	  rfact = 1
	else
	  do i = n, 2, -1
	    rfact = rfact*i
	  enddo
	endif
	return
	end

	real*8 function rfact1(n,m)
	rfact1 = 1
	if (m .le. 0) then
	  print *, ' rfact1(n,m) m must be .ge. 1'
	  stop
	if (n .eq. m) rfact1 = 1
	else if(n .lt. m) then
	  print *, ' rfact1(n,m): n must be .ge. than m'
	  print *, n, m
	  stop
	else
	  do i = n, n-m+1, -1
	    rfact1 = rfact1*i
	  enddo
	endif
	return
	stop
	end