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

feat(v4.0.0): Q.SIN and (Q.REDUCE), bit for bit with v3

Browse Source 
(Q.REDUCE) is v3's q48_reduce_angle: the angle as one cell in [-pi, pi].
x / 2pi does not fit a cell, but only the remainder is needed, and
|x| mod 2pi is two UM/MOD steps (high cell first, its remainder leading
the low cell); then the sign of x and one step of 2pi back into range.

Q.SIN is v3's q48_sin_approx on the reduced angle: the odd Taylor series
to n = 11, each term the last times x^2 over n(n-1), stopping below 10
ulp. After the reduction every value fits one cell, so the 6-cell
variable (QT) holds single cells. The loop is laid out for Q.COS to jump
into, so the pair share it without an extra call level.

Checked bit for bit against v3 (both functions ported into the test) on
33 edge angles -- around +-pi/2, +-pi, +-2pi, the 32-bit seam and both
64-bit extremes -- and 3000 pseudo-random angles of every size and sign,
at 32- and 64-bit cells, optimised and ASan+UBSan (`make test`,
`make sanitize`). Mutating either range boundary, the series length or
the sign alternation fails. Moving the 10-ulp stop to 11 fails nothing,
and a scan of all 205888 reduced angles in C shows why: no input's
result depends on it, for sine or cosine.

Headroom (data under arg / return): (Q.REDUCE) 5/5, Q.SIN 3/2.

Test results on the amd64 host only. This is a development check, not
acceptance (JUSTIFICATION.md section 16).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
rajames
2026-10-03 10:44:42 -04:00
co-authored by Claude Opus 5.5
parent e117333a87
commit 1ce705783c
2 changed files with 180 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
v4/tests/test_foundation.c
+174 -2
View File
@@ -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);