plait/gauss.awk

# generate a normally distributed random number

function boxmuller(x1, x2, w)
{
  if (gaussindex==1)
    {
      gaussindex = 0;
      return gauss1;
    }
  else
    {
      do
	{
	  x1 = 2.0 * rand() - 1.0;
	  x2 = 2.0 * rand() - 1.0;
	  w = x1 * x1 + x2 * x2;
	} while (w >= 1.0);
      
      w = sqrt((-2.0 * log(w)) / w );
      gauss1 = x1 * w;
      gaussindex = 1;
      return (x2 * w);
    }
}

function boxmuller_old(x1, x2, w, y1, y2)
{
  do 
    {
      x1 = 2.0 * rand() - 1.0;
      x2 = 2.0 * rand() - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
  
  w = sqrt( (-2.0 * log( w ) ) / w );
  y1 = x1 * w;
  y2 = x2 * w;
  
  return y1;
}

BEGIN {
  srand();
  max = 0;
  min = 100;
  total = 0;
  for (j=-20; j<20; j++)
    {
      freq[j] = 0;
    }
  for (i=1; i<1000; i++)
    {
      foo = boxmuller();
      if (foo > max) max = foo;
      if (foo < min) min = foo;
      total += foo;
      
      freq[int((foo+4)*7)]++
#      for (j=-25; j<25; j++)
#	{
#	  if (foo*7 > (j-1) && foo*7 < j)
#	    {
#	      freq[j]++;
#	    }
#	}
    }      
  print "min=" min;
  print "max=" max;
  print "avg=" total / 1000;
      
  for (j=1; j<55; j++)
    {
      printf "%3d ", j;
      for (k=1; k<freq[j]; k++)
	{
	  printf "x";
	}
      printf "\n";
    }
      
}