Disable multithreaded gr_mpoly division functions on arm64

[fredrik-johansson]

The multithreaded gr_mpoly division fails sporadically on M1; see #2815 (comment).

The first patch in this PR is a plausible but untested fix generated by Claude Opus 5.

The second patch conservatively disables the multithreaded gr_mpoly division code when FLINT_KNOW_STRONG_ORDER == 0. This could be rolled back once someone has stress tested the multithreaded division on M1 or similar.

[oscarbenjamin]

Is there something that I can help with testing here?

I have tested on Mac M3 and the original reproducer from #2815 (comment) fails about 5% of the time with the exact hardcoded polynomials but I guess it would be less with randomly generated polynomials.

I asked codex to come up with a more reliable test and it came up with this that failed 29/30 attempts on this machine (not with this PR):

Test patch

diff --git a/src/gr_mpoly/test/t-divides_heap_threaded.c b/src/gr_mpoly/test/t-divides_heap_threaded.c
index c76040987..3243c635d 100644
--- a/src/gr_mpoly/test/t-divides_heap_threaded.c
+++ b/src/gr_mpoly/test/t-divides_heap_threaded.c
@@ -12,6 +12,90 @@
 #include "test_helpers.h"
 #include "gr_mpoly.h"

+static void
+gr_mpoly_divides_heap_threaded_check_large_nmod(flint_rand_t state)
+{
+#if FLINT_USES_PTHREAD && FLINT_BITS == 64
+    const ordering_t orderings[] = {ORD_DEGLEX, ORD_DEGREVLEX};
+    const char * ordering_names[] = {"deglex", "degrevlex"};
+    gr_ctx_t cctx;
+    slong attempt, nvars, ord_index, num_threads;
+
+    gr_ctx_init_nmod(cctx, UWORD(9139524932520349247));
+
+    /* Check f*g/g = f for larger sparse examples over nmod. */
+    for (attempt = 0; attempt < 300; attempt++)
+    {
+        for (nvars = 19; nvars <= 20; nvars++)
+        {
+            for (ord_index = 0; ord_index < 2; ord_index++)
+            {
+                gr_mpoly_ctx_t ctx;
+                gr_mpoly_t f, g, h, q1, q2;
+                int status;
+
+                gr_mpoly_ctx_init(ctx, cctx, nvars, orderings[ord_index]);
+                gr_mpoly_init(f, ctx);
+                gr_mpoly_init(g, ctx);
+                gr_mpoly_init(h, ctx);
+                gr_mpoly_init(q1, ctx);
+                gr_mpoly_init(q2, ctx);
+
+                status = GR_SUCCESS;
+                status |= gr_mpoly_randtest_bits(f, state, 27, 5, ctx);
+                status |= gr_mpoly_randtest_bits(g, state, 23, 5, ctx);
+
+                if (gr_mpoly_is_zero(g, ctx) != T_FALSE)
+                    status |= gr_mpoly_one(g, ctx);
+
+                status |= gr_mpoly_mul(h, f, g, ctx);
+
+                if (status == GR_SUCCESS)
+                {
+                    for (num_threads = 2; num_threads <= 4; num_threads++)
+                    {
+                        int dstatus1, dstatus2;
+
+                        flint_set_num_threads(num_threads);
+                        dstatus1 = gr_mpoly_divides_heap(q1, h, g, ctx);
+                        dstatus2 = gr_mpoly_divides_heap_threaded(q2, h, g, ctx);
+
+                        if (dstatus1 == GR_SUCCESS && dstatus2 == GR_DOMAIN)
+                        {
+                            flint_printf("FAIL: exact division reported as GR_DOMAIN\n");
+                            flint_printf("attempt = %wd, nvars = %wd, ordering = %s, threads = %wd\n",
+                                    attempt, nvars, ordering_names[ord_index], num_threads);
+                            flint_printf("heap = %d, threaded = %d\n", dstatus1, dstatus2);
+                            fflush(stdout);
+                            flint_abort();
+                        }
+
+                        if (dstatus1 == GR_SUCCESS && dstatus2 == GR_SUCCESS &&
+                                gr_mpoly_equal(q1, q2, ctx) == T_FALSE)
+                        {
+                            flint_printf("FAIL: threaded quotient != heap quotient\n");
+                            flint_printf("attempt = %wd, nvars = %wd, ordering = %s, threads = %wd\n",
+                                    attempt, nvars, ordering_names[ord_index], num_threads);
+                            fflush(stdout);
+                            flint_abort();
+                        }
+                    }
+                }
+
+                gr_mpoly_clear(f, ctx);
+                gr_mpoly_clear(g, ctx);
+                gr_mpoly_clear(h, ctx);
+                gr_mpoly_clear(q1, ctx);
+                gr_mpoly_clear(q2, ctx);
+                gr_mpoly_ctx_clear(ctx);
+            }
+        }
+    }
+
+    gr_ctx_clear(cctx);
+#endif
+}
+
 TEST_FUNCTION_START(gr_mpoly_divides_heap_threaded, state)
 {
     slong i, j;
@@ -19,6 +103,8 @@ TEST_FUNCTION_START(gr_mpoly_divides_heap_threaded, state)
     gr_ctx_t ctx1;
     gr_ctx_init_nmod(ctx1, 17);

+    gr_mpoly_divides_heap_threaded_check_large_nmod(state);
+
     /* Check that the threaded quotient matches the serial heap quotient */
     for (i = 0; i < 1000 * flint_test_multiplier(); i++)
     {

With this PR when I run that test it passes with multiplier 1000:

$ FLINT_TEST_MULTIPLIER=1000 build/gr_mpoly/test/main gr_mpoly_divides_heap_threaded
gr_mpoly_divides_heap_threaded...
gr_mpoly_divides_heap_threaded                   87.79   (PASS)

On main it failed 3 times in a row with multiplier 0.01

$ FLINT_TEST_MULTIPLIER=0.01 build/gr_mpoly/test/main gr_mpoly_divides_heap_threaded
gr_mpoly_divides_heap_threaded...
FAIL: exact division reported as GR_DOMAIN
attempt = 38, nvars = 19, ordering = degrevlex, threads = 3
heap = 0, threaded = 1
zsh: abort      FLINT_TEST_MULTIPLIER=0.01 build/gr_mpoly/test/main

The first commit in this PR adds fences while the second just disables threaded division. Using only the first commit the test still passes

$ FLINT_TEST_MULTIPLIER=1000 build/gr_mpoly/test/main gr_mpoly_divides_heap_threaded
gr_mpoly_divides_heap_threaded...
gr_mpoly_divides_heap_threaded                  109.42   (PASS)

Maybe it is worth adding something like that test in a separate PR (with no fixes) to check that CI can fail reliably.

[oscarbenjamin]

I also asked codex to review the first commit:

Findings

High: src/gr_mpoly/divides_heap_threaded.c and src/gr_mpoly/divrem_ideal_heap_threaded.c

The commit uses atomic_thread_fence() around ordinary volatile fields (producer, Q->length). That may work empirically on ARM with the current compiler, but it is not a valid C11 synchronization scheme. A release/acquire fence pair only establishes inter-thread synchronization through atomic operations; here the flag and length loads/stores are still non-atomic, so the code still has data races in the C memory model.

The same issue exists in the ideal variant.

So: yes, the change appears to fix the observed Apple Silicon failure, but as a patch I would be uncomfortable approving it as-is because it relies on volatile + fences rather than actual atomic publication.

The more defensible version would make the synchronization variable atomic, e.g. producer as _Atomic int, with:

atomic_store_explicit(&next->producer, 1, memory_order_release);

and:

if (atomic_load_explicit(&L->producer, memory_order_acquire) == 1)

Then separately, gr_mpoly_ts_struct.length should probably also be atomic if readers use it to decide which coefficient/exponent entries are visible. The existing gr_mpoly_ts_append() release fence before A->length = newlength has the same problem if length remains just volatile slong.

Medium: H->cur is also shared and read outside the mutex in worker_loop(). The new release fence is placed before H->cur = next and next->producer = 1, but there is no matching acquire load of H->cur. This is mostly pre-existing design debt, but if the patch is presented as making the weak-memory threading protocol correct, H->cur is part of that protocol too.

Bottom Line

The experiment strongly suggests the fence commit fixes the demonstrated failure in practice on M3. But as a code-review conclusion: I would not present this as a fully sound atomics fix unless the shared publication fields are changed from volatile to actual C atomics, or the relevant accesses are put under the existing mutex.

[fredrik-johansson]

If the first patch reliably makes stress testing pass, then we can just apply that alone; no need to outright disable multithreading.

The codex suggestion probably makes sense, but that would indeed be a more substantial patch.

[oscarbenjamin]

I think it would be good to first make sure that the test suite does actually stress this in CI. I can open a PR with the test shown to check that.

Also I can look into this a bit more but from a few runs it seems that the first commit is slower than the second commit when running the test suite. It is possible that the fenced threading is slower than just not using threads but I guess that depends a lot on the coefficient ring, number of terms etc. Maybe the test suite is not testing particularly important cases though from a performance perspective.

[fredrik-johansson]

It's normal that the tests are slower with threading than without. The test code exercises the threaded versions even in regimes where they are slower than the serial versions.

OTOH, we should also verify that any changes to atomics handling don't slow down actual large threaded multiplications on x86-64.

[oscarbenjamin]

I adjusted the parameters of the test code in gh-2821:

• Choose big (more variables) rings 50% of the time instead of 10%.
• Always use at least 2 threads.
• Use 5x more iterations

With those I see the test fail on M3 about 95% of the time at multiplier 1. The last 5 CI runs in gh-2821 show that 4 out of 5 failed so about 80% failure rate on M1 with 0.5 multiplier.

Maybe there are more reliable ways to stress the contention but from the failure cases I've seen there isn't an obvious way to increase failure probability dramatically without a runtime increase that corresponds more or less to just increasing the number of iterations.