Casino game 30 and 40

###  Exact P(31-40 in k cards drawn from 6 deck shoe)

require(compiler)
enableJIT(3)

max.cards = 28
p = rep(0,28)
for (k in 4:max.cards) {
  m = combn(10+k-1,k)   # hands with k cards -> columns of m in * and | form
  hands.tot = 0
  hands.gt.30 = 0
  for (j in 1:choose(10+k-1,k) ) {                  # for each hand
    cards = rle(cumsum(c(m[1,j],diff(m[,j]) - 1)))  # convert cards to run length form
    total = sum(cards$values * cards$lengths)       # sum of cards
    hands = 1
    if ( any(cards$values < 10) )
      hands = hands*prod(choose(24,cards$lengths[cards$values < 10]))
    if ( any(cards$values == 10) )
      hands = hands*choose(96,cards$lengths[cards$values == 10])
    if (total > 30) hands.gt.30 = hands.gt.30 + hands
      hands.tot = hands.tot + hands
  }
  p[k] = hands.gt.30/hands.tot - sum(p[3:(k-1)])
  cat('\ncards =',k, '  prob =',p[k])
  flush.console()
  remove(m)
}

###  Simulation of P(31-40 in k cards drawn from 6 deck shoe)

require(compiler)
enableJIT(3)

shoe = c(rep(1:9,24), rep(10,96))
count = rep(0,28)
sims = 0
while(1) {
  cards = sample(shoe,28,replace=FALSE)
  num = which(cumsum(cards) > 30)[1]
  count[num] = count[num] + 1
  sims = sims + 1
}
p = count/sims
error = 3.29*sqrt(p*(1-p)/sims)
p
error
sims

### Infinite deck approximation of P(31-40 in k cards drawn from 6 deck shoe)

conv = function(f,g) convolve(f,rev(g),type="o")
p = rep(0,28)
shoe = c(rep(24/312,9), 96/312)
pmf = conv(shoe,shoe)
pmf = conv(pmf,shoe)
for (k in 4:28) {
  pmf = conv(pmf,shoe)
  p[k] = sum(pmf[(31-k+1):length(pmf)]) - sum(p[3:(k-1)])
  cat('\ncards =',k, '  prob =',p[k])
  flush.console()
}

### Exact P(31-40 in k cards drawn from 6 deck shoe) after Either

require(compiler)
enableJIT(3)

n = c(rep(24,9),96)   # number of each card
Q = rep(0,28)         # holds final probabilities
for (K in 4:28) {     # for each possible number of cards K
  for (k in 21:30) {  # for each sum of first K-1 cards k
    for (k10 in 0:3) {
      k.v = k10
      for (k9 in 0:max(0,min((k-10*k10)%/%9,K-1-k10))) {
        k.v = c(k9,k.v)
        for (k8 in 0:max(0,min((k-sum((9:10)*k.v))%/%8, K-1-sum(k.v)))) {
          k.v = c(k8,k.v)
          for (k7 in 0:max(0,min((k-sum((8:10)*k.v))%/%7, K-1-sum(k.v)))) {
            k.v = c(k7,k.v)
            for (k6 in 0:max(0,min((k-sum((7:10)*k.v))%/%6, K-1-sum(k.v)))) {
              k.v = c(k6,k.v)
              for (k5 in 0:max(0,min((k-sum((6:10)*k.v))%/%5, K-1-sum(k.v)))) {
                k.v = c(k5,k.v)
                for (k4 in 0:max(0,min((k-sum((5:10)*k.v))%/%4, K-1-sum(k.v)))) {
                  k.v = c(k4,k.v)
                  for (k3 in 0:max(0,min((k-sum((4:10)*k.v))%/%3, K-1-sum(k.v)))) {
                    k.v = c(k3,k.v)
                    for (k2 in 0:max(0,min((k-sum((3:10)*k.v))%/%2, K-1-sum(k.v)))) {
                      k.v = c(k2,k.v)
                      k1 = K-1-sum(k.v)
                      k.v = c(k1,k.v)
                      if ( sum(k.v) + 1 == K & sum(k.v*(1:10)) == k ) {
                        Q[K] = Q[K] + prod(choose(24,k.v[1:9]))*choose(96,k.v[10])*
                                       sum(n[(31-k):10]-k.v[(31-k):10]) / (choose(312,K)*K)
                      }
                      k.v = k.v[-(1:2)]    # Removes k1 and k2
                    } 
                    k.v = k.v[-1]          # Removes k3
                  } 
                  k.v = k.v[-1]            # Removes k4
                } 
                k.v = k.v[-1]              # Removes k5
              } 
              k.v = k.v[-1]                # Removes k6
            } 
            k.v = k.v[-1]                  # Removes k7
          } 
          k.v = k.v[-1]                    # Removes k8
        } 
        k.v = k.v[-1]                      # Removes k9
      } 
      k.v = k.v[-1]                        # Removes k10
    }
  }
  cat('\ncards =',K,'  prob =',Q[K])
  flush.console()
}