Re: circular random number distribution

[Jochen Küpper <jochen at fhi-berlin dot mpg dot de>]

Heiko Bauke <Heiko.Bauke@Physik.Uni-Magdeburg.DE> writes:

> | do {
> |     x = gsl_ran_flat(rng, -1, 1);
> |     y = gsl_ran_flat(rng, -1, 1);
> | } while(sqrt(x*x + y*y) > 1.);

Oops, of course! (Stupid me.)

Ok, I wrote a "real" benchmark, see attached code, generating 1e8
random pairs in a 1/1-circle using gsl_rng_taus and the different
methods discussed (see attached code for details).

Intel(R) Xeon(TM) CPU 2.80GHz machine with gcc-3.3.3:
,----
| reject: 15.59 s
| polar:  27.25 s
| dir_2d: 24.38 s
| trig2d: 27.68 s
`----
AMD Opteron(tm) Processor 850 2.4GHz with gcc-3.3.3:
,----
| reject: 4.48 s
| polar:  16.25 s
| dir_2d: 8.29 s
| trig2d: 16.87 s
`----

Greetings,
Jochen
-- 
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
    Liberté, Égalité, Fraternité                GnuPG key: CC1B0B4D
        (Part 3 you find in my messages before fall 2003.)

#include <cmath>
#include <ctime>
#include <iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

using namespace std;

const unsigned number = 100000000;


inline double time()
{
    return double(clock()) / CLOCKS_PER_SEC;
}


int main()
{
    gsl_rng_env_setup();
    const gsl_rng_type *rng_type = gsl_rng_taus;
    gsl_rng *rng = gsl_rng_alloc(rng_type);
    double x, y;
    // reject method
    double starttime = time();
    for(unsigned i=0; i<number; ++i) {
        do {
            x = gsl_ran_flat(rng, -1, 1);
            y = gsl_ran_flat(rng, -1, 1);
        } while((x*x + y*y) > 1.);
    }
    cout << "reject: " << time()-starttime << " s" << endl;
    // polar coordinates
    starttime = time();
    for(unsigned i=0; i<number; ++i) {
        double phi = gsl_ran_flat(rng, 0, 2.0*M_PI);
        double r = gsl_ran_flat(rng, 0, 1);
        x = sqrt(r)*cos(phi);
        y = sqrt(r)*sin(phi);
    }
    cout << "polar:  " << time()-starttime << " s" << endl;
    // Spherical Vector Distributions
    starttime = time();
    for(unsigned i=0; i<number; ++i) {
        gsl_ran_dir_2d(rng, &x, &y);
        double r = gsl_ran_flat(rng, 0, 1);
        x *= sqrt(r);
        y *= sqrt(r);
    }
    cout << "dir_2d: " << time()-starttime << " s" << endl;
    starttime = time();
    for(unsigned i=0; i<number; ++i) {
        gsl_ran_dir_2d_trig_method(rng, &x, &y);
        double r = gsl_ran_flat(rng, 0, 1);
        x *= sqrt(r);
        y *= sqrt(r);
    }
    cout << "trig2d: " << time()-starttime << " s" << endl;
}

// Local Variables:
// c-file-style: "Stroustrup"
// fill-column: 80
// End:

Attachment: pgp00000.pgp
Description: PGP signature