/* Box-Muller method       */
/* (c) 2025 cepstrum.co.jp */

clear;

SEED='a5a5a5a5';
//SEED='5a5a5a5a';
//SEED='0f0f0f0f';
//SEED='f0f0f0f0';

grand('setgen', 'mt');                    // use Mersenne-Twister 
grand('setsd',  hex2dec(SEED));           // set seed

DATLEN =2^20;                             // random data length (2^20=1048576)
HIST_DIV=100;                             // histragem parameter

x1=grand(1, DATLEN, 'def');               // range "[0, 1)"
x2=grand(1, DATLEN, 'def');               // range "[0, 1)"
x1=1.0-x1;                                // range "(0, 1]"
x2=1.0-x2;                                // range "(0, 1]"

y=sqrt(-2.0*log(x1)).*cos(2.0*%pi*x2);    // Box-Muller method

// calculate histgram
//ymax=1.05*max(abs(y));
ymax=6.0;
hist_range=2.0*ymax/(2*HIST_DIV+1);
hist=zeros(1:2*HIST_DIV+1);
for i=1:DATLEN
  ix=floor((y(i)+ymax)/hist_range)+1;
  hist(ix)=hist(ix)+1;
end

// plot histgram
scf(1);
//clf;           // clear screen
gcf().axes_size = [800, 600];
plot2d(((1:2*HIST_DIV+1)-HIST_DIV-1)*hist_range, hist, axesflag=1, style=2, rect=[-ymax, 0, ymax, 28000]);
plot2d([0, 0], [0, 28000],style=35);
title fontsize 2
title("cepstrum.co.jp");
xlabel "output histgram (Gaussian)"