deviates.h

   1 #ifndef DEVIATES_H
   2 #define DEVIATES_H
   3
   4 #include "ran.h"
   5
   6 struct Normaldev : Ran {
   7     Doub mu,sig;
   8     Normaldev(Doub mmu, Doub ssig, Ullong i)
   9     : Ran(i), mu(mmu), sig(ssig){}
  10     // Constructor arguments are mu, sigma, and a random sequence seed.
  11     Doub dev() {
  12     // Return a normal deviate.
  13         Doub u,v,x,y,q;
  14         do {
  15             u = doub();
  16             v = 1.7156*(doub()-0.5);
  17             x = u - 0.449871;
  18             y = abs(v) + 0.386595;
  19             q = SQR(x) + y*(0.19600*y-0.25472*x);
  20         } while (q > 0.27597
  21             && (q > 0.27846 || SQR(v) > -4.*log(u)*SQR(u)));
  22         return mu + sig*v/u;
  23     }
  24 };
  25
  26 #endif // DEVIATES_H