Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

[nt,kf,xa,xb,xc,xd]=airyzo(nt,kf,xa,xb,xc,xd);
```function [nt,kf,xa,xb,xc,xd]=airyzo(nt,kf,xa,xb,xc,xd);

%                Ai(x) and Ai'(x), a and a', and the associated
%                values of Ai(a') and Ai'(a); and the first NT
%                zeros of Airy functions Bi(x) and Bi'(x), b and
%                b', and the associated values of Bi(b') and
%                Bi'(b)
%       Input :  NT    --- Total number of zeros
%                KF    --- Function code
%                          KF=1 for Ai(x) and Ai'(x)
%                          KF=2 for Bi(x) and Bi'(x)
%       Output:  XA(m) --- a, the m-th zero of Ai(x) or
%                          b, the m-th zero of Bi(x)
%                XB(m) --- a', the m-th zero of Ai'(x) or
%                          b', the m-th zero of Bi'(x)
%                XC(m) --- Ai(a') or Bi(b')
%                XD(m) --- Ai'(a) or Bi'(b)
%                          ( m --- Serial number of zeros )
%       Routine called: AIRYB for computing Airy functions and
%                       their derivatives
%       =======================================================

pi=3.141592653589793d0;
for  i=1:nt;
if (kf == 1) ;
u=3.0.*pi.*(4.0.*i-1)./8.0d0;
u1=1./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((-15.5902.*u1+.929844).*u1-.138889).*u1+.10416667d0).*u1+1.0d0);
elseif (kf == 2);
if (i == 1) ;
rt0=-1.17371;
else;
u=3.0.*pi.*(4.0.*i-3.0)./8.0;
u1=1.0d0./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((-15.5902.*u1+.929844).*u1-.138889).*u1+.10416667).*u1+1.0);
end;
end;
rt=1.0e300;
while (abs((rt-rt0)./rt) > 1.d-9);
x=rt0;
if (kf == 1) rt=rt0-ai./ad; end;
if (kf == 2) rt=rt0-bi./bd; end;
if (abs((rt-rt0)./rt) > 1.d-9) ;
rt0=rt;
end;
end;
xa(i)=rt;
if (kf == 1) xd(i)=ad; end;
if (kf == 2) xd(i)=bd; end;
end;
for  i=1:nt;
if (kf == 1) ;
if (i == 1) ;
rt0=-1.01879;
else;
u=3.0.*pi.*(4.0.*i-3.0)./8.0;
u1=1./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((15.0168.*u1-.873954).*u1+.121528).*u1-.145833d0).*u1+1.0d0);
end;
elseif (kf == 2);
if (i == 1) ;
rt0=-2.29444;
else;
u=3.0.*pi.*(4.0.*i-1.0)./8.0;
u1=1.0./(u.*u);
rt0=-(u.*u).^(1.0./3.0).*((((15.0168.*u1-.873954).*u1+.121528).*u1-.145833).*u1+1.0);
end;
end;
rt=1.0e300;
while (abs((rt-rt0)./rt) > 1.0d-9);
x=rt0;
if (kf == 1) rt=rt0-ad./(ai.*x); end;
if (kf == 2) rt=rt0-bd./(bi.*x); end;
if (abs((rt-rt0)./rt) > 1.0d-9) ;
rt0=rt;
end;
end;
xb(i)=rt;
if (kf == 1) xc(i)=ai; end;
if (kf == 2) xc(i)=bi; end;
end;
return;

%       Input:   x  --- Argument of Airy function
%       Output:  AI --- Ai(x)
%                BI --- Bi(x)
%                BD --- Bi'(x)
%       =======================================================

ck=zeros(41,1);dk=zeros(41,1);

eps=1.0d-15;
pi=3.141592653589793d0;
c1=0.355028053887817d0;
c2=0.258819403792807d0;
sr3=1.732050807568877d0;
xa=abs(x);
xq=sqrt(xa);
if (x > 0.0d0) xm=5.0; end;
if (x <= 0.0d0) xm=8.0; end;
if (x == 0.0d0) ;
ai=c1;
bi=sr3.*c1;
bd=sr3.*c2;
return;
end;
if (xa <= xm) ;
fx=1.0d0;
r=1.0d0;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k-1.0d0).*x;
fx=fx+r;
if (abs(r./fx) < eps) break; end;
end;
gx=x;
r=x;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k+1.0d0).*x;
gx=gx+r;
if (abs(r./gx) < eps) break; end;
end;
ai=c1.*fx-c2.*gx;
bi=sr3.*(c1.*fx+c2.*gx);
df=.5d0.*x.*x;
r=df;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k+2.0d0).*x;
df=df+r;
if (abs(r./df) < eps) break; end;
end;
dg=1.0d0;
r=1.0d0;
for  k=1:40;
r=r.*x./(3.0d0.*k).*x./(3.0d0.*k-2.0d0).*x;
dg=dg+r;
if (abs(r./dg) < eps) break; end;
end;
bd=sr3.*(c1.*df+c2.*dg);
else;
xe=xa.*xq./1.5d0;
xr1=1.0d0./xe;
xar=1.0d0./xq;
xf=sqrt(xar);
rp=.5641895835477563d0;
r=1.0d0;
for  k=1:40;
r=r.*(6.0d0.*k-1.0d0)./216.0d0.*(6.0d0.*k-3.0d0)./k.*(6.0d0.*k-5.0d0)./(2.0d0.*k-1.0d0);
ck(k)=r;
dk(k)=-(6.0d0.*k+1.0d0)./(6.0d0.*k-1.0d0).*ck(k);
end;
km=fix(24.5-xa);
if (xa < 6.0) km=14; end;
if (xa > 15.0) km=10; end;
if (x > 0.0d0) ;
sai=1.0d0;
r=1.0d0;
for  k=1:km;
r=-r.*xr1;
sai=sai+ck(k).*r;
end;
sbi=1.0d0;
sbd=1.0d0;
r=1.0d0;
for  k=1:km;
r=r.*xr1;
sbi=sbi+ck(k).*r;
sbd=sbd+dk(k).*r;
end;
xp1=exp(-xe);
ai=.5d0.*rp.*xf.*xp1.*sai;
bi=rp.*xf./xp1.*sbi;
bd=rp./xf./xp1.*sbd;
else;
xcs=cos(xe+pi./4.0d0);
xss=sin(xe+pi./4.0d0);
ssa=1.0d0;
sda=1.0d0;
r=1.0d0;
xr2=1.0d0./(xe.*xe);
for  k=1:km;
r=-r.*xr2;
ssa=ssa+ck(2.*k).*r;
sda=sda+dk(2.*k).*r;
end;
ssb=ck(1).*xr1;
sdb=dk(1).*xr1;
r=xr1;
for  k=1:km;
r=-r.*xr2;
ssb=ssb+ck(2.*k+1).*r;
sdb=sdb+dk(2.*k+1).*r;
end;
ai=rp.*xf.*(xss.*ssa-xcs.*ssb);
bi=rp.*xf.*(xcs.*ssa+xss.*ssb);