Reorganize source tree: kernel/, v3/, v4/ split and board infrastructure
Source tree reorganization: - Move StarForth v3 engine to v3/ (src/, include/, Makefile) - Move kernel to kernel/ (src/, include/, linker/, Makefile) - Create v4/ skeleton for F18-ISA golden model (DECOMPOSITION.md, JUSTIFICATION.md) - Move FABRIC-0..4.md to docs/fabric/ - Move ONTOLOGY.md and ROADMAP.md to docs/ Board infrastructure: - Add boards/ser5/, boards/raspi/, boards/milkv/, boards/zynq7020/ - Each board has board.mk (ISA, CPU flags, boot recipe) and README.md - Root Makefile becomes thin dispatcher: boot_image, all, clean, docs take TARGET - make boot_image TARGET=SER5|RASPI|MILKV builds one GPT/MBR image per board - ZYNQ7020 target exists but stops with clear error (ARMv7 port not built yet) - scripts/mkdiskimage.sh builds disk images for all boards Docs pipeline: - docs/book/ with LaTeX master (main.tex) and Makefile - pandoc converts Markdown to LaTeX at build time - Two Lua filters: table-widths.lua (wide tables wrap), code-breaks.lua (inline code breaks) - make docs builds single PDF (754 pages, 0 missing characters) - make docs TARGET=<board> adds board appendix - build/docs/<book|board>/meta.tex stamps git commit into PDF Bug fixes: - 42 include paths that only worked by accident now use correct relative paths - clang-18 hardcode replaced with configurable CC variable (fixed aarch64 build) - Pi 5: kernel_2712.img linked at 0x80000, .bss zeroed, memory reserved - Doxyfile, .clang-tidy, README.md, Kconfig paths updated Verified: - Hosted v3 build passes 1012 tests, 0 failures - SER5 image boots in QEMU (OVMF), POST passes, K exact (65536 = Q48_ONE) - Milk-V image boots in QEMU (OpenSBI + U-Boot + bootefi), POST passes - make clean TARGET=<board> removes only that board and its ISA objects - make all builds all boards, hosted v3, and docs in one run Co-authored-by: Junie <junie@jetbrains.com>
This commit is contained in:
rajames
@@ -0,0 +1,263 @@
|
||||
/*
|
||||
StarForth — Steady-State Virtual Machine Runtime
|
||||
|
||||
Copyright (c) 2023–2025 Robert A. James
|
||||
All rights reserved.
|
||||
|
||||
This file is part of the StarForth project.
|
||||
|
||||
Licensed under the StarForth License, Version 1.0 (the "License");
|
||||
you may not use this file except in compliance with the License.
|
||||
|
||||
You may obtain a copy of the License at:
|
||||
https://github.com/star.4th@proton.me/StarForth/LICENSE.txt
|
||||
|
||||
This software is provided "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
express or implied, including but not limited to the warranties of
|
||||
merchantability, fitness for a particular purpose, and noninfringement.
|
||||
|
||||
See the License for the specific language governing permissions and
|
||||
limitations under the License.
|
||||
|
||||
StarForth — Steady-State Virtual Machine Runtime
|
||||
Copyright (c) 2023–2025 Robert A. James
|
||||
All rights reserved.
|
||||
|
||||
This file is part of the StarForth project.
|
||||
|
||||
Licensed under the StarForth License, Version 1.0 (the "License");
|
||||
you may not use this file except in compliance with the License.
|
||||
|
||||
You may obtain a copy of the License at:
|
||||
https://github.com/star.4th@proton.me/StarForth/LICENSE.txt
|
||||
|
||||
This software is provided "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
express or implied, including but not limited to the warranties of
|
||||
merchantability, fitness for a particular purpose, and noninfringement.
|
||||
|
||||
See the License for the specific language governing permissions and
|
||||
limitations under the License.
|
||||
|
||||
*/
|
||||
|
||||
/*
|
||||
* math_portable.c
|
||||
*
|
||||
* Pure 64-bit fixed-point Q48.16 math functions.
|
||||
* Suitable for embedded environments with no libm or floating-point
|
||||
* (HISTORICAL: also L4Re, through mid-2026; removed as an active target).
|
||||
* All arithmetic is integer-only for maximum production performance.
|
||||
*
|
||||
* Q48.16 Format:
|
||||
* - 48-bit signed integer + 16-bit fractional part
|
||||
* - 1.0 = 65536 (2^16)
|
||||
* - Precision: 2^-16 ≈ 0.0000153
|
||||
* - Range: ±140 trillion
|
||||
*/
|
||||
|
||||
#include <stdint.h>
|
||||
#include "math_portable.h"
|
||||
|
||||
/**
|
||||
* @brief Compute exp(@p x) via Taylor series in Q48.16 fixed-point.
|
||||
*
|
||||
* Evaluates the classic Taylor expansion:
|
||||
* @code
|
||||
* exp(x) = 1 + x + x²/2! + x³/3! + … (up to 40 terms)
|
||||
* @endcode
|
||||
* entirely in 64-bit integer arithmetic using the Q48.16 representation
|
||||
* (1.0 = 65536 = 2^16).
|
||||
*
|
||||
* The function is intentionally scoped to the range the physics engine
|
||||
* actually uses — the argument to @c erf_q48() is @c -x² times a ratio
|
||||
* near 1, yielding values in roughly the @c [-20, 0] interval. Outside
|
||||
* this range the implementation provides coarse-grained clamping:
|
||||
* - @p x == 0 → returns @c Q48 (1.0 exactly)
|
||||
* - @p x > 20.0 → saturates at @c Q48 << 8 (approx e^20, avoids overflow)
|
||||
* - @p x < -100.0 → underflows to 0 (result is negligible)
|
||||
*
|
||||
* The term loop exits early on the first term smaller than 1 ULP (the
|
||||
* series has converged) and guards against intermediate overflow by
|
||||
* clamping @c result at @c INT64_MAX-1.
|
||||
*
|
||||
* @note This is a file-local (@c static) helper; callers are
|
||||
* @c erf_q48() only.
|
||||
*
|
||||
* @param x Input value in Q48.16 format (1.0 = 65536). Intended range:
|
||||
* approximately [-100×65536, +20×65536].
|
||||
* @return exp(@p x) in Q48.16 format; 0 on underflow; saturated at
|
||||
* @c INT64_MAX-1 on overflow.
|
||||
*/
|
||||
static int64_t exp_q48(int64_t x) {
|
||||
const int64_t Q48 = 65536LL; /* 2^16 */
|
||||
|
||||
/* Handle extremes */
|
||||
if (x == 0) return Q48;
|
||||
if (x > (20LL << 16)) return Q48 << 8; /* Approximate e^20 */
|
||||
if (x < -(100LL << 16)) return 0; /* Underflow to zero */
|
||||
|
||||
/* Taylor series: exp(x) = Σ x^n / n! (in fixed-point) */
|
||||
int64_t result = Q48; /* Start with 1 */
|
||||
int64_t term = Q48;
|
||||
int64_t x_power = x; /* x^1 */
|
||||
|
||||
for (int i = 1; i < 40; i++) {
|
||||
/* term = x^i / i! */
|
||||
term = (x_power / (int64_t)i) / Q48; /* Scale down to avoid overflow */
|
||||
if (term > -1 && term < 1) break; /* Negligible */
|
||||
|
||||
result += term;
|
||||
x_power = (x_power * x) / Q48; /* Keep x^n scaled */
|
||||
|
||||
/* Guard against overflow */
|
||||
if (result < 0 && result < (1LL << 62)) break;
|
||||
if (result > 0 && result > (1LL << 62)) {
|
||||
result = 9223372036854775806LL; /* Saturate at INT64_MAX-1 */
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Compute sqrt(@p x) via Newton-Raphson in Q48.16 fixed-point.
|
||||
*
|
||||
* Applies the Newton-Raphson recurrence:
|
||||
* @code
|
||||
* r_{n+1} = (r_n² + x) / (2 × r_n)
|
||||
* @endcode
|
||||
* for up to 20 iterations, with all arithmetic scaled to the Q48.16
|
||||
* representation (1.0 = 65536). Intermediate products use 8-bit right
|
||||
* shifts on both operands before squaring to prevent 64-bit overflow:
|
||||
* @code
|
||||
* guess_squared = (guess >> 8) * (guess >> 8)
|
||||
* @endcode
|
||||
* This sacrifices 16 bits of precision in the squared term, which is
|
||||
* acceptable because the division by @c 2*guess restores the correct
|
||||
* order of magnitude. Convergence is declared when successive iterates
|
||||
* differ by at most 1 in Q48.16 units (sub-nanosecond for timing uses).
|
||||
*
|
||||
* Special cases handled without iteration:
|
||||
* - @p x ≤ 0 → returns 0 (no imaginary results)
|
||||
* - @p x == 65536 (1.0) → returns 65536 exactly
|
||||
*
|
||||
* Initial guess: @c x/2 for @p x > 1.0, else @c 1.0. Values above 1.0
|
||||
* converge from above; values in (0, 1.0) converge from the 1.0 initial
|
||||
* guess.
|
||||
*
|
||||
* @param x Non-negative value in Q48.16 format (1.0 = 65536).
|
||||
* Negative values are treated as zero.
|
||||
* @return sqrt(@p x) in Q48.16 format; 0 for @p x ≤ 0.
|
||||
*/
|
||||
int64_t sqrt_q48(int64_t x) {
|
||||
const int64_t Q48 = 65536LL;
|
||||
|
||||
/* Handle special cases */
|
||||
if (x < 0) return 0;
|
||||
if (x == 0) return 0;
|
||||
if (x == Q48) return Q48; /* sqrt(1.0) = 1.0 */
|
||||
|
||||
/* Initial guess: average of x and 1 */
|
||||
int64_t guess = (x > Q48) ? (x >> 1) : Q48;
|
||||
|
||||
/* Newton-Raphson iterations */
|
||||
for (int i = 0; i < 20; i++) {
|
||||
/* next = (guess^2 + x) / (2 * guess)
|
||||
* In Q48.16: (guess^2) >> 16 + x, then divide by 2*guess
|
||||
*/
|
||||
int64_t guess_squared = (guess >> 8) * (guess >> 8); /* Avoid overflow */
|
||||
int64_t next = (guess_squared + (x >> 16)) / (2 * guess);
|
||||
|
||||
/* Check convergence */
|
||||
int64_t diff = next - guess;
|
||||
if (diff > -2 && diff < 2) return next;
|
||||
|
||||
guess = next;
|
||||
}
|
||||
|
||||
return guess;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Compute erf(@p x) via Abramowitz & Stegun 7.1.26 in Q48.16 fixed-point.
|
||||
*
|
||||
* Implements the rational approximation from Abramowitz & Stegun, formula 7.1.26:
|
||||
* @code
|
||||
* erf(x) ≈ sign(x) × sqrt(1 − exp(−x² × (4/π + a·x²) / (1 + a·x²)))
|
||||
* @endcode
|
||||
* where @c a ≈ 0.147 (in Q48.16: 9633). The maximum error of the true
|
||||
* mathematical approximation is @c |ε| < 1.5×10⁻⁷; the fixed-point
|
||||
* encoding introduces additional rounding of approximately 1 ULP in Q48.16.
|
||||
*
|
||||
* All arithmetic is performed in 64-bit integer Q48.16 format — no
|
||||
* floating-point operations. The computation sequence:
|
||||
* 1. Extract sign; work with @c |x|.
|
||||
* 2. Compute @c x² with 8-bit prescaling to prevent overflow:
|
||||
* @c x_sq = (x >> 8) * (x >> 8).
|
||||
* 3. Evaluate the A&S numerator @c (4/π + a·x²) and denominator @c (1 + a·x²).
|
||||
* 4. Compute @c ratio = numerator / denominator (shift to preserve Q48.16 scale).
|
||||
* 5. Compute @c exp_arg = −x² × ratio and evaluate @c exp_q48(exp_arg).
|
||||
* 6. Compute @c base = 1 − exp_val, clamp to [0, 1].
|
||||
* 7. Return @c sign × sqrt_q48(base).
|
||||
*
|
||||
* Q48.16 constants used:
|
||||
* - @c a_q48 = 9633 (0.147 × 65536)
|
||||
* - @c four_over_pi_q48 = 83328 (4/π × 65536 ≈ 1.2732 × 65536)
|
||||
*
|
||||
* Consumed by @c hotwords_speedup_estimate() to compute P(speedup > 1.1×)
|
||||
* and P(speedup > 2.0×) via the normal CDF identity:
|
||||
* @code
|
||||
* P(Z ≤ z) = 0.5 × (1 + erf(z / √2))
|
||||
* @endcode
|
||||
* where @c √2 in Q48.16 = 92681.
|
||||
*
|
||||
* @param x Input value in Q48.16 format (dimensionless Z-score;
|
||||
* 1.0 = 65536). May be negative; sign determines the sign
|
||||
* of the result.
|
||||
* @return erf(@p x) in Q48.16 format, in the range [−65536, +65536]
|
||||
* corresponding to (−1.0, +1.0).
|
||||
*/
|
||||
int64_t erf_q48(int64_t x) {
|
||||
const int64_t Q48 = 65536LL;
|
||||
|
||||
if (x == 0) return 0;
|
||||
|
||||
/* Compute sign and work with absolute value */
|
||||
int sign = (x < 0) ? -1 : 1;
|
||||
if (x < 0) x = -x;
|
||||
|
||||
/* Abramowitz & Stegun coefficient a ≈ 0.147 (in Q48.16) */
|
||||
const int64_t a_q48 = 9633LL; /* 0.147 * 65536 */
|
||||
const int64_t four_over_pi_q48 = 83328LL; /* 4/π * 65536 */
|
||||
|
||||
/* Compute x^2 in Q48.16 */
|
||||
int64_t x_sq = (x >> 8) * (x >> 8); /* Scaled to avoid overflow */
|
||||
|
||||
/* Compute denominator: 1 + a*x^2 */
|
||||
int64_t denom = Q48 + ((a_q48 * x_sq) >> 16);
|
||||
if (denom == 0) denom = 1; /* Guard divide-by-zero */
|
||||
|
||||
/* Compute numerator: 4/π + a*x^2 */
|
||||
int64_t numer = four_over_pi_q48 + ((a_q48 * x_sq) >> 16);
|
||||
|
||||
/* Compute ratio = numerator / denominator */
|
||||
int64_t ratio = (numer << 16) / denom;
|
||||
|
||||
/* Compute exp argument: -x^2 * ratio */
|
||||
int64_t exp_arg = -(x_sq * ratio / Q48);
|
||||
|
||||
/* Compute exp(-x^2 * ratio) */
|
||||
int64_t exp_val = exp_q48(exp_arg);
|
||||
|
||||
/* Compute base: 1 - exp(exp_arg) */
|
||||
int64_t base = Q48 - exp_val;
|
||||
if (base < 0) base = 0;
|
||||
if (base > Q48) base = Q48;
|
||||
|
||||
/* Compute sqrt(base) */
|
||||
int64_t sqrt_base = sqrt_q48(base);
|
||||
|
||||
/* Apply sign and return */
|
||||
return sign * sqrt_base;
|
||||
}
|
||||
Reference in New Issue
Block a user