update rhalton sampling to double precision

src/random_ray/random_ray.cpp

@@ -180,12 +180,12 @@ float exponentialG2(float tau)
 // Implementation of the Fisher-Yates shuffle algorithm.
 // Algorithm adapted from:
 //    https://en.cppreference.com/w/cpp/algorithm/random_shuffle#Version_3
-void fisher_yates_shuffle(vector<int>& arr, uint64_t* seed)
+void fisher_yates_shuffle(vector<int64_t>& arr, uint64_t* seed)
 {
   // Loop over the array from the last element down to the second
-  for (size_t i = arr.size() - 1; i > 0; --i) {
+  for (int i = arr.size() - 1; i > 0; --i) {
     // Generate a random index in the range [0, i]
-    size_t j = uniform_int_distribution(0, i, seed);
+    int j = uniform_int_distribution(0, i, seed);
     // Swap arr[i] with arr[j]
     std::swap(arr[i], arr[j]);
   }
@@ -196,27 +196,27 @@ void fisher_yates_shuffle(vector<int>& arr, uint64_t* seed)
 // Algorithm adapted from:
 //      A. B. Owen. A randomized halton algorithm in r. Arxiv, 6 2017.
 //      URL https://arxiv.org/abs/1706.02808
-vector<vector<float>> rhalton(int N, int dim, uint64_t* seed, int64_t skip = 0)
+vector<vector<double>> rhalton(int64_t N, int dim, uint64_t* seed, int64_t skip = 0)
 {
-  int b;
+  int64_t b;
   double b2r;
   vector<double> ans(N);
-  vector<int> ind(N);
-  vector<int> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
-  vector<vector<float>> halton(N, vector<float>(dim, 0.0));
+  vector<int64_t> ind(N);
+  vector<int64_t> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
+  vector<vector<double>> halton(N, vector<double>(dim, 0.0));
 
   std::iota(ind.begin(), ind.end(), skip);
 
   for (int D = 0; D < dim; ++D) {
     b = primes[D];
     b2r = 1.0 / b;
-    vector<int> res(ind);
+    vector<int64_t> res(ind);
     std::fill(ans.begin(), ans.end(), 0.0);
 
     while ((1.0 - b2r) < 1.0) {
-      vector<int> dig(N);
+      vector<int64_t> dig(N);
       // randomaly permute a sequence from skip to skip+N
-      vector<int> perm(b);
+      vector<int64_t> perm(b);
       std::iota(perm.begin(), perm.end(), 0);
       fisher_yates_shuffle(perm, seed);
 
@@ -638,7 +638,7 @@ SourceSite RandomRay::sample_halton()
   stream() = STREAM_TRACKING;
 
   // Calculate next samples in LDS
-  vector<vector<float>> samples = rhalton(1, 5, current_seed(), skip = skip);
+  vector<vector<double>> samples = rhalton(1, 5, current_seed(), skip = skip);
 
   // Get spatial box of ray_source_
   SpatialBox* sb = dynamic_cast<SpatialBox*>(
@@ -651,10 +651,10 @@ SourceSite RandomRay::sample_halton()
            xi * ((sb->upper_right() - shift) - (sb->lower_left() + shift));
 
   // Sample Polar cosine and azimuthal angles
-  float mu = 2.0 * samples[0][3] - 1.0;
-  float azi = 2.0 * PI * samples[0][4];
+  double mu = 2.0 * samples[0][3] - 1.0;
+  double azi = 2.0 * PI * samples[0][4];
   // Convert to Cartesian coordinates
-  float c = std::sqrt(1.0 - mu * mu);
+  double c = std::sqrt(1.0 - mu * mu);
   site.u.x = mu;
   site.u.y = std::cos(azi) * c;
   site.u.z = std::sin(azi) * c;

tests/regression_tests/random_ray_halton_samples/results_true.dat

@@ -1,62 +1,62 @@
 k-combined:
-8.388050E-01 7.383283E-03
+8.388051E-01 7.383265E-03
 tally 1:
-5.033306E+00
-5.072159E+00
+5.033308E+00
+5.072162E+00
 1.917335E+00
 7.360725E-01
 4.666410E+00
 4.360038E+00
-2.851811E+00
-1.629361E+00
-4.365591E-01
-3.818885E-02
+2.851812E+00
+1.629362E+00
+4.365590E-01
+3.818884E-02
 1.062497E+00
-2.262072E-01
-1.697620E+00
-5.829331E-01
+2.262071E-01
+1.697621E+00
+5.829333E-01
 5.639912E-02
-6.427569E-04
+6.427568E-04
 1.372642E-01
 3.807294E-03
-2.376682E+00
+2.376683E+00
 1.151027E+00
-8.060903E-02
+8.060902E-02
 1.323179E-03
-1.961863E-01
-7.837694E-03
-7.145451E+00
+1.961862E-01
+7.837693E-03
+7.145452E+00
 1.037540E+01
-8.551805E-02
-1.486270E-03
-2.081364E-01
-8.803959E-03
+8.551803E-02
+1.486269E-03
+2.081363E-01
+8.803955E-03
 2.053205E+01
-8.469503E+01
-3.235620E-02
-2.102893E-04
-8.006316E-02
-1.287560E-03
-1.326546E+01
-3.519488E+01
-1.867472E-01
-6.975145E-03
-5.194280E-01
-5.396293E-02
+8.469498E+01
+3.235618E-02
+2.102891E-04
+8.006311E-02
+1.287559E-03
+1.326545E+01
+3.519484E+01
+1.867471E-01
+6.975133E-03
+5.194275E-01
+5.396284E-02
 7.558115E+00
 1.142535E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-3.386210E+00
+3.386211E+00
 2.294414E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
 1.827274E+00
-6.782303E-01
+6.782305E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -74,7 +74,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 1.828685E+01
-6.719608E+01
+6.719606E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -86,7 +86,7 @@ tally 1:
 0.000000E+00
 0.000000E+00
 4.590713E+00
-4.224106E+00
+4.224107E+00
 1.705847E+00
 5.831967E-01
 4.151691E+00
@@ -98,8 +98,8 @@ tally 1:
 9.911973E-01
 1.970084E-01
 1.664732E+00
-5.598667E-01
-5.385644E-02
+5.598668E-01
+5.385645E-02
 5.858353E-04
 1.310758E-01
 3.470126E-03
@@ -111,61 +111,61 @@ tally 1:
 7.080943E-03
 7.105765E+00
 1.025960E+01
-8.287513E-02
+8.287512E-02
 1.396474E-03
 2.017039E-01
-8.272054E-03
+8.272053E-03
 2.099024E+01
-8.854252E+01
+8.854251E+01
 3.191885E-02
-2.048369E-04
-7.898097E-02
-1.254176E-03
+2.048368E-04
+7.898095E-02
+1.254175E-03
 1.355820E+01
-3.676863E+01
+3.676862E+01
 1.815102E-01
-6.590730E-03
-5.048615E-01
-5.098892E-02
+6.590727E-03
+5.048614E-01
+5.098890E-02
 5.093659E+00
 5.192360E+00
-1.874631E+00
-7.031791E-01
-4.562477E+00
-4.165198E+00
+1.874632E+00
+7.031793E-01
+4.562478E+00
+4.165199E+00
 2.870213E+00
 1.650126E+00
 4.244068E-01
-3.608744E-02
+3.608745E-02
 1.032921E+00
-2.137597E-01
+2.137598E-01
 1.703400E+00
-5.873028E-01
-5.464556E-02
-6.042852E-04
+5.873029E-01
+5.464557E-02
+6.042855E-04
 1.329964E-01
-3.579412E-03
+3.579413E-03
 2.389118E+00
-1.163673E+00
-7.832235E-02
-1.251253E-03
-1.906209E-01
-7.411650E-03
+1.163674E+00
+7.832237E-02
+1.251254E-03
+1.906210E-01
+7.411654E-03
 7.162706E+00
 1.042515E+01
-8.273829E-02
+8.273831E-02
 1.391799E-03
 2.013709E-01
-8.244357E-03
-2.043146E+01
-8.383560E+01
+8.244359E-03
+2.043145E+01
+8.383557E+01
 3.096158E-02
 1.924116E-04
-7.661227E-02
-1.178099E-03
+7.661226E-02
+1.178098E-03
 1.314148E+01
-3.454146E+01
-1.771733E-01
-6.279891E-03
-4.927985E-01
-4.858413E-02
+3.454143E+01
+1.771732E-01
+6.279887E-03
+4.927984E-01
+4.858410E-02