Take the 2-minute tour ×
Programming Puzzles & Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It's 100% free, no registration required.

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
share|improve this question
    
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 Nov 24 '12 at 17:28
show 2 more comments

5 Answers

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], 
TableHeadings->{t[k,{k,0,3.9,.1}],t[k,{k,0,.09,.01}]}]

enter image description here

share|improve this answer
    
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. –  David Carraher Nov 25 '12 at 16:25
add comment

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.

share|improve this answer
add comment

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)

share|improve this answer
    
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
add comment

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).

share|improve this answer
add comment

k

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

+`q`p!(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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.