Fix M/MOD: hand-rolled 128/64 bit-serial division (FABRIC-3.md §XII.4)

mixed_math_word_m_slash_mod() had the identical bug class already fixed
in M* (commit 9a09949): it reconstructed the dividend from
`(dhigh << 32) | (dlow & 0xFFFFFFFF)`, the same wrong "32-bit halves of a
64-bit value" assumption. Latent for small inputs (fits in 32 bits, so the
reconstruction coincidentally worked), confirmed genuinely broken for a
true wide double -- feeding M*'s own correct 10^24 output into it gave a
quotient/remainder wrong by many orders of magnitude and the wrong sign.

Unlike M*, this one can't just switch to __int128 -- __int128's own `/`/`%`
need libgcc's __udivti3/__umodti3 for 128-bit division, unavailable in
this freestanding, -nostdlib build (the exact constraint
src/starkernel/arch/amd64/timer.c's own doc comment already flagged:
__int128 multiply/shift-by-constant/compare/subtract all compile clean,
only division doesn't). Fixed instead with a hand-rolled unsigned
128-by-64-bit bit-serial (restoring) long division -- 128 iterations of
shift-by-1/compare/subtract only, all in the safe set. Operates on
magnitudes via unsigned negation from 0 (well-defined even for the
extreme negative edge); signs reapplied afterward matching the same C99
truncating-toward-zero convention the previous, narrower implementation
already used, unchanged.

Verified: rebuilt all three architectures, confirmed clean link with no
__udivti3/__umodti3 undefined-symbol errors. Booted and tested all three,
identical results: 1000000000000 1000000000000 M* SWAP 1000000 M/MOD ->
1000000000000000000 remainder 0 (exact division, the case that was wrong
by orders of magnitude before); three sign-combination cases all correct
(-+, +-, --), confirming sign handling survived the magnitude-only
rewrite. T15 (original small-input case) unaffected. Full 24-case
exerciser reran clean on all three, no regressions.

All bugs found by the std79 exerciser campaign, including this one found
while fixing another, are now closed.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EXieurDfDSsDFdnSyusuWo

src/word_source/mixed_arithmetic_words.c

@@ -172,11 +172,35 @@ void mixed_math_word_m_star(VM *vm) {
  * @brief M/MOD ( d n -- rem quot )
  *
  * Divides the double-cell dividend @c d by the single-cell divisor @c n,
- * pushing the remainder deeper and the quotient on TOS. On 64-bit builds,
- * the double is reconstructed as a @c long long by shifting the high 32 bits
- * left; the UBSan-safe path uses an unsigned shift before casting to signed to
- * avoid undefined behaviour on negative values. On 32-bit builds, only handles
- * the case where the high cell is zero; sets @c vm->error = 1 otherwise.
+ * pushing the remainder deeper and the quotient on TOS. A "double" here is
+ * two full @c cell_t-width cells (128 bits total on a 64-bit build) -- same
+ * convention as M* (see that word's own doc comment for the full story).
+ * The previous implementation reconstructed the dividend as a plain 64-bit
+ * @c long long from `(dhigh << 32) | (dlow & 0xFFFFFFFF)`, the same wrong
+ * "32-bit halves of a 64-bit value" assumption M* had -- silently discarding
+ * the upper 32 bits of @c dlow and losing @c dhigh entirely once it wasn't
+ * just 0 or -1. Latent for any dividend whose true magnitude fit in 32 bits
+ * (never exercised otherwise until M*'s own fix started producing genuine
+ * wide doubles to feed it -- FABRIC-3.md §XII.4).
+ *
+ * Fixed via a hand-rolled unsigned 128-bit-by-64-bit bit-serial (restoring)
+ * long division, not @c __int128's own `/`/`%` operators: those need
+ * libgcc's `__udivti3`/`__umodti3`, unavailable in this freestanding,
+ * `-nostdlib` build (confirmed in `src/starkernel/arch/amd64/timer.c`'s own
+ * doc comment -- `__int128` multiply/shift-by-constant/compare/subtract all
+ * compile cleanly with zero undefined symbols on all three toolchains, only
+ * division doesn't). Shift-by-1, compare, and subtract on `unsigned
+ * __int128` are all in that safe set, so the classic 128-iteration
+ * shift-and-subtract algorithm works without any libgcc dependency.
+ * Operates on magnitudes (converted via unsigned negation from 0, which is
+ * well-defined even for the extreme negative edge, unlike signed negation);
+ * signs are reapplied afterward matching C99 truncating-toward-zero
+ * division (quotient sign = XOR of operand signs, remainder sign = the
+ * dividend's own sign) -- the same convention the previous, narrower
+ * implementation already used via plain C `/`/`%`, unchanged here.
+ *
+ * On 32-bit builds, only handles the case where the high cell is zero; sets
+ * @c vm->error = 1 otherwise.
  *
  * Stack effect: ( d_high d_low n -- remainder quotient )   TOS = quotient
  *
@@ -198,12 +222,28 @@ void mixed_math_word_m_slash_mod(VM *vm) {
     }
 
     if (sizeof(cell_t) == 8) {
-        /* UBSan fix: Perform shift on unsigned, then cast to signed (2025-12-09)
-         * Left-shifting negative signed values is undefined behavior */
-        long long dividend = (long long)(((unsigned long long) dhigh << 32) |
-                                         ((unsigned long long) dlow & 0xFFFFFFFFLL));
-        cell_t quotient = (cell_t)(dividend / n);
-        cell_t remainder = (cell_t)(dividend % n);
+        int dividend_neg = (dhigh < 0);
+        int divisor_neg = (n < 0);
+
+        unsigned __int128 udividend = ((unsigned __int128) (ucell_t) dhigh << 64) |
+                                       (unsigned __int128) (ucell_t) dlow;
+        if (dividend_neg) udividend = (unsigned __int128) 0 - udividend;
+
+        ucell_t udivisor = (ucell_t) n;
+        if (divisor_neg) udivisor = (ucell_t) 0 - udivisor;
+
+        unsigned __int128 uquot = 0, urem = 0;
+        for (int i = 127; i >= 0; i--) {
+            urem = (urem << 1) | ((udividend >> i) & 1);
+            if (urem >= (unsigned __int128) udivisor) {
+                urem -= (unsigned __int128) udivisor;
+                uquot |= ((unsigned __int128) 1 << i);
+            }
+        }
+
+        cell_t quotient = (dividend_neg ^ divisor_neg) ? -(cell_t) uquot : (cell_t) uquot;
+        cell_t remainder = dividend_neg ? -(cell_t) urem : (cell_t) urem;
+
         vm_push(vm, remainder); // remainder first (deeper)
         vm_push(vm, quotient); // quotient last (TOS)
     } else {