admin/LithosAnanake
1
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
a041b401ea34fbf528a2ce9925ccbc4f89a8dc49
LithosAnanake/v4/tests/test_foundation.c
T
rajamesandClaude Opus 5.5 481d484e93 feat(v4.0.0): the mixed and double leftovers 
M- M* M/MOD MOD */ */MOD, D0< D2* D2/ 2ROT, 2DROP and 2>R 2R@ 2R>,
beside UM* and SM/REM in test_foundation.c.  Executed on the golden
model at both cell widths against C and results recorded from the v3
binary.

M- widens n before negating it, so the most negative n is right.
*/MOD goes through a full double product.  D2/ is one +* step.
M- and M/MOD take the double in the standard order ( lo hi ), as M+
does; v3 took its low cell on top.

The foundation test's node is now full: 958 of the 960 words below its
variables.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-03 18:37:53 -04:00

3135 lines
142 KiB
C
Raw Blame History

/* test_foundation.c -- the first DECOMPOSITION.md definitions, executed.
*
* Section 4 says of its colon definitions that each "has been traced by hand,
* but none has been executed". This file assembles the helpers, the sign and
* zero tests, U< and UM* exactly as written there (plus 2DUP and - from
* section 5, which U< needs) and runs them on the golden model against the C
* operation each one stands for.
*
* UM* is the full-range version that section 4 gives under D-3, checked
* against the reference v4_umul over every pair of the edge vectors and 20000
* pseudo-random pairs at each cell width. UM/MOD is section 4's call-free
* version, checked against q*d + r = uhi:ulo, r < d over the edge vectors and
* 20000 pseudo-random cases with uhi < ud. SM/REM and DNEGATE are the
* call-free versions in sections 4 and 5.7; SM/REM is checked on dividends
* built as q*n + r. D+ is the call-free version in 5.7, checked against C
* over every quadruple of the edge vectors and pseudo-random pairs. /MOD,
* U>, ABS, S>D, DABS, M+, D-, D0= and D= are as written in section 5 and
* checked against C. Q.+, Q.-, Q.ABS and Q.NEG (section 5.26: the D+, D-,
* DABS and DNEGATE words) are checked against v3's 64-bit q48_add, q48_sub,
* q48_abs and 0 - q, with Q values two cells wide at both widths (D-10).
* Q.FROM-INT and Q.TO-INT are checked against n * 2^16 and q >> 16 in C and
* against v3's q48_from_u64 and q48_to_u64 where v4 keeps v3's meaning.
* <, =, D<, 2OVER, the call-free (D<), DMAX, DMIN and 2SWAP, and the Q
* comparisons on them, are checked against C, signed (D-8), and against v3
* where v3's unsigned Q comparisons agree. Q.* is checked against an independent limb-by-limb
* 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.
* 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).
* Q.LOG is checked bit for bit against v3's q48_log_approx for x > 0, and
* for 0 with NODE-ERROR at or below zero.
* (Q.REDUCE), Q.SIN and Q.COS are checked bit for bit against v3's
* q48_reduce_angle, q48_sin_approx and q48_cos_approx on angles of every
* size and sign. The in-line constants Q.1, Q.0 and Q.SCALE are checked
* against v3's values and in use. Every word is
* also probed for how much of the 10- and 9-deep circular stacks (D-2) it
* leaves to its caller.
*
* Every call is made with a canary under the arguments, and the canary must
* still be directly under the results afterwards: a definition that leaves
* the right answer but an unbalanced stack is wrong.
*/
#include "v4/asm.h"
#include "v4/testcode.h"
#include "v4/umul.h"
#include <stdint.h>
#include <stdio.h>
static int failures = 0, checks = 0;
#define CHECK(c,...) do{checks++; if(!(c)){failures++; printf("FAIL %s:%d: ",__FILE__,__LINE__); printf(__VA_ARGS__); printf("\n");}}while(0)
#define CANARY ((v4_cell)0x0C0FFEE5)
/* NODE-ERROR (section 7); D-4 leaves its address open, so the test puts it
* in the word next to the halt address. */
#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))
/* (QR), Q.SQRT's 5-cell variable: q, x (two cells each), rounds left. */
#define QR ((v4_cell)(V4_NODE_WORDS - 24u))
/* (QL), Q.LOG's 7-cell variable: m, y, k, rounds left, term, sum, n*1.0. */
#define QL ((v4_cell)(V4_NODE_WORDS - 32u))
/* (QT), Q.SIN / Q.COS's 6-cell variable: sign, x^2, term, sum, n, subtract. */
#define QT ((v4_cell)(V4_NODE_WORDS - 40u))
/* Scratch words for the byte-access tests. */
#define BYTES ((v4_cell)(V4_NODE_WORDS - 48u))
#define MAXU ((v4_ucell)~(v4_ucell)0)
#define FLAG(c) ((c) ? V4_ALL_ONES : (v4_cell)0)
static v4_node n;
static v4_exec_state es;
static v4_heat h;
static v4_asm as;
static v4_cell w_nip, w_swap, w_or, w_negate, w_rot, w_zless, w_zequal,
w_2dup, w_minus, w_uless, w_umstar, w_ummod,
w_ugreater, w_abs, w_s2d, w_dplus, w_dnegate, w_dabs, w_smrem,
w_slashmod, w_star, w_slash, w_mminus, w_mstar, w_mslashmod, w_mod, w_starslashmod,
w_starslash, w_d0less, w_d2star, w_d2slash, w_2rot, w_2drop, t_2r, t_2rorder, 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_qsqrt, w_qlog,
w_qreduce, w_qsin, l_trig, w_qcos,
w_q1, w_q0, w_qscale, w_q1_times, w_q0_plus, w_q1_toint, w_one_fromint,
w_lshift, w_rshift, w_cfetch, w_cstore;
#define O(name) v4_asm_op(&as, V4_OP_##name)
#define LIT(v) v4_asm_lit(&as, (v4_cell)(v))
#define CALL(w) v4_asm_branch(&as, V4_OP_CALL, (w))
/* Section 2's capsule IF ... THEN: the flag is dropped on both paths. */
static v4_asm_ref if_(void)
{
v4_asm_ref r = v4_asm_branch_fwd(&as, V4_OP_IF);
v4_asm_op(&as, V4_OP_DROP);
return r;
}
static void then_(v4_asm_ref r)
{
v4_asm_ref j = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, r, v4_asm_label(&as));
v4_asm_op(&as, V4_OP_DROP);
v4_asm_resolve(&as, j, v4_asm_label(&as));
}
static void build(void)
{
v4_asm_ref ref;
v4_node_reset(&n);
v4_asm_begin(&as, &n, 16);
/* : NIP push drop pop ; */
w_nip = v4_asm_label(&as);
O(PUSH); O(DROP); O(RPOP); O(SEMI);
/* : SWAP over push push drop pop pop ; */
w_swap = v4_asm_label(&as);
O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); O(SEMI);
/* : OR over inv and xor ; */
w_or = v4_asm_label(&as);
O(OVER); O(INV); O(AND); O(XOR); O(SEMI);
/* : NEGATE inv 1 + ; */
w_negate = v4_asm_label(&as);
O(INV); LIT(1); O(ADD); O(SEMI);
/* : ROT push SWAP pop SWAP ; */
w_rot = v4_asm_label(&as);
O(PUSH); CALL(w_swap); O(RPOP); CALL(w_swap); O(SEMI);
/* : 0< -if L1 drop -1 ; L1: drop 0 ; */
w_zless = v4_asm_label(&as);
ref = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(-1); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); LIT(0); O(SEMI);
/* : 0= if L1 drop 0 ; L1: drop -1 ; */
w_zequal = v4_asm_label(&as);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); LIT(0); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); LIT(-1); O(SEMI);
/* : 2DUP over over ; (section 5.1) */
w_2dup = v4_asm_label(&as);
O(OVER); O(OVER); O(SEMI);
/* : - NEGATE + ; (section 5.4), NEGATE in line */
w_minus = v4_asm_label(&as);
O(INV); LIT(1); O(ADD); O(ADD); O(SEMI);
/* : U< 2DUP xor 0< IF NIP 0< ELSE - 0< THEN ;
* with section 2's IF: the flag is dropped on both arms. */
w_uless = v4_asm_label(&as);
CALL(w_2dup); O(XOR); CALL(w_zless);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); CALL(w_nip); CALL(w_zless); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); CALL(w_minus); CALL(w_zless); O(SEMI);
/* : UM* ( u1 u2 -- ulo uhi ) section 4, D-3
* over 1 and inv 1 + over 2/ and u1 u2 t0 t0 = u2 2/ if u1 odd
* over a! 31 push R: loop count (FOR's push, early)
* push over 2/ pop u1 u2 s t0 A: u2, s = u1 2/
* L: +* unext u1 u2 s hi A: lo
* push drop pop u1 u2 hi
* 2* a -if L0 drop 1 + jump L1 L0: drop L1: push R: 2hi + top bit of lo
* over -if L2 drop dup jump L3 L2: drop 0 L3: push R: + (u1<0 ? u2 : 0)
* dup -if L4 drop over 1 and jump L5 L4: drop 0 L5: u1 u2 (u2<0 ? u1&1 : 0)
* pop + pop + push u1 u2 R: hi'
* and 1 and a 2* + pop ; lo' hi'
* u1 and u2 wait on the data stack under the loop (+* touches only T, S
* and A), so the corrections are made afterwards and the return stack
* holds at most two temporaries. 31 is the cell width less one, pushed
* before s is made so the data stack never holds more than four; the loop
* body is the start of its own word so that unext restarts it. */
{
v4_asm_ref l0, l1, l2, l3, l4, l5;
w_umstar = v4_asm_label(&as);
O(OVER); LIT(1); O(AND); O(INV); LIT(1); O(ADD); O(OVER); O(TWO_SLASH); O(AND);
O(OVER); O(BANG_A); LIT(V4_CELL_BITS - 1); O(PUSH);
O(PUSH); O(OVER); O(TWO_SLASH); O(RPOP);
(void)v4_asm_label(&as);
O(MUL_STEP); O(UNEXT);
O(PUSH); O(DROP); O(RPOP);
O(TWO_STAR); O(PUSH_A);
l0 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(1); O(ADD);
l1 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l0, v4_asm_label(&as));
O(DROP);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
O(PUSH);
O(OVER);
l2 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(DUP);
l3 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l2, v4_asm_label(&as));
O(DROP); LIT(0);
v4_asm_resolve(&as, l3, v4_asm_label(&as));
O(PUSH);
O(DUP);
l4 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(OVER); LIT(1); O(AND);
l5 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l4, v4_asm_label(&as));
O(DROP); LIT(0);
v4_asm_resolve(&as, l5, v4_asm_label(&as));
O(RPOP); O(ADD); O(RPOP); O(ADD); O(PUSH);
O(AND); LIT(1); O(AND); O(PUSH_A); O(TWO_STAR); O(ADD); O(RPOP);
O(SEMI);
}
/* : UM/MOD ( ulo uhi ud -- urem uquot ) section 4
* a! 31 FOR
* -if L0
* 2* over -if L1 drop 1 + jump L2 L1: drop L2: push 2* pop
* jump SUB
* L0:
* 2* over -if L3 drop 1 + jump L4 L3: drop L4: push 2* pop
* dup a xor -if L5
* drop -if NOSUB jump SUB
* L5: drop dup inv a + inv -if L6
* drop jump NOSUB
* L6: push drop pop jump SETBIT
* SUB: inv a + inv
* SETBIT: push 1 + pop
* NOSUB:
* NEXT
* over push push drop pop pop ;
* No calls inside the loop, and SWAP is in line. */
{
v4_cell loop, l_sub, l_setbit, l_nosub;
v4_asm_ref r0, r1, r2, r3, r4, r5, r6, to_sub1, to_sub2, to_nosub1,
to_nosub2, to_setbit;
w_ummod = v4_asm_label(&as);
O(BANG_A); LIT(V4_CELL_BITS - 1); O(PUSH);
loop = v4_asm_label(&as);
r0 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
/* hi's top bit was set: shift, then subtract regardless */
O(TWO_STAR); O(OVER);
r1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(1); O(ADD);
r2 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, r1, v4_asm_label(&as));
O(DROP);
v4_asm_resolve(&as, r2, v4_asm_label(&as));
O(PUSH); O(TWO_STAR); O(RPOP);
to_sub1 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
/* L0: top bit clear: shift, then compare hi' with d */
v4_asm_resolve(&as, r0, v4_asm_label(&as));
O(TWO_STAR); O(OVER);
r3 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(1); O(ADD);
r4 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, r3, v4_asm_label(&as));
O(DROP);
v4_asm_resolve(&as, r4, v4_asm_label(&as));
O(PUSH); O(TWO_STAR); O(RPOP);
O(DUP); O(PUSH_A); O(XOR);
r5 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP);
to_nosub1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
to_sub2 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, r5, v4_asm_label(&as)); /* L5 */
O(DROP); O(DUP); O(INV); O(PUSH_A); O(ADD); O(INV);
r6 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP);
to_nosub2 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, r6, v4_asm_label(&as)); /* L6 */
O(PUSH); O(DROP); O(RPOP);
to_setbit = v4_asm_branch_fwd(&as, V4_OP_JUMP);
l_sub = v4_asm_label(&as);
O(INV); O(PUSH_A); O(ADD); O(INV);
l_setbit = v4_asm_label(&as);
O(PUSH); LIT(1); O(ADD); O(RPOP);
l_nosub = v4_asm_label(&as);
v4_asm_branch(&as, V4_OP_NEXT, loop);
O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); O(SEMI);
v4_asm_resolve(&as, to_sub1, l_sub);
v4_asm_resolve(&as, to_sub2, l_sub);
v4_asm_resolve(&as, to_setbit, l_setbit);
v4_asm_resolve(&as, to_nosub1, l_nosub);
v4_asm_resolve(&as, to_nosub2, l_nosub);
}
/* ---- section 5 words that SM/REM and /MOD rest on, as written there. */
/* : U> SWAP U< ; 5.5 */
w_ugreater = v4_asm_label(&as);
CALL(w_swap); CALL(w_uless); O(SEMI);
/* : ABS dup 0< IF NEGATE THEN ; 5.4 */
w_abs = v4_asm_label(&as);
O(DUP); CALL(w_zless); ref = if_(); CALL(w_negate); then_(ref); O(SEMI);
/* : S>D dup 0< ; 5.7 */
w_s2d = v4_asm_label(&as);
O(DUP); CALL(w_zless); O(SEMI);
/* : D+ ( d1 d2 -- d3 ) 5.7, call-free
* push over push push drop pop al bl R: bh ah
* over over xor -if L1 top bits of al, bl differ:
* drop + -if C1 jump C0 carry iff sum's top bit clear
* L1: drop over -if L2 same: carry iff both set
* drop + jump C1
* L2: drop +
* C0: pop pop + ;
* C1: pop pop + 1 + ; */
{
v4_asm_ref l1, l2, c1a, c1b, c0a;
v4_cell c0, c1;
w_dplus = v4_asm_label(&as);
O(PUSH); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP);
O(OVER); O(OVER); O(XOR);
l1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(ADD);
c1a = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
c0a = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
O(DROP); O(OVER);
l2 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(ADD);
c1b = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l2, v4_asm_label(&as));
O(DROP); O(ADD);
c0 = v4_asm_label(&as);
O(RPOP); O(RPOP); O(ADD); O(SEMI);
c1 = v4_asm_label(&as);
O(RPOP); O(RPOP); O(ADD); LIT(1); O(ADD); O(SEMI);
v4_asm_resolve(&as, c0a, c0);
v4_asm_resolve(&as, c1a, c1);
v4_asm_resolve(&as, c1b, c1);
}
/* : DNEGATE ( d -- -d ) 5.7, call-free
* inv over if L1 drop push inv 1 + pop ;
* L1: drop 1 + ;
* -d = ~d + 1: the + 1 carries into hi exactly when lo = 0. */
w_dnegate = v4_asm_label(&as);
O(INV); O(OVER); ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(PUSH); O(INV); LIT(1); O(ADD); O(RPOP); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); LIT(1); O(ADD); O(SEMI);
/* : DABS dup 0< IF DNEGATE THEN ; 5.7 */
w_dabs = v4_asm_label(&as);
O(DUP); CALL(w_zless); ref = if_(); CALL(w_dnegate); then_(ref); O(SEMI);
/* : SM/REM ( d n -- rem quot ) section 4
* over over xor push R: quotient sign (top bit)
* over push R: + remainder sign (of d)
* -if L0 inv 1 + L0: push R: + |n|
* -if L1 DNEGATE L1: |d|
* pop UM/MOD urem uquot
* pop -if L2 drop push inv 1 + pop jump L3 L2: drop L3:
* pop -if L4 drop inv 1 + ; L4: drop ;
* Sign tests are native -if, as in 0<; NEGATE is in line. */
{
v4_asm_ref l0, l1, l2, l3, l4;
w_smrem = v4_asm_label(&as);
O(OVER); O(OVER); O(XOR); O(PUSH); O(OVER); O(PUSH);
l0 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(INV); LIT(1); O(ADD);
v4_asm_resolve(&as, l0, v4_asm_label(&as));
O(PUSH);
l1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
CALL(w_dnegate);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
O(RPOP); CALL(w_ummod);
O(RPOP);
l2 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(PUSH); O(INV); LIT(1); O(ADD); O(RPOP);
l3 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l2, v4_asm_label(&as));
O(DROP);
v4_asm_resolve(&as, l3, v4_asm_label(&as));
O(RPOP);
l4 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(INV); LIT(1); O(ADD); O(SEMI);
v4_asm_resolve(&as, l4, v4_asm_label(&as));
O(DROP); O(SEMI);
}
/* : /MOD push S>D pop SM/REM ; 5.6 */
w_slashmod = v4_asm_label(&as);
O(PUSH); CALL(w_s2d); O(RPOP); CALL(w_smrem); O(SEMI);
/* : * UM* drop ; 5.4 */
w_star = v4_asm_label(&as);
CALL(w_umstar); O(DROP); O(SEMI);
/* : / /MOD NIP ; 5.4
* with /MOD's body and NIP in line (push S>D pop SM/REM push drop pop),
* so it runs no deeper than /MOD. */
w_slash = v4_asm_label(&as);
O(PUSH); CALL(w_s2d); O(RPOP); CALL(w_smrem); O(PUSH); O(DROP); O(RPOP); O(SEMI);
/* : M+ S>D D+ ; 5.6 */
w_mplus = v4_asm_label(&as);
CALL(w_s2d); CALL(w_dplus); O(SEMI);
/* : D- DNEGATE D+ ; 5.7 */
w_dminus = v4_asm_label(&as);
CALL(w_dnegate); CALL(w_dplus); O(SEMI);
/* : D0= OR 0= ; 5.7 */
w_d0equal = v4_asm_label(&as);
CALL(w_or); CALL(w_zequal); O(SEMI);
/* : D= D- D0= ; 5.7 */
w_dequal = v4_asm_label(&as);
CALL(w_dminus); CALL(w_d0equal); O(SEMI);
/* +* with S = 0 never adds: each step is an exact arithmetic right shift
* of the double T:A by one bit. Both words below are built on that. */
/* : Q.FROM-INT ( n -- q ) 5.26
* push 0 a! 0 pop 15 FOR +* UNEXT push drop a pop ;
* T:A = n:0 is n * 2^N; shifted right N-16 bits it is n * 2^16. 15 is
* N-17 at 32-bit cells (47 at 64). Clobbers A. */
w_qfromint = v4_asm_label(&as);
O(PUSH); LIT(0); O(BANG_A); LIT(0); O(RPOP);
LIT(V4_CELL_BITS - 17); O(PUSH);
(void)v4_asm_label(&as);
O(MUL_STEP); O(UNEXT);
O(PUSH); O(DROP); O(PUSH_A); O(RPOP); O(SEMI);
/* : Q.TO-INT ( q -- n ) 5.26
* push a! 0 pop 15 FOR +* UNEXT drop drop a ;
* T:A = q shifted right 16 bits; the low cell is A. 15 at every cell
* width. Clobbers A. */
w_qtoint = v4_asm_label(&as);
O(PUSH); O(BANG_A); LIT(0); O(RPOP);
LIT(15); O(PUSH);
(void)v4_asm_label(&as);
O(MUL_STEP); O(UNEXT);
O(DROP); O(DROP); O(PUSH_A); O(SEMI);
/* : = xor 0= ; 5.5 */
w_equal = v4_asm_label(&as);
O(XOR); CALL(w_zequal); O(SEMI);
/* : < 2DUP xor 0< IF drop 0< ELSE - 0< THEN ; 5.5 */
w_less = v4_asm_label(&as);
O(OVER); O(OVER); O(XOR); CALL(w_zless);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(DROP); CALL(w_zless);
{
v4_asm_ref j = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); CALL(w_minus); CALL(w_zless);
v4_asm_resolve(&as, j, v4_asm_label(&as));
}
O(SEMI);
/* : D< ROT 2DUP = IF 2DROP U< ELSE SWAP < NIP NIP THEN ; 5.7 */
w_dless = v4_asm_label(&as);
CALL(w_rot); O(OVER); O(OVER); CALL(w_equal);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(DROP); O(DROP); CALL(w_uless);
{
v4_asm_ref j = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); CALL(w_swap); CALL(w_less); CALL(w_nip); CALL(w_nip);
v4_asm_resolve(&as, j, v4_asm_label(&as));
}
O(SEMI);
/* : 2SWAP ROT push ROT pop ; 5.7, call-free
* with ROT (push SWAP pop SWAP) and SWAP (over push push drop pop pop)
* written in line:
* push over push push drop pop pop pop over push push drop pop pop
* push
* push over push push drop pop pop pop over push push drop pop pop
* pop ; */
#define ROT_INLINE() do { O(PUSH); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); \
O(RPOP); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0)
w_2swap = v4_asm_label(&as);
ROT_INLINE(); O(PUSH); ROT_INLINE(); O(RPOP); O(SEMI);
#undef ROT_INLINE
/* : 2OVER push push 2DUP pop pop 2SWAP ; 5.7 */
w_2over = v4_asm_label(&as);
O(PUSH); O(PUSH); O(OVER); O(OVER); O(RPOP); O(RPOP); CALL(w_2swap); O(SEMI);
/* : (D<) ( d1 d2 -- d1 d2 flag ) 5.7, call-free
* Signed d1 < d2, leaving both doubles. Copies of ah and bh go on top;
* if they differ, flag comes from them, else from al U< bl. In every
* case the tested value has its top bit set exactly when d1 < d2:
* highs, signs differ: ah (d1 < d2 iff ah < 0)
* highs, signs agree: ah - bh (cannot overflow)
* lows, top bits differ: bl (al U< bl iff bl's is set)
* lows, top bits agree: al - bl
* x - y with y on top is `push inv pop + inv`.
*
* dup push push over pop al ah bl ah bh R: bh
* over over xor if TIE
* drop over over xor -if HS
* drop drop jump S1 al ah bl ah
* HS: drop push inv pop + inv al ah bl ah-bh
* S1: -if N1 drop -1 jump D1 N1: drop 0
* D1: pop over push push drop pop pop ; al ah bl bh f
* TIE: drop drop drop al ah bl R: bh
* dup push push over pop al ah al bl R: bh bl
* over over xor -if LS
* drop push drop pop jump S2 al ah bl
* LS: drop push inv pop + inv al ah al-bl
* S2: -if N2 drop -1 jump D2 N2: drop 0
* D2: pop pop al ah f bl bh
* push over push push drop pop pop pop over push push drop pop pop ; */
{
v4_asm_ref tie, hs, s1a, n1, d1, ls, s2a, n2, d2;
w_dltkeep = v4_asm_label(&as);
O(DUP); O(PUSH); O(PUSH); O(OVER); O(RPOP);
O(OVER); O(OVER); O(XOR);
tie = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(OVER); O(OVER); O(XOR);
hs = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(DROP);
s1a = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, hs, v4_asm_label(&as));
O(DROP); O(PUSH); O(INV); O(RPOP); O(ADD); O(INV);
v4_asm_resolve(&as, s1a, v4_asm_label(&as));
n1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(-1);
d1 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, n1, v4_asm_label(&as));
O(DROP); LIT(0);
v4_asm_resolve(&as, d1, v4_asm_label(&as));
O(RPOP); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); O(SEMI);
v4_asm_resolve(&as, tie, v4_asm_label(&as));
O(DROP); O(DROP); O(DROP);
O(DUP); O(PUSH); O(PUSH); O(OVER); O(RPOP);
O(OVER); O(OVER); O(XOR);
ls = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(PUSH); O(DROP); O(RPOP);
s2a = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, ls, v4_asm_label(&as));
O(DROP); O(PUSH); O(INV); O(RPOP); O(ADD); O(INV);
v4_asm_resolve(&as, s2a, v4_asm_label(&as));
n2 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); LIT(-1);
d2 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, n2, v4_asm_label(&as));
O(DROP); LIT(0);
v4_asm_resolve(&as, d2, v4_asm_label(&as));
O(RPOP); O(RPOP);
O(PUSH); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP);
O(RPOP); O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); O(SEMI);
}
/* : DMAX ( d1 d2 -- d3 ) (D<) if L drop push push drop drop pop pop ;
* L: drop drop drop ; 5.7 */
w_dmax = v4_asm_label(&as);
CALL(w_dltkeep);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(PUSH); O(PUSH); O(DROP); O(DROP); O(RPOP); O(RPOP); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); O(DROP); O(DROP); O(SEMI);
/* : DMIN ( d1 d2 -- d3 ) (D<) if L drop drop drop ;
* L: drop push push drop drop pop pop ; 5.7 */
w_dmin = v4_asm_label(&as);
CALL(w_dltkeep);
ref = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(DROP); O(DROP); O(SEMI);
v4_asm_resolve(&as, ref, v4_asm_label(&as));
O(DROP); O(PUSH); O(PUSH); O(DROP); O(DROP); O(RPOP); O(RPOP); O(SEMI);
/* Q.> as 5.26 gives it, SWAP D<, and with 2SWAP. */
w_qgt_doc = v4_asm_label(&as);
CALL(w_swap); CALL(w_dless); O(SEMI);
w_qgt = v4_asm_label(&as);
CALL(w_2swap); CALL(w_dless); O(SEMI);
/* : Q.* ( a b -- c ) 5.26
* a = a0 a1, b = b0 b1. Cells 0..2 of the signed product a*b, then
* shifted right 16. The cell products are unsigned; reading a1 and b1
* as signed takes (a1<0 ? b0 : 0) + (b1<0 ? a0 : 0) off cell 2, the
* only cell above the unsigned product's that the result reaches.
* Products in the order a1*b1 (low cell), a0*b1, a1*b0, a0*b0, each
* input dropped after its last use; b0 waits on the return stack.
* SWAP push over over UM* drop a0 a1 b1 t3 R: b0
* push push over pop a0 a1 a0 b1 R: b0 t3
* -if L1 SWAP jump L2 L1: SWAP drop 0 L2: a0 a1 b1 c2
* pop SWAP - push a0 a1 b1 R: b0 t3-c2
* push over pop UM* pop + a0 a1 x0 x1 R: b0
* ROT pop dup push a0 x0 x1 a1 b0
* over -if L3 drop dup jump L4 L3: drop 0 L4: ... a1 b0 c1
* push UM* pop - D+ a0 X0 X1
* ROT pop UM* X0 X1 l00 h00
* SWAP push 0 D+ pop c1 c2 c0
* push over pop SWAP Q.TO-INT push Q.TO-INT pop SWAP ;
* SWAP and ROT in line. Clobbers A. */
#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)
{
v4_asm_ref l1, l2, l3, l4;
w_qstar = v4_asm_label(&as);
SWAP_INLINE(); O(PUSH); O(OVER); O(OVER); CALL(w_umstar); O(DROP);
O(PUSH); O(PUSH); O(OVER); O(RPOP);
l1 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
SWAP_INLINE();
l2 = v4_asm_branch_fwd(&as, V4_OP_JUMP);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
SWAP_INLINE(); O(DROP); LIT(0);
v4_asm_resolve(&as, l2, v4_asm_label(&as));
O(RPOP); SWAP_INLINE(); CALL(w_minus); O(PUSH);
O(PUSH); O(OVER); O(RPOP); CALL(w_umstar); O(RPOP); O(ADD);
ROT_INLINE(); O(RPOP); O(DUP); O(PUSH);
O(OVER);
l3 = v4_asm_branch_fwd(&as, V4_OP_MINUS_IF);
O(DROP); O(DUP);
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); CALL(w_umstar); O(RPOP); CALL(w_minus); CALL(w_dplus);
ROT_INLINE(); O(RPOP); CALL(w_umstar);
SWAP_INLINE(); O(PUSH); LIT(0); CALL(w_dplus); O(RPOP);
O(PUSH); O(OVER); O(RPOP); SWAP_INLINE(); CALL(w_qtoint);
O(PUSH); CALL(w_qtoint); O(RPOP); SWAP_INLINE();
O(SEMI);
}
#undef ROT_INLINE
#undef SWAP_INLINE
/* : 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. 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);
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(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(OVER); O(ADD); O(ADD);
O(PUSH); O(TWO_STAR); O(PUSH_A); O(ADD); O(RPOP); O(RPOP);
O(SEMI);
}
/* : (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
* 2^(2N-1) (Q./ checks both first). Restoring division, 2N+16 steps,
* on one shifting register: quotient q (starts as a) above remainder r.
* The bit leaving q's top enters r; each step's quotient bit enters q at
* the next step's shift, and a final shift takes in the last one. With
* no overflow the bits leaving q in the last 16 steps are zeros.
* q and b live in the scratch variable (Q/) = QS: q0 q1 b0 b1, so the
* 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
* 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
* drop drop -if NT1 jump TAKE NT1: jump NOTAKE
* SAMEH: drop over SWAP push inv pop + inv r0 r1 r1-b1
* -if T2 drop jump NOTAKE T2: drop jump TAKE
* EQH: drop drop over QS 2 + b! @b r0 r1 r0 b0
* over over xor -if SAMEL
* drop drop -if NT3 drop jump TAKEL NT3: drop jump NOTAKEL
* SAMEL: drop push inv pop + inv r0 r1 r0-b0
* -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: 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
* 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)
#define ROT_INLINE() do { O(PUSH); SWAP_INLINE(); O(RPOP); SWAP_INLINE(); } while (0)
#define SUB_INLINE() do { O(PUSH); O(INV); O(RPOP); O(ADD); O(INV); } while (0)
#define FWD(op) v4_asm_branch_fwd(&as, V4_OP_##op)
#define HERE_(r) v4_asm_resolve(&as, (r), v4_asm_label(&as))
{
v4_cell loop, l_take, l_notake, l_takel, l_notakel, l_end;
v4_asm_ref cmp, eqh, sameh, nt1, t2, samel, nt3, t4;
v4_asm_ref j_take[4], j_notake[3], j_takel[3], j_notakel[3], j_end[3];
unsigned nt = 0, nn = 0, ntl = 0, nnl = 0, ne = 0, k;
w_uqdiv = v4_asm_label(&as);
LIT(QS + 2); O(BANG_A); SWAP_INLINE(); O(STORE_INC); O(STORE_A);
LIT(QS); O(BANG_A); SWAP_INLINE(); O(STORE_INC); O(STORE_A);
LIT(0); LIT(0); LIT(0);
LIT(2 * V4_CELL_BITS + 15); O(PUSH);
loop = v4_asm_label(&as);
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);
HERE_(cmp); /* CMP */
O(DROP); LIT(QS + 3); O(BANG_B); O(FETCH_B);
O(OVER); O(OVER); O(XOR); eqh = FWD(IF);
sameh = FWD(MINUS_IF);
O(DROP); O(DROP); nt1 = FWD(MINUS_IF);
j_take[nt++] = FWD(JUMP);
HERE_(nt1); j_notake[nn++] = FWD(JUMP);
HERE_(sameh); /* SAMEH */
O(DROP); O(OVER); SWAP_INLINE(); SUB_INLINE(); t2 = FWD(MINUS_IF);
O(DROP); j_notake[nn++] = FWD(JUMP);
HERE_(t2); O(DROP); j_take[nt++] = FWD(JUMP);
HERE_(eqh); /* EQH */
O(DROP); O(DROP); O(OVER); LIT(QS + 2); O(BANG_B); O(FETCH_B);
O(OVER); O(OVER); O(XOR); samel = FWD(MINUS_IF);
O(DROP); O(DROP); nt3 = FWD(MINUS_IF);
O(DROP); j_takel[ntl++] = FWD(JUMP);
HERE_(nt3); O(DROP); j_notakel[nnl++] = FWD(JUMP);
HERE_(samel); /* SAMEL */
O(DROP); SUB_INLINE(); t4 = FWD(MINUS_IF);
O(DROP); j_notakel[nnl++] = FWD(JUMP);
HERE_(t4); O(DROP); j_takel[ntl++] = FWD(JUMP);
l_takel = v4_asm_label(&as); /* TAKEL */
O(DROP); LIT(QS + 2); O(BANG_B); O(FETCH_B); SUB_INLINE(); LIT(0); LIT(1);
j_end[ne++] = FWD(JUMP);
l_notakel = v4_asm_label(&as); /* NOTAKEL */
LIT(0); j_end[ne++] = FWD(JUMP);
l_take = v4_asm_label(&as); /* TAKE */
{
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);
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);
for (k = 0; k < nn; k++) v4_asm_resolve(&as, j_notake[k], l_notake);
for (k = 0; k < ntl; k++) v4_asm_resolve(&as, j_takel[k], l_takel);
for (k = 0; k < nnl; k++) v4_asm_resolve(&as, j_notakel[k], l_notakel);
for (k = 0; k < ne; k++) v4_asm_resolve(&as, j_end[k], l_end);
}
#undef HERE_
#undef FWD
#undef SUB_INLINE
#undef ROT_INLINE
#undef SWAP_INLINE
/* : Q./ ( a b -- q ) 5.26, D-11
* over over OR if ZERO drop
* push over pop SWAP over xor QS 4 + b! !b a0 a1 b0 b1 (Q/)+4: sign
* -if L1 DNEGATE L1: push push -if L2 DNEGATE L2: pop pop |a| |b|
* dup if CHK drop jump DIV b1 != 0: no overflow
* CHK: drop over -131072 and if CHK2 drop jump DIV b0 >= 2^17: none
* 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 a0 a1 b0 0
* OV: drop drop drop pop pop drop drop
* QS 4 + b! @b -if OVP drop 0 MSB ; OVP: drop -1 MAXHI ;
* DIV: (UQ/) QS 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 ;
* MSB and MAXHI are the high cells of Q min and Q max. Clobbers A and B. */
#define SWAP_INLINE() do { O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0)
#define SUB_INLINE() do { O(PUSH); O(INV); O(RPOP); O(ADD); O(INV); } while (0)
#define FWD(op) v4_asm_branch_fwd(&as, V4_OP_##op)
#define HERE_(r) v4_asm_resolve(&as, (r), v4_asm_label(&as))
{
v4_asm_ref zero, l1, l2, chk, chk2, same, noov1, ov1, ov2, ovp, l3, z0, zp;
v4_asm_ref d1, d2, d3;
v4_cell l_div, l_noov, l_ov;
w_qslash = v4_asm_label(&as);
O(OVER); O(OVER); CALL(w_or); zero = FWD(IF);
O(DROP);
O(PUSH); O(OVER); O(RPOP); SWAP_INLINE(); O(OVER); O(XOR);
LIT(QS + 4); O(BANG_B); O(STORE_B);
l1 = FWD(MINUS_IF); CALL(w_dnegate); HERE_(l1);
O(PUSH); O(PUSH);
l2 = FWD(MINUS_IF); CALL(w_dnegate); HERE_(l2);
O(RPOP); O(RPOP);
O(DUP); chk = FWD(IF);
O(DROP); d1 = FWD(JUMP);
HERE_(chk);
O(DROP); O(OVER); LIT(-131072); O(AND); chk2 = FWD(IF);
O(DROP); d2 = FWD(JUMP);
HERE_(chk2);
O(DROP); O(PUSH); O(PUSH); O(DUP); O(RPOP); O(DUP); O(PUSH);
LIT(V4_CELL_BITS - 18); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_STAR); O(UNEXT);
O(OVER); O(OVER); O(XOR); same = FWD(MINUS_IF);
O(DROP); O(DROP); noov1 = FWD(MINUS_IF);
ov1 = FWD(JUMP);
HERE_(same);
O(DROP); SUB_INLINE(); ov2 = FWD(MINUS_IF);
l_noov = v4_asm_label(&as);
O(DROP); O(RPOP); O(RPOP); d3 = FWD(JUMP);
l_ov = v4_asm_label(&as);
O(DROP); O(DROP); O(DROP); O(RPOP); O(RPOP); O(DROP); O(DROP);
LIT(QS + 4); O(BANG_B); O(FETCH_B); ovp = FWD(MINUS_IF);
O(DROP); LIT(0); LIT((v4_cell)V4_MSB); O(SEMI);
HERE_(ovp); O(DROP); LIT(-1); LIT((v4_cell)(MAXU >> 1)); O(SEMI);
l_div = v4_asm_label(&as);
CALL(w_uqdiv); LIT(QS + 4); O(BANG_B); O(FETCH_B); l3 = FWD(MINUS_IF);
O(DROP); CALL(w_dnegate); O(SEMI);
HERE_(l3); O(DROP); O(SEMI);
HERE_(zero);
O(DROP); O(DROP); O(DROP);
LIT(NODE_ERROR); O(BANG_B); LIT(-1); O(STORE_B);
O(OVER); O(OVER); CALL(w_or); z0 = FWD(IF);
O(DROP); zp = FWD(MINUS_IF);
O(DROP); O(DROP); LIT(0); LIT((v4_cell)V4_MSB); O(SEMI);
HERE_(zp); O(DROP); O(DROP); LIT(-1); LIT((v4_cell)(MAXU >> 1)); O(SEMI);
HERE_(z0); O(DROP); O(SEMI);
v4_asm_resolve(&as, d1, l_div);
v4_asm_resolve(&as, d2, l_div);
v4_asm_resolve(&as, d3, l_div);
v4_asm_resolve(&as, noov1, l_noov);
v4_asm_resolve(&as, ov1, l_ov);
v4_asm_resolve(&as, ov2, l_ov);
}
#undef HERE_
#undef FWD
#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
/* : 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
/* : Q.LOG ( x -- ln x ) 5.26, v3's q48_log_approx
* x = 2^k * m with 1.0 <= m < 2.0; ln m by up to 6 Newton rounds on
* e^y = m from y = m - 1.0; result y + k * ln 2 (ln 2 = 45426). e^y is
* Q.EXP's Taylor loop written here in line, so Q.LOG calls only Q.*, Q./
* and UM*. After the reduction m, y, the Taylor term and sum all fit one
* cell, so (QL) holds single cells: m y k rounds term sum n*1.0.
* x <= 0 gives 0 and sets NODE-ERROR (D-12). D2/ is the logical double
* shift right, as in Q.SQRT.
* dup -if POS drop jump ERR
* POS: drop over over OR if ZERO drop 0 QL 2 + b! !b
* DOWN: dup if CHKLO drop jump SHR
* CHKLO: drop over -131072 and if DOWNDONE drop
* SHR: D2/ QL 2 + b! @b 1 + !b jump DOWN x >= 2.0: halve, k+1
* DOWNDONE: drop drop m
* UP: dup -65536 and if UPSHIFT drop jump NORM
* UPSHIFT: drop 2* QL 2 + b! @b -1 + !b jump UP m < 1.0: double, k-1
* NORM: dup QL b! !b -65536 + QL 1 + b! !b 6 QL 3 + b! !b
* NEWTON: QL 1 + b! @b dup QL 4 + b! !b 65536 + QL 5 + b! !b
* 131072 QL 6 + b! !b
* EXPL: QL 4 + b! @b 0 QL 1 + b! @b 0 Q.* drop
* 0 QL 6 + b! @b 0 Q./ drop dup QL 4 + b! !b
* dup QL 5 + b! @b + !b
* -50 + -if EXPC drop jump EXPD
* EXPC: drop QL 6 + b! @b 65536 + dup !b -720896 + -if EXPD2 drop jump EXPL
* EXPD2: drop
* EXPD: QL b! @b QL 5 + b! @b push inv pop + inv delta = m - e^y
* dup QL 4 + b! !b delta waits in the term's cell
* if DZ -if DPOS
* inv 1 + 0 QL 5 + b! @b 0 Q./ drop corr = -delta / e^y
* QL 1 + b! @b SWAP push inv pop + inv -if KEEP drop 0 KEEP: QL 1 + b! !b
* jump TEST y = max(y - corr, 0)
* DPOS: 0 QL 5 + b! @b 0 Q./ drop QL 1 + b! @b + !b jump TEST
* DZ: drop
* TEST: QL 4 + b! @b -if AP inv 1 + AP: -100 + -if CONT drop jump FIN
* CONT: drop QL 3 + b! @b -1 + dup !b if FIN0 drop jump NEWTON
* FIN0: drop
* FIN: QL 1 + b! @b QL 2 + b! @b
* -if KP inv 1 + 45426 UM* drop inv 1 + jump KS KP: 45426 UM* drop
* KS: + dup -if SP drop -1 ; SP: drop 0 ;
* ZERO: drop ERR: drop drop NODE-ERROR b! -1 !b 0 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 JUMP_TO(l) v4_asm_branch(&as, V4_OP_JUMP, (l))
#define VGET(a) do { LIT(a); O(BANG_B); O(FETCH_B); } while (0)
#define VSET(a) do { LIT(a); O(BANG_B); O(STORE_B); } while (0)
#define SWAP_INLINE() do { O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0)
#define SUB_INLINE() do { O(PUSH); O(INV); O(RPOP); O(ADD); O(INV); } 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, j_err, zero, chklo, j_shr, downdone, upshift, j_norm, expc, j_expd,
expd2, dzero, dpos, keep, j_test, ap, contn, j_fin, fin0, kp, j_ks, sp;
v4_cell l_down, l_up, l_newton, l_expl, l_test, l_err;
w_qlog = v4_asm_label(&as);
O(DUP); pos = FWD(MINUS_IF);
O(DROP); j_err = FWD(JUMP);
HERE_(pos);
O(DROP); O(OVER); O(OVER); CALL(w_or); zero = FWD(IF);
O(DROP); LIT(0); VSET(QL + 2);
l_down = v4_asm_label(&as);
O(DUP); chklo = FWD(IF);
O(DROP); j_shr = FWD(JUMP);
HERE_(chklo);
O(DROP); O(OVER); LIT(-131072); O(AND); downdone = FWD(IF);
O(DROP);
HERE_(j_shr);
D2SLASH(); VGET(QL + 2); LIT(1); O(ADD); O(STORE_B); JUMP_TO(l_down);
HERE_(downdone);
O(DROP); O(DROP);
l_up = v4_asm_label(&as);
O(DUP); LIT(-65536); O(AND); upshift = FWD(IF);
O(DROP); j_norm = FWD(JUMP);
HERE_(upshift);
O(DROP); O(TWO_STAR); VGET(QL + 2); LIT(-1); O(ADD); O(STORE_B); JUMP_TO(l_up);
HERE_(j_norm);
O(DUP); VSET(QL); LIT(-65536); O(ADD); VSET(QL + 1); LIT(6); VSET(QL + 3);
l_newton = v4_asm_label(&as);
VGET(QL + 1); O(DUP); VSET(QL + 4); LIT(65536); O(ADD); VSET(QL + 5);
LIT(131072); VSET(QL + 6);
l_expl = v4_asm_label(&as);
VGET(QL + 4); LIT(0); VGET(QL + 1); LIT(0); CALL(w_qstar); O(DROP);
LIT(0); VGET(QL + 6); LIT(0); CALL(w_qslash); O(DROP);
O(DUP); VSET(QL + 4);
O(DUP); VGET(QL + 5); O(ADD); O(STORE_B);
LIT(-50); O(ADD); expc = FWD(MINUS_IF);
O(DROP); j_expd = FWD(JUMP);
HERE_(expc);
O(DROP); VGET(QL + 6); LIT(65536); O(ADD); O(DUP); O(STORE_B);
LIT(-720896); O(ADD); expd2 = FWD(MINUS_IF);
O(DROP); JUMP_TO(l_expl);
HERE_(expd2);
O(DROP);
HERE_(j_expd);
VGET(QL); VGET(QL + 5); SUB_INLINE(); O(DUP); VSET(QL + 4);
dzero = FWD(IF);
dpos = FWD(MINUS_IF);
O(INV); LIT(1); O(ADD); LIT(0); VGET(QL + 5); LIT(0); CALL(w_qslash); O(DROP);
VGET(QL + 1); SWAP_INLINE(); SUB_INLINE();
keep = FWD(MINUS_IF); O(DROP); LIT(0); HERE_(keep);
VSET(QL + 1); j_test = FWD(JUMP);
HERE_(dpos);
LIT(0); VGET(QL + 5); LIT(0); CALL(w_qslash); O(DROP);
VGET(QL + 1); O(ADD); O(STORE_B);
{
v4_asm_ref j_test2 = FWD(JUMP);
HERE_(dzero);
O(DROP);
l_test = v4_asm_label(&as);
v4_asm_resolve(&as, j_test2, l_test);
}
VGET(QL + 4);
ap = FWD(MINUS_IF); O(INV); LIT(1); O(ADD); HERE_(ap);
LIT(-100); O(ADD); contn = FWD(MINUS_IF);
O(DROP); j_fin = FWD(JUMP);
HERE_(contn);
O(DROP); VGET(QL + 3); LIT(-1); O(ADD); O(DUP); O(STORE_B);
fin0 = FWD(IF);
O(DROP); JUMP_TO(l_newton);
HERE_(fin0);
O(DROP);
HERE_(j_fin);
VGET(QL + 1); VGET(QL + 2);
kp = FWD(MINUS_IF);
O(INV); LIT(1); O(ADD); LIT(45426); CALL(w_umstar); O(DROP); O(INV); LIT(1); O(ADD);
j_ks = FWD(JUMP);
HERE_(kp);
LIT(45426); CALL(w_umstar); O(DROP);
HERE_(j_ks);
O(ADD); O(DUP); sp = FWD(MINUS_IF);
O(DROP); LIT(-1); O(SEMI);
HERE_(sp); O(DROP); LIT(0); O(SEMI);
HERE_(zero);
O(DROP);
l_err = v4_asm_label(&as);
O(DROP); O(DROP); LIT(NODE_ERROR); O(BANG_B); LIT(-1); O(STORE_B); LIT(0); LIT(0); O(SEMI);
v4_asm_resolve(&as, j_err, l_err);
v4_asm_resolve(&as, j_test, l_test);
}
#undef D2SLASH
#undef SUB_INLINE
#undef SWAP_INLINE
#undef VSET
#undef VGET
#undef JUMP_TO
#undef HERE_
#undef FWD
#define FWD(op) v4_asm_branch_fwd(&as, V4_OP_##op)
#define HERE_(r) v4_asm_resolve(&as, (r), v4_asm_label(&as))
#define JUMP_TO(l) v4_asm_branch(&as, V4_OP_JUMP, (l))
#define VGET(a) do { LIT(a); O(BANG_B); O(FETCH_B); } while (0)
#define VSET(a) do { LIT(a); O(BANG_B); O(STORE_B); } while (0)
/* : (Q.REDUCE) ( x -- r ) 5.26, v3's q48_reduce_angle
* The angle x reduced to one cell r, -pi <= r <= pi (pi = 205887,
* 2pi = 411774): the remainder of x / 2pi with the sign of x, then one
* step of 2pi back into range. x / 2pi does not fit a cell, but only
* the remainder is wanted: |x| mod 2pi is two UM/MOD steps, the high
* cell first, its remainder then leading the low cell.
* dup QT b! !b -if L1 DNEGATE L1: u0 u1 (QT): sign of x
* 0 411774 UM/MOD drop 411774 UM/MOD drop |x| mod 2pi
* QT b! @b -if P1 drop inv 1 + jump J1 P1: drop J1: r
* dup -205888 + -if BIG drop
* dup 205887 + -if OK drop 411774 + ;
* OK: drop ;
* BIG: drop -411774 + ;
* Clobbers A and B. */
{
v4_asm_ref l1, p1, j1, big, ok;
w_qreduce = v4_asm_label(&as);
O(DUP); VSET(QT);
l1 = FWD(MINUS_IF); CALL(w_dnegate); HERE_(l1);
LIT(0); LIT(411774); CALL(w_ummod); O(DROP);
LIT(411774); CALL(w_ummod); O(DROP);
VGET(QT); p1 = FWD(MINUS_IF);
O(DROP); O(INV); LIT(1); O(ADD); j1 = FWD(JUMP);
HERE_(p1); O(DROP);
HERE_(j1);
O(DUP); LIT(-205888); O(ADD); big = FWD(MINUS_IF);
O(DROP);
O(DUP); LIT(205887); O(ADD); ok = FWD(MINUS_IF);
O(DROP); LIT(411774); O(ADD); O(SEMI);
HERE_(ok); O(DROP); O(SEMI);
HERE_(big); O(DROP); LIT(-411774); O(ADD); O(SEMI);
}
/* : Q.SIN ( x -- sin x ) 5.26, v3's q48_sin_approx
* On r = (Q.REDUCE) x and xu = |r|: xu - xu^3/3! + xu^5/5! - ... to
* n = 11, each term the last times xu^2 over n(n-1), stopping once a
* term is below 10 ulp; negated if r < 0. After the reduction every
* value fits one cell, so (QT) holds single cells: sign, x^2, term,
* sum, n, subtract flag. The loop is TRIG; Q.COS jumps into it.
* (Q.REDUCE) dup QT b! !b -if A1 inv 1 + A1: xu
* dup QT 2 + b! !b dup QT 3 + b! !b term = sum = xu
* 0 over 0 Q.* drop QT 1 + b! !b x^2
* 3 QT 4 + b! !b -1 QT 5 + b! !b
* TRIG: QT 2 + b! @b 0 QT 1 + b! @b 0 Q.* drop
* 0 QT 4 + b! @b dup -1 + UM* drop 65536 UM* drop 0 Q./ drop
* dup QT 2 + b! !b term
* QT 5 + b! @b if ADDT drop inv 1 + jump ACC ADDT: drop
* ACC: QT 3 + b! @b + !b QT 5 + b! @b inv !b sum, flip the flag
* QT 2 + b! @b -10 + -if CONT drop jump END
* CONT: drop QT 4 + b! @b 2 + dup !b -12 + -if END2 drop jump TRIG
* END2: drop
* END: QT 3 + b! @b QT b! @b -if PR drop inv 1 + jump SD PR: drop
* SD: dup -if SP drop -1 ; SP: drop 0 ;
* Clobbers A and B. */
{
v4_asm_ref a1, addt, j_acc, contn, j_end, end2, pr, j_sd, sp;
w_qsin = v4_asm_label(&as);
CALL(w_qreduce); O(DUP); VSET(QT);
a1 = FWD(MINUS_IF); O(INV); LIT(1); O(ADD); HERE_(a1);
O(DUP); VSET(QT + 2); O(DUP); VSET(QT + 3);
LIT(0); O(OVER); LIT(0); CALL(w_qstar); O(DROP); VSET(QT + 1);
LIT(3); VSET(QT + 4); LIT(-1); VSET(QT + 5);
l_trig = v4_asm_label(&as);
VGET(QT + 2); LIT(0); VGET(QT + 1); LIT(0); CALL(w_qstar); O(DROP);
LIT(0); VGET(QT + 4); O(DUP); LIT(-1); O(ADD); CALL(w_umstar); O(DROP);
LIT(65536); CALL(w_umstar); O(DROP); LIT(0); CALL(w_qslash); O(DROP);
O(DUP); VSET(QT + 2);
VGET(QT + 5); addt = FWD(IF);
O(DROP); O(INV); LIT(1); O(ADD); j_acc = FWD(JUMP);
HERE_(addt); O(DROP);
HERE_(j_acc);
VGET(QT + 3); O(ADD); O(STORE_B);
VGET(QT + 5); O(INV); O(STORE_B);
VGET(QT + 2); LIT(-10); O(ADD); contn = FWD(MINUS_IF);
O(DROP); j_end = FWD(JUMP);
HERE_(contn);
O(DROP); VGET(QT + 4); LIT(2); O(ADD); O(DUP); O(STORE_B);
LIT(-12); O(ADD); end2 = FWD(MINUS_IF);
O(DROP); JUMP_TO(l_trig);
HERE_(end2); O(DROP);
HERE_(j_end);
VGET(QT + 3); VGET(QT); pr = FWD(MINUS_IF);
O(DROP); O(INV); LIT(1); O(ADD); j_sd = FWD(JUMP);
HERE_(pr); O(DROP);
HERE_(j_sd);
O(DUP); sp = FWD(MINUS_IF);
O(DROP); LIT(-1); O(SEMI);
HERE_(sp); O(DROP); LIT(0); O(SEMI);
}
/* : Q.COS ( x -- cos x ) 5.26, v3's q48_cos_approx
* On xu = |(Q.REDUCE) x|: 1 - xu^2/2! + xu^4/4! - ... to n = 10, by
* Q.SIN's loop: term and sum start at 1.0, n at 2, and the sign cell is
* 0 so nothing is negated at the end.
* (Q.REDUCE) -if A1 inv 1 + A1: xu
* 0 QT b! !b 0 over 0 Q.* drop QT 1 + b! !b x^2
* 65536 QT 2 + b! !b 65536 QT 3 + b! !b
* 2 QT 4 + b! !b -1 QT 5 + b! !b jump TRIG
* Clobbers A and B. */
{
v4_asm_ref a1;
w_qcos = v4_asm_label(&as);
CALL(w_qreduce);
a1 = FWD(MINUS_IF); O(INV); LIT(1); O(ADD); HERE_(a1);
LIT(0); VSET(QT);
LIT(0); O(OVER); LIT(0); CALL(w_qstar); O(DROP); VSET(QT + 1);
LIT(65536); VSET(QT + 2); LIT(65536); VSET(QT + 3);
LIT(2); VSET(QT + 4); LIT(-1); VSET(QT + 5);
JUMP_TO(l_trig);
}
#undef VSET
#undef VGET
#undef JUMP_TO
#undef HERE_
#undef FWD
/* Q.1, Q.0, Q.SCALE 5.26, fate IN
* In-line double-cell constants: the compiler places two literals, low
* cell first. Q.1 and Q.SCALE are `65536 0` (1.0), Q.0 is `0 0`; v3's
* are 65536, 0 and 65536. Each expansion is assembled here followed by
* `;` so that it can be run, and then in the middle of a definition, as
* it would be used:
* Q.1 Q.* ( q -- q ) Q.0 Q.+ ( q -- q )
* Q.1 Q.TO-INT ( -- 1 ) 1 Q.FROM-INT ( -- q ), which must be Q.1 */
#define Q_ONE() do { LIT(65536); LIT(0); } while (0)
#define Q_ZERO() do { LIT(0); LIT(0); } while (0)
#define Q_SCALE() do { LIT(65536); LIT(0); } while (0)
w_q1 = v4_asm_label(&as); Q_ONE(); O(SEMI);
w_q0 = v4_asm_label(&as); Q_ZERO(); O(SEMI);
w_qscale = v4_asm_label(&as); Q_SCALE(); O(SEMI);
w_q1_times = v4_asm_label(&as); Q_ONE(); CALL(w_qstar); O(SEMI);
w_q0_plus = v4_asm_label(&as); Q_ZERO(); CALL(w_dplus); O(SEMI);
w_q1_toint = v4_asm_label(&as); Q_ONE(); CALL(w_qtoint); O(SEMI);
w_one_fromint = v4_asm_label(&as); LIT(1); CALL(w_qfromint); O(SEMI);
#undef Q_SCALE
#undef Q_ZERO
#undef Q_ONE
/* : LSHIFT ( x n -- x<<n ) BEGIN dup WHILE 1- SWAP 2* SWAP REPEAT drop ; 5.5
* : RSHIFT ( x n -- x>>n ) BEGIN dup WHILE 1- SWAP 2/ MSB inv and SWAP REPEAT drop ;
* As written. Section 2 expands BEGIN ... WHILE ... REPEAT as it does IF:
* L0: dup if L1 drop <body> jump L0 L1: drop
* 1- is in line (-1 +); SWAP is a call. */
{
v4_cell l0;
v4_asm_ref l1;
w_lshift = v4_asm_label(&as);
l0 = v4_asm_label(&as);
O(DUP); l1 = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); LIT(-1); O(ADD); CALL(w_swap); O(TWO_STAR); CALL(w_swap);
v4_asm_branch(&as, V4_OP_JUMP, l0);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
O(DROP); O(DROP); O(SEMI);
w_rshift = v4_asm_label(&as);
l0 = v4_asm_label(&as);
O(DUP); l1 = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); LIT(-1); O(ADD); CALL(w_swap);
O(TWO_SLASH); LIT((v4_cell)V4_MSB); O(INV); O(AND); CALL(w_swap);
v4_asm_branch(&as, V4_OP_JUMP, l0);
v4_asm_resolve(&as, l1, v4_asm_label(&as));
O(DROP); O(DROP); O(SEMI);
}
/* Byte access on a word-addressed node (D-1), 5.3: four
* bytes to a cell, little-endian; a byte address is 4 * word address +
* byte index. `@` and `!` here are the opcodes, addressing through A.
* : C@ ( baddr -- c ) call-free
* dup 2/ 2/ a! 3 and k A: word address
* if K0 -1 + if K1 -1 + if K2
* drop @ 23 FOR 2/ UNEXT 255 and ; byte 3
* K2: drop @ 15 FOR 2/ UNEXT 255 and ; byte 2
* K1: drop @ 7 FOR 2/ UNEXT 255 and ; byte 1
* K0: drop @ 255 and ; byte 0
* : C! ( c baddr -- ) call-free
* dup 2/ 2/ a! 3 and push 255 and pop c' k A: word address
* if K0 -1 + if K1 -1 + if K2
* drop 23 FOR 2* UNEXT @ 4278190080 inv and + ! ; byte 3
* K2: drop 15 FOR 2* UNEXT @ -16711681 and + ! ; byte 2
* K1: drop 7 FOR 2* UNEXT @ -65281 and + ! ; byte 1
* K0: drop @ -256 and + ! ; byte 0
* One case per byte position: the byte is shifted up, that byte of the
* cell cleared with a constant mask, and the two added. */
w_cfetch = v4_asm_label(&as);
{
v4_asm_ref f0, f1, f2;
O(DUP); O(TWO_SLASH); O(TWO_SLASH); O(BANG_A); LIT(3); O(AND);
f0 = v4_asm_branch_fwd(&as, V4_OP_IF);
LIT(-1); O(ADD); f1 = v4_asm_branch_fwd(&as, V4_OP_IF);
LIT(-1); O(ADD); f2 = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); O(FETCH_A); LIT(23); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_SLASH); O(UNEXT);
LIT(255); O(AND); O(SEMI);
v4_asm_resolve(&as, f2, v4_asm_label(&as));
O(DROP); O(FETCH_A); LIT(15); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_SLASH); O(UNEXT);
LIT(255); O(AND); O(SEMI);
v4_asm_resolve(&as, f1, v4_asm_label(&as));
O(DROP); O(FETCH_A); LIT(7); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_SLASH); O(UNEXT);
LIT(255); O(AND); O(SEMI);
v4_asm_resolve(&as, f0, v4_asm_label(&as));
O(DROP); O(FETCH_A); LIT(255); O(AND); O(SEMI);
}
{
v4_asm_ref k0, k1, k2;
w_cstore = v4_asm_label(&as);
O(DUP); O(TWO_SLASH); O(TWO_SLASH); O(BANG_A);
LIT(3); O(AND); O(PUSH); LIT(255); O(AND); O(RPOP);
k0 = v4_asm_branch_fwd(&as, V4_OP_IF);
LIT(-1); O(ADD); k1 = v4_asm_branch_fwd(&as, V4_OP_IF);
LIT(-1); O(ADD); k2 = v4_asm_branch_fwd(&as, V4_OP_IF);
O(DROP); LIT(23); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_STAR); O(UNEXT);
O(FETCH_A); LIT((v4_cell)(v4_ucell)0xFF000000u); O(INV); O(AND); O(ADD); O(STORE_A); O(SEMI);
v4_asm_resolve(&as, k2, v4_asm_label(&as));
O(DROP); LIT(15); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_STAR); O(UNEXT);
O(FETCH_A); LIT(-16711681); O(AND); O(ADD); O(STORE_A); O(SEMI);
v4_asm_resolve(&as, k1, v4_asm_label(&as));
O(DROP); LIT(7); O(PUSH);
(void)v4_asm_label(&as);
O(TWO_STAR); O(UNEXT);
O(FETCH_A); LIT(-65281); O(AND); O(ADD); O(STORE_A); O(SEMI);
v4_asm_resolve(&as, k0, v4_asm_label(&as));
O(DROP); O(FETCH_A); LIT(-256); O(AND); O(ADD); O(STORE_A); O(SEMI);
}
/* ---- the rest of 5.6 and 5.7 ---- */
#define FWD(op) v4_asm_branch_fwd(&as, V4_OP_##op)
#define HERE_(r) v4_asm_resolve(&as, (r), v4_asm_label(&as))
#define SWAP_INLINE() do { O(OVER); O(PUSH); O(PUSH); O(DROP); O(RPOP); O(RPOP); } while (0)
{
v4_asm_ref a, b;
/* : M- ( d n -- d ) S>D DNEGATE jump D+
* S>D in line: dup -if P drop -1 jump J P: drop 0 J:
* d - n as d + (-n) with n widened first, so the most negative n is
* subtracted correctly (NEGATE M+ would add it). */
w_mminus = v4_asm_label(&as);
O(DUP); a = FWD(MINUS_IF); O(DROP); LIT(-1); b = FWD(JUMP);
HERE_(a); O(DROP); LIT(0);
HERE_(b);
CALL(w_dnegate); v4_asm_branch(&as, V4_OP_JUMP, w_dplus);
/* : M* ( n1 n2 -- d )
* over over xor push R: sign of the product
* -if A inv 1 + A: push -if B inv 1 + B: pop |n1| |n2|
* UM* pop -if P drop jump DNEGATE P: drop ; */
w_mstar = v4_asm_label(&as);
O(OVER); O(OVER); O(XOR); O(PUSH);
a = FWD(MINUS_IF); O(INV); LIT(1); O(ADD); HERE_(a);
O(PUSH);
a = FWD(MINUS_IF); O(INV); LIT(1); O(ADD); HERE_(a);
O(RPOP);
CALL(w_umstar); O(RPOP); a = FWD(MINUS_IF);
O(DROP); v4_asm_branch(&as, V4_OP_JUMP, w_dnegate);
HERE_(a); O(DROP); O(SEMI);
/* : M/MOD ( d n -- rem quot ) jump SM/REM */
w_mslashmod = v4_asm_label(&as);
v4_asm_branch(&as, V4_OP_JUMP, w_smrem);
/* : MOD ( n1 n2 -- rem ) /MOD drop, /MOD's body in line */
w_mod = v4_asm_label(&as);
O(PUSH); CALL(w_s2d); O(RPOP); CALL(w_smrem); O(DROP); O(SEMI);
// : */MOD ( n1 n2 n3 -- rem quot ) push M* pop jump SM/REM
w_starslashmod = v4_asm_label(&as);
O(PUSH); CALL(w_mstar); O(RPOP); v4_asm_branch(&as, V4_OP_JUMP, w_smrem);
// : */ ( n1 n2 n3 -- quot ) push M* pop SM/REM push drop pop ; that is, */MOD NIP
w_starslash = v4_asm_label(&as);
O(PUSH); CALL(w_mstar); O(RPOP); CALL(w_smrem); O(PUSH); O(DROP); O(RPOP); O(SEMI);
/* : D0< ( d -- flag ) push drop pop jump 0< NIP 0< */
w_d0less = v4_asm_label(&as);
O(PUSH); O(DROP); O(RPOP); v4_asm_branch(&as, V4_OP_JUMP, w_zless);
/* : D2* ( d -- 2d )
* 2* over -if P drop 1 + jump J P: drop J: push 2* pop ;
* the low cell's top bit enters the high cell. */
w_d2star = v4_asm_label(&as);
O(TWO_STAR); O(OVER); a = FWD(MINUS_IF); O(DROP); LIT(1); O(ADD); b = FWD(JUMP);
HERE_(a); O(DROP);
HERE_(b);
O(PUSH); O(TWO_STAR); O(RPOP); O(SEMI);
/* : D2/ ( d -- d/2 ) push a! 0 pop +* push drop a pop ;
* one +* with S = 0 is an exact arithmetic right shift of T:A (as in
* Q.FROM-INT). Clobbers A. */
w_d2slash = v4_asm_label(&as);
O(PUSH); O(BANG_A); LIT(0); O(RPOP); O(MUL_STEP); O(PUSH); O(DROP); O(PUSH_A); O(RPOP); O(SEMI);
/* : 2ROT ( d1 d2 d3 -- d2 d3 d1 ) push push 2SWAP pop pop jump 2SWAP */
w_2rot = v4_asm_label(&as);
O(PUSH); O(PUSH); CALL(w_2swap); O(RPOP); O(RPOP); v4_asm_branch(&as, V4_OP_JUMP, w_2swap);
/* 2DROP is drop drop */
w_2drop = v4_asm_label(&as);
O(DROP); O(DROP); O(SEMI);
/* ( d -- d d ) 2>R 2R@ 2R> in line:
* 2>R is SWAP push push
* 2R@ is pop pop 2DUP push push SWAP
* 2R> is pop pop SWAP */
t_2r = v4_asm_label(&as);
SWAP_INLINE(); O(PUSH); O(PUSH);
O(RPOP); O(RPOP); O(OVER); O(OVER); O(PUSH); O(PUSH); SWAP_INLINE();
O(RPOP); O(RPOP); SWAP_INLINE();
O(SEMI);
/* ( lo hi -- hi lo ) 2>R R> R> the high cell is on top of R, as in v3 */
t_2rorder = v4_asm_label(&as);
SWAP_INLINE(); O(PUSH); O(PUSH); O(RPOP); O(RPOP); O(SEMI);
}
#undef SWAP_INLINE
#undef HERE_
#undef FWD
CHECK(v4_asm_ok(&as), "foundation words assemble");
CHECK(v4_asm_label(&as) <= BYTES, "code stays below the variables");
printf(" code: %ld of %u words\n", (long)v4_asm_label(&as) - 16, (unsigned)V4_NODE_WORDS);
}
/* Call `word` with a canary and up to three arguments on fresh stacks. */
static int call(v4_cell word, unsigned argc, v4_cell a, v4_cell b, v4_cell c)
{
v4_dstack_reset(&n.ds);
v4_rstack_reset(&n.rs);
v4_exec_reset(&es);
v4_heat_reset(&h);
v4_dstack_push(&n.ds, CANARY);
if (argc > 0) v4_dstack_push(&n.ds, a);
if (argc > 1) v4_dstack_push(&n.ds, b);
if (argc > 2) v4_dstack_push(&n.ds, c);
return v4_test_call(&n, &es, &h, word, 4000000) > 0;
}
/* The same with four arguments. */
static int call4(v4_cell word, v4_cell a, v4_cell b, v4_cell c, v4_cell d)
{
v4_dstack_reset(&n.ds);
v4_rstack_reset(&n.rs);
v4_exec_reset(&es);
v4_dstack_push(&n.ds, CANARY);
v4_dstack_push(&n.ds, a);
v4_dstack_push(&n.ds, b);
v4_dstack_push(&n.ds, c);
v4_dstack_push(&n.ds, d);
return v4_test_call(&n, &es, &h, word, 4000000) > 0;
}
/* The results, top first, then the canary. */
static int left1(v4_cell t)
{
return n.ds.t == t && n.ds.s == CANARY;
}
static int left2(v4_cell s, v4_cell t)
{
if (n.ds.t != t || n.ds.s != s) return 0;
(void)v4_dstack_pop(&n.ds);
return n.ds.s == CANARY;
}
static int left3(v4_cell third, v4_cell s, v4_cell t)
{
if (n.ds.t != t) return 0;
(void)v4_dstack_pop(&n.ds);
return left2(third, s);
}
static const v4_cell vec[] = {
0, 1, 2, 3, -1, -2, 12345, -12345,
(v4_cell)(V4_MSB - 1u), (v4_cell)V4_MSB, (v4_cell)(V4_MSB + 1u),
(v4_cell)(V4_MSB >> 1), (v4_cell)((V4_MSB >> 1) - 1u),
(v4_cell)(MAXU / 3u), (v4_cell)(MAXU / 3u * 2u)
};
#define NVEC (sizeof vec / sizeof vec[0])
/* Does UM* as written give the true product? */
static int umstar_exact(v4_ucell u1, v4_ucell u2)
{
v4_ucell lo, hi;
v4_umul(u1, u2, &lo, &hi);
return call(w_umstar, 2, (v4_cell)u1, (v4_cell)u2, 0) && left2((v4_cell)lo, (v4_cell)hi);
}
/* UM/MOD against the identity it must satisfy: q*d + r = uhi:ulo, r < d.
* Checked through v4_umul (itself tested in test_umul.c) rather than against
* a second division routine, so no C divider has to be trusted. */
static int ummod_exact(v4_ucell lo, v4_ucell hi, v4_ucell d)
{
v4_ucell r, q, plo, phi, slo;
if (!call(w_ummod, 3, (v4_cell)lo, (v4_cell)hi, (v4_cell)d)) return 0;
q = (v4_ucell)n.ds.t;
r = (v4_ucell)n.ds.s;
if (!left2((v4_cell)r, (v4_cell)q)) return 0;
v4_umul(q, d, &plo, &phi);
slo = plo + r;
phi += (slo < plo);
return r < d && slo == lo && phi == hi;
}
/* Doubles in C, as (lo, hi) with hi the high cell, for the section 5.7 and
* SM/REM checks. */
static void dadd(v4_ucell al, v4_ucell ah, v4_ucell bl, v4_ucell bh,
v4_ucell *rl, v4_ucell *rh)
{
*rl = al + bl;
*rh = ah + bh + (*rl < al);
}
static void dneg(v4_ucell l, v4_ucell h, v4_ucell *rl, v4_ucell *rh)
{
dadd(~l, ~h, 1u, 0u, rl, rh);
}
/* Signed q * n as a signed double. */
static void smul(v4_cell q, v4_cell nn, v4_ucell *rl, v4_ucell *rh)
{
v4_ucell uq = (v4_ucell)q, un = (v4_ucell)nn;
if (q < 0) uq = 0u - uq;
if (nn < 0) un = 0u - un;
v4_umul(uq, un, rl, rh);
if ((q < 0) != (nn < 0)) dneg(*rl, *rh, rl, rh);
}
/* SM/REM on the dividend d = q*n + r, built so that (r, q) is the
* truncating answer: |r| < |n|, and r is zero or has the sign of d. */
static int smrem_exact(v4_cell q, v4_cell nn, v4_cell r)
{
v4_ucell pl, ph, dl, dh;
smul(q, nn, &pl, &ph);
dadd(pl, ph, (v4_ucell)r, r < 0 ? MAXU : 0u, &dl, &dh);
return call(w_smrem, 3, (v4_cell)dl, (v4_cell)dh, nn) && left2(r, q);
}
/* Six arguments, for 2ROT. */
static int call6(v4_cell word, const v4_cell *a, unsigned dfill, unsigned rfill)
{
unsigned i;
v4_dstack_reset(&n.ds);
v4_rstack_reset(&n.rs);
v4_exec_reset(&es);
for (i = 0; i < dfill; i++) v4_dstack_push(&n.ds, (v4_cell)(0x5A000000 + i));
v4_dstack_push(&n.ds, CANARY);
for (i = 0; i < 6; i++) v4_dstack_push(&n.ds, a[i]);
for (i = 0; i < rfill; i++) v4_rstack_push(&n.rs, (v4_cell)(0x6B000000 + i));
return v4_test_call(&n, &es, &h, word, 4000000) > 0;
}
static int rot6_ok(const v4_cell *a, unsigned dfill, unsigned rfill)
{
static const unsigned from[6] = { 2, 3, 4, 5, 0, 1 }; /* d2 d3 d1 */
unsigned i;
if (!call6(w_2rot, a, dfill, rfill)) return 0;
for (i = 6; i-- > 0; ) if (v4_dstack_pop(&n.ds) != a[from[i]]) return 0;
if (v4_dstack_pop(&n.ds) != CANARY) return 0;
for (i = dfill; i-- > 0; ) if (v4_dstack_pop(&n.ds) != (v4_cell)(0x5A000000 + i)) return 0;
for (i = rfill; i-- > 0; ) if (v4_rstack_pop(&n.rs) != (v4_cell)(0x6B000000 + i)) return 0;
return 1;
}
/* n1 * n2 / n3 through a double product. Returns 1 when the truncated
* quotient fits a cell and (r, q) is it: q*n3 + r = n1*n2, |r| < |n3|, r zero
* or of the product's sign; 0 when it fits and (r, q) is wrong; -1 when it
* does not fit (unspecified, as for SM/REM). */
static int starslash_ref(v4_cell n1, v4_cell n2, v4_cell n3, v4_cell r, v4_cell q)
{
v4_ucell pl, ph, al, ah, ql, qh, sl, sh, un3 = n3 < 0 ? 0u - (v4_ucell)n3 : (v4_ucell)n3, ur;
int pneg;
smul(n1, n2, &pl, &ph);
pneg = (ph & V4_MSB) != 0;
al = pl; ah = ph;
if (pneg) dneg(pl, ph, &al, &ah);
if (ah >= (V4_MSB >> 1)) return -1;
if (((ah << 1) | (al >> (V4_CELL_BITS - 1))) >= un3) return -1; /* |p| >= |n3| * 2^(N-1) */
smul(q, n3, &ql, &qh);
dadd(ql, qh, (v4_ucell)r, r < 0 ? MAXU : 0u, &sl, &sh);
ur = r < 0 ? 0u - (v4_ucell)r : (v4_ucell)r;
return sl == pl && sh == ph && ur < un3 && (r == 0 || (r < 0) == pneg);
}
/* Q48.16 (section 5.26). v3's Q.+ and Q.- are uint64_t a + b and a - b,
* wrapping (v3/include/q48_16.h). Here a Q value is a double: at 32-bit
* cells its two halves, at 64-bit cells the value sign-extended (D-8). */
static void qcells(uint64_t q, v4_cell *lo, v4_cell *hi)
{
#if V4_CELL_BITS == 32
*lo = (v4_cell)(v4_ucell)(q & 0xFFFFFFFFu);
*hi = (v4_cell)(v4_ucell)(q >> 32);
#else
*lo = (v4_cell)(v4_ucell)q;
*hi = (q >> 63) ? V4_ALL_ONES : 0;
#endif
}
/* Run Q.+ (the D+ word) or Q.- (the D- word) on a and b and compare with
* v3's 64-bit wrapping result. At 32-bit cells the double must equal v3's
* value bit for bit. At 64-bit cells the low cell must equal v3's value and
* the high cell must be the true sign of the unwrapped result. */
static int q_matches_v3(v4_cell word, uint64_t a, uint64_t b, uint64_t v3)
{
v4_cell al, ah, bl, bh, rl, rh;
qcells(a, &al, &ah);
qcells(b, &bl, &bh);
if (!call4(word, al, ah, bl, bh)) return 0;
#if V4_CELL_BITS == 32
qcells(v3, &rl, &rh);
return left2(rl, rh);
#else
{
v4_ucell xl, xh, nl, nh;
if (word == w_dminus) { dneg((v4_ucell)bl, (v4_ucell)bh, &nl, &nh); }
else { nl = (v4_ucell)bl; nh = (v4_ucell)bh; }
dadd((v4_ucell)al, (v4_ucell)ah, nl, nh, &xl, &xh);
rl = (v4_cell)(v4_ucell)v3;
rh = (v4_cell)xh;
return (v4_ucell)rl == xl && left2(rl, rh);
}
#endif
}
/* The same for the one-operand Q words: Q.ABS (the DABS word) and Q.NEG
* (the DNEGATE word). `wl`, `wh` is the exact double result, which at
* 64-bit cells may differ from v3 in the high cell only (D-10). */
static int q1_matches_v3(v4_cell word, uint64_t a, uint64_t v3)
{
v4_cell al, ah, rl, rh;
v4_ucell xl, xh;
qcells(a, &al, &ah);
if (!call(word, 2, al, ah, 0)) return 0;
#if V4_CELL_BITS == 32
(void)xl; (void)xh;
qcells(v3, &rl, &rh);
return left2(rl, rh);
#else
if (word == w_dnegate || ah < 0) dneg((v4_ucell)al, (v4_ucell)ah, &xl, &xh);
else { xl = (v4_ucell)al; xh = (v4_ucell)ah; }
rl = (v4_cell)(v4_ucell)v3;
rh = (v4_cell)xh;
return (v4_ucell)rl == xl && left2(rl, rh);
#endif
}
/* x shifted right `k` bits, arithmetic, without C's implementation-defined
* right shift of a negative value. */
static v4_ucell asr_u(v4_ucell x, unsigned k)
{
v4_ucell r = x >> k;
if (k && (x & V4_MSB)) r |= ~(MAXU >> k);
return r;
}
static uint64_t asr64(uint64_t x, unsigned k)
{
uint64_t r = x >> k;
if (k && (x >> 63)) r |= ~((~(uint64_t)0) >> k);
return r;
}
/* The k results, given bottom first, then the canary under them. */
static int leftk(const v4_cell *want, unsigned k)
{
while (k-- > 0)
if (v4_dstack_pop(&n.ds) != want[k]) return 0;
return v4_dstack_pop(&n.ds) == CANARY;
}
/* Signed double a < b, in C. */
static int dlt(v4_ucell al, v4_ucell ah, v4_ucell bl, v4_ucell bh)
{
if (ah != bh) return (v4_cell)ah < (v4_cell)bh;
return al < bl;
}
/* Q48.16 as v4 reads it (D-8: signed). */
static int64_t qs(uint64_t q)
{
return (q >> 63) ? -(int64_t)(~q) - 1 : (int64_t)q;
}
/* Reference Q.* for any cell width, by a different method from the FORTH:
* the Q values (two cells each) as signed integers in 16-bit limbs, the
* magnitudes multiplied schoolbook, the product negated if the signs differ,
* then shifted right 16 (one limb) and cut to two cells. Shifting the
* two's-complement product rounds toward minus infinity. */
#define QLIMBS (2 * V4_CELL_BITS / 16)
static void q_to_limbs(v4_ucell lo, v4_ucell hi, uint32_t *l, int *neg)
{
unsigned i;
*neg = (hi & V4_MSB) != 0;
if (*neg) dneg(lo, hi, &lo, &hi);
for (i = 0; i < QLIMBS / 2; i++) {
l[i] = (uint32_t)((lo >> (16 * i)) & 0xFFFFu);
l[i + QLIMBS / 2] = (uint32_t)((hi >> (16 * i)) & 0xFFFFu);
}
}
static void qmul_ref(v4_ucell al, v4_ucell ah, v4_ucell bl, v4_ucell bh,
v4_ucell *rl, v4_ucell *rh)
{
uint32_t a[QLIMBS], b[QLIMBS], p[2 * QLIMBS];
int an, bn;
unsigned i, j;
q_to_limbs(al, ah, a, &an);
q_to_limbs(bl, bh, b, &bn);
for (i = 0; i < 2 * QLIMBS; i++) p[i] = 0;
for (i = 0; i < QLIMBS; i++) {
uint32_t carry = 0;
for (j = 0; j < QLIMBS; j++) {
uint32_t t = a[i] * b[j] + p[i + j] + carry; /* < 2^32 */
p[i + j] = t & 0xFFFFu;
carry = t >> 16;
}
p[i + QLIMBS] = carry;
}
if (an != bn) { /* negate, 2's complement */
uint32_t carry = 1;
for (i = 0; i < 2 * QLIMBS; i++) {
uint32_t t = (~p[i] & 0xFFFFu) + carry;
p[i] = t & 0xFFFFu;
carry = t >> 16;
}
}
*rl = 0; *rh = 0; /* limbs 1 .. QLIMBS */
for (i = 0; i < QLIMBS / 2; i++) {
*rl |= (v4_ucell)p[1 + i] << (16 * i);
*rh |= (v4_ucell)p[1 + QLIMBS / 2 + i] << (16 * i);
}
}
/* v3's q48_mul: unsigned (a * b) >> 16, cut to 64 bits, as its __int128
* branch computes it on 64-bit gcc/clang hosts. This is the portable
* branch of v3/src/word_source/q48_16_words.c with its last line corrected:
* v3 returns (result_hi << 16) | (result_lo >> 16), which is wrong whenever
* the product reaches past bit 64 (the high cell must move up 48 bits, not
* 16). v3 uses that branch only on hosts without __int128. */
static uint64_t v3_q48_mul(uint64_t a, uint64_t b)
{
uint64_t a_hi = a >> 32;
uint64_t a_lo = a & 0xFFFFFFFFULL;
uint64_t b_hi = b >> 32;
uint64_t b_lo = b & 0xFFFFFFFFULL;
uint64_t p_ll = a_lo * b_lo;
uint64_t p_lh = a_lo * b_hi;
uint64_t p_hl = a_hi * b_lo;
uint64_t p_hh = a_hi * b_hi;
uint64_t carry = 0;
uint64_t mid = p_lh + p_hl;
uint64_t result_hi, result_lo;
if (mid < p_lh) carry++;
result_hi = p_hh + (mid >> 32) + (carry << 32);
result_lo = p_ll + ((mid & 0xFFFFFFFFULL) << 32);
if (result_lo < p_ll) result_hi++;
return (result_hi << 48) | (result_lo >> 16);
}
/* Reference Q./ (D-11), by a different method from the FORTH: magnitudes
* in 16-bit limbs, |a| * 2^16 divided by |b| one bit at a time, then the
* sign applied (rounding toward zero) and the result saturated to Q max or
* Q min. Division by zero gives Q max / Q min by the dividend's sign, 0 for
* 0 / 0, and *err = 1. */
static void qdiv_ref(v4_ucell al, v4_ucell ah, v4_ucell bl, v4_ucell bh,
v4_ucell *rl, v4_ucell *rh, int *err)
{
uint32_t a[QLIMBS], b[QLIMBS], num[QLIMBS + 1], quo[QLIMBS + 1], rem[QLIMBS + 1];
int an, bn, neg, big = 0;
unsigned i, k;
*err = 0;
if (bl == 0 && bh == 0) {
*err = 1;
if (al == 0 && ah == 0) { *rl = 0; *rh = 0; }
else if (ah & V4_MSB) { *rl = 0; *rh = V4_MSB; }
else { *rl = MAXU; *rh = MAXU >> 1; }
return;
}
q_to_limbs(al, ah, a, &an);
q_to_limbs(bl, bh, b, &bn);
neg = an != bn;
num[0] = 0; /* |a| * 2^16 */
for (i = 0; i < QLIMBS; i++) num[i + 1] = a[i];
for (i = 0; i <= QLIMBS; i++) { quo[i] = 0; rem[i] = 0; }
for (k = 16 * (QLIMBS + 1); k-- > 0; ) {
uint32_t bit = (num[k / 16] >> (k % 16)) & 1u, c = bit;
int ge = 1;
for (i = 0; i <= QLIMBS; i++) { /* rem = 2 rem + bit */
uint32_t t = (rem[i] << 1) | c;
rem[i] = t & 0xFFFFu;
c = t >> 16;
}
for (i = QLIMBS + 1; i-- > 0; ) { /* rem >= |b| ? */
uint32_t bv = i < QLIMBS ? b[i] : 0;
if (rem[i] != bv) { ge = rem[i] > bv; break; }
}
if (ge) {
uint32_t br = 0;
for (i = 0; i <= QLIMBS; i++) {
uint32_t bv = i < QLIMBS ? b[i] : 0;
uint32_t t = rem[i] - bv - br;
br = (t >> 16) & 1u;
rem[i] = t & 0xFFFFu;
}
quo[k / 16] |= 1u << (k % 16);
}
}
if (quo[QLIMBS] || (quo[QLIMBS - 1] & 0x8000u)) big = 1; /* >= 2^(2N-1) */
if (big) {
if (neg) { *rl = 0; *rh = V4_MSB; }
else { *rl = MAXU; *rh = MAXU >> 1; }
return;
}
*rl = 0; *rh = 0;
for (i = 0; i < QLIMBS / 2; i++) {
*rl |= (v4_ucell)quo[i] << (16 * i);
*rh |= (v4_ucell)quo[QLIMBS / 2 + i] << (16 * i);
}
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;
}
/* 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;
}
/* v3's q48_log_approx, verbatim but for the names. */
static uint64_t v3_q48_log(uint64_t x)
{
const uint64_t LN2_Q48 = 45426;
int k = 0, iter;
uint64_t m = x, y, result;
if (x == 0) return 0;
if (x == 65536) return 0;
if (m >= 131072) {
while (m >= 131072) { m >>= 1; k++; }
} else if (m < 65536) {
while (m < 65536) { m <<= 1; k--; }
}
y = (m > 65536) ? (m - 65536) : 0;
for (iter = 0; iter < 6; iter++) {
uint64_t exp_y = v3_q48_exp(y), correction;
int64_t delta_signed;
if (exp_y == 0) break;
delta_signed = (int64_t)m - (int64_t)exp_y;
if (delta_signed > 0) {
correction = v3_q48_div((uint64_t)delta_signed, exp_y);
y = y + correction;
} else if (delta_signed < 0) {
correction = v3_q48_div((uint64_t)(-delta_signed), exp_y);
if (y > correction) y = y - correction;
else y = 0;
}
if (delta_signed < 100 && delta_signed > -100) break;
}
result = y;
if (k > 0) result = result + v3_q48_mul((uint64_t)k << 16, LN2_Q48);
else if (k < 0) result = result - v3_q48_mul((uint64_t)(-k) << 16, LN2_Q48);
return result;
}
/* v3's q48_reduce_angle and q48_sin_approx, verbatim but for the names. */
static int64_t v3_q48_reduce(int64_t x)
{
const int64_t PI_Q48 = 205887;
const int64_t TWO_PI_Q48 = 2 * PI_Q48;
int64_t k = x / TWO_PI_Q48;
x -= k * TWO_PI_Q48;
while (x > PI_Q48) x -= TWO_PI_Q48;
while (x < -PI_Q48) x += TWO_PI_Q48;
return x;
}
static uint64_t v3_q48_sin(uint64_t q)
{
int64_t x = v3_q48_reduce((int64_t)q);
int negative = (x < 0), subtract = 1, nn;
uint64_t xu, x2, term, result;
if (negative) x = -x;
xu = (uint64_t)x;
x2 = v3_q48_mul(xu, xu);
term = xu;
result = xu;
for (nn = 3; nn <= 11; nn += 2) {
term = v3_q48_mul(term, x2);
term = v3_q48_div(term, (uint64_t)(nn * (nn - 1)) << 16);
result = subtract ? result - term : result + term;
subtract = !subtract;
if (term < 10) break;
}
return negative ? (uint64_t)0 - result : result;
}
/* v3's q48_cos_approx, verbatim but for the names. */
static uint64_t v3_q48_cos(uint64_t q)
{
int64_t x = v3_q48_reduce((int64_t)q);
uint64_t xu = (uint64_t)(x < 0 ? -x : x);
uint64_t x2 = v3_q48_mul(xu, xu);
uint64_t term = 65536, result = term;
int subtract = 1, nn;
for (nn = 2; nn <= 10; nn += 2) {
term = v3_q48_mul(term, x2);
term = v3_q48_div(term, (uint64_t)(nn * (nn - 1)) << 16);
result = subtract ? result - term : result + term;
subtract = !subtract;
if (term < 10) break;
}
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
* are circular (D-2), so a word that needs more depth than is free does not
* fault: it silently overwrites the oldest cells, and this is how that shows. */
#define DFILL(i) ((v4_cell)(0x5A000000 + (i)))
#define RFILL(i) ((v4_cell)(0x6B000000 + (i)))
static int fits(v4_cell word, unsigned argc, const v4_cell *arg,
unsigned nres, unsigned dfill, unsigned rfill)
{
v4_cell want[V4_DATA_DEPTH];
unsigned i;
v4_dstack_reset(&n.ds);
v4_rstack_reset(&n.rs);
v4_exec_reset(&es);
for (i = 0; i < argc; i++) v4_dstack_push(&n.ds, arg[i]);
if (v4_test_call(&n, &es, &h, word, 4000000) <= 0) return 0;
if (nres > V4_DATA_DEPTH) return 0;
for (i = 0; i < nres; i++) want[i] = v4_dstack_pop(&n.ds);
v4_dstack_reset(&n.ds);
v4_rstack_reset(&n.rs);
v4_exec_reset(&es);
for (i = 0; i < dfill; i++) v4_dstack_push(&n.ds, DFILL(i));
v4_dstack_push(&n.ds, CANARY);
for (i = 0; i < argc; i++) v4_dstack_push(&n.ds, arg[i]);
for (i = 0; i < rfill; i++) v4_rstack_push(&n.rs, RFILL(i));
if (v4_test_call(&n, &es, &h, word, 4000000) <= 0) return 0;
for (i = 0; i < nres; i++) if (v4_dstack_pop(&n.ds) != want[i]) return 0;
if (v4_dstack_pop(&n.ds) != CANARY) return 0;
for (i = dfill; i-- > 0; ) if (v4_dstack_pop(&n.ds) != DFILL(i)) return 0;
for (i = rfill; i-- > 0; ) if (v4_rstack_pop(&n.rs) != RFILL(i)) return 0;
return 1;
}
/* The most marked cells that survive on each stack, over a set of argument
* tuples chosen to take every branch. The data figure counts cells below the
* canary, so the caller may hold canary + headroom cells under the arguments. */
static void headroom(const char *name, v4_cell word, unsigned argc,
const v4_cell (*args)[4], unsigned nargs, unsigned nres,
int *dh, int *rh)
{
int d, r;
unsigned k;
for (d = 0; d < V4_DATA_DEPTH; d++) {
for (k = 0; k < nargs; k++) if (!fits(word, argc, args[k], nres, (unsigned)d + 1u, 0)) break;
if (k < nargs) break;
}
for (r = 0; r < V4_RET_DEPTH; r++) {
for (k = 0; k < nargs; k++) if (!fits(word, argc, args[k], nres, 0, (unsigned)r + 1u)) break;
if (k < nargs) break;
}
*dh = d; *rh = r;
printf(" %s headroom: data %d below canary, return %d below its return address\n",
name, d, r);
}
int main(void)
{
printf("v4 foundation tests: V4_CELL_BITS=%d\n", V4_CELL_BITS);
build();
for (unsigned i = 0; i < NVEC; i++) {
v4_cell a = vec[i];
v4_ucell ua = (v4_ucell)a;
CHECK(call(w_negate, 1, a, 0, 0) && left1((v4_cell)(0u - ua)), "NEGATE [%u]", i);
CHECK(call(w_zless, 1, a, 0, 0) && left1(FLAG(ua & V4_MSB)), "0< [%u]", i);
CHECK(call(w_zequal, 1, a, 0, 0) && left1(FLAG(a == 0)), "0= [%u]", i);
for (unsigned j = 0; j < NVEC; j++) {
v4_cell b = vec[j];
v4_ucell ub = (v4_ucell)b;
CHECK(call(w_nip, 2, a, b, 0) && left1(b), "NIP [%u,%u]", i, j);
CHECK(call(w_swap, 2, a, b, 0) && left2(b, a), "SWAP [%u,%u]", i, j);
CHECK(call(w_or, 2, a, b, 0) && left1((v4_cell)(ua | ub)), "OR [%u,%u]", i, j);
CHECK(call(w_2dup, 2, a, b, 0) && n.ds.t == b && n.ds.s == a
&& (v4_dstack_pop(&n.ds), v4_dstack_pop(&n.ds), left2(a, b)),
"2DUP [%u,%u]", i, j);
CHECK(call(w_minus, 2, a, b, 0) && left1((v4_cell)(ua - ub)), "- [%u,%u]", i, j);
CHECK(call(w_uless, 2, a, b, 0) && left1(FLAG(ua < ub)), "U< [%u,%u]", i, j);
for (unsigned k = 0; k < NVEC; k++) {
v4_cell c = vec[k];
CHECK(call(w_rot, 3, a, b, c) && left3(b, c, a), "ROT [%u,%u,%u]", i, j, k);
}
CHECK(umstar_exact(ua, ub), "UM* [%u,%u]", i, j);
}
}
/* UM* over the full range: products of pseudo-random operands, and of
* operands chosen near the edges D-3 makes dangerous (top bit set, all
* ones, just under a power of two). */
{
v4_ucell x = (v4_ucell)0x9E3779B9u;
for (unsigned i = 0; i < 20000; i++) {
v4_ucell u1, u2;
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
u1 = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
u2 = x;
switch (i & 3u) {
case 1: u1 |= V4_MSB; break;
case 2: u1 |= V4_MSB; u2 |= V4_MSB; break;
case 3: u1 = MAXU - (u1 & 7u); break;
default: break;
}
CHECK(umstar_exact(u1, u2), "UM* random [%u]", i);
}
}
/* The two cases that broke UM* as first written in section 4. */
CHECK(umstar_exact(MAXU, MAXU), "UM* MAX*MAX: carry out of T (D-3)");
CHECK(umstar_exact(V4_MSB - 1u, 3u), "UM* (2^(n-1)-1)*3");
/* UM/MOD over its defined range, uhi < ud. */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++)
for (unsigned k = 0; k < NVEC; k++) {
v4_ucell lo = (v4_ucell)vec[i], hi = (v4_ucell)vec[j], d = (v4_ucell)vec[k];
if (hi < d)
CHECK(ummod_exact(lo, hi, d), "UM/MOD [%u,%u,%u]", i, j, k);
}
{
v4_ucell x = (v4_ucell)0x2545F491u;
for (unsigned i = 0; i < 20000; i++) {
v4_ucell lo, hi, d;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; lo = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; d = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; hi = x;
switch (i & 3u) {
case 1: d |= V4_MSB; break; /* top bit of d set */
case 2: d >>= (unsigned)(x & 31u); break; /* small divisors */
case 3: d = MAXU - (d & 7u); break; /* d just under 2^n */
default: break;
}
if (d == 0) d = 1;
hi %= d; /* uhi < ud */
CHECK(ummod_exact(lo, hi, d), "UM/MOD random [%u]", i);
}
}
/* Section 5 words under SM/REM and /MOD, against their C meaning. */
for (unsigned i = 0; i < NVEC; i++) {
v4_cell a = vec[i];
v4_ucell ua = (v4_ucell)a;
CHECK(call(w_abs, 1, a, 0, 0) && left1(a < 0 ? (v4_cell)(0u - ua) : a), "ABS [%u]", i);
CHECK(call(w_s2d, 1, a, 0, 0) && left2(a, FLAG(a < 0)), "S>D [%u]", i);
for (unsigned j = 0; j < NVEC; j++) {
v4_cell b = vec[j];
v4_ucell ub = (v4_ucell)b, rl, rh;
CHECK(call(w_ugreater, 2, a, b, 0) && left1(FLAG(ua > ub)), "U> [%u,%u]", i, j);
dneg(ua, ub, &rl, &rh);
CHECK(call(w_dnegate, 2, a, b, 0) && left2((v4_cell)rl, (v4_cell)rh), "DNEGATE [%u,%u]", i, j);
if (b >= 0) { rl = ua; rh = ub; }
CHECK(call(w_dabs, 2, a, b, 0) && left2((v4_cell)rl, (v4_cell)rh), "DABS [%u,%u]", i, j);
for (unsigned k = 0; k < NVEC; k++)
for (unsigned l = 0; l < NVEC; l++) {
v4_ucell sl, sh;
dadd(ua, ub, (v4_ucell)vec[k], (v4_ucell)vec[l], &sl, &sh);
CHECK(call4(w_dplus, a, b, vec[k], vec[l])
&& left2((v4_cell)sl, (v4_cell)sh), "D+ [%u,%u,%u,%u]", i, j, k, l);
}
if (b != 0 && !(a == (v4_cell)V4_MSB && b == -1))
CHECK(call(w_slashmod, 2, a, b, 0) && left2(a % b, a / b), "/MOD [%u,%u]", i, j);
if (b != 0 && !(a == (v4_cell)V4_MSB && b == -1))
CHECK(call(w_slash, 2, a, b, 0) && left1(a / b), "/ [%u,%u]", i, j);
CHECK(call(w_star, 2, a, b, 0) && left1((v4_cell)(ua * ub)), "* [%u,%u]", i, j);
}
}
/* The rest of 5.6 and 5.7. "v3:" marks results recorded from the v3
* binary on 2026-10-03; v3's M- and M/MOD took the double with its low
* cell on top, so their arguments are in v4's order here. */
CHECK(call(w_mstar, 2, 6, 7, 0) && left2(42, 0), "v3: 6 7 M*");
CHECK(call(w_mstar, 2, -6, 7, 0) && left2(-42, -1), "v3: -6 7 M*");
CHECK(call(w_mod, 2, 17, 5, 0) && left1(2), "v3: 17 5 MOD");
CHECK(call(w_mod, 2, -17, 5, 0) && left1(-2), "v3: -17 5 MOD");
CHECK(call(w_mod, 2, 17, -5, 0) && left1(2), "v3: 17 -5 MOD");
CHECK(call(w_starslash, 3, 7, 3, 2) && left1(10), "v3: 7 3 2 */");
CHECK(call(w_starslash, 3, -7, 3, 2) && left1(-10), "v3: -7 3 2 */");
CHECK(call(w_starslashmod, 3, 7, 3, 2) && left2(1, 10), "v3: 7 3 2 */MOD");
CHECK(call(w_starslashmod, 3, -7, 3, 2) && left2(-1, -10), "v3: -7 3 2 */MOD");
CHECK(call(w_d0less, 2, 5, 0, 0) && left1(0), "v3: 5 0 D0<");
CHECK(call(w_d0less, 2, 5, -1, 0) && left1(-1), "v3: 5 -1 D0<");
CHECK(call(w_d2star, 2, 3, 0, 0) && left2(6, 0), "v3: 3 0 D2*");
CHECK(call(w_d2star, 2, -1, 0, 0) && left2(-2, 1), "v3: -1 0 D2*");
CHECK(call(w_d2slash, 2, 6, 0, 0) && left2(3, 0), "v3: 6 0 D2/");
CHECK(call(w_d2slash, 2, 1, 1, 0) && left2((v4_cell)V4_MSB, 0), "v3: 1 1 D2/");
CHECK(call(w_d2slash, 2, -4, -1, 0) && left2(-2, -1), "v3: -4 -1 D2/");
CHECK(call(w_mslashmod, 3, 100, 0, 7) && left2(2, 14), "v3: 100 7 M/MOD");
CHECK(call(w_mslashmod, 3, -100, -1, 7) && left2(-2, -14), "v3: -100 7 M/MOD");
CHECK(call(w_mminus, 3, 100, 0, 7) && left2(93, 0), "v3: 100 7 M-");
CHECK(call(w_mminus, 3, 100, 0, -7) && left2(107, 0), "v3: 100 -7 M-");
CHECK(call(w_2drop, 3, 1, 2, 3) && left1(1), "v3: 1 2 3 2DROP");
CHECK(call(t_2r, 2, 1, 2, 0) && v4_dstack_pop(&n.ds) == 2 && v4_dstack_pop(&n.ds) == 1 && left2(1, 2),
"v3: 1 2 2>R 2R@ 2R>");
CHECK(call(t_2rorder, 2, 1, 2, 0) && left2(2, 1), "2>R leaves the high cell on top of R");
{
static const v4_cell six[6] = { 1, 2, 3, 4, 5, 6 };
CHECK(rot6_ok(six, 0, 0), "v3: 1 2 3 4 5 6 2ROT");
}
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++) {
v4_cell a = vec[i], b = vec[j];
v4_ucell ua = (v4_ucell)a, ub = (v4_ucell)b, pl, ph;
smul(a, b, &pl, &ph);
CHECK(call(w_mstar, 2, a, b, 0) && left2((v4_cell)pl, (v4_cell)ph), "M* [%u,%u]", i, j);
CHECK(call(w_d0less, 2, a, b, 0) && left1(FLAG(b < 0)), "D0< [%u,%u]", i, j);
CHECK(call(w_d2star, 2, a, b, 0)
&& left2((v4_cell)(ua << 1), (v4_cell)((ub << 1) | (ua >> (V4_CELL_BITS - 1)))), "D2* [%u,%u]", i, j);
CHECK(call(w_d2slash, 2, a, b, 0)
&& left2((v4_cell)((ua >> 1) | (ub << (V4_CELL_BITS - 1))), (v4_cell)((ub >> 1) | (ub & V4_MSB))),
"D2/ [%u,%u]", i, j);
CHECK(call(w_2drop, 3, 77, a, b) && left1(77), "2DROP [%u,%u]", i, j);
CHECK(call(t_2r, 2, a, b, 0) && v4_dstack_pop(&n.ds) == b && v4_dstack_pop(&n.ds) == a && left2(a, b),
"2>R 2R@ 2R> [%u,%u]", i, j);
if (b != 0 && !(a == (v4_cell)V4_MSB && b == -1)) {
CHECK(call(w_mod, 2, a, b, 0) && left1(a % b), "MOD [%u,%u]", i, j);
CHECK(call(w_mslashmod, 3, a, a < 0 ? -1 : 0, b) && left2(a % b, a / b), "M/MOD [%u,%u]", i, j);
}
for (unsigned k = 0; k < NVEC; k++) {
v4_cell c = vec[k], r, q;
int ok;
{
v4_cell six[6];
six[0] = a; six[1] = b; six[2] = c; six[3] = vec[(i + 5) % NVEC]; six[4] = vec[(j + 7) % NVEC]; six[5] = vec[(k + 3) % NVEC];
CHECK(rot6_ok(six, 0, 0), "2ROT [%u,%u,%u]", i, j, k);
}
if (c == 0) continue;
if (!call(w_starslashmod, 3, a, b, c)) { CHECK(0, "*/MOD returns [%u,%u,%u]", i, j, k); continue; }
q = n.ds.t; r = n.ds.s;
ok = starslash_ref(a, b, c, r, q);
CHECK(ok != 0, "*/MOD [%u,%u,%u]", i, j, k);
if (ok == 1) {
CHECK(left2(r, q), "*/MOD leaves two cells [%u,%u,%u]", i, j, k);
CHECK(call(w_starslash, 3, a, b, c) && left1(q), "*/ [%u,%u,%u]", i, j, k);
}
}
}
/* M/MOD is SM/REM: a double that does not fit a cell */
CHECK(call(w_mslashmod, 3, 0, 5, 10) && n.ds.s == 0
&& (v4_ucell)n.ds.t == (v4_ucell)1 << (V4_CELL_BITS - 1), "M/MOD of 5 * 2^N by 10");
/* SM/REM: every quotient, divisor and in-range remainder sign drawn from
* the edge vectors, then pseudo-random. */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++) {
v4_cell q = vec[i], nn = vec[j];
v4_ucell un = nn < 0 ? 0u - (v4_ucell)nn : (v4_ucell)nn;
int neg = (q < 0) != (nn < 0);
if (nn == 0) continue;
CHECK(smrem_exact(q, nn, 0), "SM/REM r=0 [%u,%u]", i, j);
if (un > 1u) {
v4_cell r = (v4_cell)(un - 1u);
if (q == 0) {
CHECK(smrem_exact(q, nn, r), "SM/REM r>0 q=0 [%u,%u]", i, j);
CHECK(smrem_exact(q, nn, (v4_cell)(0u - (v4_ucell)r)), "SM/REM r<0 q=0 [%u,%u]", i, j);
} else {
CHECK(smrem_exact(q, nn, neg ? (v4_cell)(0u - (v4_ucell)r) : r),
"SM/REM r=max [%u,%u]", i, j);
}
}
}
{
v4_ucell x = (v4_ucell)0x6C078965u;
for (unsigned i = 0; i < 20000; i++) {
v4_cell q, nn, r;
v4_ucell un, ur;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; q = (v4_cell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; nn = (v4_cell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; ur = x;
if (i & 1u) nn = (v4_cell)((v4_ucell)nn >> (x & 31u)); /* small divisors */
if (i & 2u) q = (v4_cell)((v4_ucell)q >> (x & 31u)); /* small quotients */
if (nn == 0) nn = 1;
un = nn < 0 ? 0u - (v4_ucell)nn : (v4_ucell)nn;
r = (v4_cell)(ur % un);
if (q == 0 ? (x & 4u) : ((q < 0) != (nn < 0))) r = (v4_cell)(0u - (v4_ucell)r);
CHECK(smrem_exact(q, nn, r), "SM/REM random [%u]", i);
}
}
/* D+ on pseudo-random doubles, weighted toward low cells whose top bits
* decide the carry. */
{
v4_ucell x = (v4_ucell)0x41C64E6Du, v[4], sl, sh;
for (unsigned i = 0; i < 20000; i++) {
for (unsigned k = 0; k < 4; k++) { x ^= x << 13; x ^= x >> 7; x ^= x << 17; v[k] = x; }
if (i & 1u) v[0] |= V4_MSB;
if (i & 2u) v[2] |= V4_MSB;
if (i & 4u) v[2] = 0u - v[0]; /* low sum wraps to 0 */
dadd(v[0], v[1], v[2], v[3], &sl, &sh);
CHECK(call4(w_dplus, (v4_cell)v[0], (v4_cell)v[1], (v4_cell)v[2], (v4_cell)v[3])
&& left2((v4_cell)sl, (v4_cell)sh), "D+ random [%u]", i);
}
}
/* M+, D-, D0= and D= against C: every edge-vector combination, then
* pseudo-random doubles. */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++) {
v4_ucell ua = (v4_ucell)vec[i], ub = (v4_ucell)vec[j], rl, rh;
CHECK(call(w_d0equal, 2, vec[i], vec[j], 0) && left1(FLAG(ua == 0 && ub == 0)),
"D0= [%u,%u]", i, j);
for (unsigned k = 0; k < NVEC; k++) {
v4_cell c = vec[k];
dadd(ua, ub, (v4_ucell)c, c < 0 ? MAXU : 0u, &rl, &rh);
CHECK(call(w_mplus, 3, vec[i], vec[j], c) && left2((v4_cell)rl, (v4_cell)rh),
"M+ [%u,%u,%u]", i, j, k);
{
v4_ucell nl, nh, dl, dh;
dneg((v4_ucell)c, c < 0 ? MAXU : 0u, &nl, &nh);
dadd(ua, ub, nl, nh, &dl, &dh);
CHECK(call(w_mminus, 3, vec[i], vec[j], c) && left2((v4_cell)dl, (v4_cell)dh),
"M- [%u,%u,%u]", i, j, k);
}
for (unsigned l = 0; l < NVEC; l++) {
v4_ucell nl, nh;
dneg((v4_ucell)vec[k], (v4_ucell)vec[l], &nl, &nh);
dadd(ua, ub, nl, nh, &rl, &rh);
CHECK(call4(w_dminus, vec[i], vec[j], vec[k], vec[l])
&& left2((v4_cell)rl, (v4_cell)rh), "D- [%u,%u,%u,%u]", i, j, k, l);
CHECK(call4(w_dequal, vec[i], vec[j], vec[k], vec[l])
&& left1(FLAG(i == k && j == l)), "D= [%u,%u,%u,%u]", i, j, k, l);
}
}
}
{
v4_ucell x = (v4_ucell)0x5851F42Du, v[4], rl, rh, nl, nh;
for (unsigned i = 0; i < 20000; i++) {
for (unsigned k = 0; k < 4; k++) { x ^= x << 13; x ^= x >> 7; x ^= x << 17; v[k] = x; }
if (i & 1u) v[2] = v[0]; /* equal low cells */
if (i & 2u) v[3] = v[1]; /* equal high cells */
if (i & 4u) v[2] = 0u - v[0]; /* low sum wraps */
dadd(v[0], v[1], v[2], (v4_cell)v[2] < 0 ? MAXU : 0u, &rl, &rh);
CHECK(call(w_mplus, 3, (v4_cell)v[0], (v4_cell)v[1], (v4_cell)v[2])
&& left2((v4_cell)rl, (v4_cell)rh), "M+ random [%u]", i);
dneg(v[2], v[3], &nl, &nh);
dadd(v[0], v[1], nl, nh, &rl, &rh);
CHECK(call4(w_dminus, (v4_cell)v[0], (v4_cell)v[1], (v4_cell)v[2], (v4_cell)v[3])
&& left2((v4_cell)rl, (v4_cell)rh), "D- random [%u]", i);
CHECK(call4(w_dequal, (v4_cell)v[0], (v4_cell)v[1], (v4_cell)v[2], (v4_cell)v[3])
&& left1(FLAG(v[0] == v[2] && v[1] == v[3])), "D= random [%u]", i);
CHECK(call(w_d0equal, 2, (v4_cell)(v[0] & v[2]), (v4_cell)(v[1] & v[3]), 0)
&& left1(FLAG((v[0] & v[2]) == 0 && (v[1] & v[3]) == 0)), "D0= random [%u]", i);
}
}
/* Q.+ and Q.- (section 5.26: the D+ and D- words) against v3. */
{
static const uint64_t qv[] = {
0, 1, 0x8000u, 0x10000u, 0x18000u, /* 0, ulp, 0.5, 1.0, 1.5 */
(uint64_t)0 - 0x10000u, (uint64_t)0 - 1u, /* -1.0, -ulp */
0xFFFFu, 0x10001u, 0xFFFFFFFFu, 0x100000000u, /* across the 32-bit seam */
((uint64_t)1 << 63) - 1u, (uint64_t)1 << 63, /* Q max, Q min */
(uint64_t)12345 << 16, (uint64_t)0 - ((uint64_t)12345 << 16)
};
const unsigned nq = sizeof qv / sizeof qv[0];
uint64_t x = 0x9E3779B97F4A7C15u;
unsigned over = 0;
for (unsigned i = 0; i < nq; i++)
for (unsigned j = 0; j < nq; j++) {
CHECK(q_matches_v3(w_dplus, qv[i], qv[j], qv[i] + qv[j]), "Q.+ [%u,%u]", i, j);
CHECK(q_matches_v3(w_dminus, qv[i], qv[j], qv[i] - qv[j]), "Q.- [%u,%u]", i, j);
}
for (unsigned i = 0; i < 20000; i++) {
uint64_t a, b;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; a = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; b = x;
if (i & 1u) a >>= (unsigned)(x & 63u); /* mixed magnitudes */
if (i & 2u) b = (uint64_t)0 - b;
if (((a + b) ^ a) & ((a + b) ^ b) & ((uint64_t)1 << 63)) over++;
CHECK(q_matches_v3(w_dplus, a, b, a + b), "Q.+ random [%u]", i);
CHECK(q_matches_v3(w_dminus, a, b, a - b), "Q.- random [%u]", i);
}
printf(" Q.+ random cases that overflow Q48.16: %u of 20000\n", over);
/* Q.FROM-INT: n * 2^16 as a signed double. v3 clamps n < 0 to 0
* (section 5.26 retires that) and wraps at 64 bits, so v3's value is
* checked only for n >= 0, and at 64-bit cells only in the low cell
* (D-10). */
for (unsigned i = 0; i < NVEC + 20000u; i++) {
v4_cell nn, el, eh;
if (i < NVEC) nn = vec[i];
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
nn = (v4_cell)(v4_ucell)x;
if (i & 1u) nn = (v4_cell)asr_u((v4_ucell)nn, (unsigned)(x >> 58) % V4_CELL_BITS);
}
el = (v4_cell)((v4_ucell)nn << 16);
eh = (v4_cell)asr_u((v4_ucell)nn, V4_CELL_BITS - 16);
CHECK(call(w_qfromint, 1, nn, 0, 0) && left2(el, eh), "Q.FROM-INT [%u]", i);
if (nn >= 0) {
uint64_t v3 = (uint64_t)(v4_ucell)nn << 16;
v4_cell vl, vh;
qcells(v3, &vl, &vh);
CHECK(call(w_qfromint, 1, nn, 0, 0) && n.ds.s == vl
&& (V4_CELL_BITS == 64 || n.ds.t == vh), "Q.FROM-INT v3 [%u]", i);
}
/* and back */
CHECK(call(w_qfromint, 1, nn, 0, 0)
&& (v4_test_call(&n, &es, &h, w_qtoint, 100000) > 0) && left1(nn),
"Q.TO-INT Q.FROM-INT [%u]", i);
}
/* Q.TO-INT: v3's (int64_t)q >> 16, arithmetic. At 32-bit cells the
* result is its low cell; at 64-bit cells it is v3's value exactly. */
for (unsigned i = 0; i < nq + 20000u; i++) {
uint64_t q;
v4_cell ql, qh;
if (i < nq) q = qv[i];
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17; q = x;
if (i & 1u) q = asr64(q, (unsigned)(x >> 58));
}
qcells(q, &ql, &qh);
CHECK(call(w_qtoint, 2, ql, qh, 0) && left1((v4_cell)(v4_ucell)asr64(q, 16)),
"Q.TO-INT v3 [%u]", i);
}
/* and on any double, not only sign-extended ones */
for (unsigned i = 0; i < 20000; i++) {
v4_ucell lo, hi;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; lo = (v4_ucell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; hi = (v4_ucell)x;
CHECK(call(w_qtoint, 2, (v4_cell)lo, (v4_cell)hi, 0)
&& left1((v4_cell)((lo >> 16) | (hi << (V4_CELL_BITS - 16)))),
"Q.TO-INT double [%u]", i);
}
/* Q.ABS and Q.NEG against v3's q48_abs and 0 - q. */
for (unsigned i = 0; i < nq; i++) {
uint64_t q = qv[i];
CHECK(q1_matches_v3(w_dabs, q, (q < ((uint64_t)1 << 63)) ? q : (uint64_t)0 - q),
"Q.ABS [%u]", i);
CHECK(q1_matches_v3(w_dnegate, q, (uint64_t)0 - q), "Q.NEG [%u]", i);
}
for (unsigned i = 0; i < 20000; i++) {
uint64_t q;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; q = x;
if (i & 1u) q >>= (unsigned)(x & 63u);
if (i & 2u) q = (uint64_t)0 - q;
CHECK(q1_matches_v3(w_dabs, q, (q < ((uint64_t)1 << 63)) ? q : (uint64_t)0 - q),
"Q.ABS random [%u]", i);
CHECK(q1_matches_v3(w_dnegate, q, (uint64_t)0 - q), "Q.NEG random [%u]", i);
}
}
/* <, =, D<, 2SWAP, 2OVER, DMAX and DMIN as written in section 5, against C. */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++) {
v4_cell a = vec[i], b = vec[j];
CHECK(call(w_less, 2, a, b, 0) && left1(FLAG(a < b)), "< [%u,%u]", i, j);
CHECK(call(w_equal, 2, a, b, 0) && left1(FLAG(a == b)), "= [%u,%u]", i, j);
for (unsigned k = 0; k < NVEC; k++)
for (unsigned l = 0; l < NVEC; l++) {
v4_cell c = vec[k], d = vec[l];
int lt = dlt((v4_ucell)a, (v4_ucell)b, (v4_ucell)c, (v4_ucell)d);
v4_cell sw[4], ov[6], mx[2], mn[2];
sw[0] = c; sw[1] = d; sw[2] = a; sw[3] = b;
ov[0] = a; ov[1] = b; ov[2] = c; ov[3] = d; ov[4] = a; ov[5] = b;
mx[0] = lt ? c : a; mx[1] = lt ? d : b;
mn[0] = lt ? a : c; mn[1] = lt ? b : d;
CHECK(call4(w_dless, a, b, c, d) && left1(FLAG(lt)), "D< [%u,%u,%u,%u]", i, j, k, l);
{
v4_cell kp[5];
kp[0] = a; kp[1] = b; kp[2] = c; kp[3] = d; kp[4] = FLAG(lt);
CHECK(call4(w_dltkeep, a, b, c, d) && leftk(kp, 5), "(D<) [%u,%u,%u,%u]", i, j, k, l);
}
CHECK(call4(w_2swap, a, b, c, d) && leftk(sw, 4), "2SWAP [%u,%u,%u,%u]", i, j, k, l);
CHECK(call4(w_2over, a, b, c, d) && leftk(ov, 6), "2OVER [%u,%u,%u,%u]", i, j, k, l);
CHECK(call4(w_dmax, a, b, c, d) && leftk(mx, 2), "DMAX [%u,%u,%u,%u]", i, j, k, l);
CHECK(call4(w_dmin, a, b, c, d) && leftk(mn, 2), "DMIN [%u,%u,%u,%u]", i, j, k, l);
}
}
/* Q.=, Q.<, Q.>, Q.0=, Q.MAX, Q.MIN (section 5.26: D=, D<, 2SWAP D<, D0=,
* DMAX, DMIN) on Q values, signed (D-8). v3 compares Q values unsigned,
* so v3 agrees only when both have the same sign; that is checked too. */
{
static const uint64_t qc[] = {
0, 1, 0x8000u, 0x10000u, 0x18000u,
(uint64_t)0 - 0x10000u, (uint64_t)0 - 1u, (uint64_t)0 - 0x8000u,
0xFFFFFFFFu, 0x100000000u, (uint64_t)0 - 0x100000000u,
((uint64_t)1 << 63) - 1u, (uint64_t)1 << 63
};
const unsigned nc = sizeof qc / sizeof qc[0];
uint64_t x = 0x2545F4914F6CDD1Du;
unsigned doc_wrong = 0, doc_runs = 0;
for (unsigned i = 0; i < nc * nc + 20000u; i++) {
uint64_t a, b;
v4_cell al, ah, bl, bh, w[2];
int lt, gt, eq;
if (i < nc * nc) { a = qc[i / nc]; b = qc[i % nc]; }
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17; a = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; b = x;
if (i & 1u) b = a; /* ties */
if (i & 2u) b ^= (uint64_t)1 << (x >> 58); /* near */
}
qcells(a, &al, &ah);
qcells(b, &bl, &bh);
lt = qs(a) < qs(b); gt = qs(a) > qs(b); eq = a == b;
CHECK(call4(w_dequal, al, ah, bl, bh) && left1(FLAG(eq)), "Q.= [%u]", i);
CHECK(call4(w_dless, al, ah, bl, bh) && left1(FLAG(lt)), "Q.< [%u]", i);
CHECK(call4(w_qgt, al, ah, bl, bh) && left1(FLAG(gt)), "Q.> [%u]", i);
CHECK(call(w_d0equal, 2, al, ah, 0) && left1(FLAG(a == 0)), "Q.0= [%u]", i);
w[0] = lt ? bl : al; w[1] = lt ? bh : ah;
CHECK(call4(w_dmax, al, ah, bl, bh) && leftk(w, 2), "Q.MAX [%u]", i);
w[0] = lt ? al : bl; w[1] = lt ? ah : bh;
CHECK(call4(w_dmin, al, ah, bl, bh) && leftk(w, 2), "Q.MIN [%u]", i);
if ((a >> 63) == (b >> 63)) { /* v3 agrees here */
CHECK(lt == (a < b) && gt == (a > b), "Q compare v3 [%u]", i);
}
/* Q.> as section 5.26 writes it, SWAP D<: count wrong answers. */
doc_runs++;
if (!(call4(w_qgt_doc, al, ah, bl, bh) && left1(FLAG(gt)))) doc_wrong++;
}
printf(" Q.> as written (SWAP D<): %u wrong of %u\n", doc_wrong, doc_runs);
}
/* Q.* against the limb reference (signed, D-8), on edge Q values and
* pseudo-random ones, and against v3's q48_mul for non-negative
* operands, where v3's unsigned product agrees. */
{
static const uint64_t qm[] = {
0, 1, 0x8000u, 0x10000u, 0x18000u, 0x20000u, 0xFFFFu,
(uint64_t)0 - 0x10000u, (uint64_t)0 - 1u, (uint64_t)0 - 0x18000u,
0xFFFFFFFFu, 0x100000000u, 0x7FFFFFFFu, (uint64_t)0 - 0x100000000u,
((uint64_t)1 << 63) - 1u, (uint64_t)1 << 63, (uint64_t)1 << 47
};
const unsigned nm = sizeof qm / sizeof qm[0];
uint64_t x = 0xD1B54A32D192ED03u;
for (unsigned i = 0; i < nm * nm + 20000u; i++) {
uint64_t a, b;
v4_cell al, ah, bl, bh;
v4_ucell rl, rh;
if (i < nm * nm) { a = qm[i / nm]; b = qm[i % nm]; }
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17; a = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; b = x;
if (i & 1u) a = asr64(a, (unsigned)(x >> 58));
if (i & 2u) b = asr64(b, (unsigned)(x >> 52) & 63u);
}
qcells(a, &al, &ah);
qcells(b, &bl, &bh);
qmul_ref((v4_ucell)al, (v4_ucell)ah, (v4_ucell)bl, (v4_ucell)bh, &rl, &rh);
CHECK(call4(w_qstar, al, ah, bl, bh) && left2((v4_cell)rl, (v4_cell)rh),
"Q.* [%u]", i);
if (!(a >> 63) && !(b >> 63)) {
v4_cell vl, vh;
qcells(v3_q48_mul(a, b), &vl, &vh);
CHECK(call4(w_qstar, al, ah, bl, bh) && n.ds.s == vl
&& (V4_CELL_BITS == 64 || n.ds.t == vh), "Q.* v3 [%u]", i);
}
}
}
/* (UQ/) against the Q./ reference, for non-negative a and b where the
* quotient fits (|a| * 2^16 < |b| * 2^(2N-1), i.e. not (b < 2^17 and
* a1 >= b0 << (N-17))). */
{
uint64_t x = 0x94D049BB133111EBu;
unsigned ran = 0;
for (unsigned i = 0; i < 5000; i++) {
v4_ucell al, ah, bl, bh, rl, rh;
int err;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; al = (v4_ucell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; ah = (v4_ucell)x >> 1;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; bl = (v4_ucell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; bh = (v4_ucell)x >> 1;
switch (i & 7u) {
case 1: bh = 0; break;
case 2: bh = 0; bl >>= (unsigned)(x >> 59); break;
case 3: ah >>= (unsigned)(x >> 59) % V4_CELL_BITS; break;
case 4: bh >>= (unsigned)(x >> 58) % V4_CELL_BITS; break;
case 5: ah = 0; bh = 0; break;
case 6: al = bl; ah = bh; break; /* a = b */
default: break;
}
if (bl == 0 && bh == 0) bl = 1;
if (bh == 0 && (bl >> 17) == 0 && ah >= (bl << (V4_CELL_BITS - 17))) continue;
qdiv_ref(al, ah, bl, bh, &rl, &rh, &err);
ran++;
CHECK(call4(w_uqdiv, (v4_cell)al, (v4_cell)ah, (v4_cell)bl, (v4_cell)bh)
&& left2((v4_cell)rl, (v4_cell)rh), "(UQ/) [%u]", i);
}
/* Divisors in [2^(2N-2), 2^(2N-1)), where the remainder's and the
* divisor's top bits can differ. (Q./ passes magnitudes, so b never
* exceeds 2^(2N-1); that one value is covered through Q./.) */
for (unsigned i = 0; i < 4000; i++) {
v4_ucell al, ah, bl, bh, rl, rh;
int err;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; al = (v4_ucell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; ah = (v4_ucell)x >> 1;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; bl = (v4_ucell)x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; bh = (v4_ucell)x;
bh = (bh >> 2) | (V4_MSB >> 1); /* b in [2^(2N-2), 2^(2N-1)) */
if (i & 1u) bh |= (V4_MSB >> 1) - 1u; /* near the top of it */
qdiv_ref(al, ah, bl, bh, &rl, &rh, &err);
ran++;
CHECK(call4(w_uqdiv, (v4_cell)al, (v4_cell)ah, (v4_cell)bl, (v4_cell)bh)
&& left2((v4_cell)rl, (v4_cell)rh), "(UQ/) large b [%u]", i);
}
printf(" (UQ/) cases run: %u\n", ran);
}
/* D2*C against C. */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned j = 0; j < NVEC; j++)
for (unsigned c = 0; c < 2; c++) {
v4_ucell lo = (v4_ucell)vec[i], hi = (v4_ucell)vec[j];
v4_cell w[3];
w[0] = (v4_cell)((lo << 1) | c);
w[1] = (v4_cell)((hi << 1) | (lo >> (V4_CELL_BITS - 1)));
w[2] = (v4_cell)(hi >> (V4_CELL_BITS - 1));
CHECK(call(w_d2starc, 3, vec[i], vec[j], (v4_cell)c) && leftk(w, 3),
"D2*C [%u,%u,%u]", i, j, c);
}
/* Q./ against the reference (D-11), and NODE-ERROR set exactly on
* division by zero; against v3's q48_div for non-negative a < 2^48 and
* b != 0 where v3's quotient is below 2^63. */
{
static const uint64_t qd[] = {
0, 1, 3, 0x8000u, 0x10000u, 0x18000u, 0x30000u, 0xFFFFu,
(uint64_t)0 - 0x10000u, (uint64_t)0 - 1u, (uint64_t)0 - 0x30000u,
0xFFFFFFFFu, 0x100000000u, (uint64_t)0 - 0x100000000u,
((uint64_t)1 << 47) - 1u, (uint64_t)1 << 47,
((uint64_t)1 << 63) - 1u, (uint64_t)1 << 63
};
const unsigned nd = sizeof qd / sizeof qd[0];
uint64_t x = 0xBF58476D1CE4E5B9u;
for (unsigned i = 0; i < nd * nd + 5000u; i++) {
uint64_t a, b;
v4_cell al, ah, bl, bh;
v4_ucell rl, rh;
int err;
if (i < nd * nd) { a = qd[i / nd]; b = qd[i % nd]; }
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17; a = x;
x ^= x << 13; x ^= x >> 7; x ^= x << 17; b = x;
if (i & 1u) a = asr64(a, (unsigned)(x >> 58));
if (i & 2u) b = asr64(b, (unsigned)(x >> 52) & 63u);
if ((i & 12u) == 12u) b = 0;
}
qcells(a, &al, &ah);
qcells(b, &bl, &bh);
qdiv_ref((v4_ucell)al, (v4_ucell)ah, (v4_ucell)bl, (v4_ucell)bh, &rl, &rh, &err);
v4_node_store(&n, NODE_ERROR, 0);
CHECK(call4(w_qslash, al, ah, bl, bh) && left2((v4_cell)rl, (v4_cell)rh),
"Q./ [%u]", i);
CHECK(v4_node_load(&n, NODE_ERROR) == (err ? -1 : 0), "Q./ NODE-ERROR [%u]", i);
if (!(a >> 63) && !(b >> 63) && b != 0 && a < ((uint64_t)1 << 48)) {
uint64_t v3 = (a << 16) / b;
if (!(v3 >> 63)) {
v4_cell vl, vh;
qcells(v3, &vl, &vh);
CHECK(call4(w_qslash, al, ah, bl, bh) && n.ds.s == vl
&& (V4_CELL_BITS == 64 || n.ds.t == vh), "Q./ v3 [%u]", i);
}
}
}
}
/* Q.* and Q./ on edge values built from cells, so that they are the true
* extremes at either cell width: 0, ulp, -ulp, 0.5, 1.0, -1.0, 2.0, Q max,
* Q max - ulp, Q min, Q min + ulp, and the largest Q below 1.0 * 2^(N-16). */
{
v4_ucell ce[13][2];
unsigned k = 0, ne;
ce[k][0] = 0; ce[k][1] = 0; k++;
ce[k][0] = 1; ce[k][1] = 0; k++;
ce[k][0] = MAXU; ce[k][1] = MAXU; k++;
ce[k][0] = 0x8000u; ce[k][1] = 0; k++;
ce[k][0] = 0x10000u; ce[k][1] = 0; k++;
ce[k][0] = (v4_ucell)0 - 0x10000u; ce[k][1] = MAXU; k++;
ce[k][0] = 0x20000u; ce[k][1] = 0; k++;
ce[k][0] = MAXU; ce[k][1] = MAXU >> 1; k++;
ce[k][0] = MAXU - 1u; ce[k][1] = MAXU >> 1; k++;
ce[k][0] = 0; ce[k][1] = V4_MSB; k++;
ce[k][0] = 1; ce[k][1] = V4_MSB; k++;
ce[k][0] = 0; ce[k][1] = 1; k++;
ce[k][0] = MAXU; ce[k][1] = 0; k++;
ne = k;
for (unsigned i = 0; i < ne; i++)
for (unsigned j = 0; j < ne; j++) {
v4_ucell rl, rh;
int err;
qmul_ref(ce[i][0], ce[i][1], ce[j][0], ce[j][1], &rl, &rh);
CHECK(call4(w_qstar, (v4_cell)ce[i][0], (v4_cell)ce[i][1],
(v4_cell)ce[j][0], (v4_cell)ce[j][1])
&& left2((v4_cell)rl, (v4_cell)rh), "Q.* cell edge [%u,%u]", i, j);
qdiv_ref(ce[i][0], ce[i][1], ce[j][0], ce[j][1], &rl, &rh, &err);
v4_node_store(&n, NODE_ERROR, 0);
CHECK(call4(w_qslash, (v4_cell)ce[i][0], (v4_cell)ce[i][1],
(v4_cell)ce[j][0], (v4_cell)ce[j][1])
&& left2((v4_cell)rl, (v4_cell)rh)
&& v4_node_load(&n, NODE_ERROR) == (err ? -1 : 0),
"Q./ cell edge [%u,%u]", i, j);
}
}
/* 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.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.LOG against v3's q48_log_approx for x > 0, bit for bit; x <= 0
* gives 0 and NODE-ERROR (D-12). */
{
static const int64_t ql[] = {
1, 2, 3, 0x7FFF, 0x8000, 0xFFFF, 0x10000, 0x10001, 0x18000, 0x1FFFF, 0x20000,
0x2B7E1, 0x30000, 0xA0000, 0x640000, 0x7FFFFFFF, 0x80000000, 0xFFFFFFFF,
0x100000000, ((int64_t)1 << 47), ((int64_t)1 << 62) + 12345,
(int64_t)(((uint64_t)1 << 63) - 1u), 0, -1, -0x10000
};
uint64_t x = 0xA0761D6478BD642Fu;
for (unsigned i = 0; i < sizeof ql / sizeof ql[0] + 2000u; i++) {
int64_t q;
v4_cell xl, xh, el, eh;
if (i < sizeof ql / sizeof ql[0]) q = ql[i];
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
q = (int64_t)((x >> 1) >> (x & 63u)); /* every magnitude */
if (q == 0) q = 1;
}
qcells((uint64_t)q, &xl, &xh);
if (q <= 0) { el = 0; eh = 0; }
else qcells(v3_q48_log((uint64_t)q), &el, &eh);
v4_node_store(&n, NODE_ERROR, 0);
CHECK(call(w_qlog, 2, xl, xh, 0) && left2(el, eh)
&& v4_node_load(&n, NODE_ERROR) == (q <= 0 ? -1 : 0),
"Q.LOG [%u] x=%lld", i, (long long)q);
}
/* The reduced m takes only 65536 values, 1.0 <= m < 2.0; sweep every
* 17th directly (k = 0), which covers the Newton iteration itself. */
for (int64_t m = 65536; m < 131072; m += 17) {
v4_cell xl, xh, el, eh;
qcells((uint64_t)m, &xl, &xh);
qcells(v3_q48_log((uint64_t)m), &el, &eh);
CHECK(call(w_qlog, 2, xl, xh, 0) && left2(el, eh), "Q.LOG m=%lld", (long long)m);
}
}
/* (Q.REDUCE), Q.SIN and Q.COS against v3's q48_reduce_angle,
* q48_sin_approx and q48_cos_approx, bit for bit, on angles of every size
* and sign. */
{
static const int64_t qa[] = {
0, 1, -1, 0x8000, 0x10000, -0x10000, 102943, 102944, -102943, -102944,
205886, 205887, 205888, -205886, -205887, -205888,
411773, 411774, 411775, -411773, -411774, -411775, 617661, -617661,
0x7FFFFFFF, -0x7FFFFFFF, 0x80000000, 0x100000000, -0x100000000,
((int64_t)1 << 40) + 7, -(((int64_t)1 << 40) + 7),
(int64_t)(((uint64_t)1 << 63) - 1u), (int64_t)((uint64_t)1 << 63)
};
uint64_t x = 0xE7037ED1A0B428DBu;
for (unsigned i = 0; i < sizeof qa / sizeof qa[0] + 3000u; i++) {
int64_t q;
v4_cell xl, xh, el, eh;
if (i < sizeof qa / sizeof qa[0]) q = qa[i];
else {
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
q = (int64_t)asr64(x, (unsigned)(x & 63u)); /* every size, both signs */
if (i & 1u) q %= 1000000; /* and a few turns of the circle */
}
qcells((uint64_t)q, &xl, &xh);
CHECK(call(w_qreduce, 2, xl, xh, 0) && left1((v4_cell)v3_q48_reduce(q)),
"(Q.REDUCE) [%u] x=%lld", i, (long long)q);
qcells(v3_q48_sin((uint64_t)q), &el, &eh);
CHECK(call(w_qsin, 2, xl, xh, 0) && left2(el, eh), "Q.SIN [%u] x=%lld", i, (long long)q);
qcells(v3_q48_cos((uint64_t)q), &el, &eh);
CHECK(call(w_qcos, 2, xl, xh, 0) && left2(el, eh), "Q.COS [%u] x=%lld", i, (long long)q);
}
}
/* Q.1, Q.0 and Q.SCALE: v3's values, and how they behave in use. */
{
v4_cell ol, oh, zl, zh;
uint64_t x = 0x8EBC6AF09C88C6E3u;
qcells(65536u, &ol, &oh); /* v3 Q.1 and Q.SCALE */
qcells(0u, &zl, &zh); /* v3 Q.0 */
CHECK(call(w_q1, 0, 0, 0, 0) && left2(ol, oh), "Q.1 is v3's 65536");
CHECK(call(w_q0, 0, 0, 0, 0) && left2(zl, zh), "Q.0 is v3's 0");
CHECK(call(w_qscale, 0, 0, 0, 0) && left2(ol, oh), "Q.SCALE is v3's 65536");
CHECK(call(w_q1_toint, 0, 0, 0, 0) && left1(1), "Q.1 Q.TO-INT is 1");
CHECK(call(w_one_fromint, 0, 0, 0, 0) && left2(ol, oh), "1 Q.FROM-INT is Q.1");
for (unsigned i = 0; i < 2000; i++) {
uint64_t q;
v4_cell ql, qh;
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
q = asr64(x, (unsigned)(x & 63u));
if (i == 0) q = 0;
if (i == 1) q = ((uint64_t)1 << 63) - 1u;
if (i == 2) q = (uint64_t)1 << 63;
qcells(q, &ql, &qh);
CHECK(call(w_q1_times, 2, ql, qh, 0) && left2(ql, qh), "q Q.1 Q.* is q [%u]", i);
CHECK(call(w_q0_plus, 2, ql, qh, 0) && left2(ql, qh), "q Q.0 Q.+ is q [%u]", i);
}
}
/* LSHIFT and RSHIFT as written, against C, for every count 0 .. N (a
* count of N must give 0). */
for (unsigned i = 0; i < NVEC; i++)
for (unsigned c = 0; c <= V4_CELL_BITS; c++) {
v4_ucell u = (v4_ucell)vec[i];
v4_ucell l = c < V4_CELL_BITS ? u << c : 0u;
v4_ucell r = c < V4_CELL_BITS ? u >> c : 0u;
CHECK(call(w_lshift, 2, vec[i], (v4_cell)c, 0) && left1((v4_cell)l), "LSHIFT [%u,%u]", i, c);
CHECK(call(w_rshift, 2, vec[i], (v4_cell)c, 0) && left1((v4_cell)r), "RSHIFT [%u,%u]", i, c);
}
/* C@ and C!: every byte of two adjacent words, several
* values, the other bytes and the neighbouring words left alone. */
{
static const v4_ucell cv[] = { 0, 1, 0x41, 0x7F, 0x80, 0xFF, 0x1FF, 0xABCD };
const v4_ucell pat = (v4_ucell)0xA5C3F00Fu | ((MAXU >> 16 >> 16) << 16 << 16);
for (unsigned b = 0; b < 8; b++)
for (unsigned k = 0; k < sizeof cv / sizeof cv[0]; k++) {
v4_cell w = BYTES + 1 + (v4_cell)(b / 4u);
unsigned sh = 8u * (b % 4u);
v4_ucell want = (pat & ~((v4_ucell)0xFFu << sh)) | ((cv[k] & 0xFFu) << sh);
v4_cell baddr = (BYTES + 1) * 4 + (v4_cell)b;
for (unsigned j = 0; j < 4; j++) v4_node_store(&n, BYTES + (v4_cell)j, (v4_cell)pat);
CHECK(call(w_cstore, 2, (v4_cell)cv[k], baddr, 0) && n.ds.t == CANARY,
"C! returns [%u,%u]", b, k);
CHECK((v4_ucell)v4_node_load(&n, w) == want, "C! stores [%u,%u]", b, k);
for (unsigned j = 0; j < 4; j++)
if (BYTES + (v4_cell)j != w)
CHECK((v4_ucell)v4_node_load(&n, BYTES + (v4_cell)j) == pat, "C! neighbours [%u,%u,%u]", b, k, j);
CHECK(call(w_cfetch, 1, baddr, 0, 0) && left1((v4_cell)(cv[k] & 0xFFu)), "C@ [%u,%u]", b, k);
}
}
/* Q./ at the overflow boundary: |a| * 2^16 against |b| * 2^(2N-1), for
* divisors around 2^16 and 2^17, every sign. */
{
static const v4_ucell bv[] = { 0x10000u, 0x10001u, 0x18000u, 0x1FFFFu, 0x20000u, 0x20001u, 0xFFFFu };
for (unsigned i = 0; i < sizeof bv / sizeof bv[0]; i++)
for (unsigned j = 0; j < 4; j++)
for (unsigned sg = 0; sg < 4; sg++) {
v4_ucell thr = bv[i] << (V4_CELL_BITS - 17), al, ah, bl = bv[i], bh = 0, rl, rh;
int err;
switch (j) {
case 0: al = 0; ah = thr; break; /* at the threshold */
case 1: al = MAXU; ah = thr - 1u; break; /* just below */
case 2: al = 1; ah = thr; break; /* just above */
default: al = MAXU; ah = thr; break;
}
if (ah & V4_MSB) continue; /* |a| must be a Q value */
if (sg & 1u) dneg(al, ah, &al, &ah);
if (sg & 2u) dneg(bl, bh, &bl, &bh);
qdiv_ref(al, ah, bl, bh, &rl, &rh, &err);
CHECK(call4(w_qslash, (v4_cell)al, (v4_cell)ah, (v4_cell)bl, (v4_cell)bh)
&& left2((v4_cell)rl, (v4_cell)rh), "Q./ boundary [%u,%u,%u]", i, j, sg);
}
}
/* Stack headroom of the two longest definitions. */
{
static const v4_cell umstar_args[][4] = {
{ 0, 0, 0 }, { -1, -1, 0 }, { 3, -1, 0 }, { -1, 3, 0 }, { 12345, -12345, 0 }
};
static const v4_cell ummod_args[][4] = {
{ 0, 0, 1 }, { -1, -2, -1 }, { 12345, 0, 7 }, { -1, 0x7FFF, 0x8000 },
{ 0, (v4_cell)(V4_MSB - 1u), (v4_cell)V4_MSB }
};
int dh, rh;
headroom("UM*", w_umstar, 2, umstar_args, 5, 2, &dh, &rh);
CHECK(dh >= 0 && rh >= 0, "UM* runs at all");
headroom("UM/MOD", w_ummod, 3, ummod_args, 5, 2, &dh, &rh);
CHECK(dh >= 0 && rh >= 0, "UM/MOD runs at all");
{
static const v4_cell one_args[][4] = { { 5, 0, 0 }, { -5, 0, 0 }, { 0, 0, 0 } };
static const v4_cell two_args[][4] = {
{ 5, 0, 0 }, { 0, -1, 0 }, { -5, -1, 0 }, { -1, 5, 0 }, { 0, 0, 0 }
};
static const v4_cell smrem_args[][4] = {
{ 7, 0, 2 }, { -7, -1, 2 }, { 7, 0, -2 }, { -7, -1, -2 }, { 0, 0, -3 }
};
static const v4_cell slashmod_args[][4] = {
{ 7, 2, 0 }, { -7, 2, 0 }, { 7, -2, 0 }, { -7, -2, 0 }
};
headroom("ABS", w_abs, 1, one_args, 3, 1, &dh, &rh);
headroom("U>", w_ugreater, 2, two_args, 5, 1, &dh, &rh);
static const v4_cell dplus_args[][4] = {
{ -1, 0, 1, 0 }, { 5, 7, 9, 11 }, { -1, -1, -1, -1 },
{ (v4_cell)V4_MSB, 0, 1, 0 }, { (v4_cell)V4_MSB, 0, (v4_cell)V4_MSB, 0 },
{ 1, 0, 2, 0 }
};
headroom("D+", w_dplus, 4, dplus_args, 6, 2, &dh, &rh);
CHECK(rh >= 5, "D+ leaves 5 return entries");
headroom("DNEGATE", w_dnegate, 2, two_args, 5, 2, &dh, &rh);
headroom("DABS", w_dabs, 2, two_args, 5, 2, &dh, &rh);
static const v4_cell mplus_args[][4] = {
{ -1, 0, 1, 0 }, { 5, 7, -9, 0 }, { 0, 0, -1, 0 }, { (v4_cell)V4_MSB, 0, (v4_cell)V4_MSB, 0 }
};
static const v4_cell dsub_args[][4] = {
{ 0, 0, 1, 0 }, { 5, 7, 9, 11 }, { -1, -1, -1, -1 }, { 0, 1, 0, 1 },
{ (v4_cell)V4_MSB, 0, (v4_cell)V4_MSB, 0 }, { 1, 0, 2, 0 }
};
headroom("M+", w_mplus, 3, mplus_args, 4, 2, &dh, &rh);
headroom("M-", w_mminus, 3, mplus_args, 4, 2, &dh, &rh);
CHECK(dh >= 4 && rh >= 4, "M- leaves room");
{
static const v4_cell six[6] = { 1, 2, 3, 4, 5, 6 };
static const v4_cell d1_args[][4] = { { 5, -1, 0, 0 }, { -1, 0, 0, 0 }, { 1, 1, 0, 0 } };
static const v4_cell ms_args[][4] = { { 7, 3, 0, 0 }, { -7, 3, 0, 0 }, { -7, -3, 0, 0 } };
static const v4_cell ss_args[][4] = { { 7, 3, 2, 0 }, { -7, 3, 2, 0 }, { 7, -3, -2, 0 } };
int d, r;
headroom("M*", w_mstar, 2, ms_args, 3, 2, &dh, &rh);
CHECK(dh >= 4 && rh >= 4, "M* leaves room");
headroom("M/MOD", w_mslashmod, 3, smrem_args, 5, 2, &dh, &rh);
headroom("MOD", w_mod, 2, slashmod_args, 4, 1, &dh, &rh);
CHECK(dh >= 3 && rh >= 2, "MOD leaves room");
headroom("*/MOD", w_starslashmod, 3, ss_args, 3, 2, &dh, &rh);
CHECK(dh >= 3 && rh >= 2, "*/MOD leaves room");
headroom("*/", w_starslash, 3, ss_args, 3, 1, &dh, &rh);
CHECK(dh >= 3 && rh >= 2, "*/ leaves room");
headroom("D0<", w_d0less, 2, d1_args, 3, 1, &dh, &rh);
headroom("D2*", w_d2star, 2, d1_args, 3, 2, &dh, &rh);
headroom("D2/", w_d2slash, 2, d1_args, 3, 2, &dh, &rh);
CHECK(dh >= 5 && rh >= 6, "D2/ leaves room");
for (d = 0; d < V4_DATA_DEPTH; d++) if (!rot6_ok(six, (unsigned)d + 1u, 0)) break;
for (r = 0; r < V4_RET_DEPTH; r++) if (!rot6_ok(six, 0, (unsigned)r + 1u)) break;
printf(" 2ROT headroom: data %d below canary, return %d below its return address\n", d, r);
CHECK(d >= 2 && r >= 1, "2ROT leaves room");
}
headroom("D-", w_dminus, 4, dsub_args, 6, 2, &dh, &rh);
headroom("D0=", w_d0equal, 2, two_args, 5, 1, &dh, &rh);
headroom("D=", w_dequal, 4, dsub_args, 6, 1, &dh, &rh);
CHECK(dh >= 0 && rh >= 0, "D= runs at all");
static const v4_cell cmp_args[][4] = {
{ 0, 0, 0, 0 }, { 1, 0, 2, 0 }, { 2, 0, 1, 0 }, { 0, -1, 0, 1 }, { 0, 1, 0, -1 },
{ 0, 5, 0, 7 }, { 0, 7, 0, 5 }, { -1, 0, 1, 0 }, { 1, 0, -1, 0 }
};
headroom("<", w_less, 2, two_args, 5, 1, &dh, &rh);
headroom("=", w_equal, 2, two_args, 5, 1, &dh, &rh);
headroom("D<", w_dless, 4, cmp_args, 9, 1, &dh, &rh);
headroom("(D<)", w_dltkeep, 4, cmp_args, 9, 5, &dh, &rh);
headroom("2SWAP", w_2swap, 4, cmp_args, 9, 4, &dh, &rh);
headroom("2OVER", w_2over, 4, cmp_args, 9, 6, &dh, &rh);
headroom("DMAX", w_dmax, 4, cmp_args, 9, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "DMAX leaves room");
headroom("DMIN", w_dmin, 4, cmp_args, 9, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "DMIN leaves room");
headroom("Q.> (2SWAP D<)", w_qgt, 4, cmp_args, 9, 1, &dh, &rh);
static const v4_cell qstar_args[][4] = {
{ 0, 1, 0, 1 }, { -1, -1, -1, -1 }, { 0, -1, 0, 1 }, { 5, 7, 9, 11 },
{ (v4_cell)V4_MSB, -1, 3, (v4_cell)V4_MSB }
};
headroom("Q.*", w_qstar, 4, qstar_args, 5, 2, &dh, &rh);
CHECK(dh >= 0 && rh >= 0, "Q.* runs at all");
static const v4_cell qdiv_args[][4] = {
{ 0, 1, 0, 1 }, { 0, 3, 0, 1 }, { 0, -1, 0, 1 }, { 5, 7, 9, 11 },
{ 0, 1, 0, 0 }, { 0, 0x7FFF, 1, 0 }, { -1, 0x7FFF, 0, -3 }
};
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");
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");
static const v4_cell qlog_args[][4] = {
{ 0x10000, 0, 0, 0 }, { 0x2B7E1, 0, 0, 0 }, { 3, 0, 0, 0 }, { -1, 0x7FFF, 0, 0 },
{ 0, 0, 0, 0 }, { 5, -1, 0, 0 }
};
headroom("Q.LOG", w_qlog, 2, qlog_args, 6, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "Q.LOG leaves room");
static const v4_cell qtrig_args[][4] = {
{ 0, 0, 0, 0 }, { 0x10000, 0, 0, 0 }, { -0x10000, -1, 0, 0 }, { 205887, 0, 0, 0 },
{ -1, 0x7FFF, 0, 0 }, { 5, -3, 0, 0 }
};
headroom("(Q.REDUCE)", w_qreduce, 2, qtrig_args, 6, 1, &dh, &rh);
headroom("Q.SIN", w_qsin, 2, qtrig_args, 6, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "Q.SIN leaves room");
headroom("Q.COS", w_qcos, 2, qtrig_args, 6, 2, &dh, &rh);
CHECK(dh >= 2 && rh >= 2, "Q.COS leaves room");
static const v4_cell shift_args[][4] = { { -1, 0, 0, 0 }, { -1, 5, 0, 0 }, { 12345, 31, 0, 0 } };
headroom("LSHIFT", w_lshift, 2, shift_args, 3, 1, &dh, &rh);
headroom("RSHIFT", w_rshift, 2, shift_args, 3, 1, &dh, &rh);
static const v4_cell cf_args[][4] = { { (V4_NODE_WORDS - 47) * 4 + 3, 0, 0, 0 }, { (V4_NODE_WORDS - 47) * 4, 0, 0, 0 } };
static const v4_cell cs_args[][4] = { { 0x41, (V4_NODE_WORDS - 47) * 4 + 3, 0, 0 }, { 0xFF, (V4_NODE_WORDS - 47) * 4, 0, 0 } };
headroom("C@", w_cfetch, 1, cf_args, 2, 1, &dh, &rh);
CHECK(dh >= 7 && rh >= 7, "C@ is call-free");
headroom("C!", w_cstore, 2, cs_args, 2, 0, &dh, &rh);
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);
headroom("/", w_slash, 2, slashmod_args, 4, 1, &dh, &rh);
CHECK(dh >= 3 && rh >= 2, "/ runs no deeper than /MOD");
headroom("*", w_star, 2, slashmod_args, 4, 1, &dh, &rh);
CHECK(dh >= 4 && rh >= 5, "* leaves room");
CHECK(dh >= 0 && rh >= 0, "/MOD runs at all");
}
}
CHECK(v4_node_guards_intact(&n), "guards intact");
printf(" %d checks, %d failures\n", checks, failures);
return failures ? 1 : 0;
}
Reference in New Issue View Git Blame Copy Permalink