/* CONVERTED FROM WIKI FORTRAN */
Normal_Random: procedure options (main);
   declare (array(1000), pi, temp,
            mean initial (1.0), sd initial (0.5)) float (18);
   declare (i, n) fixed binary;

   n = hbound(array, 1);
   pi = 4.0*ATAN(1.0);
   array = random(); /* Uniform distribution */
   /* Now convert to normal distribution */
   DO i = 1 to n-1 by 2;
      temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean;
      array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean;
      array(i) = temp;
   END;
   /* Check mean and standard deviation */
   mean = SUM(array)/n;
   sd = SQRT(SUM((array - mean)**2)/n);
   put skip edit ( "Mean = ", mean ) (a, F(18,16) );
   put skip edit ( "Standard Deviation = ", sd) (a, F(18,16));
END Normal_Random;