From c7212595564e2477955ee4f2c9c7ba627184aa83 Mon Sep 17 00:00:00 2001
From: Zaher Salman <zaher.salman@gmail.com>
Date: Fri, 19 Jun 2026 17:14:53 +0200
Subject: [PATCH] Forgottent seeded pcg32

---
 fortran/Makefile    |  2 +-
 fortran/rng_pcg32.c | 35 ++++++++++++++++++++++++++++-------
 2 files changed, 29 insertions(+), 8 deletions(-)

diff --git a/fortran/Makefile b/fortran/Makefile
index 4442afb..612a5cd 100644
--- a/fortran/Makefile
+++ b/fortran/Makefile
@@ -5,7 +5,7 @@
 # 4-Feb-2013 ZS, clean up makefile and include random number generator
 #            in the code. Use gfortran which is now a standard part of gcc.
 # 18-Jun-2026 ZS, add explicit release/debug/profile build modes.
-# 19-Jun-2026 ZA, add optional random number generators to improve speed
+# 19-Jun-2026 ZS, add optional random number generators to improve speed
 
 FC = gfortran
 PROGRAM = trimspNL
diff --git a/fortran/rng_pcg32.c b/fortran/rng_pcg32.c
index 3f3889d..2048ca6 100644
--- a/fortran/rng_pcg32.c
+++ b/fortran/rng_pcg32.c
@@ -9,8 +9,13 @@
  * which gfortran resolves to pcg32_random_uniform_().
  *
  * The generator is deterministic and self-contained.  It returns REAL*4
- * values in (0,1) to avoid exact zero in legacy transport code.  The high 24 bits of each 32-bit PCG output are used,
- * matching the precision of a single-precision uniform variate.
+ * values in (0,1) to avoid exact zero in legacy transport code.
+ * The high 24 bits of each 32-bit PCG output are used, matching the
+ * precision of a single-precision uniform variate.
+ *
+ * trimspNL calls pcg32_seed_() with the three existing input seeds
+ * RI, RI2 and RI3.  This keeps PCG32 runs reproducible from the same
+ * input deck without changing the file format.
  */
 
 static uint64_t pcg32_state = 0x853c49e6748fea9bULL;
@@ -61,12 +66,28 @@ pcg32_random_uniform_(float *values, int *n_values)
     }
 }
 
-/* Optional explicit seed hook for later experiments.  Not used by trimspNL yet. */
-void
-pcg32_seed_(int *seed1, int *seed2)
+static uint64_t
+pcg32_mix(uint64_t x)
 {
-    uint64_t initstate = (uint32_t)(*seed1);
-    uint64_t initseq   = (uint32_t)(*seed2);
+    x += 0x9e3779b97f4a7c15ULL;
+    x = (x ^ (x >> 30u)) * 0xbf58476d1ce4e5b9ULL;
+    x = (x ^ (x >> 27u)) * 0x94d049bb133111ebULL;
+    return x ^ (x >> 31u);
+}
+
+void
+pcg32_seed_(int *seed1, int *seed2, int *seed3)
+{
+    uint64_t s1 = (uint32_t)(*seed1);
+    uint64_t s2 = (uint32_t)(*seed2);
+    uint64_t s3 = (uint32_t)(*seed3);
+    uint64_t initstate;
+    uint64_t initseq;
+
+    initstate = pcg32_mix(s1 ^ (s2 << 21u) ^ (s3 << 42u)
+        ^ 0x4d595df4d0f33173ULL);
+    initseq = pcg32_mix(s3 ^ (s2 << 17u) ^ (s1 << 34u)
+        ^ 0x14057b7ef767814fULL);
 
     pcg32_state = 0U;
     pcg32_inc = (initseq << 1u) | 1u;