Who created the first standard normal table? I'm about to introduce the standard normal table in my introductory statistics class, and that got me wondering: who created the first standard normal table? How did they do it before computers came along? I shudder to think of someone brute-force computing a thousand Riemann sums by hand.
 A: According to H.A. David [1] Laplace recognized the need for tables of the normal distribution "as early as 1783" and the first normal table was produced by Kramp in 1799. 
Laplace suggested two series approximations, one for the integral from $0$ to $x$ of $e^{-t^2}$ (which is proportional to a normal distribution with variance $\frac{_1}{^2}$) and one for the upper tail.
However, Kramp didn't use these series of Laplace, since there was a gap in the intervals for which they could be usefully applied.
In effect he starts with the integral for the tail area from 0 and then applies a Taylor expansion about the last calculated integral -- that is, as he calculates new values in the table he shifts the $x$ of his Taylor expansion of $G(x+h)$ (where $G$ is the integral giving the upper tail area).
To be specific, quoting the relevant couple of sentences:

he simply uses a Taylor series expansion of $G(x + h)$ about $G(x)$, with $h = .01$, up to the term in $h^3$. This enables him to proceed step by step from $x = 0$ to $x = h, 2h, 3h,...$, upon multiplying $he^{-x^2}$ by $$1-hx+ \frac13(2x^2 - 1)h^2 - \frac16(2x^3 - 3x)h^3.$$ Thus, at $x = 0$ this product reduces to $$.01 (1 - \frac13 \times .0001 ) = .00999967,\qquad\qquad (4)$$ so that at $G(.01) = .88622692 - .00999967 = .87622725$. The next term on the left of (4) can be shown to be $10^{-9}$, so that its omission is justified. 

David indicates that the tables were widely used.
So rather than thousands of Riemann sums it was hundreds of Taylor expansions.

On a smaller note, in a pinch (stuck with only a calculator and a few remembered values from the normal table) I have quite successfully applied Simpson's rule (and related rules for numerical integration) to get a good approximation at other values; it's not all that tedious to produce an abbreviated table* to a few figures of accuracy. [To produce tables of the scale and accuracy of Kramp's would be a fairly large task, though, even using a cleverer method, as he did.]
* By an abbreviated table, I mean one where you can basically get away with interpolation in between tabulated values without losing too much accuracy. If you only want say around 3 figure accuracy you really don't need to compute all that many values. I have effectively used polynomial interpolation (more precisely, applied finite difference techniques), which allows for a table with fewer values than linear interpolation -- if somewhat more effort at the interpolation step -- and also have done interpolation with a logit transformation, which makes linear interpolation considerably more effective, but is only much use if you have a good calculator). 
[1] Herbert A. David (2005),
"Tables Related to the Normal Distribution: A Short History"
The American Statistician, Vol. 59, No. 4 (Nov.), pp. 309-311 
[2] Kramp (1799),
Analyse des Réfractions Astronomiques et Terrestres,
Leipzig: Schwikkert  
A: Interesting issue! I think the first idea did not come through the integration of complex formula; rather, the result of applying the asymptotics in combinatorics. Pen and paper method may take several weeks; not so tough for Karl Gauss compared to calculation of pie for his predecessors. I think Gauss's idea was courageous; calculation was easy for him.
Example of creating standard z table from scratch-
1. Take a population of n (say n is 20) numbers and list all possible samples of size r (say r is 5) from that.
2. calculate the sample means. You get nCr sample means (here, 20c5=15504 means).
3. Their mean is the same as population mean. Find the stdev of sample means.
4. Find z scores of sample means using those pop mean and stdev of sample means.
5. Sort z in ascending order and find the probability of z being in a range in your nCr z values.
6. Compare values with normal tables. Smaller n is good for hand calculations. Larger n will produce closer approximates of the normal table values.  
The following code is in r:
n <- 20  
r <- 5  

p <- sample(1:40,n)  # Don't be misled!! Here, 'sample' is an r function  
                     used to produce n random numbers between 1 and 40.  
                     You can take any 20 numbers, possibly all different.  

c <- combn(p, r)     # all the nCr samples listed  
cmean <- array(0)  

for(i in 1:choose(n,r)) {  
    cmean[i] <- mean(c[,i])  
                }  

z <- array(0)  
for(i in 1:choose(n,r)) {  
    z[i] <- (cmean[i]-mean(c))/sd(cmean)  
                }  

ascend <- sort(z, decreasing = FALSE)  

Probability of z falling between 0 and positive value q below; compare with a known table. Manipulate q below between 0 and 3.5 to compare.  
q <- 1  
probability <- (length(ascend[ascend<q])-length(ascend[ascend<0]))/choose(n,r)   
probability   # For example, if you use n=30 and r=5, then for q=1, you  
              will get probability is 0.3413; for q=2, prob is 0.4773

