starkernel: item 4.3.3a -- Q48.16 trigonometry (Q.SIN/Q.COS)
Adds q48_reduce_angle() (range-reduce a signed Q48.16 angle into [-PI_Q48, PI_Q48] via one integer division plus a bounded fix-up loop) and q48_sin_approx/q48_cos_approx (Taylor series, terms n=3,5,7,9,11 for sin and n=2,4,6,8,10 for cos, early exit below 10). Q.SIN/Q.COS registered as FORTH words in q48_words.c, same pattern as Q.LOG/Q.EXP/Q.SQRT. Found mid-implementation: this codebase has two independent Q48.16 implementations -- src/word_source/q48_16_words.c (hosted/vendored) and src/starkernel/math/q48_16.c (kernel-only; the kernel build does not compile the former at all). The hosted build linked fine after the first pass; the kernel build failed with undefined references until the same two functions were added to both .c files and both q48_16.h headers (include/q48_16.h and include/starkernel/q48_16.h). Not fixed at the root -- Q.LOG/Q.EXP/Q.SQRT already had this same four-file duplication, unremarked until now -- just navigated correctly for this item. Verified live on amd64 via serial injection: sin/cos at 0, +-pi/2, pi, and 3pi (range-reduction across multiple turns) all match expected values within Taylor-series truncation error (<0.2%). All three architectures boot clean to ok> with the DoE completing; dict_hash identical across all three (0x291a660b05fa7b52). FABRIC.md item 4.3.3a marked done with full acceptance evidence.
This commit is contained in:
@@ -258,3 +258,85 @@ q48_16_t q48_sqrt_approx(q48_16_t q)
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
/* ============================================================================
|
||||
* Approximation: Sine / Cosine (Taylor Series, Integer-Only)
|
||||
* ============================================================================
|
||||
*
|
||||
* FABRIC.md item 4.3.3a -- needed by the Console drawing fabric's
|
||||
* CIRCLE/ARC/ELLIPSE (item 4.3.3b). Radian input.
|
||||
*
|
||||
* PI_Q48 = 205887 (pi * 65536, rounded). TWO_PI_Q48 is derived as
|
||||
* 2 * PI_Q48 rather than independently rounded, so the range-reduction
|
||||
* boundary at +-PI_Q48 is self-consistent (no seam).
|
||||
*/
|
||||
|
||||
/* Reduce a signed Q48.16 angle into [-PI_Q48, PI_Q48]. A single integer
|
||||
* division on the raw Q48.16 representations gives the correct (scale-
|
||||
* independent) quotient of how many full 2*pi turns to remove, then a
|
||||
* bounded fix-up loop (at most one or two iterations) handles the
|
||||
* remainder landing just outside the target interval. */
|
||||
static int64_t q48_reduce_angle(int64_t x)
|
||||
{
|
||||
const int64_t PI_Q48 = 205887; /* pi * 65536, rounded */
|
||||
const int64_t TWO_PI_Q48 = 2 * PI_Q48; /* derived, not independently rounded */
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
q48_16_t q48_sin_approx(q48_16_t q)
|
||||
{
|
||||
int64_t x = q48_reduce_angle((int64_t)q);
|
||||
|
||||
int negative = (x < 0);
|
||||
if (negative) x = -x;
|
||||
|
||||
q48_16_t xu = (q48_16_t)x;
|
||||
q48_16_t x2 = q48_mul(xu, xu);
|
||||
q48_16_t term = xu;
|
||||
q48_16_t result = xu;
|
||||
int subtract = 1;
|
||||
|
||||
/* sin(x) = x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - x^11/11! ... */
|
||||
for (int n = 3; n <= 11; n += 2) {
|
||||
term = q48_mul(term, x2);
|
||||
term = q48_div(term, q48_from_u64((uint64_t)(n * (n - 1))));
|
||||
|
||||
result = subtract ? q48_sub(result, term) : q48_add(result, term);
|
||||
subtract = !subtract;
|
||||
|
||||
if (term < 10) break;
|
||||
}
|
||||
|
||||
return negative ? (q48_16_t)(0ULL - result) : result;
|
||||
}
|
||||
|
||||
q48_16_t q48_cos_approx(q48_16_t q)
|
||||
{
|
||||
int64_t x = q48_reduce_angle((int64_t)q);
|
||||
|
||||
q48_16_t xu = (q48_16_t)(x < 0 ? -x : x);
|
||||
q48_16_t x2 = q48_mul(xu, xu);
|
||||
q48_16_t term = Q48_ONE;
|
||||
q48_16_t result = term;
|
||||
int subtract = 1;
|
||||
|
||||
/* cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + x^8/8! - x^10/10! ... */
|
||||
for (int n = 2; n <= 10; n += 2) {
|
||||
term = q48_mul(term, x2);
|
||||
term = q48_div(term, q48_from_u64((uint64_t)(n * (n - 1))));
|
||||
|
||||
result = subtract ? q48_sub(result, term) : q48_add(result, term);
|
||||
subtract = !subtract;
|
||||
|
||||
if (term < 10) break;
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user