Add uniform range function: uniform(low, high). int64 not supported

include/openrand/base_state.h

@@ -87,27 +87,51 @@ class BaseRNG {
       if constexpr (std::is_integral_v<T>)
         return static_cast<T>(x);
       else
-        return uniform<float, uint32_t>(x);
+        return u01<float, uint32_t>(x);
     } else {
       const uint64_t x = gen().template draw<uint64_t>();
       if constexpr (std::is_integral_v<T>)
         return static_cast<T>(x);
       else
-        return uniform<double, uint64_t>(x);
+        return u01<double, uint64_t>(x);
     }
   }
 
+  /**
+   * @brief Generates a number from a uniform distribution between a and b.
+   *
+   * @Note For integer type, this method is slightly biased towards lower numbers.
+   *
+   * @tparam T Data type to be returned. Can be float or double.
+   * double.
+   * @param low lower bound of the uniform distribution
+   * @param high upper bound of the uniform distribution
+   * @return T random number from a uniform distribution between a and b
+   */
   template <typename T = float>
-  DEVICE void fill_random(T *array, const int N) {
-    for (int i = 0; i < N; i++) array[i] = rand<T>();
+  DEVICE T uniform(const T low, const T high) {
+    // TODO: Allow 64 bit integers
+    static_assert(!(std::is_integral_v<T> && sizeof(T) > sizeof(int32_t)),
+                  "64 bit int not yet supported");
+
+    T range = high - low;
+
+    if constexpr (std::is_floating_point_v<T>) {
+      return low + range * rand<T>();
+    } else if constexpr (std::is_integral_v<T>) {
+      // Thanks to (https://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-modulo-reduction/)
+      return low + T(((uint64_t)rand<uint32_t>() * (uint64_t)range) >> 32);
+    }
   }
 
-  template <typename Ftype, typename Utype>
-  inline DEVICE Ftype uniform(const Utype in) const {
-    constexpr Ftype factor =
-        Ftype(1.) / (Ftype(~static_cast<Utype>(0)) + Ftype(1.));
-    constexpr Ftype halffactor = Ftype(0.5) * factor;
-    return Utype(in) * factor + halffactor;
+  /*
+   * @brief Fills an array with random numbers from a uniform distribution [0, 1)
+   *
+   * @Note This array is filled serially. `N` ideally should not be large.
+   */
+  template <typename T = float>
+  DEVICE void fill_random(T *array, const int N) {
+    for (int i = 0; i < N; i++) array[i] = rand<T>();
   }
 
   /**
@@ -199,13 +223,26 @@ class BaseRNG {
     }
   }
 
+ protected:
+  /*
+   * @brief Converts a random number integer to a floating point number between [0., 1.)
+   */
+  template <typename Ftype, typename Utype>
+  inline DEVICE Ftype u01(const Utype in) const {
+    constexpr Ftype factor =
+        Ftype(1.) / (Ftype(~static_cast<Utype>(0)) + Ftype(1.));
+    constexpr Ftype halffactor = Ftype(0.5) * factor;
+    return Utype(in) * factor + halffactor;
+  }
+
  private:
   /**
    * @brief Returns a reference to the random number generator.
    */
   DEVICE __inline__ RNG &gen() {
     return *static_cast<RNG *>(this);
   }
+
 };  // class BaseRNG
 
 }  // namespace openrand

include/openrand/phillox.h

@@ -78,8 +78,8 @@ class Phillox : public BaseRNG<Phillox> {
   openrand::float4 draw_float4() {
     generate();
     return openrand::float4{
-        uniform<float, uint32_t>(_out[0]), uniform<float, uint32_t>(_out[1]),
-        uniform<float, uint32_t>(_out[2]), uniform<float, uint32_t>(_out[3])};
+        u01<float, uint32_t>(_out[0]), u01<float, uint32_t>(_out[1]),
+        u01<float, uint32_t>(_out[2]), u01<float, uint32_t>(_out[3])};
   }
 
  private:

tests/test_uniform.cpp

@@ -73,6 +73,43 @@ TEST(RNG, basic) {
   test_basic<openrand::Squares>();
 }
 
+template <typename RNG>
+void test_range(){
+  RNG rng (1234567, 1234567);
+
+  double mean = 0;
+  for(int i = 0; i < 1000; i++){
+    auto x = rng.uniform(10.0, 20.0);
+    mean += x;
+    EXPECT_TRUE((x >= 10.0) && (x <= 20.0));
+  }
+  EXPECT_NEAR(mean / 1000.0, 15.0, 0.36);  // 99.99 % confidence
+
+  mean = 0;
+  for(int i = 0; i < 1000; i++){
+    auto x = rng.uniform(-20.0, -10.0);
+    mean += x;
+    EXPECT_TRUE((x >= -20.0) && (x <= -10.0));
+  }
+  EXPECT_NEAR(mean / 1000.0, -15.0, 0.36);
+
+  // For integer type, this method is slightly biased towards lower numbers.
+  mean = 0;
+  for(int i = 0; i < 1000; i++){
+    auto x = rng.template uniform<int>(10,20);
+    mean += (float)x;
+    EXPECT_TRUE((x >= 10) && (x <= 20));
+  }
+  EXPECT_NEAR(mean / 1000.0, 15.0, 1.0);
+}
+
+TEST(RNG, range){
+  test_range<openrand::Phillox>();
+  test_range<openrand::Tyche>();
+  test_range<openrand::Threefry>();
+  test_range<openrand::Squares>();
+}
+
 template <typename RNG>
 void test_mean() {
   RNG rng(0, 0);