# Normal distribution values

Write the shortest program possible which outputs a table of Normal distribution values

phi(z) = Integral[t from -infty to z of e^(-t^2/2) /sqrt(2 pi) dt]


for z from 0.00 to 3.99, giving values accurate to 4 decimal places.

The canonical format for such a table appears to be:

    0.00    0.01    0.02    0.03    0.04    0.05    0.06    0.07    0.08    0.09
0.0 0.5000  0.5040  0.5080  0.5120  0.5160  0.5199  0.5239  0.5279  0.5319  0.5359
0.1 0.5398  0.5438  0.5478  0.5517  0.5557  0.5596  0.5636  0.5675  0.5714  0.5753
0.2 0.5793  0.5832  0.5871  0.5910  0.5948  0.5987  0.6026  0.6064  0.6103  0.6141
0.3 0.6179  0.6217  0.6255  0.6293  0.6331  0.6368  0.6406  0.6443  0.6480  0.6517
0.4 0.6554  0.6591  0.6628  0.6664  0.6700  0.6736  0.6772  0.6808  0.6844  0.6879
0.5 0.6915  0.6950  0.6985  0.7019  0.7054  0.7088  0.7123  0.7157  0.7190  0.7224
0.6 0.7257  0.7291  0.7324  0.7357  0.7389  0.7422  0.7454  0.7486  0.7517  0.7549
0.7 0.7580  0.7611  0.7642  0.7673  0.7704  0.7734  0.7764  0.7794  0.7823  0.7852
0.8 0.7881  0.7910  0.7939  0.7967  0.7995  0.8023  0.8051  0.8078  0.8106  0.8133
0.9 0.8159  0.8186  0.8212  0.8238  0.8264  0.8289  0.8315  0.8340  0.8365  0.8389
1.0 0.8413  0.8438  0.8461  0.8485  0.8508  0.8531  0.8554  0.8577  0.8599  0.8621
1.1 0.8643  0.8665  0.8686  0.8708  0.8729  0.8749  0.8770  0.8790  0.8810  0.8830
1.2 0.8849  0.8869  0.8888  0.8907  0.8925  0.8944  0.8962  0.8980  0.8997  0.9015
1.3 0.9032  0.9049  0.9066  0.9082  0.9099  0.9115  0.9131  0.9147  0.9162  0.9177
1.4 0.9192  0.9207  0.9222  0.9236  0.9251  0.9265  0.9279  0.9292  0.9306  0.9319
1.5 0.9332  0.9345  0.9357  0.9370  0.9382  0.9394  0.9406  0.9418  0.9429  0.9441
1.6 0.9452  0.9463  0.9474  0.9484  0.9495  0.9505  0.9515  0.9525  0.9535  0.9545
1.7 0.9554  0.9564  0.9573  0.9582  0.9591  0.9599  0.9608  0.9616  0.9625  0.9633
1.8 0.9641  0.9649  0.9656  0.9664  0.9671  0.9678  0.9686  0.9693  0.9699  0.9706
1.9 0.9713  0.9719  0.9726  0.9732  0.9738  0.9744  0.9750  0.9756  0.9761  0.9767
2.0 0.9772  0.9778  0.9783  0.9788  0.9793  0.9798  0.9803  0.9808  0.9812  0.9817
2.1 0.9821  0.9826  0.9830  0.9834  0.9838  0.9842  0.9846  0.9850  0.9854  0.9857
2.2 0.9861  0.9864  0.9868  0.9871  0.9875  0.9878  0.9881  0.9884  0.9887  0.9890
2.3 0.9893  0.9896  0.9898  0.9901  0.9904  0.9906  0.9909  0.9911  0.9913  0.9916
2.4 0.9918  0.9920  0.9922  0.9925  0.9927  0.9929  0.9931  0.9932  0.9934  0.9936
2.5 0.9938  0.9940  0.9941  0.9943  0.9945  0.9946  0.9948  0.9949  0.9951  0.9952
2.6 0.9953  0.9955  0.9956  0.9957  0.9959  0.9960  0.9961  0.9962  0.9963  0.9964
2.7 0.9965  0.9966  0.9967  0.9968  0.9969  0.9970  0.9971  0.9972  0.9973  0.9974
2.8 0.9974  0.9975  0.9976  0.9977  0.9977  0.9978  0.9979  0.9979  0.9980  0.9981
2.9 0.9981  0.9982  0.9982  0.9983  0.9984  0.9984  0.9985  0.9985  0.9986  0.9986
3.0 0.9987  0.9987  0.9987  0.9988  0.9988  0.9989  0.9989  0.9989  0.9990  0.9990
3.1 0.9990  0.9991  0.9991  0.9991  0.9992  0.9992  0.9992  0.9992  0.9993  0.9993
3.2 0.9993  0.9993  0.9994  0.9994  0.9994  0.9994  0.9994  0.9995  0.9995  0.9995
3.3 0.9995  0.9995  0.9995  0.9996  0.9996  0.9996  0.9996  0.9996  0.9996  0.9997
3.4 0.9997  0.9997  0.9997  0.9997  0.9997  0.9997  0.9997  0.9997  0.9997  0.9998
3.5 0.9998  0.9998  0.9998  0.9998  0.9998  0.9998  0.9998  0.9998  0.9998  0.9998
3.6 0.9998  0.9998  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999
3.7 0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999
3.8 0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999  0.9999
3.9 1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000

• do you want to exclude the inevitable boring answers in languages which have a built-in function to do this? – Peter Taylor Nov 24 '12 at 13:37
• Yes. No build-in functions. – student Nov 24 '12 at 13:38
• Also, do you want to place a limit on running time or are you willing to accept Monte-Carlo approaches (probably based on Box-Müller)? – Peter Taylor Nov 24 '12 at 13:39
• I haven't thought that but it would be nice that the program will halt in reasonable time. But it is hard to say what would be reasonable. – student Nov 24 '12 at 13:42
• @PeterTaylor Nice idea. As far as run time, MC solutions would only need to run on order of 10^9 total throws to achieve the expected precision. Maybe 10^11 to be confident. On a multi-GHz, superscaler machine that halts in a few minutes. Maybe better if you are able to take advantage of of the SIMD unit (thought is is not clear that you can). – dmckee --- ex-moderator kitten Nov 24 '12 at 17:28

# Mathematica 132 136 127

If Erfc is acceptable, here's my attempt, accurate to 6 places:

t=Table;TableForm[Partition[t[0.5 Erfc[-x/Sqrt[2]],{x,0,3.99,.01}],10],


• Someone needs to tell Wolfram what a significant figure is. 0.5 isn't the same as 0.500000. – Peter Taylor Nov 25 '12 at 15:21
• Interesting point. If one uses ScientificForm[number,4], the output will display, at most, four digits (times a power of 10), but it still truncates zeros to the right. So ScientificForm[0.5 ,3] returns 5. * 10^(-1). EngineeringForm[.5,3] shows similar inconsistencies. – DavidC Nov 25 '12 at 16:25
• (I know this answer is old but) If you want to be explicit about the precision, use 0.51 or 0.56. 0.5 and 0.500000 both have machine precision. – user202729 Nov 19 '18 at 5:15

## GolfScript (124 chars)

10:&,~]'    0.0'*40,{n\:I&/'.'I&%&,{'   '\I&*+&8?:^*100/:x^54,-2%{:i*x.{*^/}:$~$i).)*/^\-}/$3989423$15000500+5<(49-'.'@}%}/


Note: the only whitespace in there are two tab characters. MarkDown is turning them into triple-spaces.

The interesting part (i.e. the calculation) is done in x.8 fixpoint using identity 26.2.10 of Abramowitz and Stegun (1972). The program runs in under 2 seconds on my machine.

## Perl - 137 133 bytes

print"\t0.0$_"for 0..9;for($f=2.506628;$i++<80;$f*=-++$i){printf' %.1f',$x,s;;+\$x**$i/$i/$f;;printf"\t%.4f",eval,$x+=.01for(.5.$_)x10}


(both \t should be replaced with a literal tab character.)

Uses a Taylor Series expansion[1], adding one term per row, which seems to be good enough. The script produces the output exactly as expected, although each row does have a trailing whitespace. Runs in about 60ms on my lappy.

Edit as per comment by Peter Taylor.

[1] http://mathworld.wolfram.com/NormalDistributionFunction.html (formula 10)

• Hmm - MarkDown probably clobbered the tabs that were in the table when I pasted it. If you replace the spaces in the first format string by a single tab character and change the last format string into <TAB>%.4f then you get the correct output and save 4 chars. – Peter Taylor Nov 24 '12 at 17:15

## c -- 372

Just for interest this is an implementation of Peter's Monte Carlo suggestion from the comments

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define P printf
#define F(c,n) for(c=0;c<n;++c)
long c=9e9,i,j,k;double a[400],b=.01,t;
double B(){return sqrt(-2*log(drand48()))*cos(2*M_PI*drand48());}main()
{srand48(time(0));F(i,c){t=B();F(k,400)a[k] += t<b*k;}P("%4c",' ');F(j,10)
P("%4.2f    ",b*j);F(i,40){P("\n%3.1f ",10*i*b);F(j,10)
P("%6.4f  ",a[10*i+j]/c);}}


Lots of repetitive code. A very straight forward implementation, ungolfed it looks like

#include "stdio.h"
#include "stdlib.h"
#include "math.h"

long c=9e9, i, j, k;  /* demands longs longer than 32 bits */

double a[400], bw=.01, t;

double B/*oxMuller*/(){
/* normal distribution, centered on 0, std dev 1 */
return sqrt(-2*log(drand48())) * cos(2*M_PI*drand48());
}

main(){
srand48(time(0));
for(i=0; i<c; ++i) { /* compute the probabilities by MC integration */
t=B/*oxMuller*/();
for(k=0; k<400; ++k) a[k] += t<bw*k;
}
/* Print the table */
printf("%4c",' ');
for(j=0; j<10; ++j) printf("%4.2f    ",bw*j);
for(i=0; i<40; ++i) {
printf("\n%3.1f ",10*i*bw);
for(j=0; j<10; ++j) printf("%6.4f  ",a[10*i+j]/c);
}
printf("\n"); /* make the end pretty */
}


It is slow, about 60 seconds to do 100,000,000 throws with absolutely all the optimization bells and whistles turned on (-m64 -Ofast) on my 2.2 GHz Core 2 Duo laptop. Worse, it needs circa 50 times that many to provide real confidence that you have the right value in each and every bin. As written I'm asking for even more.

That seems too long to me, but I suspect that I underestimate how much work drand48 is doing (and I have sacrificed about a factor of 2 in the random number generation at the alter of small code).

It also demands that long is a bit bigger than 32 bits (thus the -m64 above).

• drand48() is not the problem; you're being hit by branch misprediction. I got a >3X performance improvement when I replaced for(k=0; k<400; ++k) a[k] += t<bw*k; with for( k=t<0?0:ceil(t/bw); k<400; ++k) a[k]+=1; – ceilingcat Apr 3 '18 at 4:28

# k

This is based on the Zelen and Severo algorithm here. 162 characters. Can improve the character count significantly by reducing precision.

+qp!(q;{abs(x>0)-(exp[-.5*x*x]%sqrt 2*acos -1)*t*.31938153+t*-.356563782+t*1.781477937+t*-1.821255978+1.330274429*t:1%1+.2316419*abs x}(q:0.1*!391)+\:0.01*!100)
`

Runs in about 15ms on a Pentium E2160.