/* To see where a sample mean ranks among samples.
  Normal distribution */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define n 25

#define ss 100000000

typedef double R;

static R v(int j){return j<2 ? j+1 : 3.1415926535897933*v(j-2)*2/j;}

// Each invocation of rnd returns a new random normal deviate.
R rnd(){static int f = 0; R static const scale = 0x1.p-30;
  static R m;
  if(f) {f = 0; return m;}
  while (1){R x = scale*random()-1, y = scale*random()-1,
    s = x*x +y*y; 
    if(s<1) {R r = sqrt(- 2*log(s)/s);
    f = 1; m = x*r; return y*r;}}}

int main(){int j = 10000;
  if(RAND_MAX != 0x7fffffff) {printf("futz\n"); exit(0);}
  printf("%d samples for the %d-simplex\n", ss, n-1);
  if(0) while(j--) printf("%20.14f\n", rnd());
if(0){R s1=0, s2=0, s3=0;
  while (j--){R s = rnd(); s1 += s; s2 += s*s; s3 += s*s*s;}
  printf("%20.14f %20.14f %20.14f\n", s1, s2, s3);}
{int j = ss; int h[n+1];
  {int k=n+1; while(k--) h[k]=0;}
  while(j--){R a[n]; R s = 0;
    {int k=n; while(k--) {a[k] = rnd(); s += a[k];}}
    s /= n;
    {int k=n, m=0; while(k--) if(s<a[k]) ++m; ++h[m];}}
  if(h[0] || h[n]) {printf("Scrungee!!\n"); exit(9);}
  {int j; printf("j  h[j]\n");
    for(j=0; j<=n; ++j) printf("%d %d\n", j, h[j]);
    printf("\n");}
    {float A[n+1]; {int j=n+1; while(j--) A[j]=0;}
       {int j = n; while(j--){
        int C = 1, k=j+1;
        while(k--){
          printf("A[%d,%d] += %d*h[%d]\n", n-1,k,C,j);
          A[k] += (R)C*h[j];
          C = C*k/(j+1-k);}
        printf("\n");
        }}
      printf(" k  A[%d,k]\n", n-1);
      {int k, C=1;
      for(k=0; k<=n; ++k) {R a = A[k]/((R)C*ss);
          printf("%2d %8.5f, %8.5f %d\n", k, a, a*k*v(k), C);
          C = C*(n-k)/(k+1);}}}
return 0;}}