diff --git a/docs/v4.0.0/DECOMPOSITION.md b/docs/v4.0.0/DECOMPOSITION.md index 48ccf306..895cbd90 100644 --- a/docs/v4.0.0/DECOMPOSITION.md +++ b/docs/v4.0.0/DECOMPOSITION.md @@ -147,8 +147,8 @@ IF body1 ELSE body2 THEN → if L1 drop body1 jump L2 times, as on the F18. `FOR ... UNEXT` is the same but the body must fit in one instruction word. **Register conventions.** `A` and `B` are caller-saved. A word that uses them says so. Words in this -document that clobber `A`: `@ ! +! -! 2@ 2! C@ C! UM* * UM/MOD Q.FROM-INT Q.TO-INT Q.* Q./ Q.EXP Q.SQRT Q.LOG SEND RECV`. Words that clobber `B`: -`Q./ Q.EXP Q.SQRT Q.LOG SEND RECV`. +document that clobber `A`: `@ ! +! -! 2@ 2! C@ C! UM* * UM/MOD Q.FROM-INT Q.TO-INT Q.* Q./ Q.EXP Q.SQRT Q.LOG Q.SIN SEND RECV`. Words that clobber `B`: +`Q./ Q.EXP Q.SQRT Q.LOG Q.SIN SEND RECV`. **Return-stack words** (`>R R> R@ 2>R 2R> 2R@ I J UNLOOP` and the loop runtimes) are always IN. @@ -664,7 +664,10 @@ word. Per D-10 it occupies two cells on a 64-bit node too. Per D-8, v4 Q values | `D2*C` | CAP | `( lo hi cin -- lo' hi' cout )`: the double shifted left one bit, `cin` entering at the bottom and `cout` the bit leaving the top (both 0 or 1). `a! dup -if L1 drop 1 jump L2 L1: drop 0 L2: push over -if L3 drop 1 jump L4 L3: drop 0 L4: over + + push 2* a + pop pop` — `cin` waits in `A`, and `2hi + m` is `over + +`, so it needs one return entry. Executed on the golden model (2026-10-02). Clobbers `A`. | | `(Q/)` | CAP | Variable, 5 cells: `Q./`'s quotient register, divisor and result sign. Like `BASE` and the hold buffer, it lives in node memory. | | `Q.ABS` `Q.NEG` | CAP | `DABS`, `DNEGATE` — executed on the golden model (2026-10-02) against v3's `q48_abs` and `0 - q`: bit-for-bit at 32-bit cells, including v3's wrap of Q min to itself; at 64-bit cells the low cell is v3's and the high cell is the true sign, per D-10. | -| `Q.SIN` `Q.COS` | CAP | Algorithms ported from `q48_words.c`; the hosted C versions are the golden model. | +| `Q.COS` | CAP | Algorithm ported from `q48_words.c`; the hosted C version is the golden model. | +| `Q.SIN` | CAP | 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 the variable `(QT)` holds single cells and the only calls are `(Q.REDUCE)`, `Q.*`, `Q./` and `UM*`. The Taylor loop is written so that `Q.COS` can jump into it. Executed on the golden model (2026-10-03): bit for bit v3's result on angles of every size and sign. Leaves its caller 3 data cells under its argument and 2 return entries. Clobbers `A` and `B`. | +| `(Q.REDUCE)` | CAP | `( x -- r )`, v3's `q48_reduce_angle`: the angle reduced to one cell, `-pi <= r <= pi` with `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 and its remainder then leading the low cell. `dup (QT) b! !b -if L1 DNEGATE L1: 0 411774 UM/MOD drop 411774 UM/MOD drop (QT) b! @b -if P1 drop inv 1 + jump J1 P1: drop J1: dup -205888 + -if BIG drop dup 205887 + -if OK drop 411774 + ; OK: drop ; BIG: drop -411774 + ;` Executed on the golden model (2026-10-03), bit for bit v3's. Clobbers `A` and `B`. | +| `(QT)` | CAP | Variable, 6 cells: `Q.SIN` / `Q.COS`'s sign, `x^2`, term, sum, `n` and subtract flag. | | `Q.LOG` | CAP | 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`, stopping once `|m - e^y| < 100` ulp; result `y + k * ln 2` with `ln 2 = 45426`. `e^y` is `Q.EXP`'s Taylor loop **written in line** (flattened, ruled 2026-10-02), so `Q.LOG` calls only `Q.*`, `Q./` and `UM*` and the call chain is one level shorter than through `Q.EXP`. After the reduction `m`, `y`, the Taylor term and its sum each fit one cell at either width, so the variable `(QL)` holds single cells. `x <= 0` returns 0 and sets `NODE-ERROR` (D-12). Executed on the golden model (2026-10-03): bit for bit v3's result for every `x > 0` tested, including a sweep of the reduced `m`. v3's clamp `y = max(y - corr, 0)` is kept for fidelity but cannot fire: over all 65536 values of `m` it is never reached, and the iteration never needs more than 4 rounds. Leaves its caller 3 data cells under its argument and 2 return entries. Clobbers `A` and `B`. | | `(QL)` | CAP | Variable, 7 cells: `Q.LOG`'s `m`, `y`, `k`, rounds left, Taylor term (then `delta`), sum, and `n * 1.0`. | | `Q.SQRT` | CAP | v3's `q48_sqrt_approx`: Newton from `x0 = q/2 + 0.25`, up to 8 rounds of `x' = (x + q/x) / 2`, returning `x` as soon as `|x' - x| < 10` ulp; `sqrt(0) = 0`, `sqrt(1.0) = 1.0`. `q < 0` returns 0 and sets `NODE-ERROR` (D-12). `q`, `x` and the rounds left live in the variable `(QR)`. The halving is a logical double shift, `push a! 0 pop +* MAXHI and push drop a pop`. Executed on the golden model (2026-10-02): bit for bit v3's result for `0 <= q < 2^48`; above that v3's `q48_div` saturates and v4 divides correctly, so they part. Leaves its caller 3 data cells under its argument and 2 return entries. Clobbers `A` and `B`. | diff --git a/v4/tests/test_foundation.c b/v4/tests/test_foundation.c index b0b0eaeb..197ba5d1 100644 --- a/v4/tests/test_foundation.c +++ b/v4/tests/test_foundation.c @@ -31,7 +31,9 @@ * 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. Every word is + * for 0 with NODE-ERROR at or below zero. + * (Q.REDUCE) and Q.SIN are checked bit for bit against v3's + * q48_reduce_angle and q48_sin_approx on angles of every size and sign. Every word is * also probed for how much of the 10- and 9-deep circular stacks (D-2) it * leaves to its caller. * @@ -60,6 +62,8 @@ static int failures = 0, checks = 0; #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)) #define MAXU ((v4_ucell)~(v4_ucell)0) #define FLAG(c) ((c) ? V4_ALL_ONES : (v4_cell)0) @@ -73,7 +77,8 @@ static v4_cell w_nip, w_swap, w_or, w_negate, w_rot, w_zless, w_zequal, w_ugreater, w_abs, w_s2d, w_dplus, w_dnegate, w_dabs, w_smrem, w_slashmod, w_mplus, w_dminus, w_d0equal, w_dequal, w_qfromint, w_qtoint, w_less, w_equal, w_dless, w_2swap, w_2over, - w_dmax, w_dmin, w_qgt_doc, w_qgt, w_dltkeep, w_qstar, w_d2starc, w_uqdiv, w_qslash, w_qexp, w_qsqrt, w_qlog; + 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; #define O(name) v4_asm_op(&as, V4_OP_##name) #define LIT(v) v4_asm_lit(&as, (v4_cell)(v)) @@ -1172,6 +1177,106 @@ static void build(void) #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); + } +#undef VSET +#undef VGET +#undef JUMP_TO +#undef HERE_ #undef FWD CHECK(v4_asm_ok(&as), "foundation words assemble"); @@ -1612,6 +1717,37 @@ static uint64_t v3_q48_log(uint64_t x) 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; +} + /* 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 @@ -2332,6 +2468,35 @@ int main(void) } } + /* (Q.REDUCE) and Q.SIN against v3's q48_reduce_angle and q48_sin_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); + } + } + /* Q./ at the overflow boundary: |a| * 2^16 against |b| * 2^(2N-1), for * divisors around 2^16 and 2^17, every sign. */ { @@ -2451,6 +2616,13 @@ int main(void) }; 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("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);