Probabilistic linear algebra may fail when rand() == 0

[wegank]

In #343 (comment), I found an issue where on Windows,

$ ./msolve.exe -f input_files/eco10-31.ms -o out.ms --random-seed 1779878689 -P 2 -d 0 -L 0 -l 44 -t 1 -v 2

writes nothing to out.ms. Actually, when I look at the first line where the results differ between macOS and Windows, I see that it says 0 new 9 zero instead of 1 new 8 zero.

I initially thought that the bug must be Windows-specific, but that's not the case. In fact, the bug is PRNG-specific: on Windows, rand() returns 0 just before the results start to vary! One may say that rand() == 0 isn't too rare there, as RAND_MAX is 32767 on Windows. However, if I force rand() to return 0 at the same call on macOS, it fails there too!

To reproduce the issue, apply the following patch to hijack the PRNG:

diff --git a/src/neogb/la_ff_32.c b/src/neogb/la_ff_32.c
index 90953f2..f351ef2 100644
--- a/src/neogb/la_ff_32.c
+++ b/src/neogb/la_ff_32.c
@@ -28,6 +28,17 @@
 #include <arm_neon.h>
 #endif
 
+static inline int32_t fake_rand(void)
+{
+    static uint64_t call_idx = 0;
+    int32_t rv = rand();
+    if (call_idx == 25789) {
+        rv = 0;
+    }
+    call_idx++;
+    return rv;
+}
+
 static inline cf32_t *normalize_dense_matrix_row_ff_32(
         cf32_t *row,
         const hm_t len,
@@ -2382,7 +2393,7 @@ static void probabilistic_sparse_reduced_echelon_form_ff_32(
 
                 /* fill random value array */
                 for (j = 0; j < nrbl; ++j) {
-                    mull[j] = (int64_t)rand() & mask;
+                    mull[j] = (int64_t)fake_rand() & mask;
                 }
                 /* generate one dense row as random linear combination
                  * of the rows of the block */

Then, build msolve with the patch applied and run the following command. Observe that out.ms is empty.

$ ./msolve -f input_files/eco10-31.ms -o out.ms --random-seed 1779878689 -P 2 -d 0 -L 0 -l 44 -t 1 -v 2

[ederc]

Well, we could pack this into a do while like

do {
  mull[j] = (int64_t)rand() & mask;
} while (mull[j] == 0);

Okay, this might run into an infinite loop, but adding 1 puts a bias on 1 thus is also not a good solution.