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

feat(v4.0.0): Q.EXP bit for bit with v3; lighter D2*C and (UQ/)

Browse Source 
Q.EXP is v3's q48_exp_approx: the Taylor series to 10 terms on |q|,
stopping below 50 ulp, 1.0 Q./ e^|q| for q < 0, e^0 = 1.0, and for
|q| >= 16.0 either 0 or Q max (v3 returned all ones, which is -ulp read
signed under D-8). x, term, sum, sign and term index live in a new
8-cell variable (QE). Checked bit for bit against v3's q48_exp_approx
(ported into the test) on 23 edge values and 3000 pseudo-random q in
(-16.0, 16.0), at both cell widths, optimised and ASan+UBSan; dropping a
term, moving the 50-ulp stop or skipping the reciprocal all fail.

Q.EXP calls Q./, which calls (UQ/), which calls D2*C, and at first it
left its caller no return entries. Lightened along the chain:
- D2*C takes cin in A and builds 2hi+m with `over + +`: no ROT, no SWAP.
- (UQ/) parks the quotient bit instead of using ROT, stores q without
  SWAP, and does its full trial subtraction (with borrow) in line
  instead of calling DNEGATE D+.
- Q.EXP keeps its term index in (QE) instead of a FOR count.

Headroom (data under args / return): (UQ/) 4/3 -> 4/4, Q./ 3/2 -> 3/3,
Q.EXP 3/2.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
rajames
2026-10-03 00:26:36 -04:00
co-authored by Claude Opus 5.5
parent 9e03bdab19
commit 29cf2688d2
2 changed files with 214 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files
v4/tests/test_foundation.c
+206 -25
View File
@@ -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);