diff --git a/docs/v4.0.0/DECOMPOSITION.md b/docs/v4.0.0/DECOMPOSITION.md index a858a7a5..74df7f25 100644 --- a/docs/v4.0.0/DECOMPOSITION.md +++ b/docs/v4.0.0/DECOMPOSITION.md @@ -147,8 +147,8 @@ IF body1 ELSE body2 THEN → if L1 drop body1 jump L2 times, as on the F18. `FOR ... UNEXT` is the same but the body must fit in one instruction word. **Register conventions.** `A` and `B` are caller-saved. A word that uses them says so. Words in this -document that clobber `A`: `@ ! +! -! 2@ 2! C@ C! UM* * UM/MOD Q.FROM-INT Q.TO-INT Q.* Q./ SEND RECV`. Words that clobber `B`: -`Q./ SEND RECV`. +document that clobber `A`: `@ ! +! -! 2@ 2! C@ C! UM* * UM/MOD Q.FROM-INT Q.TO-INT Q.* Q./ Q.EXP SEND RECV`. Words that clobber `B`: +`Q./ Q.EXP SEND RECV`. **Return-stack words** (`>R R> R@ 2>R 2R> 2R@ I J UNLOOP` and the loop runtimes) are always IN. @@ -658,12 +658,14 @@ word. Per D-10 it occupies two cells on a 64-bit node too. Per D-8, v4 Q values | --- | --- | --- | | `Q.+` `Q.-` | CAP | `D+`, `D-` — executed on the golden model (2026-10-02) against v3's `q48_add`/`q48_sub` (`uint64_t` wrapping): bit-for-bit at 32-bit cells; at 64-bit cells the low cell is v3's, per D-10. | | `Q.*` | CAP | `SWAP push over over UM* drop push push over pop -if L1 SWAP jump L2 L1: SWAP drop 0 L2: pop SWAP - push push over pop UM* pop + ROT pop dup push over -if L3 drop dup jump L4 L3: drop 0 L4: push UM* pop - D+ ROT pop UM* SWAP push 0 D+ pop push over pop SWAP Q.TO-INT push Q.TO-INT pop SWAP` (`SWAP` and `ROT` in line) — `floor(a*b / 2^16)`, signed (D-8), cut to two cells. Cells 0–2 of the product come from the unsigned cell products `a1*b1` (low cell), `a0*b1`, `a1*b0`, `a0*b0`, in that order, each input dropped after its last use and `b0` waiting on the return stack. Reading `a1` and `b1` as signed takes `(b1<0 ? a0 : 0)` and `(a1<0 ? b0 : 0)` off cell 2; each is folded in as soon as its operands are adjacent. The result is cells 0–2 shifted right 16 with `Q.TO-INT` twice. Executed on the golden model (2026-10-02) against an independent limb-by-limb reference, and against v3's `q48_mul` for non-negative operands. Leaves its caller 3 data cells under its arguments and 3 return entries. Clobbers `A`. | -| `Q./` | CAP | `over over OR if ZERO drop push over pop SWAP over xor (Q/) 4 + b! !b -if L1 DNEGATE L1: push push -if L2 DNEGATE L2: pop pop dup if CHK drop jump DIV CHK: drop over -131072 and if CHK2 drop jump DIV CHK2: drop push push dup pop dup push N-18 FOR 2* UNEXT over over xor -if SAME drop drop -if NOOV jump OV SAME: drop push inv pop + inv -if OV NOOV: drop pop pop jump DIV OV: drop drop drop pop pop drop drop (Q/) 4 + b! @b -if OVP drop 0 MSB ; OVP: drop -1 MAXHI ; DIV: (UQ/) (Q/) 4 + b! @b -if L3 drop DNEGATE ; L3: drop ; ZERO: drop drop drop NODE-ERROR b! -1 !b over over OR if Z0 drop -if ZP drop drop 0 MSB ; ZP: drop drop -1 MAXHI ; Z0: drop ;` — D-11: `a * 2^16 / b`, rounded toward zero; saturated to Q max / Q min on overflow; on division by zero, Q max / Q min by the dividend's sign (0 for `0 / 0`) and `NODE-ERROR` set. The sign of the result waits in `(Q/)`. Overflow needs no division: the quotient reaches `2^(2N-1)` exactly when `|a| * 2^16 >= |b| * 2^(2N-1)`, possible only for `|b| < 2^17`, where it is `a1 >= b0 << (N-17)` unsigned. `MSB` and `MAXHI` are the high cells of Q min and Q max. Executed on the golden model (2026-10-02) against an independent limb-by-limb reference (including `NODE-ERROR`), on width-native edge values, and against v3's `q48_div` for non-negative `a < 2^48`. Leaves its caller 3 data cells under its arguments and 2 return entries. Clobbers `A` and `B`. | -| `(UQ/)` | CAP | Internal to `Q./`: unsigned `floor(a * 2^16 / b)` for `b != 0` and no overflow. Restoring division, `2N+16` steps, on one shifting register — quotient (starting as `a`) above remainder — with the quotient and divisor in `(Q/)` so the stacks carry only the remainder, the quotient bit and the loop count. Each step's quotient bit enters at the next step's shift; with no overflow the bits leaving the quotient in the last 16 steps are zeros. Trial subtraction by in-line sign tests (`x - y` with `y` on top is `push inv pop + inv`). Executed on the golden model (2026-10-02). Clobbers `A` and `B`. | -| `D2*C` | CAP | `( lo hi cin -- lo' hi' cout )`: the double shifted left one bit, `cin` entering at the bottom and `cout` the bit leaving the top (both 0 or 1). `ROT dup -if L1 drop 1 jump L2 L1: drop 0 L2: push 2* + SWAP dup -if L3 drop 1 jump L4 L3: drop 0 L4: push 2* pop pop ROT + SWAP`, `ROT` and `SWAP` in line. Executed on the golden model (2026-10-02). | +| `Q./` | CAP | `over over OR if ZERO drop push over pop SWAP over xor (Q/) 4 + b! !b -if L1 DNEGATE L1: push push -if L2 DNEGATE L2: pop pop dup if CHK drop jump DIV CHK: drop over -131072 and if CHK2 drop jump DIV CHK2: drop push push dup pop dup push N-18 FOR 2* UNEXT over over xor -if SAME drop drop -if NOOV jump OV SAME: drop push inv pop + inv -if OV NOOV: drop pop pop jump DIV OV: drop drop drop pop pop drop drop (Q/) 4 + b! @b -if OVP drop 0 MSB ; OVP: drop -1 MAXHI ; DIV: (UQ/) (Q/) 4 + b! @b -if L3 drop DNEGATE ; L3: drop ; ZERO: drop drop drop NODE-ERROR b! -1 !b over over OR if Z0 drop -if ZP drop drop 0 MSB ; ZP: drop drop -1 MAXHI ; Z0: drop ;` — D-11: `a * 2^16 / b`, rounded toward zero; saturated to Q max / Q min on overflow; on division by zero, Q max / Q min by the dividend's sign (0 for `0 / 0`) and `NODE-ERROR` set. The sign of the result waits in `(Q/)`. Overflow needs no division: the quotient reaches `2^(2N-1)` exactly when `|a| * 2^16 >= |b| * 2^(2N-1)`, possible only for `|b| < 2^17`, where it is `a1 >= b0 << (N-17)` unsigned. `MSB` and `MAXHI` are the high cells of Q min and Q max. Executed on the golden model (2026-10-02) against an independent limb-by-limb reference (including `NODE-ERROR`), on width-native edge values, and against v3's `q48_div` for non-negative `a < 2^48`. Leaves its caller 3 data cells under its arguments and 3 return entries. Clobbers `A` and `B`. | +| `(UQ/)` | CAP | Internal to `Q./`: unsigned `floor(a * 2^16 / b)` for `b != 0` and no overflow. Restoring division, `2N+16` steps, on one shifting register — quotient (starting as `a`) above remainder — with the quotient and divisor in `(Q/)` so the stacks carry only the remainder, the quotient bit and the loop count. Each step's quotient bit enters at the next step's shift; with no overflow the bits leaving the quotient in the last 16 steps are zeros. Trial subtraction by in-line sign tests (`x - y` with `y` on top is `push inv pop + inv`), including the full double subtraction with its borrow, so the loop calls nothing but `D2*C`. Executed on the golden model (2026-10-02). Clobbers `A` and `B`. | +| `D2*C` | CAP | `( lo hi cin -- lo' hi' cout )`: the double shifted left one bit, `cin` entering at the bottom and `cout` the bit leaving the top (both 0 or 1). `a! dup -if L1 drop 1 jump L2 L1: drop 0 L2: push over -if L3 drop 1 jump L4 L3: drop 0 L4: over + + push 2* a + pop pop` — `cin` waits in `A`, and `2hi + m` is `over + +`, so it needs one return entry. Executed on the golden model (2026-10-02). Clobbers `A`. | | `(Q/)` | CAP | Variable, 5 cells: `Q./`'s quotient register, divisor and result sign. Like `BASE` and the hold buffer, it lives in node memory. | | `Q.ABS` `Q.NEG` | CAP | `DABS`, `DNEGATE` — executed on the golden model (2026-10-02) against v3's `q48_abs` and `0 - q`: bit-for-bit at 32-bit cells, including v3's wrap of Q min to itself; at 64-bit cells the low cell is v3's and the high cell is the true sign, per D-10. | -| `Q.LOG` `Q.EXP` `Q.SQRT` `Q.SIN` `Q.COS` | CAP | Algorithms ported from `q48_words.c`; the hosted C versions are the golden model. | +| `Q.LOG` `Q.SQRT` `Q.SIN` `Q.COS` | CAP | Algorithms ported from `q48_words.c`; the hosted C versions are the golden model. | +| `Q.EXP` | CAP | v3's `q48_exp_approx`: `1 + x + x^2/2! + …` to 10 terms on `x = |q|`, stopping once a term is below 50 ulp; `1.0 Q./ e^|q|` for `q < 0`; `e^0 = 1.0`; `|q| >= 16.0` gives 0 (`q < 0`) or **Q max** (v3 returned all ones, which reads as −ulp under D-8). `x`, the term, the sum, the sign and the term index live in the variable `(QE)`, so only `Q.*` and `Q./` arguments sit on the stacks. Executed on the golden model (2026-10-02): bit for bit v3's result for `|q| < 16.0` (at 64-bit cells as the same value, high cell 0). Leaves its caller 3 data cells under its argument and 2 return entries. Clobbers `A` and `B`. | +| `(QE)` | CAP | Variable, 8 cells: `Q.EXP`'s `x`, term, sum, sign and term index. | | `Q.FROM-INT` | CAP | `push 0 a! 0 pop 15 FOR +* UNEXT push drop a pop` — `n * 2^16` as a signed double. `+*` with `S = 0` never adds, so each step is an exact arithmetic right shift of `T:A`; starting from `T:A = n:0` (that is, `n * 2^N`) and shifting `N-16` bits leaves `n * 2^16`. The count is `N-17`: 15 at 32-bit cells, 47 at 64. Negative values are no longer clamped to 0. Executed on the golden model (2026-10-02): matches v3's `q48_from_u64` for `n >= 0` (at 64-bit cells in the low cell, where v3 wraps for `n >= 2^47`, per D-10). Clobbers `A`. | | `Q.TO-INT` | CAP | `push a! 0 pop 15 FOR +* UNEXT drop drop a` — the same `+*` shift, 16 bits at every cell width, leaving the low cell in `A`. Rounds toward minus infinity, as v3's arithmetic shift does (−1.5 gives −2). Executed on the golden model (2026-10-02): v3's `q48_to_u64` exactly at 64-bit cells, its low cell at 32. Clobbers `A`. | | `Q.1` `Q.0` `Q.SCALE` | IN | Double-cell constants. | diff --git a/v4/tests/test_foundation.c b/v4/tests/test_foundation.c index 5d9f94bb..3d229643 100644 --- a/v4/tests/test_foundation.c +++ b/v4/tests/test_foundation.c @@ -26,7 +26,8 @@ * reference, signed, and against v3's q48_mul for non-negative operands. * Q./ (D-11) and its core (UQ/) are checked against an independent * limb-by-limb long division, including NODE-ERROR, on width-native edge - * values, and against v3's q48_div where v4 keeps v3's meaning. Every word is + * values, and against v3's q48_div where v4 keeps v3's meaning. + * Q.EXP is checked bit for bit against v3's q48_exp_approx for |q| < 16.0. Every word is * also probed for how much of the 10- and 9-deep circular stacks (D-2) it * leaves to its caller. * @@ -49,6 +50,8 @@ static int failures = 0, checks = 0; #define NODE_ERROR ((v4_cell)(V4_NODE_WORDS - 2u)) /* (Q/), Q./'s 5-cell scratch variable: q0 q1 b0 b1 sign. */ #define QS ((v4_cell)(V4_NODE_WORDS - 8u)) +/* (QE), Q.EXP's 8-cell variable: x, term, sum (two cells each), sign, n. */ +#define QE ((v4_cell)(V4_NODE_WORDS - 16u)) #define MAXU ((v4_ucell)~(v4_ucell)0) #define FLAG(c) ((c) ? V4_ALL_ONES : (v4_cell)0) @@ -62,7 +65,7 @@ static v4_cell w_nip, w_swap, w_or, w_negate, w_rot, w_zless, w_zequal, w_ugreater, w_abs, w_s2d, w_dplus, w_dnegate, w_dabs, w_smrem, w_slashmod, w_mplus, w_dminus, w_d0equal, w_dequal, w_qfromint, w_qtoint, w_less, w_equal, w_dless, w_2swap, w_2over, - w_dmax, w_dmin, w_qgt_doc, w_qgt, w_dltkeep, w_qstar, w_d2starc, w_uqdiv, w_qslash; + w_dmax, w_dmin, w_qgt_doc, w_qgt, w_dltkeep, w_qstar, w_d2starc, w_uqdiv, w_qslash, w_qexp; #define O(name) v4_asm_op(&as, V4_OP_##name) #define LIT(v) v4_asm_lit(&as, (v4_cell)(v)) @@ -605,37 +608,34 @@ static void build(void) /* : D2*C ( lo hi cin -- lo' hi' cout ) 5.26, for Q./ * The double shifted left one bit, cin (0 or 1) entering at the bottom, - * cout (0 or 1) the bit that left the top. - * ROT dup -if L1 drop 1 jump L2 L1: drop 0 L2: push hi cin lo R: m - * 2* + SWAP lo' hi - * dup -if L3 drop 1 jump L4 L3: drop 0 L4: push R: m cout - * 2* pop pop ROT + SWAP ; lo' hi' cout - * ROT and SWAP in line. */ -#define SWAP_INLINE() do { O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0) -#define ROT_INLINE() do { O(PUSH); SWAP_INLINE(); O(RPOP); SWAP_INLINE(); } while (0) + * cout (0 or 1) the bit that left the top. cin waits in A. + * a! dup -if L1 drop 1 jump L2 L1: drop 0 L2: push lo hi R: cout + * over -if L3 drop 1 jump L4 L3: drop 0 L4: lo hi m m: lo's top bit + * over + + lo 2hi+m + * push 2* a + pop pop ; lo' hi' cout + * Clobbers A. */ { v4_asm_ref l1, l2, l3, l4; w_d2starc = v4_asm_label(&as); - ROT_INLINE(); O(DUP); + O(BANG_A); O(DUP); l1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF); O(DROP); LIT(1); l2 = v4_asm_branch_fwd(&as, V4_OP_JUMP); v4_asm_resolve(&as, l1, v4_asm_label(&as)); O(DROP); LIT(0); v4_asm_resolve(&as, l2, v4_asm_label(&as)); - O(PUSH); O(TWO_STAR); O(ADD); SWAP_INLINE(); - O(DUP); + O(PUSH); O(OVER); l3 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF); O(DROP); LIT(1); l4 = v4_asm_branch_fwd(&as, V4_OP_JUMP); v4_asm_resolve(&as, l3, v4_asm_label(&as)); O(DROP); LIT(0); v4_asm_resolve(&as, l4, v4_asm_label(&as)); - O(PUSH); O(TWO_STAR); O(RPOP); O(RPOP); ROT_INLINE(); O(ADD); SWAP_INLINE(); + O(OVER); O(ADD); O(ADD); + O(PUSH); O(TWO_STAR); O(PUSH_A); O(ADD); O(RPOP); O(RPOP); O(SEMI); } -#undef ROT_INLINE -#undef SWAP_INLINE + /* : (UQ/) ( a0 a1 b0 b1 -- q0 q1 ) 5.26, core of Q./ * Unsigned floor(a * 2^16 / b), for b != 0 and a quotient below @@ -648,8 +648,8 @@ static void build(void) * stacks hold only r, the quotient bit and the loop count. * QS 2 + a! SWAP !+ ! QS a! SWAP !+ ! 0 0 0 r0 r1 qb * 2N+15 FOR - * QS a! @+ @ ROT D2*C r0 r1 q0 q1 c - * push SWAP QS a! !+ ! pop D2*C r0 r1 t + * push QS a! @+ @ pop D2*C r0 r1 q0 q1 c + * push QS 1 + a! ! QS a! ! pop D2*C r0 r1 t * if CMP drop jump TAKE * CMP: drop QS 3 + b! @b r0 r1 b1 * over over xor if EQH -if SAMEH @@ -663,11 +663,15 @@ static void build(void) * -if T4 drop jump NOTAKEL T4: drop jump TAKEL * TAKEL: drop QS 2 + b! @b push inv pop + inv 0 1 jump END * NOTAKEL: 0 jump END - * TAKE: QS 2 + a! @+ @ DNEGATE D+ 1 jump END + * TAKE: push QS 2 + b! @b over over push inv pop + inv push r0 b0 R: r1 d0 + * over over xor -if SB drop NIP -if NB drop 1 jump BD NB: drop 0 jump BD + * SB: drop push inv pop + inv -if NB2 drop 1 jump BD NB2: drop 0 + * BD: pop pop QS 3 + b! @b push inv pop + inv borrow d0 r1-b1 + * push SWAP pop SWAP - 1 jump END d0 r1-b1-borrow 1 * NOTAKE: 0 * END: * NEXT - * QS a! @+ @ ROT D2*C drop push push drop drop pop pop ; + * push QS a! @+ @ pop D2*C drop push push drop drop pop pop ; * x - y with y on top is push inv pop + inv. SWAP, ROT in line. * Clobbers A and B. */ #define SWAP_INLINE() do { O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0) @@ -687,8 +691,8 @@ static void build(void) LIT(0); LIT(0); LIT(0); LIT(2 * V4_CELL_BITS + 15); O(PUSH); loop = v4_asm_label(&as); - LIT(QS); O(BANG_A); O(FETCH_INC); O(FETCH_A); ROT_INLINE(); CALL(w_d2starc); - O(PUSH); SWAP_INLINE(); LIT(QS); O(BANG_A); O(STORE_INC); O(STORE_A); O(RPOP); + O(PUSH); LIT(QS); O(BANG_A); O(FETCH_INC); O(FETCH_A); O(RPOP); CALL(w_d2starc); + O(PUSH); LIT(QS + 1); O(BANG_A); O(STORE_A); LIT(QS); O(BANG_A); O(STORE_A); O(RPOP); CALL(w_d2starc); cmp = FWD(IF); O(DROP); j_take[nt++] = FWD(JUMP); @@ -719,13 +723,28 @@ static void build(void) l_notakel = v4_asm_label(&as); /* NOTAKEL */ LIT(0); j_end[ne++] = FWD(JUMP); l_take = v4_asm_label(&as); /* TAKE */ - LIT(QS + 2); O(BANG_A); O(FETCH_INC); O(FETCH_A); - CALL(w_dnegate); CALL(w_dplus); LIT(1); j_end[ne++] = FWD(JUMP); + { + v4_asm_ref sb, nb, bd, nb2, bd2; + O(PUSH); LIT(QS + 2); O(BANG_B); O(FETCH_B); + O(OVER); O(OVER); SUB_INLINE(); O(PUSH); + O(OVER); O(OVER); O(XOR); sb = FWD(MINUS_IF); + O(DROP); O(PUSH); O(DROP); O(RPOP); nb = FWD(MINUS_IF); + O(DROP); LIT(1); bd = FWD(JUMP); + HERE_(nb); O(DROP); LIT(0); bd2 = FWD(JUMP); + HERE_(sb); + O(DROP); SUB_INLINE(); nb2 = FWD(MINUS_IF); + O(DROP); LIT(1); j_end[ne] = FWD(JUMP); + HERE_(nb2); O(DROP); LIT(0); + HERE_(bd); HERE_(bd2); v4_asm_resolve(&as, j_end[ne], v4_asm_label(&as)); + O(RPOP); O(RPOP); LIT(QS + 3); O(BANG_B); O(FETCH_B); SUB_INLINE(); + O(PUSH); SWAP_INLINE(); O(RPOP); SWAP_INLINE(); CALL(w_minus); + LIT(1); j_end[ne++] = FWD(JUMP); + } l_notake = v4_asm_label(&as); /* NOTAKE */ LIT(0); l_end = v4_asm_label(&as); /* END */ v4_asm_branch(&as, V4_OP_NEXT, loop); - LIT(QS); O(BANG_A); O(FETCH_INC); O(FETCH_A); ROT_INLINE(); CALL(w_d2starc); O(DROP); + O(PUSH); LIT(QS); O(BANG_A); O(FETCH_INC); O(FETCH_A); O(RPOP); CALL(w_d2starc); O(DROP); O(PUSH); O(PUSH); O(DROP); O(DROP); O(RPOP); O(RPOP); O(SEMI); for (k = 0; k < nt; k++) v4_asm_resolve(&as, j_take[k], l_take); @@ -822,6 +841,101 @@ static void build(void) #undef SUB_INLINE #undef SWAP_INLINE + /* : Q.EXP ( q -- e ) 5.26, v3's q48_exp_approx + * e^q by the Taylor series 1 + x + x^2/2! + ... to 10 terms, stopping + * once a term is below 50 ulp, on x = |q|; for q < 0 the result is + * 1.0 Q./ e^|q|. e^0 = 1.0; |q| >= 16.0 gives 0 (q < 0) or Q max (v3 + * returned all ones, which is -ulp read signed, D-8). x, the term and + * the sum live in (QE), so only Q.* and Q./ arguments sit on the stack. + * over over OR if ONE drop + * dup QE 6 + b! !b -if L1 DNEGATE L1: x (QE)+6: sign + * over over -1048576 -1 D+ push drop pop -if BIG drop + * over over QE 2 + a! SWAP !+ ! over over QE a! SWAP !+ ! + * 65536 0 D+ QE 4 + a! SWAP !+ ! term = x, sum = 1 + x + * 2 QE 7 + b! !b (QE)+7: n + * LOOP: + * QE 2 + a! @+ @ QE a! @+ @ Q.* + * QE 7 + b! @b Q.FROM-INT Q./ term*x / n + * over over QE 2 + a! SWAP !+ ! + * QE 4 + a! @+ @ D+ QE 4 + a! SWAP !+ ! sum += term + * QE 2 + a! @+ @ if HI0 drop drop jump CONT + * HI0: drop -if SM drop jump CONT + * SM: -50 + -if CONT1 drop jump DONE term < 50: stop + * CONT1: drop + * CONT: QE 7 + b! @b 1 + dup !b -11 + -if DONE1 drop jump LOOP + * DONE1: drop after n = 10 + * DONE: QE 4 + a! @+ @ QE 6 + b! @b -if POS drop 65536 0 2SWAP Q./ ; + * POS: drop ; + * BIG: drop drop drop QE 6 + b! @b -if BIGP drop 0 0 ; BIGP: drop -1 MAXHI ; + * ONE: drop drop drop 65536 0 ; + * Clobbers A and B. */ +#define FWD(op) v4_asm_branch_fwd(&as, V4_OP_##op) +#define HERE_(r) v4_asm_resolve(&as, (r), v4_asm_label(&as)) +#define STORE2(addr) do { LIT(addr); O(BANG_A); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); \ + O(STORE_INC); O(STORE_A); } while (0) +#define FETCH2(addr) do { LIT(addr); O(BANG_A); O(FETCH_INC); O(FETCH_A); } while (0) + { + v4_asm_ref one, l1, big, hi0, sm, cont1, j_cont1, j_cont2, j_done, pos, bigp; + v4_cell loop, l_cont, l_done; + + w_qexp = v4_asm_label(&as); + O(OVER); O(OVER); CALL(w_or); one = FWD(IF); + O(DROP); + O(DUP); LIT(QE + 6); O(BANG_B); O(STORE_B); + l1 = FWD(MINUS_IF); CALL(w_dnegate); HERE_(l1); + O(OVER); O(OVER); LIT(-1048576); LIT(-1); CALL(w_dplus); + O(PUSH); O(DROP); O(RPOP); + big = FWD(MINUS_IF); + O(DROP); + O(OVER); O(OVER); STORE2(QE + 2); + O(OVER); O(OVER); STORE2(QE); + LIT(65536); LIT(0); CALL(w_dplus); STORE2(QE + 4); + LIT(2); LIT(QE + 7); O(BANG_B); O(STORE_B); + loop = v4_asm_label(&as); + FETCH2(QE + 2); FETCH2(QE); CALL(w_qstar); + LIT(QE + 7); O(BANG_B); O(FETCH_B); CALL(w_qfromint); CALL(w_qslash); + O(OVER); O(OVER); STORE2(QE + 2); + FETCH2(QE + 4); CALL(w_dplus); STORE2(QE + 4); + FETCH2(QE + 2); hi0 = FWD(IF); + O(DROP); O(DROP); j_cont1 = FWD(JUMP); + HERE_(hi0); + O(DROP); sm = FWD(MINUS_IF); + O(DROP); j_cont2 = FWD(JUMP); + HERE_(sm); + LIT(-50); O(ADD); cont1 = FWD(MINUS_IF); + O(DROP); j_done = FWD(JUMP); + HERE_(cont1); + O(DROP); + l_cont = v4_asm_label(&as); + LIT(QE + 7); O(BANG_B); O(FETCH_B); LIT(1); O(ADD); O(DUP); O(STORE_B); + LIT(-11); O(ADD); + { + v4_asm_ref fin = FWD(MINUS_IF); + O(DROP); v4_asm_branch(&as, V4_OP_JUMP, loop); + HERE_(fin); O(DROP); + } + l_done = v4_asm_label(&as); + FETCH2(QE + 4); LIT(QE + 6); O(BANG_B); O(FETCH_B); + pos = FWD(MINUS_IF); + O(DROP); LIT(65536); LIT(0); CALL(w_2swap); CALL(w_qslash); O(SEMI); + HERE_(pos); O(DROP); O(SEMI); + HERE_(big); + O(DROP); O(DROP); O(DROP); LIT(QE + 6); O(BANG_B); O(FETCH_B); + bigp = FWD(MINUS_IF); + O(DROP); LIT(0); LIT(0); O(SEMI); + HERE_(bigp); O(DROP); LIT(-1); LIT((v4_cell)(MAXU >> 1)); O(SEMI); + HERE_(one); + O(DROP); O(DROP); O(DROP); LIT(65536); LIT(0); O(SEMI); + + v4_asm_resolve(&as, j_cont1, l_cont); + v4_asm_resolve(&as, j_cont2, l_cont); + v4_asm_resolve(&as, j_done, l_done); + } +#undef FETCH2 +#undef STORE2 +#undef HERE_ +#undef FWD + CHECK(v4_asm_ok(&as), "foundation words assemble"); } @@ -1172,6 +1286,40 @@ static void qdiv_ref(v4_ucell al, v4_ucell ah, v4_ucell bl, v4_ucell bh, if (neg) dneg(*rl, *rh, rl, rh); } +/* v3's q48_div and q48_exp_approx (v3/src/word_source/q48_16_words.c), + * verbatim but for the names, as the parity reference for Q.EXP. They use + * v3_q48_mul above, which is v3's __int128 result. */ +static uint64_t v3_q48_div(uint64_t a, uint64_t b) +{ + if (b == 0) return 0; + if (a > 0x0000FFFFFFFFFFFFULL) return 0xFFFFFFFFFFFFFFFFULL; + return (a << 16) / b; +} +static uint64_t v3_q48_exp(uint64_t q) +{ + int64_t q_signed; + int is_negative, nn; + uint64_t x, result, term; + + if (q == 0) return 65536; + q_signed = (int64_t)q; + if (q_signed >= 1048576) return 0xFFFFFFFFFFFFFFFFULL; + if (q_signed <= -1048576) return 0; + is_negative = (q_signed < 0); + x = is_negative ? (uint64_t)(-q_signed) : q; + result = 65536; + term = x; + result = result + term; + for (nn = 2; nn <= 10; nn++) { + term = v3_q48_mul(term, x); + term = v3_q48_div(term, (uint64_t)nn << 16); + result = result + term; + if (term < 50) break; + } + if (is_negative) result = v3_q48_div(65536, result); + return result; +} + /* Stack headroom. Run `word` with `dfill` marked cells under the canary and * `rfill` marked cells under its return address, and report whether the * results, the canary and every marked cell come back intact. The F18 stacks @@ -1800,6 +1948,33 @@ int main(void) } } + /* Q.EXP against v3's q48_exp_approx for |q| < 16.0 (bit for bit at + * 32-bit cells; the same value, high cell 0, at 64), and the limits. */ + { + static const int64_t qx[] = { + 0, 1, -1, 0x8000, -0x8000, 0x10000, -0x10000, 0x18000, 0x20000, -0x20000, + 0x2B7E1, 0x50000, -0x50000, 0xA0000, -0xA0000, 0xFFFFF, -0xFFFFF, + 0x100000, -0x100000, 0x100001, -0x100001, 0x7FFFFFFF, -0x7FFFFFFF + }; + uint64_t x = 0x9FB21C651E98DF25u; + for (unsigned i = 0; i < sizeof qx / sizeof qx[0] + 3000u; i++) { + int64_t q; + uint64_t uq; + v4_cell ql, qh, el, eh; + if (i < sizeof qx / sizeof qx[0]) q = qx[i]; + else { + x ^= x << 13; x ^= x >> 7; x ^= x << 17; + q = (int64_t)(x % 0x200000u) - 0x100000; /* (-16.0, 16.0) */ + } + uq = (uint64_t)q; + qcells(uq, &ql, &qh); + if (q >= 0x100000) { el = V4_ALL_ONES; eh = (v4_cell)(MAXU >> 1); } + else if (q <= -0x100000) { el = 0; eh = 0; } + else qcells(v3_q48_exp(uq), &el, &eh); + CHECK(call(w_qexp, 2, ql, qh, 0) && left2(el, eh), "Q.EXP [%u] q=%lld", i, (long long)q); + } + } + /* Q./ at the overflow boundary: |a| * 2^16 against |b| * 2^(2N-1), for * divisors around 2^16 and 2^17, every sign. */ { @@ -1901,6 +2076,12 @@ int main(void) headroom("(UQ/)", w_uqdiv, 4, qdiv_args, 2, 2, &dh, &rh); headroom("Q./", w_qslash, 4, qdiv_args, 7, 2, &dh, &rh); CHECK(dh >= 2 && rh >= 2, "Q./ leaves room"); + static const v4_cell qexp_args[][4] = { + { 0, 0, 0, 0 }, { 0x10000, 0, 0, 0 }, { -0x10000, -1, 0, 0 }, + { 0x50000, 0, 0, 0 }, { 0x200000, 0, 0, 0 }, { -0x200000, -1, 0, 0 } + }; + headroom("Q.EXP", w_qexp, 2, qexp_args, 6, 2, &dh, &rh); + CHECK(dh >= 2 && rh >= 2, "Q.EXP leaves room"); headroom("SM/REM", w_smrem, 3, smrem_args, 5, 2, &dh, &rh); CHECK(rh >= 1, "SM/REM leaves room for /MOD's return address"); headroom("/MOD", w_slashmod, 2, slashmod_args, 4, 2, &dh, &rh);