      subroutine rmastable(xr,a) 
c alexei bazavov feb 4 2003, bernd berg oct 27 2006.
c generate random numbers for symmetric stable distributions. pdf 
c is given by characteristic function \phi(t)=exp{-|t|^a} 0 < a =< 2.
c special cases: a=1 -- cauchy
c                a=2 -- gaussian (normal scaled by sqrt(2))
c based on j.m.chambers, c.l.mallows, b.w.stuck, 
c          j. am. stat. assoc. 71 (1976) 340.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(eps=two**-25,halfmeps=half-eps)
c phi is uniform on (-pi/2,pi/2). avoid +/-pi/2 in cos().
      call ranmar(xr1)
      phi=pi*(xr1-halfmeps)
c avoid zero and one in log (w zero or infinity):
      call ranmar(xr2)
      w=-log(xr2+eps) 
c final expression for stable variate:
c   \frac{\sin(a\phi)}{(\cos\phi)^{1/a}}\times
c   \left(\frac{\cos(1-a)\phi}{w}\right)^{1/a-1}.
      ai=one/a
      xr=(sin(a*phi)/(cos(phi))**(one/a))*
     &   (cos((one-a)*phi)/w)**(ai-one)
      return
      end
