Add Philox PRNG

Add implementation and unit tests for Philox RNG.

vtkm/CMakeLists.txt

@@ -111,3 +111,6 @@ add_subdirectory(io)
 
 #add the source folder
 add_subdirectory(source)
+
+#add Pseudo Random Number Generator folder
+add_subdirectory(prng)

vtkm/prng/CMakeLists.txt

@@ -0,0 +1,15 @@
+##============================================================================
+##  Copyright (c) Kitware, Inc.
+##  All rights reserved.
+##  See LICENSE.txt for details.
+##
+##  This software is distributed WITHOUT ANY WARRANTY; without even
+##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+##  PURPOSE.  See the above copyright notice for more information.
+##============================================================================
+
+set(headers
+  Philox.h)
+
+#-----------------------------------------------------------------------------
+add_subdirectory(testing)

vtkm/prng/Philox.h

@@ -0,0 +1,140 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#ifndef vtk_m_prng_Philox_h
+#define vtk_m_prng_Philox_h
+
+#include <vtkm/Types.h>
+
+namespace vtkm
+{
+namespace prng
+{
+namespace detail
+{
+constexpr VTKM_EXEC_CONT vtkm::Vec<vtkm::UInt32, 2> mulhilo(vtkm::UInt32 a, vtkm::UInt32 b)
+{
+  vtkm::UInt64 r = static_cast<vtkm::UInt64>(a) * b;
+  vtkm::UInt32 lo = static_cast<vtkm::UInt32>(r);
+  vtkm::UInt32 hi = r >> 32;
+  return { lo, hi };
+}
+
+// FIXME: what to do with CUDA backend?
+constexpr VTKM_EXEC_CONT vtkm::Vec<vtkm::UInt64, 2> mulhilo(vtkm::UInt64 a, vtkm::UInt64 b)
+{
+  __uint128_t r = static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b);
+  vtkm::UInt64 lo = static_cast<vtkm::UInt64>(r);
+  vtkm::UInt64 hi = r >> 64;
+  return { lo, hi };
+}
+
+template <typename UIntType, std::size_t N, UIntType... consts>
+class philox_parameters;
+
+template <typename T, T M0, T C0>
+struct philox_parameters<T, 2, M0, C0>
+{
+  static const vtkm::Vec<T, 1> multipliers;
+  static const vtkm::Vec<T, 1> round_consts;
+};
+
+template <typename T, T M0, T C0>
+const vtkm::Vec<T, 1> vtkm::prng::detail::philox_parameters<T, 2, M0, C0>::multipliers =
+  vtkm::Vec<T, 1>(M0);
+template <typename T, T M0, T C0>
+const vtkm::Vec<T, 1> vtkm::prng::detail::philox_parameters<T, 2, M0, C0>::round_consts =
+  vtkm::Vec<T, 1>(C0);
+
+// TODO: make it work with C++11
+template <typename T, T M0, T C0, T M1, T C1>
+struct philox_parameters<T, 4, M0, C0, M1, C1>
+{
+  static const vtkm::Vec<T, 2> multipliers{ M0, M1 };
+  static const vtkm::Vec<T, 2> round_consts{ C0, C1 };
+};
+
+template <typename UIntType, std::size_t N, std::size_t R, UIntType... consts>
+class philox_functor;
+
+template <typename UIntType, std::size_t R, UIntType... consts>
+class philox_functor<UIntType, 2, R, consts...>
+{
+public:
+  using counters_type = vtkm::Vec<UIntType, 2>;
+  using keys_type = vtkm::Vec<UIntType, 1>;
+
+  constexpr VTKM_EXEC_CONT counters_type operator()(counters_type counters, keys_type keys) const
+  {
+    for (std::size_t i = 0; i < R; ++i)
+    {
+      counters = round(counters, keys);
+      keys = bump_keys(keys);
+    }
+    return counters;
+  }
+
+private:
+  static constexpr VTKM_EXEC_CONT counters_type round(counters_type counters, keys_type round_keys)
+  {
+    vtkm::Vec<UIntType, 2> r =
+      mulhilo(philox_parameters<UIntType, 2, consts...>::multipliers[0], counters[0]);
+    return { r[1] ^ round_keys[0] ^ counters[1], r[0] };
+  }
+
+  static constexpr VTKM_EXEC_CONT keys_type bump_keys(keys_type keys)
+  {
+    return { keys[0] + philox_parameters<UIntType, 2, consts...>::round_consts[0] };
+  }
+};
+
+template <typename UIntType, std::size_t R, UIntType... consts>
+class philox_functor<UIntType, 4, R, consts...>
+{
+  using counters_type = vtkm::Vec<UIntType, 4>;
+  using keys_type = vtkm::Vec<UIntType, 2>;
+
+  static constexpr VTKM_EXEC_CONT counters_type round(counters_type counters, keys_type round_keys)
+  {
+    vtkm::Vec<UIntType, 2> r0 =
+      mulhilo(philox_parameters<UIntType, 4, consts...>::multipliers[0], counters[0]);
+    vtkm::Vec<UIntType, 2> r1 =
+      mulhilo(philox_parameters<UIntType, 4, consts...>::multipliers[1], counters[2]);
+    return {
+      r1[1] ^ round_keys[0] ^ counters[1], r1[0], r0[1] ^ round_keys[1] ^ counters[3], r0[0]
+    };
+  }
+
+  static constexpr VTKM_EXEC_CONT keys_type bump_key(keys_type keys)
+  {
+    keys[0] += philox_parameters<UIntType, 4, consts...>::round_consts[0];
+    keys[1] += philox_parameters<UIntType, 4, consts...>::round_consts[1];
+    return keys;
+  }
+
+public:
+  constexpr VTKM_EXEC_CONT counters_type operator()(counters_type counters, keys_type keys) const
+  {
+    for (std::size_t i = 0; i < R; ++i)
+    {
+      counters = round(counters, keys);
+      keys = bump_key(keys);
+    }
+    return counters;
+  }
+};
+
+} // namespace detail
+
+using philox_functor2x32x7 = detail::philox_functor<vtkm::UInt32, 2, 7, 0xD256D193, 0x9E3779B9>;
+using philox_functor2x32x10 = detail::philox_functor<vtkm::UInt32, 2, 10, 0xD256D193, 0x9E3779B9>;
+
+} // namespace prng
+} // namespace vtkm
+#endif //vtk_m_prng_Philox_h

vtkm/prng/testing/CMakeLists.txt

@@ -0,0 +1,15 @@
+##============================================================================
+##  Copyright (c) Kitware, Inc.
+##  All rights reserved.
+##  See LICENSE.txt for details.
+##
+##  This software is distributed WITHOUT ANY WARRANTY; without even
+##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+##  PURPOSE.  See the above copyright notice for more information.
+##============================================================================
+
+set(unit_tests
+  UnitTestPhiloxRNG.cxx)
+
+vtkm_unit_tests(SOURCES ${unit_tests}
+  ALL_BACKENDS)

vtkm/prng/testing/UnitTestPhiloxRNG.cxx

@@ -0,0 +1,49 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/prng/Philox.h>
+
+void TestPhiloxRNG2x32x7()
+{
+  vtkm::prng::philox_functor2x32x7 f2x32x7;
+  using counters_type = typename vtkm::prng::philox_functor2x32x7::counters_type;
+
+  // test cases from original Random123 implementation
+  VTKM_TEST_ASSERT(f2x32x7({}, {}) == counters_type{ 0x257a3673, 0xcd26be2a });
+  VTKM_TEST_ASSERT(f2x32x7({ 0xffffffff, 0xffffffff }, { 0xffffffff }) ==
+                   counters_type{ 0xab302c4d, 0x3dc9d239 });
+  VTKM_TEST_ASSERT(f2x32x7({ 0x243f6a88, 0x85a308d3 }, { 0x13198a2e }) ==
+                   counters_type{ 0xbedbbe6b, 0xe4c770b3 });
+}
+
+void TestPhiloxRNG2x32x10()
+{
+  vtkm::prng::philox_functor2x32x10 f2x32x10;
+  using counters_type = typename vtkm::prng::philox_functor2x32x10::counters_type;
+
+  // test cases from original Random123 implementation
+  VTKM_TEST_ASSERT(f2x32x10({}, {}) == counters_type{ 0xff1dae59, 0x6cd10df2 });
+  VTKM_TEST_ASSERT(f2x32x10({ 0xffffffff, 0xffffffff }, { 0xffffffff }) ==
+                   counters_type{ 0x2c3f628b, 0xab4fd7ad });
+  VTKM_TEST_ASSERT(f2x32x10({ 0x243f6a88, 0x85a308d3 }, { 0x13198a2e }) ==
+                   counters_type{ 0xdd7ce038, 0xf62a4c12 });
+}
+
+void TestPhiloxRNG()
+{
+  TestPhiloxRNG2x32x7();
+  TestPhiloxRNG2x32x10();
+}
+
+int UnitTestPhiloxRNG(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestPhiloxRNG, argc, argv);
+}