normal_distribution.hpp

From: Levente Farkas (lfarkas_at_[hidden])
Date: 2001-10-01 12:47:30

hi,
I'm just look into the normal_distribution.hpp and see it's use sin and
cos for calculating the result (use the Box-Muller Method), altough it
can be "faster" with the Polar Method (see eg: 
http://www.dspguru.com/howto/tech/wgn2.htm)
the strange thing that the current version are faster even it use
trigonometric functions(?) at least on a PIII with MSVC 6sp4 when
I compile release (optimezed) code.
(I append the patch which I made for the test at the end of this mail).
I can't find the reason. but another intersting question if the 
current implementation use the Box-Muller Method while don't have 
a do loop:
do {
  _r2 = _rng();
} while (_r2 == 0);
and what's the reason of the 
_cached_rho = sqrt(-2 * log(1.0-_r2));
and not the 
_cached_rho = sqrt(-2 * log(r2));
?
these are just my 2c.

---------------------
--- normal_distribution.hpp     Mon Aug 20 16:42:10 2001
+++ normal_distribution.hpp.new Mon Oct  1 19:18:55 2001
@@ -51,21 +51,23 @@
   {
 #ifndef BOOST_NO_STDC_NAMESPACE
     // allow for Koenig lookup
-    using std::sqrt; using std::log; using std::sin; using std::cos;
+    using std::sqrt; using std::log;
 #endif
     if(!_valid) {
-      _r1 = _rng();
-      _r2 = _rng();
-      _cached_rho = sqrt(-2 * log(1.0-_r2));
+      double s, r1, r2;
+      do {
+        r1 = 2 * _rng() - 1;
+        r2 = 2 * _rng() - 1;
+        s = r1 * r1 + r2 * r2;
+      } while (s >= 1.0);
+      double temp = sqrt(-2.0 * log(s) / s) * _sigma;
+      _cached = r2 * temp + _mean;
       _valid = true;
+      return r1 * temp + _mean;
     } else {
       _valid = false;
+      return _cached;
     }
-    // Can we have a boost::mathconst please?
-    const double pi = 3.14159265358979323846;
-
-    return _cached_rho * (_valid ? cos(2*pi*_r1) : sin(2*pi*_r1)) * _sigma +
-      _mean;
   }
 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
   friend bool operator==(const normal_distribution& x,
@@ -85,7 +87,7 @@
 private:
   uniform_01<base_type, RealType> _rng;
   const result_type _mean, _sigma;
-  result_type _r1, _r2, _cached_rho;
+  result_type _cached;
   bool _valid;
 };
---------------------

 
 -- Levente                               "Si vis pacem para bellum!"