test Philox against known data

pystencils_tests/test_random.py

@@ -4,6 +4,12 @@ import pystencils as ps
 from pystencils.rng import PhiloxFourFloats, PhiloxTwoDoubles
 
 
+# curand_Philox4x32_10(make_uint4(124, i, j, 0), make_uint2(0, 0))
+philox_reference = np.array([[[3576608082, 1252663339, 1987745383,  348040302],
+                              [1032407765,  970978240, 2217005168, 2424826293]],
+                             [[2958765206, 3725192638, 2623672781, 1373196132],
+                              [ 850605163, 1694561295, 3285694973, 2799652583]]])
+
 def test_philox_double():
     for target in ('cpu', 'gpu'):
         dh = ps.create_data_handling((2, 2), default_ghost_layers=0, default_target=target)
@@ -24,6 +30,12 @@ def test_philox_double():
         arr = dh.gather_array('f')
         assert np.logical_and(arr <= 1.0, arr >= 0).all()
 
+        x = philox_reference[:,:,0::2]
+        y = philox_reference[:,:,1::2]
+        z = x ^ y << (53 - 32)
+        double_reference = z * 2.**-53 + 2.**-54
+        assert(np.allclose(arr, double_reference, rtol=0, atol=np.finfo(np.float64).eps))
+
 
 def test_philox_float():
     for target in ('cpu', 'gpu'):
@@ -41,3 +53,6 @@ def test_philox_float():
         dh.all_to_cpu()
         arr = dh.gather_array('f')
         assert np.logical_and(arr <= 1.0, arr >= 0).all()
+
+        float_reference = philox_reference * 2.**-32 + 2.**-33
+        assert(np.allclose(arr, float_reference, rtol=0, atol=np.finfo(np.float32).eps))