Improved random number generation (#16)

* Implement mt19937 random number generator.
- Rendering speed improved up to 10 times.
- Need to implement thread-local random number generators.

* Fix macOS build errors and off-by-one error.

* Add checks for thread_local compiler support.
- Fix minor typo when setting compiler flags for debug target.

* Add compiler-specific thread local storage.
- It turns out that macOS and Windows runners doesn't support thread local variables from C11
- Implementation tested only with Linux so far.
- Force Github Actions to use at least 2 threads for rendering.

* Fix MSVC detection for thread local storage.

.github/workflows/cmake.yml

@@ -51,7 +51,7 @@ jobs:
     - name: Render Sample Image
       working-directory: ${{runner.workspace}}/build/bin
       shell: bash
-      run: ./ray_tracing_one_week --verbose --width 300 --height 200 -s 100 --scene metal_test sample_image.ppm
+      run: ./ray_tracing_one_week --verbose -t 2 --width 300 --height 200 -s 100 --scene metal_test sample_image.ppm
 
     - name: Upload Sample Image
       uses: actions/upload-artifact@v2

CMakeLists.txt

@@ -5,6 +5,8 @@ project(ray_tracing_one_week C)
 set(CMAKE_C_STANDARD 11)
 
 include(CheckLibraryExists)
+include(CheckIncludeFile)
+include(CheckSymbolExists)
 
 CHECK_LIBRARY_EXISTS(m ceil "" HAVE_LIB_M)
 if (HAVE_LIB_M)
@@ -34,6 +36,8 @@ add_executable(ray_tracing_one_week rt_camera.c rt_colour.c rt_aabb.c rt_perlin.
                textures/rt_texture_noise.c textures/rt_texture_image.c
                # Scenes
                scenes/rt_scenes.c
+               # Random number generation
+               random/rt_random.c random/rt_random.h
                # Threading
                ${THREADING_IMPLEMENTATION} threads/rt_thread_pool.c)
 
@@ -49,7 +53,7 @@ endif ()
 if (CMAKE_BUILD_TYPE EQUAL "DEBUG")
     # Enable inline optimization under debug configurations.
     target_compile_options(ray_tracing_one_week PRIVATE $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>,$<CXX_COMPILER_ID:GNU>>:
-                           -finline-funcitons>
+                           -finline-functions>
                            $<$<CXX_COMPILER_ID:MSVC>:
                            /Ob>)
 endif ()

main.c

@@ -133,6 +133,8 @@ int main(int argc, char const *argv[])
     bool verbose = false;
     bool render_recursive = true;
 
+    rt_random_seed(RT_RANDOM_DEFAULT_SEED);
+
     // Parse console arguments
     for (int i = 1; i < argc; ++i)
     {

random/rt_random.c

@@ -0,0 +1,110 @@
+/**
+ * Copyright (c) 2021, Evgeniy Morozov
+ * All rights reserved.
+ * This source code is licensed under the BSD-style license found in the
+ * LICENSE file in the root directory of this source tree.
+ */
+#include <stdlib.h>
+#include <assert.h>
+#include "rt_random.h"
+#include "rt_thread.h"
+
+struct rt_mt19937_s
+{
+    int w, n, m, r;
+
+    uint64_t a;
+    uint64_t u, d;
+    uint64_t s, b;
+    uint64_t t, c;
+
+    int l;
+    uint64_t f;
+
+    uint64_t upper_mask, lower_mask;
+
+    uint64_t index;
+    uint64_t *MT;
+};
+
+static RT_THREAD_LOCAL struct rt_mt19937_s gs_generator = {
+    .index = UINT64_MAX
+};
+
+static void twist(struct rt_mt19937_s *gen);
+
+void rt_random_seed(uint64_t seed)
+{
+    gs_generator.w = 64;
+    gs_generator.n = 312;
+    gs_generator.m = 156;
+    gs_generator.r = 31;
+
+    gs_generator.a = 0xB5026F5AA96619E9ULL;
+
+    gs_generator.u = 29;
+    gs_generator.d = 0x5555555555555555ULL;
+
+    gs_generator.s = 17;
+    gs_generator.b = 0x71D67FFFEDA60000;
+
+    gs_generator.t = 37;
+    gs_generator.c = 0xFFF7EEE000000000;
+
+    gs_generator.l = 43;
+
+    gs_generator.f = 6364136223846793005;
+
+    gs_generator.lower_mask = (1 << gs_generator.r) - 1;
+    gs_generator.upper_mask = ~gs_generator.lower_mask;
+
+    gs_generator.MT = malloc(sizeof(*gs_generator.MT) * gs_generator.n);
+    assert(NULL != gs_generator.MT);
+
+    gs_generator.index = gs_generator.n;
+    gs_generator.MT[0] = seed;
+
+    for (int i = 1; i < gs_generator.n; ++i)
+    {
+        gs_generator.MT[i] = gs_generator.f * (gs_generator.MT[i - 1] ^ (gs_generator.MT[i - 1] >> (gs_generator.w - 2))) + i;
+    }
+}
+
+uint64_t rt_random(void)
+{
+    if (gs_generator.index >= gs_generator.n)
+    {
+        if (gs_generator.index > gs_generator.n)
+        {
+            rt_random_seed(RT_RANDOM_DEFAULT_SEED);
+        }
+
+        twist(&gs_generator);
+    }
+
+    uint64_t y = gs_generator.MT[gs_generator.index++];
+    y ^= (y >> gs_generator.u) & gs_generator.d;
+    y ^= (y << gs_generator.s) & gs_generator.b;
+    y ^= (y << gs_generator.t) & gs_generator.c;
+    y ^= y >> gs_generator.l;
+
+    return y;
+}
+
+static void twist(struct rt_mt19937_s *gen)
+{
+    assert(NULL != gen);
+
+    for (int i = 0; i < gen->n; ++i)
+    {
+        uint64_t x = (gen->MT[i] & gen->upper_mask) + (gen->MT[(i + 1) % gen->n] & gen->lower_mask);
+        uint64_t xA = x >> 1;
+
+        if (x % 2 != 0)
+        {
+            xA ^= gen->a;
+        }
+        gen->MT[i] = gen->MT[(i + gen->m) % gen->n] ^ xA;
+    }
+    gen->index = 0;
+}

random/rt_random.h

@@ -0,0 +1,22 @@
+/**
+ * Copyright (c) 2021, Evgeniy Morozov
+ * All rights reserved.
+ * This source code is licensed under the BSD-style license found in the
+ * LICENSE file in the root directory of this source tree.
+ */
+#ifndef RAY_TRACING_ONE_WEEK_RT_RANDOM_H
+#define RAY_TRACING_ONE_WEEK_RT_RANDOM_H
+
+#include <stdint.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <math.h>
+
+#define RT_RANDOM_MAX UINT64_MAX
+#define RT_RANDOM_DEFAULT_SEED 5489
+
+void rt_random_seed(uint64_t seed);
+
+uint64_t rt_random(void);
+
+#endif // RAY_TRACING_ONE_WEEK_RT_RANDOM_H

rt_weekend.h

@@ -11,6 +11,7 @@
 #include <stdlib.h>
 #include <limits.h>
 #include <math.h>
+#include <random/rt_random.h>
 
 #ifdef M_PI
 #define PI M_PI
@@ -22,7 +23,7 @@
 
 static inline double rt_random_double(double min, double max)
 {
-    return min + (max - min) * rand() / (RAND_MAX + 1.0);
+    return min + (max - min) * rt_random() / (RT_RANDOM_MAX + 0.0);
 }
 
 static inline double rt_clamp(double x, double min, double max)

threads/rt_sync.h

@@ -10,7 +10,6 @@
 int rt_sync_get_number_of_cores(void);
 
 
-
 typedef struct rt_mutex_s rt_mutex_t;
 
 rt_mutex_t *rt_mutex_init(void);
@@ -22,7 +21,6 @@ int rt_mutex_unlock(rt_mutex_t *mutex);
 void rt_mutex_deinit(rt_mutex_t *mutex);
 
 
-
 typedef struct rt_cond_s rt_cond_t;
 
 rt_cond_t *rt_cond_init(void);

threads/rt_thread.h

@@ -8,6 +8,15 @@
 #ifndef RAY_TRACING_ONE_WEEK_RT_THREAD_H
 #define RAY_TRACING_ONE_WEEK_RT_THREAD_H
 
+// TODO: Add-in compiler version check
+#if defined(__GNUC__) || defined(__clang__)
+#define RT_THREAD_LOCAL __thread
+#elif defined(_MSC_VER)
+#define RT_THREAD_LOCAL __declspec(thread)
+#else
+#warning "Thread local storage is not supported for this compiler, there might be artifacts in the rendered image in case of multi-threaded rendering"
+#endif
+
 typedef struct rt_thread_s rt_thread_t;
 
 typedef void (*rt_thread_fn_t)(void *params);