#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define sz 2000000
#define bc 19
static double fr() {return random()/(double)~(1<<31);}
int hi[bc]; int out=0;
int main(){
for (int j=0; j<bc; ++j) hi[j]=0;
for (int j=0; j<sz; ++j) {
  double x = fr()*1.4-.7, y=fr()*1.4-.7, z=fr(), ss = x*x+y*y+z*z;
  if(ss*ss < z) ++hi[(int)(z/sqrt(ss)*bc)];
  else ++out;}
for(int j=0; j<bc; ++j) printf("%d %d %8.2f\n",
    j, hi[j], (sz-out)*(j+.5)/(bc*bc/2.));
printf ("%d out of SS.\n", out);
return 0;}
/* In the above we select vectors as described and then
record with a histogram the distribution of their cosines.
Routine fr provides a floating value uniformly between 0 and 1.
Array hi is the histogram.
Variable "out" counts the points in the box falling outside SS.
bc is the size of the histogram and sz is the sample size.
At the end we print out the histogram and also the expected
number of entries according to Lambert's formula.
0 2912  2961.06
1 8962  8883.17
2 14952 14805.28
3 20686 20727.39
4 26907 26649.50
5 32584 32571.61
6 38611 38493.72
7 44221 44415.83
8 50524 50337.94
9 56322 56260.05
10 62050 62182.16
11 68117 68104.27
12 73771 74026.39
13 79809 79948.50
14 85342 85870.61
15 92098 91792.72
16 97840 97714.83
17 103449 103636.94
18 109784 109559.05
931059 out of box.
*/