#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#define BINS 200
#define TRIALS 10000

double distribution(double arg); /* individual data distribution */

int n,m,k,l;
double x,y,z, largest, bin_x, occ;
int b_value;
double numbers[1000];            /* store 1000 randoms */
double *num_ptr;

main()
{
double averages[TRIALS];         /* store averages of the 1000 randoms */
double bins[BINS+1];               /* sorted pidgeon-holed averages */
double *bin_ptr;
double DIV;

srand((unsigned long)getpid());

  for(m=0;m<TRIALS;m++){

  z=0.0;                         /* initialize the average */  

  for(n=0;n<1000;n++){
  x=(double)(rand())/(double)RAND_MAX;
  numbers[n]=distribution(x);
  z=z+numbers[n];
  }

averages[m]=z/1000.0;            /* fill array of averages */
  }

/* set up the bins */
for(k=0;k<BINS;k++){
  bins[k]=0;}
bin_ptr=&bins[0];

/* find the largest value stored in the averages array */
largest=averages[0];
for(k=0;k<TRIALS;k++){
if(averages[k]> largest)
largest=averages[k];
}

/* set up bin size */

DIV=largest/(double)BINS;

/* now create a frequency plot of the averages */

  for (k = 0; k < TRIALS; k++)
    {
      b_value=floor(averages[k]/DIV);
      *(bin_ptr+b_value)=*(bin_ptr+b_value)+1;
    }
  for (l= 0; l <=BINS; l++)
    {
      bin_x=(float)l*(float)DIV;
      occ=(*bin_ptr)*(float)BINS/(float)TRIALS;
      printf("%f\t%f\n", bin_x, occ);
      bin_ptr++;
    }                        /*    end of data processing section */





}                             /* the end */

double distribution(double arg)
  {
    /* example; the density function for sqrt distribution */
    /* make sure that this does NOT generate large numbers */
   double y;
   y=arg*arg;
   return(y);
  }