Probability of an open-ended 6 sided die roll?

Not sure how you mean that "for each die that rolls a six, roll another" line. Does that mean that for an N of 10, you might actually roll 11 or 12?

The probability of getting a success (P) is (7-M)/6.
The mean number of successes is PN.
The standard deviation of the number of successes over numerous trials is sqrt( P*(1-P)/N )
To find the probability of S successes, use the mean and standard deviation above and look up S on a normal distribution probability chart.

If you want to say that 10 rolls could actually mean 11 or 12 rolls, use N'=7N/6 in the equations.

So, if sixes continue to be rolled, you keep rolling dice (possibly, but improbably, forever)?

Yes. If any of the dice roll a 6, add one to N. Yes, this has the infinitely improbable chance of going on forever.

However, I need a formula to solve this problem, expressed as a summation. The platform on which this will be implemented can't contain a lookup table.

Yeah, exactly. To supplement my above post, let's say a person rolls six dice, and gets the results 1-2-3-5-6-6.

The person would then roll two more dice, perhaps rolling 4-6. Since the second die rolled a six, they would roll one more die, perhaps rolling 3.

Thus, the total rolled dice would be 1-2-3-5-6-6-4-6-3 and this would contain 5 successes.

Since you mentioned a platform, I'm guessing that you need to program this. I think this recursive function solves the problem. However, as with a recursive implementation of the Fibonacci sequence, it's not very efficient for large s.

If you are going to roll N dice without the six rule, then with the six rule you will on average roll 6N/5 dice. I get that from the infinite series:

N + N/6 + N/36 + N/216... = N / (1 - 1/6) = 6N/5

Compute the probability using a standard value for N then substitute 6N/5 in for N and you should get the right average value.

The mean number of successes, Sbar = (6N/5)[(7-M)/6] = N (7-M) / 5
Expected standard deviation = sqrt [ P*(1-P)/ (6N/5) ] where P = (7-M)/6
Z score is (S - Sbar)/standard deviation

There is not an exact computation of probability from z score but you can go here to see some approximations.

The probability distribution is not normal in this case (it's not even symmetric about the mean). Z scores don't apply.

I did a trial with 10,000 iterations of N=10 and M=4. The average score was 6 as expected. The distribution was somewhat normal but skewed to the low values. The sigma was a lot more than expected (2.2 vs 1.7). Here is the histogram.

Score Count
0 14
1 93
2 294
3 765
4 1325
5 1801
6 1855
7 1491
8 1107
9 641
10 341
11 164
12 63
13 24
14 13
15 6
16 2
17 1

You may have to run a large number of trials and compute the probability from that. My program differs slightly from Mark's but I'm about to run 100,000 iterations very quickly if I'm not trying to record results. It would be pretty easy to enter S then keep track of over/under.

My program is not a simulation. If I haven't screwed up, it computes exact probabilities. Here's what I get for engineer's test (computed probabilities times 10000):

0 9.765625000000000
1 81.380208333333329
2 318.739149305555429
3 782.154224537036953
4 1353.887864101080140
5 1767.068242830503777
6 1818.271003796403193
7 1525.592847561584904
8 1073.234929051491235
9 648.584100814064755
10 344.022565672190808
11 163.211453308003172
12 70.389978852555444
13 27.978765318732467
14 10.367693444363100
15 3.615824225747254
16 1.196294374404082
17 0.377949323590987
18 0.114653000031351
19 0.033550997081045
20 0.009508158506741