Fast way to calculate difference in normal CDFs I'm running a computationally intensive method where I have to calculate the difference in Normal CDF's millions of times, such as 
pnorm(y)-pnorm(x)

I have not looked into the details of how the CDF is calculated in R but I was thinking there may be a way to speed this up, perhaps into one pnorm call.  Or maybe there a mathematical simplification I can take advantage of.
Any suggestions (maybe it's not possible)?
 A: I think that a good 'one termer' approximation to $\Phi(x)$ is given by Polya's approximation:
$$(1)\;\;\;\Phi(x)\approx\frac{1}{2}\left(1+\text{sign}(x)\sqrt{1+\exp{\frac{-2x^2}{\pi}}}\right)$$
Below I place an R code:
aa<-seq(-5,5,l=1000)
x1<-pnorm(aa)
x2<-(1+sign(aa)*sqrt(1-exp(-2/pi*aa**2)))/2
plot(aa,x1)
lines(aa,x2,col="red")
legend("topleft",lty=c(1,1),lwd=c(1,1),col=c("black","red"),legend=c("pnorm","Polya's approximation"))

and a plot:

If you are willing to rewrite (1) in a compiled language it should be 
 faster than a call to pnorm since the latter is a 
much more complicated function. 
EDIT
Just for fun I implemented Polya's approximation in c++. I got 
the following (modest) speed up over R's pnorm:
> system.time(fit1<-.C("R_fastPolya",as.integer(n),as.single(x),as.single(x)))
   user  system elapsed 
  0.340   0.184   0.523 
> system.time(pnorm(x))
   user  system elapsed 
  0.856   0.024   0.879 

e.g. a mere 40% speed up :( --though, my implementation is pretty naive--. Finally, the maximum relative error is slightly less than $\Phi(x)/100$.
