function R = mrandn(varargin);

% Implements the Marsaglia ziggurat algorithm for normal random number generation.
% This function exactly reproduces the numerical behavior 
% of the built-in (C-coded) MATLAB function RANDN (since version 5).
%
% Usage:
%
%    mrandn                      1-by-1
%    mrandn(n)                  n-by-n
%    mrandn(m,n)              m-by-n
%    s = mrandn('state')     Get state
%    mrandn('state',j)         Set state determined by j
%    mrandn('state',s)        Restore state
%
% See RANDN documentation for details.
%  See also POLARALG.

% Reference: "A fast, easily implemented method for sampling
% from decreasing or symmetric unimodal density functions",
% George Marsaglia and W. T. Tsang, SIAM J. Sci. Stat. Comput.,
% Vol. 5, No. 2, June 1984, pp. 349-359.

persistent z
if isempty(z)

   % Initialize the ziggurat and the state.

   z = randnsetup;
   randuni(0);
end

if nargin > 0 && isequal(varargin{1},'state')

   % Get or set the state.

   if nargin == 1
      R = randuni('state');
   else
      randuni(varargin{2});
   end
   return
end

% Determine the size of the output.

if nargin == 0
   m = 1; n = 1;
elseif nargin == 1
   m = varargin{1}; n = varargin{1};
else
   m = varargin{1}; n = varargin{2};
end

% Ziggurat normal random number generator.

R = zeros(m,n);
for k = 1:m*n
   [u,j] = randuni;
   rk = u*z(j+1);
   if abs(rk) < z(j)
      R(k) = rk;
   else
      R(k) = randntips(rk,j,z);
   end
end


% ----------------------------------------------------

function [u,j] = randuni(state)
% RANDUNI  Uniform random number generator used by RANDN.
% Should be initialized with
%    randuni(0) 
% Then
%    [u,j] = randuni
% generates a single random real number, -1 < u < 1,
% and an optional random integer, 1 <= j <= 64.
% The internal state is a vector s of two integers that
% can be retrieved or restored with
%    s = randuni('state')
%    randuni(s)
%    randuni(j)

   persistent icng jsr
   if nargin > 0

      % Get or set the state

      if isstr(state)
         u = [icng; jsr];
      elseif length(state) == 1
         icng = 362436069;
         jsr = state;
         if jsr == 0, jsr = 521288629; end
      else
         icng = state(1);
         jsr = state(2);
      end
      return
   end

   % The algorithm combines multiplicative congruential and
   % shift register generators.

   icng = mod(69069*icng + 1234567,2^32);
   jsr = bitxor(jsr,bitshift(jsr,13,32));
   jsr = bitxor(jsr,bitshift(jsr,-17,32));
   jsr = bitxor(jsr,bitshift(jsr,5,32));
   m = mod(icng+jsr,2^32);
   if m >= 2^31, m = m - 2^32; end
   u = pow2(m,-31);
   j = mod(m,64)+1;


% ----------------------------------------------------

function z = randnsetup
% RANDNSETUP(n)  Generate ziggurat used by RANDN.

   n = 64;
   c = sqrt(pi/2)/n;
   F = inline('x.*exp(-x.^2/2)-q','x','q');
   xn = mfzero(F, [1,5], 3*c);
   z(n-2:n+1,1) = xn;
   for k = n-3:-1:1
      z(k) = sqrt(-2*log(exp(-z(k+1).^2/2) + c/z(k+1)));
   end


% ----------------------------------------------------

function z = randntips(r,j,z)
% RANDNTIPS  Detailed calculations in the zigguart tips.

   persistent a b c c1 c2 pc xn
   if isempty(xn)

      % One time parameter computation

      f = inline('exp(-x.^2/2)');
      n = length(z)-1;
      x0 = z(1);
      b = sqrt(sqrt(pi/2)*x0*(1-sum(z(1:n)./z(2:n+1))/n)/(1-f(x0)));
      a = x0/(b*(1-f(x0)));
      c = 1+a*f(-x0);
      c1 = 0.5+n*z(n-2)*(f(0.5*(z(n-2)+z(n-3)))-f(z(n-2)))/sqrt(pi/2);
      c2 = 2-x0/b;
      pc = sqrt(pi/2)/n;
      xn = z(n+1);
   end

   x = (abs(r)-z(j))/(z(j+1)-z(j));
   y = (1 + randuni)/2;
   s = x + y;
   if s > c2
      z = sign(r)*(b-b*x);
      return
   end

   if s <= c1
      z = r;
      return
   end

   x = b - b*x;
   if y > c-a*exp(-x^2/2)
      z = sign(r)*x;
      return
   end

   if exp(-z(j+1)^2/2)+y*pc/z(j+1) <= exp(-r^2/2)
      z = r;
      return
   end

   while 1
      x = log((1 + randuni)/2)/xn;
      x2 = -2.*log((1 + randuni)/2);
      if x2 > x^2
         z = sign(r)*(xn - x);
         return
      end
   end
   
   
% ----------------------------------------------------

function b = mfzero(F,ab,varargin);

%   MFZERO  M-Code version of FZERO.
%   x = mfzero(F,[a,b]) tries to find a zero of F(x) between a and b.
%   F(a) and F(b) must have opposite signs.  mfzero returns one 
%   end point of a small subinterval of [a,b] where F changes sign.
%   Arguments beyond the first two, mfzero(F,[a,b],p1,p2,...),
%   are passed on, F(x,p1,p2,..).
%
%   Example:
%      mfzero('sin(x)',[1,4])

% Make F callable by feval.
if ischar(F) & exist(F)~=2
   F = inline(F);
elseif isa(F,'sym')
   F = inline(char(F));
end 

% Initialize.
a = ab(1);
b = ab(2);
fa = feval(F,a,varargin{:});
fb = feval(F,b,varargin{:});
if sign(fa) == sign(fb)
   error('Function must change sign on the interval')
end
c = a;
fc = fa;
d = b - c;
e = d;

% Main loop, exit from middle of the loop
while fb ~= 0

   % The three current points, a, b and c, satisfy:
   %    f(x) changes sign between a and b.
   %    abs(f(b)) <= abs(f(a)).
   %    c = previous b, so c might = a.
   % The next point is chosen from:
   %    Bisection point, (a+b)/2.
   %    Secant point determined by b and c.
   %    Inverse quadratic interpolation point determined
   %    by a, b and c if they are distinct.

   if sign(fa) == sign(fb)
      a = c;  fa = fc;
      d = b - c;  e = d;
   end
   if abs(fa) < abs(fb)
      c = b;    b = a;    a = c;
      fc = fb;  fb = fa;  fa = fc;
   end
   
   % Convergence test and possible exit
   m = 0.5*(a - b);
   tol = 2.0*eps*max(abs(b),1.0);
   if (abs(m) <= tol) | (fb == 0.0)
      break
   end
   
   % Choose bisection or interpolation
   if (abs(e) < tol) | (abs(fc) <= abs(fb))
      % Bisection
      d = m;
      e = m;
   else
      % Interpolation
      s = fb/fc;
      if (a == c)
         % Linear interpolation (secant)
         p = 2.0*m*s;
         q = 1.0 - s;
      else
         % Inverse quadratic interpolation
         q = fc/fa;
         r = fb/fa;
         p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0));
         q = (q - 1.0)*(r - 1.0)*(s - 1.0);
      end;
      if p > 0, q = -q; else p = -p; end;
      % Is interpolated point acceptable
      if (2.0*p < 3.0*m*q - abs(tol*q)) & (p < abs(0.5*e*q))
         e = d;
         d = p/q;
      else
         d = m;
         e = m;
      end;
   end
   
   % Next point
   c = b;
   fc = fb;
   if abs(d) > tol
      b = b + d;
   else
      b = b - sign(b-a)*tol;
   end
   fb = feval(F,b,varargin{:});
end