%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        This Postscript is Copyright (c) 2017, Peter J Billam          %
% Permission is granted  to any individual or institution to use, copy, %
% modify or redistribute this software, so long as it is not resold for %
% profit, and provided this notice is retained. It is provided "as is", %
% without any express or implied warranty.    http://www.pjb.com.au     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

/gauss_rand_r 0.5 1024.0 div 1024.0 div 1024.0 div def  % rand scaling factor
usertime realtime add srand

% 20210915
% https://en.wikipedia.org/wiki/Perlin_noise
% https://en.wikipedia.org/wiki/OpenSimplex_noise
% https://github.com/MarcoCiaramella/OpenSimplex2/tree/main/CPU
% https://github.com/MarcoCiaramella/OpenSimplex2/tree/main/GPU/OpenCL
%     https://github.com/KhronosGroup/OpenCL-Headers
%     https://en.wikipedia.org/wiki/OpenCL
% https://uniblock.tumblr.com/post/97868843242/noise
% https://github.com/KdotJPG/OpenSimplex2`   C# and Java

% The algorithm by  Richard P. Brent (1973), Computer Centre, ANU, Canberra
% Published (1981) as CACM Algorithm 448
% http://wwwmaths.anu.edu.au/~brent/pd/rpb023i.pdf
% was replaced in 20170510, fixing some obscure bugs ...
% http://www.design.caltech.edu/erik/Misc/Gaussian.html
% The polar form of the Box-Muller transformation is faster and more robust:
%   float x1, x2, w, y1, y2;
%   do {
%      x1 = 2.0 * ranf() - 1.0;
%      x2 = 2.0 * ranf() - 1.0;
%      w = x1 * x1 + x2 * x2;
%   } while ( w >= 1.0 );
%   w = sqrt( (-2.0 * log( w ) ) / w );
%   y1 = x1 * w;
%   y2 = x2 * w;
% where ranf() obtains a random number uniformly distributed in [0,1]
%
/gauss_rand_already false def
/gauss_rand {
  gauss_rand_already {
    /gauss_rand_already false def
    gauss_rand_y2 %  alternate between this and gauss_rand_y1  
  } {
    {
      /gauss_rand_x1 rand gauss_rand_r mul  2.0 mul 1.0 sub def
      /gauss_rand_x2 rand gauss_rand_r mul  2.0 mul 1.0 sub def
      /gauss_rand_w
        gauss_rand_x1 gauss_rand_x1 mul gauss_rand_x2 gauss_rand_x2 mul add def
      gauss_rand_w 1.0 le {exit} if
    } loop
    /gauss_rand_w  gauss_rand_w ln -2.0 mul  gauss_rand_w div  sqrt def
    /gauss_rand_y1 gauss_rand_x1 gauss_rand_w mul def
    /gauss_rand_y2 gauss_rand_x2 gauss_rand_w mul def
    /gauss_rand_already true def
    gauss_rand_y1
  } ifelse
} def

gauss_rand pop  % initialise, and discard the resulting zero

/grand {  % usage:  mean stddev grand
  gauss_rand mul add
} bind def

/irand {   % usage: n irand   ->  a random integer 0 .. n-1
  5 dict begin
  /n exch round cvi def
  n 0 eq { 0   % 20210428 defend against mod 0
  } { rand n mod
  } ifelse
  end
} bind def

/urand {   % usage: urand   ->  a random number between 0 and 1
  rand gauss_rand_r mul
} bind def

/urandom {   % random integer between 0 and 2^31 - 1
   5 dict begin
   /fn (/dev/urandom) def
   fn status {
      [ /created /referenced /bytes /pages ] { exch def } forall
      /rf fn (r) file def
      rf read pop
      256 mul rf read pop add
      256 mul rf read pop add
      128 mul rf read pop add   % remains on stack
      rf closefile
   } {
      usertime realtime add  % sad kludge
   } ifelse
} def

/randomget {   % usage:  [ an array ] randomget
  dup length  rand  exch mod  get
} bind def

/randomgetn {   % usage:  [ an array ] n randomgetn
  10 dict begin [ /n /arr_in ] { exch def } forall
  n arr_in length ge {
    arr_in
  } {
    /already n dict def [
    n {
      {
        /i  rand arr_in length mod def
        /candidate arr_in i get def
        already i known not { already i true put candidate exit } if
      } loop
    } repeat
  ] } ifelse
} bind def

/rayleigh_rand {   % usage:  sigma rayleigh_rand
  1.0 urand sub ln -2 mul sqrt mul
} bind def

/zipf_rand {   % usage:  [ an array ] zipf_rand
  20 dict begin /arr_in exch def
  /n arr_in length def
  /rel_freq n dict def
  /sum 0.0 def
  /i 0 def {
    i n ge { exit } if
    rel_freq  i  1 i 1 add div  put
    /sum  sum  rel_freq i get  add def
    /i i 1 add def
  } loop
  /switchpoints n dict def
  switchpoints 0  1 sum div   put  % first element is 1/sum
  /i 1 def {   % start loop at second element
    i  n 1 sub  ge { exit } if
    switchpoints i  switchpoints i 1 sub get  rel_freq i get sum div  add put
    /i i 1 add def
  } loop
  /r urand def
 /return_index n 1 sub def   % in case it's not found earlier
  /i 0 def {
    i  n 1 sub  ge { exit } if
    r  switchpoints i get  lt  { /return_index  i def  exit }  if
    /i i 1 add def
  } loop
  arr_in return_index get   % return
  end
} bind def