I confirmed the Wizard's numbers with an R simulation to within 0.01% with 50 million hands for 8 deck dealer stands on soft 17 (S17) and also for 1 deck S17. The 8 deck results are below, and the sim can generate the other results by changing num_decks and S17. Note that my entry for 21 includes BJs while the Wizard separates these, and my entry for 22 is for all busts 22-26. An exact calculation could be done with a recursive function.
Code:
sims = 50000000
num_decks = 8
S17 = TRUE
deck = c(rep(1:9,4*num_decks),rep(10,16*num_decks))
hit = ifelse(S17,6,7)
result = rep(0,26)
for (i in 1:sims) {
cards = sample(deck, 14, F)
j = 1
total = 0
ace = FALSE
while (total < 17) {
total = total + cards[j]
if (cards[j] == 1) {ace = TRUE}
if (ace && (total > hit) && (total < 12)) {total = total + 10}
j = j + 1
}
result[total] = result[total] + 1
}
result[22] = sum(result[22:26])
result[1:22]/sims
Output:
>result[1:22]/sims
[1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[8] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
[15] 0.0000000 0.0000000 0.1452270 0.1393228 0.1336369 0.1796657 0.1203770
[22] 0.2817706
Wizard's results (17, 18, 19, 20, 21 non-BJ, BJ, Bust):
0.145215 0.139318 0.133629 0.179709 0.072838 0.0474513 0.281839