From 8e2f918dd3cbc135816367ac9ff6a1c558b7c964 Mon Sep 17 00:00:00 2001
From: NunoSempere <nuno.sempere@protonmail.com>
Date: Sat, 13 Jan 2024 12:47:14 +0100
Subject: [PATCH] add comment about cache analysis

---
 squiggle_more.c | 36 +++++++++++++++++++++++++++++-------
 1 file changed, 29 insertions(+), 7 deletions(-)

diff --git a/squiggle_more.c b/squiggle_more.c
index ae1caeb..8b87261 100644
--- a/squiggle_more.c
+++ b/squiggle_more.c
@@ -51,7 +51,6 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
         // - xorshift can't start with 0
         // - the seeds should be reasonably separated and not correlated
         cache_box[i].seed = (uint64_t)rand() * (UINT64_MAX / RAND_MAX);
-        // printf("#%ld: %lu\n",i, *seeds[i]);
 
         // Other initializations tried:
         // *seeds[i] = 1 + i;
@@ -69,17 +68,40 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
 
             for (int j = lower_bound_inclusive; j < upper_bound_not_inclusive; j++) {
                 results[j] = sampler(&(cache_box[i].seed));
-                // In principle, these results[j] could also result in two threads competing for the same cache line.
-                // In practice, though,
-                // a) this would happen infrequently
-                // b) trying to unroll loops actually makes the code slower
-                // c) 8 results[j] are 8 doubles, which fit a cache line. If n_samples/n_threads
+                /* 
+                t starts at 0 and ends at T
+                at t=0, 
+                  thread i accesses: results[i*quotient +0], 
+                  thread i+1 acccesses: results[(i+1)*quotient +0]
+                at t=T
+                  thread i accesses: results[(i+1)*quotient -1]
+                  thread i+1 acccesses: results[(i+2)*quotient -1]
+                The results[j] that are directly adjacent are 
+                  results[(i+1)*quotient -1] (accessed by thread i at time T)
+                  results[(i+1)*quotient +0] (accessed by thread i+1 at time 0)
+                and these are themselves adjacent to
+                  results[(i+1)*quotient -2] (accessed by thread i at time T-1)
+                  results[(i+1)*quotient +1] (accessed by thread i+1 at time 2)
+                If T is large enough, which it is, two threads won't access the same
+                cache line at the same time.
+                Pictorially:
+                at t=0 ....i.........I.........
+                at t=T .............i.........I
+                and the two never overlap
+                Note that results[j] is a double, a double has 8 bytes (64 bits)
+                8 doubles fill a cache line of 64 bytes.
+                So we specifically won't get problems as long as n_samples/n_threads > 8
+                n_threads is normally 16, so n_samples > 128 
+                Note also that this is only a problem in terms of speed, if n_samples<128
+                the results are still computed, it'll just be slower
+                */
             }
         }
     }
     for (int j = divisor_multiple; j < n_samples; j++) {
         results[j] = sampler(&(cache_box[0].seed));
-        // we can just reuse a seed, this isn't problematic because we are not doing multithreading
+        // we can just reuse a seed, 
+        // this isn't problematic because we;ve now stopped doing multithreading
     }
 
     free(cache_box);