Normalny rozkład zmiennych losowych - D

Mam taki kod klasy:

class NormalDistribution : IDistribution {
public:
  	this(real mean_arg = 0, real sigma_arg = 1)
  	out {
		assert(_sigma >= 0);
  	}
  	body {
  		_mean = mean_arg;
  		_sigma = sigma_arg;
  	}

  	// compiler-generated copy ctor and assignment operator are fine

  	real mean() { return _mean; }
  	real sigma() { return _sigma; }

 	void reset() { _valid = false; }

  	real opCall(IEngine eng) 
  	in {
  		assert(eng !is null);
  	}
  	out(result) {
  		assert(result != real.nan);
  	}
  	body {
    	if(!_valid) {
      		_r1 = eng.next;
      		_r2 = eng.next;
      		assert(_r1 != 0 || _r2 != 0);
      		_cached_rho = sqrt((-2.0)*log(1.0-_r2));
      		assert(_cached_rho != real.nan);
      		_valid = true;
    	} else {
    		_valid = false;
    	}    
    	return _cached_rho*(_valid?cos(2*PI*_r1):sin(2*PI*_r1))*_sigma+_mean;
  	}
private:
  	real _mean = 0, _sigma = 1;
  	real _r1 = 0, _r2 = 0, _cached_rho = 0;
  	bool _valid = false;
}

I ma on zwracać wartość zmiennej w rozkładzie normalnym, ale ciągle mi wywala NaN/0. Dane wyświetlałem printfem(wywala 0) z stdc i Stdout(wywala "nan") z Tango.

Hej, z tego co pamiętam to najprostszy sposób na rozkład normalny to zsumować 12 liczb losowych o rozkładzie jednorodnym i to jest już dość dobrym przybliżeniem (nie pamiętam wzoru normalizacyjnego by otrzymać właściwą szerokość rozkładu).

Nie znam składni tego języka, ale ja bym to rozbił wzór na proste kawałki a następnie potraktowałbym to debuggerem.
Wydaje mi się, że powinieneś stosować jedną funkcję trygonometryczną (sinus ALBO cosinus) przynajmniej tak sugeruje algorytm Boxa-Mullera.
Mam wrażenie, że przekombinowałeś.

W c++ to wygląda tak:

double normalDistribution(double mean, double sigma) {
      double x1 = (rand()+1)/(double)RAND_MAX; // przedział (0,1>
      double x2 = rand()/(double)RAND_MAX;//  // przedział <0,1)

      double r = sqrt(-2.0 * ln(x1) );
      double ro = TwoPi * x2;

      return r*sin(ro)*sigma + mean;
      // return r*cos(ro)*sigma + mean; // druga zmienna losowa do zignorowania
}