What is the probability of dealing all $52$ cards in a standard well shuffled deck getting a single quad? I would like to know if you take a well shuffled deck of $52$ cards and then deal out all of them one at a time without replacement, what is the probability that they will be dealt such that there is exactly one quad dealt in order such as $4$ consecutive kings?  That run of $4$ has to be the only run of that length, there cannot be any others.  For example, 5, ... ,K,K,K,K,7,J,J,J,J is a "loser".  Other than computer simulation I am not sure how to solve this.

It seems that your shuffling algorithm is too slow. For the first card, you can pick a number n between 1 and 52. If it's 1, your card is where it's supposed to be. Otherwise, swap cards 1 and n. Then for card 2, pick a number between 2 and 52...  While swapping seems to double the amount of movement, it means that the card to be chosen comes from a tight array so you don't have to search for it and this will make shuffling much faster.  
For detecting 4 in a row, my idea is like punching holes in data cards so that if light shines through 4 cards, it's a match. The digital analog to this is setting a bit according to the rank of the card and ANDing four consecutive values together to see if one bit remains set.  
Processed a billion simulations in 344 seconds with 1, 2, and 3 matches occurring with frequency 2347407, 3090, and 6 respectively.
program four
   implicit none
   real harvest(51)
   integer i, j
   integer, parameter :: N(51) = [(i,i=52,2,-1)]
   integer indices(51)
   integer matches(13)
   integer jmax
   integer cards(52)
   integer, parameter :: init(52) = [((ishft(1,j),j=0,12),i=1,4)]
   integer match
   integer(selected_int_kind(18)) t0,tf,rate

   matches = 0
   jmax = 1.0e9
   call random_seed()
   call system_clock(t0)
   do j = 1, jmax
      call random_number(harvest)
      indices = harvest*N
      cards = init
      do i = 1, 51
         if(indices(i) /= 0) cards([i,indices(i)+i]) = cards([indices(i)+i,i])
      end do
      match = 0
      do i = 1, 49
         if(iall(cards(i:i+3)) /= 0) match = match+1
      end do
      if(match > 0) matches(match) = matches(match)+1
   end do
   write(*,'(a,1x,i0)') 'Number of simulations = ',jmax
   write(*,'(a)') ' N Frequency'
   do i = 1, 13
      write(*,'(i2,1x,i0)') i , matches(i)
   end do
   call system_clock(tf)
   call system_clock(count_rate=rate)
   write(*,'(a,1x,g0)') 'Time in seconds = ', real(tf-t0)/rate
end program four

