🚀 go-pugleaf

Here's the summary of responses that I received concerning the
generation of Cauchy deviates:

Pretty much everyone who replied pointed out that the Cauchy distribution
can be generated by generating two standard normal deviates and taking
their quotient:  N(0,1)/N(0,1).  Barrett P. Eynon <barry@playfair.stanford.edu>
noted that the t distribution with 1 degree of freedom is equivalent to
the Cauchy distribution.  Finally, B. Narasimhan <naras@cda.mrs.umn.edu>
noted that 

	the ratio X/Y is Cauchy any time the pair (X,Y) is radially
	distributed. See Devroye, "Non Uniform Random Variate
	Generation", for more. Thus you can even use ratio of two
	coordinates of a point generated uniformly in the unit disk in
	2D.

Code which uses this last observation was provided by Lester Ingber
<ingber@umiacs.UMD.EDU> and is included at the end of this article.

Mark Plutowski <pluto@cs.ucsd.edu> noted that the Cauchy distribution
is also called the Lorentzian distribution, and is discussed as such in
Numerical Recipies.  Both NR and Ping Chou <choup@cgrb.orst.edu> noted
that you can generate a Cauchy deviate by taking a uniform deviate
on x \in [0,1] and computing tan(PI*x).  However, in NR they later
compute the Cauchy deviates using the same method used in Ingber's code,
so presumably this is more efficient.

Robert E George <regeorge@magnus.acs.ohio-state.edu> also noted that
there is a routine in IMSL for generating Cauchy deviates.

Thanks for everyone's replies!

--Bill


/** cauchy
 **
  #                      Prof. Lester Ingber                      #
  #                   ingber@alumni.caltech.edu                   #
 ** 
 **/

#define EPS     1.e-12

/* myrand() is any uniform number generator {0,1}. */

double cauchy(void)
{
        double x, y, z;
 
        for(;;) {
                x = TWO * myrand() - ONE;
                y = TWO * myrand() - ONE;
                z = x * x + y * y;
                if (z <= ONE)
                        break;
        }
        if (fabs(x) < EPS)
                x = (x < ZERO ? x - EPS : x + EPS);
 
        return (y / x);
}