Probability to roll a sum of less than 3600 with 1000 dice throws I try to estimate the probability that the sum of the numbers when rolling a die 1000 times is less than 3600. I am confident that I can approximate the probability with a normal distribution with expectation value $\mu=3500$ and the correct standard deviation. Is $\sigma=\sqrt{\frac{35n}{12}} (n=1000)$ the correct value? If I use this $\sigma$ then $p=96.8\%$ which looks pretty high to me. How can I verify my result?
 A: Pretty easy to simulate in python.
import numpy as np
num_rolls = 1000
num_trials = 100000
cutoff = 3600 

samples = np.random.choice(range(1,7),size=(num_trials, num_rolls)).sum(axis=1)

print((samples<cutoff).mean())

My sim agrees with your calculation.
A: You should have a continuity correction in your normal approximation.  For example in R, depending whether "less than $3600$" is intended to be strict or not:
maxrolls <- 1000
faces <- 6
pnorm(3600-0.5, (faces+1)*maxrolls/2, sqrt(maxrolls*(faces^2-1)/12))
# 0.9672904
pnorm(3600+0.5, (faces+1)*maxrolls/2, sqrt(maxrolls*(faces^2-1)/12)) 
# 0.9686207

and it is possible to do a more exact calculation of the probabilities
maxrolls <- 1000
faces <- 6
probs <- numeric(faces*maxrolls+1) # offset by 1
probs[1] <- 1
for (roll in 1:maxrolls){
  matprobs <- matrix(
    rep(c(rep(0,faces+1), probs), faces)[-(1:(faces))], 
    ncol=faces)
  probs <- (rowSums(matprobs) / faces)[1:(faces*maxrolls+1)]
}
sum(probs) # should be 1
# 1
sum(probs * (0:(faces*maxrolls))) - (faces+1)*maxrolls/2 # should be 0
# 1.818989e-12
sum(probs[1:3600]) / sum(probs) # prob strictly less than 3600
# 0.9672951
sum(probs[1:3601]) / sum(probs) # prob 3600 or less
# 0.9686258

and this shows that the normal approximation with a continuity correction is very close in this case, with both calculations within $0.000005$ of the correct values
A: The mean value $\mu$ for throwing dice $\mu=\frac16(1+2+...+6)=3.5$, and the standard deviation is
$$
\begin{align}
\sigma&=\sqrt{\frac16\{(6-\mu)^2+...+(1-\mu)^2\}}\\
&=\sqrt{35\over 12}
\end{align}
$$
so your numbers for the normal approximation for the distribution after $n$ trials seem to be OK.
Numerically,
package main

import (
    "fmt"
    "math/rand"
    "time"
)

func main() {
    nTrials := 100000
    nRolls := 1000
    over := 0
    rand.Seed(time.Now().UnixNano())
    for t := 0; t < nTrials; t++ {
        total := 0
        for r := 0; r < nRolls; r++ {
            total += rand.Intn(6) + 1
        }
        if total >= 3600 {
            over++
        }
    }
    fmt.Printf("%d %d %g\n", over, nTrials,
        float64(over)/float64(nTrials))
}

gives a value of
3246 100000 0.03246

which is the same as your answer
Run online here.
so it looks good.
