%
% gauss.m
% (c) Neil Gershenfeld  3/7/01
% generate Gaussian deviates and test their cumulants
%
n = 1;
npts = 100000;
moment1 = 0;
moment2 = 0;
moment3 = 0;
for i = 1:npts
   n = rem(8121*n+28411, 134456);
   x1 = (1+n)/134456;
   n = rem(8121*n+28411, 134456);
   x2 = (1+n)/134456;
   y1 = sqrt(-2*log(x1))*sin(2*pi*x2);
   y2 = sqrt(-2*log(x1))*cos(2*pi*x2);
   moment1 = moment1 + y1 + y2;
   moment2 = moment2 + y1^2 + y2^2;
   moment3 = moment3 + y1^3 + y1^3;
   end
moment1 = moment1/(2*npts);
moment2 = moment2/(2*npts);
moment3 = moment3/(2*npts);
cum1 = moment1;
cum2 = moment2 - moment1^2;
cum3 = moment3 - 3*moment1*moment2 + 2*moment1^3;
fprintf('Cumulant 1 = %f\n',cum1)
fprintf('Cumulant 2 = %f\n',cum2)
fprintf('Cumulant 3 = %f\n',cum3)