[Help-gsl] Reg :: My Gaussian Random Number Code

[Sumit Adhikari]

Hello All,
   Following is the code I wrote for generation of gaussian random number.
Can somebody tell me if there is any issue with this.


# ifndef GAUSSIAN_RANDOM_H
# define GAUSSIAN_RANDOM_H

# include <gsl/gsl_rng.h>
# include <gsl/gsl_randist.h>

template < typename P > class gaussian_random {
protected :
   const gsl_rng_type* T;
   gsl_rng* r ;

public :
   gaussian_random(){
      gsl_rng_env_setup();
      T = gsl_rng_default ;
      r = gsl_rng_alloc(T);
   }

   ~gaussian_random(){
      gsl_rng_free(r);
   }

   P read(){
      P sigma = 1.0 ;
      P rad ;
      P x ;
      P y ;

      do {
        x =  -1 + 2*gsl_rng_uniform(r) ;
        y =  -1 + 2*gsl_rng_uniform(r) ;
        rad = x*x + y*y ;
      } while(rad > 1.0 || rad == 0);

      P gauss_rand = sigma*y*sqrt(-2.0*log(rad)/rad) ;
      return (gauss_rand);

   }


};

# endif

Regards,
-- 
Sumit Adhikari,
Institute of Computer Technology,
Faculty of Electrical Engineering,
Vienna University of Technology,
Gußhausstraße 27-29,1040 Vienna