Code covered by the BSD License  

Highlights from
Generation of Random Variates

image thumbnail

Generation of Random Variates

by

James Huntley (view profile)

 

generates random variates from over 870 univariate distributions

erfz(zz)
function f = erfz(zz)
%ERFZ  Error function for complex inputs
%   f = erfz(z) is the error function for the elements of z.
%   Z may be complex and of any size.
%   Accuracy is better than 12 significant digits.
%
%   Usage:  f = erfz(z)
%
%   Ref: Abramowitz & Stegun section 7.1
%   equations 7.1.9, 7.1.23, and 7.1.29
%
%   Tested under version 5.3.1
%
%   See also erf, erfc, erfcx, erfinc, erfcore

%   Main author Paul Godfrey <pgodfrey@intersil.com>
%   Small changes by Peter J. Acklam <jacklam@math.uio.no>
%   09-26-01

   error(nargchk(1, 1, nargin));
   
   % quick exit for empty input
   if isempty(zz)
      f = zz;
      return;
   end
   
   twopi = 2*pi;
   sqrtpi=sqrt(pi);

   f = zeros(size(zz));
   ff=f;

   az=abs(zz);
   p1=find(az<=8);
   p2=find(az> 8);

if ~isempty(p1)
   z=zz(p1);

   nn = 32;

   x = real(z);
   y = imag(z);
   k1 = 2 / pi * exp(-x.*x);
   k2 = exp(-i*2*x.*y);

   s1 = erf(x);

   s2 = zeros(size(x));
   k = x ~= 0;          % when x is non-zero
   s2(k) = k1(k) ./ (4*x(k)) .* (1 - k2(k));
   k = ~k;              % when x is zero
   s2(k) = i / pi * y(k);

   f = s1 + s2;

   k = y ~= 0;          % when y is non-zero
   xk = x(k);
   yk = y(k);

   s5 = 0;
   for n = 1 : nn
      s3 = exp(-n*n/4) ./ (n*n + 4*xk.*xk);
      s4 = 2*xk - k2(k).*(2*xk.*cosh(n*yk) - i*n*sinh(n*yk));
      s5 = s5 + s3.*s4;
   end
   s6 = k1(k) .* s5;
   f(k) = f(k) + s6;
   ff(p1)=f;
end

if ~isempty(p2)
   z=zz(p2);
   pn=find(real(z)<0);

   if ~isempty(pn)
      z(pn)=-z(pn);   
   end

   nmax=193;
   s=1;
   y=2*z.*z;
   for n=nmax:-2:1
       s=1-n.*(s./y);
   end

   f=1.0-s.*exp(-z.*z)./(sqrtpi*z);

   if ~isempty(pn)
      f(pn)=-f(pn);
   end

   pa=find(real(z)==0);
%  fix along i axis problem
   if ~isempty(pa)
      f(pa)=f(pa)-1;
   end

   ff(p2)=f;
end

f=ff;

return

   %a demo of this function is
   x = -4:0.125:4;
   y = x;
   [X, Y] = meshgrid(x,y);
   z = complex(X, Y);
   f = erfz(z);
   af = abs(f);
   %let's truncate for visibility
   p = find(af > 5);
   af(p) = 5;
   mesh(x, y, af);
   view(-70, 40);
   rotate3d on;

   return

Contact us