admin/LithosAnanake
1
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(v4.0.0): Q.SQRT bit for bit with v3; D-12 domain errors

Browse Source 
D-12 (ruled 2026-10-02): Q.SQRT of a negative value and Q.LOG of zero or
a negative value return 0 and set NODE-ERROR, as Q./ does on division by
zero.

Q.SQRT is v3's q48_sqrt_approx: Newton from x0 = q/2 + 0.25, up to 8
rounds of x' = (x + q/x) / 2, returning x once |x' - x| < 10 ulp;
sqrt(0) = 0 and sqrt(1.0) = 1.0. q, x and the rounds left live in a
5-cell variable (QR). Checked bit for bit against v3's q48_sqrt_approx
(ported into the test) on 19 edge values and 3000 pseudo-random q below
2^48, at both cell widths, optimised and ASan+UBSan. Above 2^48 v3's
q48_div saturates and v4 divides correctly, so there is no parity there.

Bug found and fixed on the way: the "< 10 ulp" test subtracted 10 from
the step's low cell as a signed number, so at 32-bit cells a step of
2^31 or more read as small and the iteration stopped early. It now
tests the low cell's top bit first, as Q.EXP does. The first test run
missed it because `make build/32-test_foundation.c` rebuilds nothing
(the Makefile's rules use absolute paths); the mutation runs, compiled
from source, exposed it. Only `make test` is used from here.

Headroom (data under arg / return): Q.SQRT 3/2.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
rajames
2026-10-03 09:07:30 -04:00
co-authored by Claude Opus 5.5
parent 29cf2688d2
commit a2d21d5cda
2 changed files with 145 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/v4.0.0/DECOMPOSITION.md
+7 -4
View File
@@ -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./ Q.EXP SEND RECV`. Words that clobber `B`:
`Q./ Q.EXP SEND RECV`.
document that clobber `A`: `@ ! +! -! 2@ 2! C@ C! UM* * UM/MOD Q.FROM-INT Q.TO-INT Q.* Q./ Q.EXP Q.SQRT SEND RECV`. Words that clobber `B`:
`Q./ Q.EXP Q.SQRT SEND RECV`.
**Return-stack words** (`>R R> R@ 2>R 2R> 2R@ I J UNLOOP` and the loop runtimes) are always IN.
@@ -173,6 +173,7 @@ definition below depends on one, it says so.
| **D-9** | Instruction word width on a 64-bit-cell host (ruled 2026-10-02). | **32 bits at every cell width.** Six 5-bit slots plus 2 spare bits, as in §1.2. On a 64-bit host the instruction word is the low 32 bits of the cell and the high half is ignored, so compiled code is identical at both widths. Only data and `@p` literals are a full cell wide. |
| **D-10** | Q48.16 width on a 64-bit-cell node (ruled 2026-10-02). | **Two cells at every cell width.** A Q value is a signed double on 32- and 64-bit nodes alike, so every Q word is the same double word at both widths. At 32-bit cells this is bit-for-bit v3's 64-bit Q. At 64-bit cells the low cell is v3's value, and where v3 wraps on Q48.16 overflow (a sum past Q max, `ABS` or `NEG` of Q min) the high cell carries the true result instead. |
| **D-11** | `Q./` semantics: division by zero, rounding, overflow (ruled 2026-10-02). | **Saturate and flag; round toward zero; saturate on overflow.** The quotient of `a * 2^16 / b` is rounded toward zero, like `SM/REM`. When it does not fit a Q value it is clamped to Q max or Q min by its sign. Division by zero returns Q max or Q min by the sign of the dividend (`0 / 0` gives 0) and sets the node's `NODE-ERROR` register (§7), which `VM-ERROR?` reads. v3 returned 0 on division by zero and saturated whenever the dividend was 2^48 or more, even when the quotient would have fit. |
| **D-12** | Q approximations outside their domain (ruled 2026-10-02). | **Return 0 and set `NODE-ERROR`.** `Q.SQRT` of a negative value and `Q.LOG` of zero or a negative value return 0 and set `NODE-ERROR` (§7), as `Q./` does on division by zero (D-11). v3 returned 0 for `ln(0)` without a flag and read negative arguments as large unsigned values. |
**Consequences of D-2 that every definition must respect.** The data stack holds 10 items and the
return stack 9, and every `call`, `FOR`, `DO` loop frame and `push` uses return-stack slots. Nesting
@@ -663,7 +664,9 @@ word. Per D-10 it occupies two cells on a 64-bit node too. Per D-8, v4 Q values
| `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.SQRT` `Q.SIN` `Q.COS` | CAP | Algorithms ported from `q48_words.c`; the hosted C versions are the golden model. |
| `Q.LOG` `Q.SIN` `Q.COS` | CAP | Algorithms ported from `q48_words.c`; the hosted C versions are the golden model. |
| `Q.SQRT` | CAP | v3's `q48_sqrt_approx`: Newton from `x0 = q/2 + 0.25`, up to 8 rounds of `x' = (x + q/x) / 2`, returning `x` as soon as `|x' - x| < 10` ulp; `sqrt(0) = 0`, `sqrt(1.0) = 1.0`. `q < 0` returns 0 and sets `NODE-ERROR` (D-12). `q`, `x` and the rounds left live in the variable `(QR)`. The halving is a logical double shift, `push a! 0 pop +* MAXHI and push drop a pop`. Executed on the golden model (2026-10-02): bit for bit v3's result for `0 <= q < 2^48`; above that v3's `q48_div` saturates and v4 divides correctly, so they part. Leaves its caller 3 data cells under its argument and 2 return entries. Clobbers `A` and `B`. |
| `(QR)` | CAP | Variable, 5 cells: `Q.SQRT`'s `q`, `x` and rounds left. |
| `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`. |
@@ -771,4 +774,4 @@ Addresses are assigned in the node memory map (D-4). Names only here.
| `GOV-*` | R/W | Governor parameters and Jacquard selector state (`L8-*`, decay rate). |
| `REC-ENABLE` | R/W | DoE recorder on/off (`HB-ON` / `HB-OFF`). |
| `PORT-STATUS` | R | Per-port ready flags, for non-blocking polls. |
| `NODE-ERROR` | R/W | Arithmetic error flag: set to −1 by `Q./` on division by zero (D-11), cleared by writing 0. `VM-ERROR?` reads it. |
| `NODE-ERROR` | R/W | Arithmetic error flag: set to −1 by `Q./` on division by zero (D-11) and by `Q.SQRT` / `Q.LOG` outside their domain (D-12), cleared by writing 0. `VM-ERROR?` reads it. |
v4/tests/test_foundation.c
+138 -2
View File
@@ -27,7 +27,9 @@
* 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.
* Q.EXP is checked bit for bit against v3's q48_exp_approx for |q| < 16.0. Every word is
* Q.EXP is checked bit for bit against v3's q48_exp_approx for |q| < 16.0.
* Q.SQRT is checked bit for bit against v3's q48_sqrt_approx for
* 0 <= q < 2^48, and for 0 with NODE-ERROR below zero (D-12). Every word is
* also probed for how much of the 10- and 9-deep circular stacks (D-2) it
* leaves to its caller.
*
@@ -52,6 +54,8 @@ static int failures = 0, checks = 0;
#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))
/* (QR), Q.SQRT's 5-cell variable: q, x (two cells each), rounds left. */
#define QR ((v4_cell)(V4_NODE_WORDS - 24u))
#define MAXU ((v4_ucell)~(v4_ucell)0)
#define FLAG(c) ((c) ? V4_ALL_ONES : (v4_cell)0)
@@ -65,7 +69,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_qexp;
w_dmax, w_dmin, w_qgt_doc, w_qgt, w_dltkeep, w_qstar, w_d2starc, w_uqdiv, w_qslash, w_qexp, w_qsqrt;
#define O(name) v4_asm_op(&as, V4_OP_##name)
#define LIT(v) v4_asm_lit(&as, (v4_cell)(v))
@@ -934,6 +938,86 @@ static void build(void)
#undef FETCH2
#undef STORE2
#undef HERE_
#undef FWD
/* : Q.SQRT ( q -- r ) 5.26, v3's q48_sqrt_approx
* Newton: x0 = q/2 + 0.25, then up to 8 rounds of x' = (x + q/x) / 2,
* returning x as soon as |x' - x| < 10 ulp. sqrt(0) = 0, sqrt(1.0) =
* 1.0; q < 0 gives 0 and sets NODE-ERROR (D-12). q, x and the rounds
* left live in (QR). D2/ below is the logical double shift right,
* push a! 0 pop +* MAXHI and push drop a pop
* (one +* step with S = 0, the sign bit cleared).
* dup -if POS drop drop drop NODE-ERROR b! -1 !b 0 0 ;
* POS: drop over over OR if ZERO drop
* over 65536 xor over OR if ONE drop
* over over QR a! SWAP !+ ! D2/ 16384 0 D+ QR 2 + a! SWAP !+ !
* 8 QR 4 + b! !b
* LOOP: QR a! @+ @ QR 2 + a! @+ @ Q./ QR 2 + a! @+ @ D+ D2/ x'
* over over QR 2 + a! @+ @ DNEGATE D+ -if L1 DNEGATE L1: x' |x'-x|
* if HI0 drop drop jump NEXT
* HI0: drop -if SM drop jump NEXT low cell >= 2^(N-1): not small
* SM: -10 + -if NEXT1 drop drop drop QR 2 + a! @+ @ ;
* NEXT1: drop
* NEXT: QR 2 + a! SWAP !+ !
* QR 4 + b! @b -1 + dup !b if DONE drop jump LOOP
* DONE: drop QR 2 + a! @+ @ ;
* ZERO: drop ; ONE: drop ;
* 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)
#define D2SLASH() do { O(PUSH); O(BANG_A); LIT(0); O(RPOP); O(MUL_STEP); \
LIT((v4_cell)(MAXU >> 1)); O(AND); O(PUSH); O(DROP); O(PUSH_A); O(RPOP); } while (0)
{
v4_asm_ref pos, zero, one, l1, hi0, next1, j_next, j_next2, done;
v4_cell loop, l_next;
w_qsqrt = v4_asm_label(&as);
O(DUP); pos = FWD(MINUS_IF);
O(DROP); O(DROP); O(DROP); LIT(NODE_ERROR); O(BANG_B); LIT(-1); O(STORE_B);
LIT(0); LIT(0); O(SEMI);
HERE_(pos);
O(DROP); O(OVER); O(OVER); CALL(w_or); zero = FWD(IF);
O(DROP); O(OVER); LIT(65536); O(XOR); O(OVER); CALL(w_or); one = FWD(IF);
O(DROP);
O(OVER); O(OVER); STORE2(QR);
D2SLASH(); LIT(16384); LIT(0); CALL(w_dplus); STORE2(QR + 2);
LIT(8); LIT(QR + 4); O(BANG_B); O(STORE_B);
loop = v4_asm_label(&as);
FETCH2(QR); FETCH2(QR + 2); CALL(w_qslash); FETCH2(QR + 2); CALL(w_dplus); D2SLASH();
O(OVER); O(OVER); FETCH2(QR + 2); CALL(w_dnegate); CALL(w_dplus);
l1 = FWD(MINUS_IF); CALL(w_dnegate); HERE_(l1);
hi0 = FWD(IF);
O(DROP); O(DROP); j_next = FWD(JUMP);
HERE_(hi0);
O(DROP);
{
v4_asm_ref sm = FWD(MINUS_IF);
O(DROP); j_next2 = FWD(JUMP);
HERE_(sm);
}
LIT(-10); O(ADD); next1 = FWD(MINUS_IF);
O(DROP); O(DROP); O(DROP); FETCH2(QR + 2); O(SEMI);
HERE_(next1);
O(DROP);
l_next = v4_asm_label(&as);
STORE2(QR + 2);
LIT(QR + 4); O(BANG_B); O(FETCH_B); LIT(-1); O(ADD); O(DUP); O(STORE_B);
done = FWD(IF);
O(DROP); v4_asm_branch(&as, V4_OP_JUMP, loop);
HERE_(done);
O(DROP); FETCH2(QR + 2); O(SEMI);
HERE_(zero); O(DROP); O(SEMI);
HERE_(one); O(DROP); O(SEMI);
v4_asm_resolve(&as, j_next, l_next);
v4_asm_resolve(&as, j_next2, l_next);
}
#undef D2SLASH
#undef FETCH2
#undef STORE2
#undef HERE_
#undef FWD
CHECK(v4_asm_ok(&as), "foundation words assemble");
@@ -1320,6 +1404,24 @@ static uint64_t v3_q48_exp(uint64_t q)
return result;
}
/* v3's q48_sqrt_approx, verbatim but for the names. */
static uint64_t v3_q48_sqrt(uint64_t q)
{
uint64_t x;
int iter;
if (q == 0) return 0;
if (q == 65536) return 65536;
x = (q >> 1) + 16384;
for (iter = 0; iter < 8; iter++) {
uint64_t q_div_x = v3_q48_div(q, x);
uint64_t x_next = (x + q_div_x) >> 1;
uint64_t delta = (x_next > x) ? (x_next - x) : (x - x_next);
if (delta < 10) break;
x = x_next;
}
return x;
}
/* 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
@@ -1975,6 +2077,34 @@ int main(void)
}
}
/* Q.SQRT against v3's q48_sqrt_approx for 0 <= q < 2^48 (where v3's
* q48_div does not saturate), bit for bit; q < 0 gives 0 and NODE-ERROR. */
{
static const int64_t qr[] = {
0, 1, 2, 9, 0x4000, 0x8000, 0xFFFF, 0x10000, 0x10001, 0x20000, 0x40000,
0x90000, 0x1000000, 0x7FFFFFFF, 0x100000000, ((int64_t)1 << 47) + 12345,
((int64_t)1 << 48) - 1, -1, -0x10000
};
uint64_t x = 0xC2B2AE3D27D4EB4Fu;
for (unsigned i = 0; i < sizeof qr / sizeof qr[0] + 3000u; i++) {
int64_t q;
v4_cell ql, qh, el, eh;
if (i < sizeof qr / sizeof qr[0]) q = qr[i];
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
q = (int64_t)(x >> (16 + (x & 31u))); /* 0 .. 2^48 */
if (q >= ((int64_t)1 << 48)) q >>= 1;
}
qcells((uint64_t)q, &ql, &qh);
if (q < 0) { el = 0; eh = 0; }
else qcells(v3_q48_sqrt((uint64_t)q), &el, &eh);
v4_node_store(&n, NODE_ERROR, 0);
CHECK(call(w_qsqrt, 2, ql, qh, 0) && left2(el, eh)
&& v4_node_load(&n, NODE_ERROR) == (q < 0 ? -1 : 0),
"Q.SQRT [%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. */
{
@@ -2082,6 +2212,12 @@ int main(void)
};
headroom("Q.EXP", w_qexp, 2, qexp_args, 6, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "Q.EXP leaves room");
static const v4_cell qsqrt_args[][4] = {
{ 0, 0, 0, 0 }, { 0x10000, 0, 0, 0 }, { 0x90000, 0, 0, 0 },
{ -1, 0x7FFF, 0, 0 }, { 5, -1, 0, 0 }
};
headroom("Q.SQRT", w_qsqrt, 2, qsqrt_args, 5, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "Q.SQRT 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);