fix(risc): DIV uses the netlist's 33-bit working value (#895)

jag_sim DIVIDE.NET (TOM and JERRY identical): an unsigned non-restoring
divider whose working value is the 32-bit dividend_hi register plus the
last_neg sign flop -- both adders carry a 33rd bit "to correctly
determine the sign, as the arguments are 32-bit unsigned values".  The
core's loop kept a 32-bit working value and took the sign from bit 31,
so every divisor above $80000000 gave garbage ($FCC00000/$FCC00000 =
$F9A9342A instead of 1, found in Xenowings' joypad reader), as did
16.16 divides whose quotient overflows.  For divisors up to $80000000
with an in-range quotient the result is bit-identical to before
(checked over 20M random cases).

- src/tom/risc_div.h: risc_divide(), the netlist algorithm, header-only
  like dsp_acc40.h.  Used by gpu_opcode_div, dsp_opcode_div and the
  dead pipelined DSP_div (DSPExecP is never defined), so no two paths
  can disagree.
- The DSP's divide-by-zero guard goes: the hardware has no special case,
  every step subtracts 0, giving $FFFFFFFF (unchanged) with the dividend
  in D_REMAIN (was 0).
- Tests: test_dsp_unit DIV cases incl. remainders (an even quotient ends
  on a negative step: remainder - divisor, as DIVIDE.NET documents);
  test_gpu_ops $FCC00000/$FCC00000.

Co-Authored-By: Claude Opus 5.5 <[email protected]>

src/jerry/dsp.c

@@ -17,6 +17,7 @@
 #include <compat/msvc.h>  /* snprintf shim for MSVC < 2015 (buildbot msvc05/10) */
 #include "dsp.h"
 #include "dsp_acc40.h"
+#include "risc_div.h"
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -3489,38 +3490,10 @@ INLINE static void dsp_opcode_abs(void)
 
 INLINE static void dsp_opcode_div(void)
 {
-	unsigned i;
-	uint32_t q;
-	uint32_t r;
-
-	/* Guard against divide-by-zero: JTRM says result is undefined but no
-	 * trap occurs.  Match the pipelined version's behavior. */
-	if (RM == 0)
-	{
-		RN = 0xFFFFFFFF;
-		dsp_remain = 0;
-		return;
-	}
-
-	/* Real algorithm, courtesy of SCPCD: NYAN! */
-	q = RN;
-	r = 0;
-
-	/* If 16.16 division, stuff top 16 bits of RN into remainder and put the
-	 * bottom 16 of RN in top 16 of quotient */
-	if (dsp_div_control & 0x01)
-		q <<= 16, r = RN >> 16;
-
-	for(i=0; i<32; i++)
-	{
-		uint32_t sign = r & 0x80000000;
-		r = (r << 1) | ((q >> 31) & 0x01);
-		r += (sign ? RM : -RM);
-		q = (q << 1) | (((~r) >> 31) & 0x01);
-	}
-
-	RN = q;
-	dsp_remain = r;
+	/* 33-bit non-restoring divide, netlist-exact (risc_div.h): correct
+	 * for divisors above $80000000; a zero divisor gives $FFFFFFFF with
+	 * the dividend in D_REMAIN, as the hardware does. */
+	RN = risc_divide(RN, RM, dsp_div_control & 0x01, &dsp_remain);
 }
 
 
@@ -3934,30 +3907,7 @@ INLINE static void DSP_cmpq(void)
 
 INLINE static void DSP_div(void)
 {
-	uint32_t _Rm = PRM, _Rn = PRN;
-
-	if (_Rm)
-	{
-		if (dsp_div_control & 1)
-		{
-			dsp_remain = (((uint64_t)_Rn) << 16) % _Rm;
-			if (dsp_remain & 0x80000000)
-				dsp_remain -= _Rm;
-			PRES = (((uint64_t)_Rn) << 16) / _Rm;
-		}
-		else
-		{
-			dsp_remain = _Rn % _Rm;
-			if (dsp_remain & 0x80000000)
-				dsp_remain -= _Rm;
-			PRES = PRN / _Rm;
-		}
-	}
-	else
-	{
-		PRES = 0xFFFFFFFF;
-		dsp_remain = 0;
-	}
+	PRES = risc_divide(PRN, PRM, dsp_div_control & 0x01, &dsp_remain);
 }
 
 INLINE static void DSP_imacn(void)

src/tom/gpu.c

@@ -33,6 +33,7 @@
 #include "jaguar.h"
 #include "m68000/m68kinterface.h"
 #include "tom.h"
+#include "risc_div.h"
 #include "event.h"
 #include "settings.h"
 #include "blit_memo.h"
@@ -3730,27 +3731,9 @@ INLINE static void gpu_opcode_abs(void)
 
 INLINE static void gpu_opcode_div(void)	// RN / RM
 {
-   unsigned i;
-   // Real algorithm, courtesy of SCPCD: NYAN!
-   uint32_t q = RN;
-   uint32_t r = 0;
-
-   // If 16.16 division, stuff top 16 bits of RN into remainder and put the
-   // bottom 16 of RN in top 16 of quotient
-   if (gpu_div_control & 0x01)
-      q <<= 16, r = RN >> 16;
-
-   for(i=0; i<32; i++)
-   {
-      uint32_t sign = r & 0x80000000;
-      r = (r << 1) | ((q >> 31) & 0x01);
-      r += (sign ? RM : -RM);
-      q = (q << 1) | (((~r) >> 31) & 0x01);
-   }
-
-   RN = q;
-   gpu_remain = r;
-
+   /* 33-bit non-restoring divide, netlist-exact (risc_div.h): correct
+    * for divisors above $80000000. */
+   RN = risc_divide(RN, RM, gpu_div_control & 0x01, &gpu_remain);
 }
 
 

src/tom/risc_div.h

@@ -0,0 +1,68 @@
+/*
+ * GPU/DSP DIV unit (TOM and JERRY share one design).
+ *
+ * jag_sim netlists/tom/DIVIDE.NET and netlists/jerry/DIVIDE.NET (logic
+ * identical; only the module name differs): an unsigned non-restoring
+ * divider, two quotient bits per tick for 16 ticks.  The working value
+ * is the 32-bit `dividend_hi` register plus the `last_neg` sign flop, so
+ * every add/subtract is 33 bits wide ("The additions / subtractions
+ * actually have to be performed at 33 bits of precision to correctly
+ * determine the sign, as the arguments are 32-bit unsigned values",
+ * DIVIDE.NET header; Adder1_33/Adder2_33 and Lnegt, lines 100-166).
+ * Step: shift the working value left one, shifting in the dividend's
+ * top bit (Dend1/Dend2), add the divisor if the last result was
+ * negative, else subtract it (Divr1/Divr2 = divisor XOR last_neg\),
+ * and shift the inverted 33rd bit into the quotient (Quosh/Quoti).
+ * div_offset (DIVCTRL bit 0) preloads the working value with the top 16
+ * dividend bits and the dividend with the low 16 bits shifted up
+ * (Divhit/Divein, lines 176-215).  The remainder register (D_REMAIN /
+ * G_REMAIN, Remrd) reads the low 32 bits of the final working value: the
+ * remainder, or the remainder minus the divisor when negative.
+ *
+ * The result is exact for every divisor, including those with bit 31
+ * set, where a 32-bit working value takes the sign from the wrong bit.
+ * A zero divisor is not special-cased in the hardware: every step
+ * subtracts 0, so an integer divide by zero returns $FFFFFFFF with the
+ * dividend in the remainder.
+ */
+
+#ifndef RISC_DIV_H
+#define RISC_DIV_H
+
+#include <stdint.h>
+
+#ifndef RISC_DIV_INLINE
+#if defined(_MSC_VER)
+#define RISC_DIV_INLINE static __inline
+#else
+#define RISC_DIV_INLINE static inline
+#endif
+#endif
+
+#define RISC_DIV_MASK33 UINT64_C(0x1FFFFFFFF)
+
+RISC_DIV_INLINE uint32_t risc_divide(uint32_t dividend, uint32_t divisor,
+                                     int offset, uint32_t *remain)
+{
+	uint64_t wk, sub;
+	uint32_t q;
+	uint32_t neg;
+	unsigned i;
+
+	wk  = offset ? (uint64_t)(dividend >> 16) : 0;
+	q   = offset ? (dividend << 16) : dividend;
+	sub = (UINT64_C(1) << 33) - (uint64_t)divisor;   /* -divisor mod 2^33 */
+	neg = 0;
+	for (i = 0; i < 32; i++)
+	{
+		wk  = ((wk << 1) | (q >> 31)) & RISC_DIV_MASK33;
+		q <<= 1;
+		wk  = (wk + (neg ? (uint64_t)divisor : sub)) & RISC_DIV_MASK33;
+		neg = (uint32_t)(wk >> 32) & 1;
+		q  |= neg ^ 1;
+	}
+	*remain = (uint32_t)wk;
+	return q;
+}
+
+#endif /* RISC_DIV_H */

test/test_dsp_unit.c

@@ -721,6 +721,49 @@ static void test_moveq(void)
    p_DSPWriteLong(DSP_CTRL_ADDR, 0, 6);
 }
 
+/* DIV is a 33-bit non-restoring divider (jag_sim DIVIDE.NET, issue #895;
+ * docs/jtrm-errata.md B26): exact for divisors above $80000000, where the
+ * old 32-bit loop took the sign from the wrong bit.  D_REMAIN reads the
+ * working value; a zero divisor yields $FFFFFFFF, remainder = dividend. */
+static void div_case(uint32_t dividend, uint32_t divisor, uint32_t ctrl,
+                     uint32_t want_q, uint32_t want_r)
+{
+   uint32_t off, q, r;
+
+   p_DSPReset();
+   for (off = 0; off < 0x2000; off += 2)
+      write_dsp_ram16(off, OP_NOP);
+   write_dsp_ram16(0x100, DSP_OP(21, 1, 2));        /* div r1,r2 */
+   p_dsp_reg_bank_0[1] = divisor;
+   p_dsp_reg_bank_0[2] = dividend;
+   p_DSPWriteLong(DSP_DIVCTRL_ADDR, ctrl, 6);
+   p_DSPWriteLong(DSP_PC_ADDR, 0xF1B100, 6);
+   p_DSPWriteLong(DSP_CTRL_ADDR, DSPGO, 6);
+   p_DSPExec(40);
+   p_DSPWriteLong(DSP_CTRL_ADDR, 0, 6);
+   q = p_dsp_reg_bank_0[2];
+   r = p_DSPReadLong(DSP_DIVCTRL_ADDR, 6);
+   p_DSPWriteLong(DSP_DIVCTRL_ADDR, 0, 6);
+   if (q == want_q && r == want_r)
+      PASS("div %08X/%08X (offset %u) = %08X rem %08X", dividend, divisor,
+           ctrl & 1, q, r);
+   else
+      FAIL("div %08X/%08X (offset %u) = %08X rem %08X, expected %08X rem %08X",
+           dividend, divisor, ctrl & 1, q, r, want_q, want_r);
+}
+
+static void test_div_33bit(void)
+{
+   printf("\n=== Test: DIV 33-bit working value (#895) ===\n");
+   /* An even quotient ends on a negative step: remainder - divisor. */
+   div_case(300, 3, 0, 100, 0xFFFFFFFDu);
+   div_case(0xFCC00000u, 0xFCC00000u, 0, 1, 0);          /* old: $F9A9342A */
+   div_case(0xFFFFFFFFu, 0x80000001u, 0, 1, 0x7FFFFFFEu);
+   div_case(0x7FFFFFFFu, 0x80000000u, 0, 0, 0x7FFFFFFFu - 0x80000000u);
+   div_case(0x00030000u, 0x00018000u, 1, 0x00020000u, 0xFFFE8000u); /* 3.0/1.5 (16.16) */
+   div_case(1234, 0, 0, 0xFFFFFFFFu, 1234);
+}
+
 /* ================================================================
  * Test 12: STORE/LOAD Round Trip
  * Write a value to DSP RAM via STORE, read it back via LOAD
@@ -1468,6 +1511,7 @@ int main(int argc, char *argv[])
    test_interrupt_return_address();
    test_dsp_execution();
    test_moveq();
+   test_div_33bit();
    test_store_load();
    test_dsp_self_stop();
    test_version_readonly();

test/test_gpu_ops.c

@@ -385,6 +385,16 @@ static void test_div(void)
    run(40);
    if (REG(1) == 14) PASS("100/7=14");
    else FAIL("100/7=%u (expected 14)", REG(1));
+
+   /* Divisor above $80000000: the 33-bit working value of DIVIDE.NET
+    * (#895).  The old 32-bit loop gave $F9A9342A here. */
+   prep();
+   gwmovei(0x100, 0xFCC00000u, 1);
+   gwmovei(0x106, 0xFCC00000u, 0);
+   gw16(0x10C, OP_DIV(0, 1));
+   run(40);
+   if (REG(1) == 1) PASS("$FCC00000/$FCC00000=1");
+   else FAIL("$FCC00000/$FCC00000=%08X (expected 1)", REG(1));
 }
 
 static void test_abs(void)