Update comparing_greedy_vs_analytical.cpp

experiments/random_unit_sampling_algorithms/comparing_greedy_vs_analytical.cpp

@@ -7,6 +7,7 @@
 #include <numeric>
 #include <random>
 #include <string>
+#include <cstring>
 #include <vector>
 using namespace std;
 using namespace chrono;
@@ -57,11 +58,34 @@ class RNG
   mt19937 _rng;
   uniform_real_distribution<double> _distributor;
 
+  static constexpr size_t DIST_CIRCLEF_ZIG_COUNT = 128; /* must be a power of two and <=2^14 */
+  static constexpr float DIST_CIRCLEF_ZIG_AREA = (0.7853981633974483/DIST_CIRCLEF_ZIG_COUNT); /* pi/4/n */
+  float zig[DIST_CIRCLEF_ZIG_COUNT*2 + 2];
+
 public:
   explicit RNG(const int seed) :
     _rng(seed),
     _distributor(0, 1)
-  { }
+  {
+    for (size_t i = 0; i < DIST_CIRCLEF_ZIG_COUNT; ++i) {
+      double A = DIST_CIRCLEF_ZIG_AREA * i;
+      double h = zig_solve(A);
+      zig[i*2+0] = cos(asin(h));
+      zig[i*2+1] = h;
+    }
+    zig[DIST_CIRCLEF_ZIG_COUNT*2+0] = 0.0;
+    zig[DIST_CIRCLEF_ZIG_COUNT*2+1] = 1.0;
+  }
+
+  static double zig_solve(double area)
+  {
+    double A, h = 0.5;
+    for (size_t i = 2; i < 64; ++i) {
+      A = h*cos(asin(h))/2 + asin(h)/2;
+      h += (A < area ? 1.0 : -1.0) / (1ull<<i);
+    }
+    return h;
+  }
 
   double num(const double a, const double b)
   {
@@ -91,6 +115,7 @@ class RNG
     }
   }
 
+
   /*== Analytical (mathematically beautiful) methods  */
 
   Vec2 analytical_in_unit_disk()
@@ -121,6 +146,35 @@ class RNG
 
     return Vec3(x, y, z);
   }
+
+  static float u32tof(uint32_t x)
+  {
+    float f;
+    memcpy(&f, &x, sizeof f);
+    return f;
+  }
+
+  Vec2 ziggurat_in_unit_disk()
+  {
+    while (true)
+    {
+      Vec2 f(num(-1, 1), num(-1, 1));
+      uint32_t rnd = _rng();
+      uint64_t idx = rnd & (DIST_CIRCLEF_ZIG_COUNT-1);
+
+      float vx = zig[idx*2+0];
+      float vy = zig[idx*2+1];
+      float ny = zig[idx*2+3];
+      f.x *= vx; /* scale to box */
+      f.y = f.y * (ny - vy) + vy;
+
+      if (f.length_squared() < 1) {
+        f.x = copysignf(f.x, u32tof(rnd));
+        f.y = copysignf(f.x, u32tof(rnd<<1));
+        return f;
+      }
+    }
+  }
 };
 
 
@@ -143,6 +197,7 @@ void run_benchmark(const int rng_seed, const int number_of_runs, const int num_p
 {
   vector<int64_t> rejection_2d_times;
   vector<int64_t> analytical_2d_times;
+  vector<int64_t> ziggurat_2d_times;
   vector<int64_t> rejection_3d_times;
   vector<int64_t> analytical_3d_times;
 
@@ -152,7 +207,7 @@ void run_benchmark(const int rng_seed, const int number_of_runs, const int num_p
   vec2_bucket.reserve(num_points_to_generate);
   vec3_bucket.reserve(num_points_to_generate);
 
-  cout << "run_number: rejection_2d_ms, analytical_2d_ms, rejection_3d_ms, analytical_3d_ms" << endl;
+  cout << "run_number: rejection_2d_ms, analytical_2d_ms, ziggurat_2d_ms, rejection_3d_ms, analytical_3d_ms" << endl;
 
   for (int i = 0; i < number_of_runs; i++)
   {
@@ -183,36 +238,49 @@ void run_benchmark(const int rng_seed, const int number_of_runs, const int num_p
     const auto a2d_end = system_clock::now();
     const auto a2d_time_ms = duration_cast<milliseconds>(a2d_end - a2d_start).count();
 
-    // Rejection 3d
+    // Ziggurat 2D
     RNG rng_3(seed);
+    const auto z2d_start = system_clock::now();
+    n = 0;
+    while (n < num_points_to_generate)
+    {
+      vec2_bucket[n++] = rng_3.ziggurat_in_unit_disk();
+    }
+    const auto z2d_end = system_clock::now();
+    const auto z2d_time_ms = duration_cast<milliseconds>(z2d_end - z2d_start).count();
+
+    // Rejection 3d
+    RNG rng_4(seed);
     const auto r3d_start = system_clock::now();
     n = 0;
     while (n < num_points_to_generate)
     {
-      vec3_bucket[n++] = rng_3.rejection_in_unit_sphere();
+      vec3_bucket[n++] = rng_4.rejection_in_unit_sphere();
     }
     const auto r3d_end = system_clock::now();
     const auto r3d_time_ms = duration_cast<milliseconds>(r3d_end - r3d_start).count();
 
     // Analytical 3d
-    RNG rng_4(seed);
+    RNG rng_5(seed);
     const auto a3d_start = system_clock::now();
     n = 0;
     while (n < num_points_to_generate)
     {
-      vec3_bucket[n++] = rng_4.analytical_in_unit_sphere();
+      vec3_bucket[n++] = rng_5.analytical_in_unit_sphere();
     }
     const auto a3d_end = system_clock::now();
     const auto a3d_time_ms = duration_cast<milliseconds>(a3d_end - a3d_start).count();
 
     cout << (i + 1) << ": "
          << r2d_time_ms << ", "
          << a2d_time_ms << ", "
+         << z2d_time_ms << ", "
          << r3d_time_ms << ", "
          << a3d_time_ms << ", " << endl;
 
     rejection_2d_times.push_back(r2d_time_ms);
     analytical_2d_times.push_back(a2d_time_ms);
+    ziggurat_2d_times.push_back(z2d_time_ms);
     rejection_3d_times.push_back(r3d_time_ms);
     analytical_3d_times.push_back(a3d_time_ms);
   }
@@ -224,11 +292,13 @@ void run_benchmark(const int rng_seed, const int number_of_runs, const int num_p
   cout << "mean: "
        << compute_mean(rejection_2d_times) << ", "
        << compute_mean(analytical_2d_times) << ", "
+       << compute_mean(ziggurat_2d_times) << ", "
        << compute_mean(rejection_3d_times) << ", "
        << compute_mean(analytical_3d_times) << endl;
   cout << "median: "
        << compute_bad_median(rejection_2d_times) << ", "
        << compute_bad_median(analytical_2d_times) << ", "
+       << compute_bad_median(ziggurat_2d_times) << ", "
        << compute_bad_median(rejection_3d_times) << ", "
        << compute_bad_median(analytical_3d_times) << endl;
 }
@@ -249,5 +319,3 @@ int main(int argc, char *argv[])
 
   return 0;
 }
-
-