441 lines
15 KiB
Plaintext
441 lines
15 KiB
Plaintext
# * Soft Math Library for Viper
|
||
# * Implements a self-contained software floating-point math runtime for float (binary32) and double (binary64),
|
||
# * without external libm dependencies. Designed for embedded, bare-metal, or custom VM environments where hardware
|
||
# * FPU math routines or system-provided standard math functions are unavailable, undesirable, or require fully
|
||
# * auditable source.
|
||
# *
|
||
# * Scope and feature set:
|
||
# * - Fundamental mathematical constants: π, e, π/2, π/4, reciprocal π, 2/π, ln(2), ln(10), 2/sqrt(π), DBL_MAX.
|
||
# * - Angle conversion: radians ↔ degrees helpers.
|
||
# * - Trigonometric functions: sin, cos, tan, sinf, cosf, tanf. Argument reduction is performed using subtraction of
|
||
# * 2π until the value falls into [-π, π]; polynomial approximations are truncated series tuned for the target interval.
|
||
# * - Inverse trigonometric functions: asin, acos, atan, atan2. asin uses a minimax polynomial suitable for [-1, 1];
|
||
# * atan relies on a Taylor core limited to |x| ≤ 1 combined with range folding via reciprocal and quadrant logic.
|
||
# * - Hyperbolic trigonometry: sinh, cosh, tanh, built atop exp and exp(-x).
|
||
# * - Exponential and logarithm: exp uses range reduction into [-ln2/2, ln2/2] followed by a Taylor polynomial for the
|
||
# * reduced residual, then scales by powers of two. Natural log implements a binary normalization pass to bring the
|
||
# * mantissa into [1, 2), then uses the rational substitution y=(s-1)/(s+1) with a convergent series.
|
||
# * - Derived logarithms and exponentials: log10, log2, exp2, expm1, log1p.
|
||
# * - General power: pow for real bases and exponents, implemented as exp(y * log(x)), with special-case handling for
|
||
# * zero and negative bases. cbrt is a wrapper over pow using the 1/3 exponent.
|
||
# * - Root and distance: Newton–Raphson sqrt / sqrtf with fixed 10 iterations; hypot for Euclidean norm.
|
||
# * - Rounding primitives: floor, ceil, round, trunc, with float variants floorf, ceilf.
|
||
# * - Modulo: fmod, fmodf following the C remainder convention where the result carries the sign of the dividend.
|
||
# * - Decomposition: modf to split a double into integer and fractional parts via a pointer output.
|
||
# * - Inverse hyperbolics: asinh, acosh, atanh expressed through log and sqrt compositions.
|
||
# * - Special functions: Gamma using Stirling’s approximation; error function erf / complementary erfc using a simple
|
||
# * rational exponential approximation.
|
||
# * - Combinatorics: factorial, binomial combination, permutation, returning double for wider dynamic range.
|
||
# * - Integer helpers: abs for int, labs for long, fabs / fabsf for floating-point absolute value.
|
||
# * - IEEE-754 predicate helpers: isnan and isinf implemented via type punning union to inspect binary64 bit patterns.
|
||
# * These avoid numeric conversion which would destroy NaN/inf payloads; they examine exponent and mantissa bits directly.
|
||
# *
|
||
# * Precision notes:
|
||
# * All polynomial cores are truncated series, not infinite expansions. Errors are bounded for the reduced interval but
|
||
# * will accumulate through composition (e.g. pow depends on log and exp; trigonometric inverses depend on asin/atan cores).
|
||
# * This implementation prioritizes determinism and source transparency over last-bit IEEE compliance. It is not a drop-in
|
||
# * replacement for a high-quality libm with accurate argument reduction, correct rounding, and exhaustive edge-case
|
||
# * handling. Input values outside the designed convergence domains will see degraded accuracy.
|
||
# *
|
||
# * Range and edge behavior:
|
||
# * exp clamps extreme inputs: values above approximately 709 return DBL_MAX, values below approximately -745 return 0.0,
|
||
# * preventing uncontrolled overflow in downstream code such as softmax. log returns NEG_1E308 for non-positive inputs
|
||
# * as a sentinel instead of NaN. Division-by-zero conditions in tan/tanf return zero as a simple fallback; other edge
|
||
# * cases propagate NaN where the underlying computation produces it. The isnan/isinf predicates use direct bit inspection
|
||
# * and respect the binary64 encoding rules for quiet/signaling NaN and signed infinities.
|
||
# *
|
||
# * Usage notes:
|
||
# * This library assumes Viper’s double maps to IEEE 754 binary64 and float maps to IEEE 754 binary32. Type punning in
|
||
# * isnan/isinf relies on a union reinterpretation of the underlying memory layout; this is portable only on targets
|
||
# * with standard little-endian or big-endian IEEE float storage and must be reviewed for exotic memory models.
|
||
# * Polynomial loops are manually folded into fixed for-loops over coefficient arrays to remove redundant source lines
|
||
# * while preserving the original floating-point evaluation order, which matters for reproducibility and error accumulation.
|
||
# *
|
||
# * License: MIT
|
||
|
||
|
||
define U_M_PI: double = 3.14159265358979323846
|
||
define U_M_E: double = 2.71828182845904523536
|
||
define U_M_PI_2: double = 1.57079632679489661923
|
||
define U_M_PI_4: double = 0.78539816339744830962
|
||
define U_M_1_PI: double = 0.31830988618379067154
|
||
define U_M_2_PI: double = 0.63661977236758134308
|
||
|
||
define U_M_LN2: double = 0.69314718055994530942
|
||
define U_M_LN10: double = 2.30258509299404568402
|
||
|
||
define U_M_2_SQRT_PI: double = 1.77245385090551602730
|
||
|
||
define DBL_MAX: double = 1.7976931348623157e308
|
||
define NEG_1E308: double = -1e308
|
||
|
||
def radians(degrees: double) -> double: return degrees * (U_M_PI / 180.0)
|
||
def degrees(radians: double) -> double: return radians * (180.0 / U_M_PI)
|
||
|
||
def sin(x: double) -> double:
|
||
while x > U_M_PI:
|
||
x -= 2.0 * U_M_PI
|
||
while x < -U_M_PI:
|
||
x += 2.0 * U_M_PI
|
||
x2: double = x * x
|
||
result: double = x
|
||
term: double = x
|
||
coeffs: const [8]double = [6.0, 20.0, 42.0, 72.0, 110.0, 156.0, 210.0, 272.0]
|
||
sign: double = -1.0
|
||
for c in coeffs:
|
||
term = term * x2 / c
|
||
result += sign * term
|
||
sign = -sign
|
||
return result
|
||
|
||
def cos(x: double) -> double:
|
||
while x > U_M_PI:
|
||
x -= 2.0 * U_M_PI
|
||
while x < -U_M_PI:
|
||
x += 2.0 * U_M_PI
|
||
x2: double = x * x
|
||
result: double = 1.0
|
||
term: double = 1.0
|
||
sign: double = -1.0
|
||
coeffs: const [8]double = [2.0, 12.0, 30.0, 56.0, 90.0, 132.0, 182.0, 240.0]
|
||
for c in coeffs:
|
||
term = term * x2 / c
|
||
result += sign * term
|
||
sign = -sign
|
||
return result
|
||
|
||
def tan(x: double) -> double: return sin(x) / cos(x)
|
||
|
||
|
||
def asin(x: double) -> double:
|
||
if x > 1.0: x = 1.0
|
||
if x < -1.0: x = -1.0
|
||
x2: double = x * x
|
||
result: double = x
|
||
result += 1.0 * x * x2 / 6.0
|
||
result += 3.0 * x * x2 * x2 / 40.0
|
||
result += 5.0 * x * x2 * x2 * x2 / 112.0
|
||
result += 35.0 * x * x2 * x2 * x2 * x2 / 1152.0
|
||
result += 63.0 * x * x2 * x2 * x2 * x2 * x2 / 2816.0
|
||
return result
|
||
|
||
|
||
def acos(x: double) -> double: return U_M_PI / 2.0 - asin(x)
|
||
|
||
|
||
def _atan_core(x: double) -> double:
|
||
x2: double = x * x
|
||
result: double = x
|
||
term: double = x
|
||
sign: double = -1.0
|
||
coeffs: const [10]double = [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0]
|
||
for c in coeffs:
|
||
term = term * x2 / c
|
||
result += sign * term
|
||
sign = -sign
|
||
return result
|
||
|
||
|
||
def atan(x: double) -> double:
|
||
ax: double = x
|
||
if ax < 0.0:
|
||
ax = -ax
|
||
if ax > 1.0:
|
||
return 0.0 - _atan_core(1.0 / ax) + U_M_PI_2 if x > 0.0 else _atan_core(1.0 / ax) - U_M_PI_2
|
||
return _atan_core(x) if x >= 0.0 else 0.0 - _atan_core(ax)
|
||
|
||
|
||
def atan2(y: double, x: double) -> double:
|
||
if x > 0.0: return atan(y / x)
|
||
elif x < 0.0:
|
||
if y >= 0.0: return atan(y / x) + U_M_PI
|
||
else: return atan(y / x) - U_M_PI
|
||
else:
|
||
if y > 0.0: return U_M_PI_2
|
||
elif y < 0.0: return 0.0 - U_M_PI_2
|
||
else: return 0.0
|
||
|
||
|
||
def sinh(x: double) -> double:
|
||
ex: double = exp(x)
|
||
e_negx: double = exp(-x)
|
||
return (ex - e_negx) / 2.0
|
||
|
||
|
||
def cosh(x: double) -> double:
|
||
ex: double = exp(x)
|
||
e_negx: double = exp(-x)
|
||
return (ex + e_negx) / 2.0
|
||
|
||
|
||
def tanh(x: double) -> double:
|
||
ex: double = exp(x)
|
||
e_negx: double = exp(-x)
|
||
return (ex - e_negx) / (ex + e_negx)
|
||
|
||
|
||
def exp(x: double) -> double:
|
||
"""e^x
|
||
Range reduction Short Taylor: Directly expanding x with 18 terms of Taylor series diverges when |x| is slightly large
|
||
(x=-10 gives 54.5, true value 4.5e-5; x=-30 gives 3.7e10), on the negative side it gives huge positive values,
|
||
which can blow up the denominator in softmax → all probabilities ≈ 0. Instead, first do
|
||
x = k*ln2 + r (|r| <= ln2/2), then calculate exp(r) and multiply by 2^k.
|
||
"""
|
||
# First do a "finite and within range" positive check: comparisons with NaN are always false and would skip
|
||
# the following two checks, and int(NaN) on x86 becomes INT_MIN -> the later 2^k scaling loop would run
|
||
# 2.1 billion times, appearing as "frozen". Therefore NaN/overflow must be blocked here.
|
||
if x <= 709.0 and x >= -745.0:
|
||
pass
|
||
else:
|
||
if x > 709.0: return DBL_MAX
|
||
if x < -745.0: return 0.0
|
||
return x # NaN: Return as is, let NaN propagate instead of crashing
|
||
# Reduction: Fold x into [-ln2/2, ln2/2], accumulate the integer part into k (inline rounding to avoid forward reference)
|
||
ln2: double = U_M_LN2
|
||
q: double = x / ln2
|
||
k: int = int(q)
|
||
frac: double = q - double(k)
|
||
if frac > 0.5: k = k + 1
|
||
elif frac < -0.5: k = k - 1
|
||
r: double = x - double(k) * ln2
|
||
# IE # exp(r), |r| <= 0.3466: after 12 terms, truncation error < 1e-17
|
||
result: double = 1.0
|
||
term: double = 1.0
|
||
n: int = 1
|
||
while n <= 12:
|
||
term = term * r / double(n)
|
||
result = result + term
|
||
n += 1
|
||
while k > 0:
|
||
result = result * 2.0
|
||
k = k - 1
|
||
while k < 0:
|
||
result = result * 0.5
|
||
k = k + 1
|
||
return result
|
||
|
||
|
||
def log(x: double) -> double:
|
||
if x <= 0.0: return NEG_1E308
|
||
exponent: int = 0
|
||
while x >= 2.0:
|
||
x /= 2.0
|
||
exponent += 1
|
||
while x < 1.0:
|
||
x *= 2.0
|
||
exponent -= 1
|
||
y: double = (x - 1.0) / (x + 1.0)
|
||
y2: double = y * y
|
||
result: double = 2.0 * y
|
||
term: double = 2.0 * y
|
||
coeffs: const [10]double = [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0]
|
||
for c in coeffs:
|
||
term *= y2
|
||
result += term / c
|
||
result += double(exponent) * U_M_LN2
|
||
return result
|
||
|
||
|
||
def sqrt(x: double) -> double:
|
||
if x < 0.0: return 0.0
|
||
if x == 0.0: return 0.0
|
||
guess: double = x / 2.0
|
||
for _ in range(10):
|
||
guess = (guess + x / guess) / double(2.0)
|
||
return guess
|
||
|
||
|
||
def abs(x: int) -> int: return x if x >= 0 else -x
|
||
def fabs(x: double) -> double: return x if x >= double(0.0) else -x
|
||
def labs(x: long) -> long: return x if x >= 0 else -x
|
||
|
||
def factorial(n: int) -> double:
|
||
if n < 0: return 0.0
|
||
if n == 0 or n == 1: return 1.0
|
||
result: double = 1.0
|
||
for i in range(2, n + 1):
|
||
result *= i
|
||
return result
|
||
|
||
|
||
def combination(n: int, k: int) -> double:
|
||
if k < 0 or k > n: return 0.0
|
||
if k == 0 or k == n: return 1.0
|
||
k = k if k < n - k else n - k
|
||
result: double = 1.0
|
||
for i in range(1, k + 1):
|
||
result = result * (n - k + i) / i
|
||
return result
|
||
|
||
|
||
def permutation(n: int, k: int) -> double:
|
||
if k < 0 or k > n: return 0.0
|
||
result: double = 1.0
|
||
for i in range(0, k):
|
||
result *= (n - i)
|
||
return result
|
||
|
||
def pow(x: double, y: double) -> double:
|
||
if x == 0.0: return 0.0 if y > 0.0 else 1.0
|
||
if x < 0.0: return exp(y * log(-x))
|
||
return exp(y * log(x))
|
||
|
||
|
||
def powf(x: double, y: double) -> double: return pow(x, y)
|
||
def cbrt(x: double ) -> double: return pow(x, 1.0 / 3.0)
|
||
def hypot(x: double, y: double) -> double: return sqrt(x * x + y * y)
|
||
|
||
|
||
def floor(x: double) -> double:
|
||
if x >= 0.0: return double(int(x))
|
||
i: int = int(x)
|
||
if x < double(i): return double(i - 1)
|
||
else: return double(i)
|
||
|
||
|
||
def ceil(x: double) -> double:
|
||
if x < 0.0: return double(int(x))
|
||
i: int = int(x)
|
||
if x > double(i): return double(i + 1)
|
||
else: return double(i)
|
||
|
||
|
||
def round(x: double) -> double: return floor(x + 0.5)
|
||
def trunc(x: double) -> double: return floor(x) if x >= 0.0 else ceil(x)
|
||
|
||
|
||
def fmod(x: double, y: double) -> double:
|
||
if y == 0.0: return 0.0
|
||
result: double = x - y * floor(x / y)
|
||
if result < 0.0: result += y
|
||
return result
|
||
|
||
|
||
def fmodf(x: float, y: float) -> float:
|
||
if y == 0.0: return 0.0
|
||
result: float = x - y * floorf(x / y)
|
||
if result < 0.0: result += y
|
||
return result
|
||
|
||
|
||
def modf(x: double, iptr: rawptr[double]) -> double:
|
||
int_part: double = floor(x)
|
||
*iptr = int_part
|
||
return x - int_part
|
||
|
||
|
||
def log10(x: double) -> double: return log(x) / U_M_LN10
|
||
def log2(x: double) -> double: return log(x) / U_M_LN2
|
||
|
||
|
||
def exp2(x: double) -> double: return pow(2.0, x)
|
||
def expm1(x: double) -> double: return exp(x) - 1.0
|
||
|
||
def log1p(x: double) -> double: return log(1.0 + x)
|
||
|
||
def asinh(x: double) -> double: return log(x + sqrt(x * x + 1.0))
|
||
def acosh(x: double) -> double: return log(x + sqrt(x * x - 1.0))
|
||
def atanh(x: double) -> double: return (log(1.0 + x) - log(1.0 - x)) / 2.0
|
||
|
||
|
||
def gamma(x: double) -> double:
|
||
if x <= 0.0: return 0.0
|
||
if x == 1.0 or x == 2.0: return 1.0
|
||
ln_gamma: double = 0.5 * log(2.0 * U_M_PI / x) + x * log(x / U_M_E)
|
||
return exp(ln_gamma)
|
||
|
||
|
||
def erf(x: double) -> double:
|
||
sign: double = 1.0 if x >= 0.0 else -1.0
|
||
abs_x: double = fabs(x)
|
||
t0: double = 4.0 * abs_x * abs_x / U_M_PI
|
||
exp_term: double = exp(-t0)
|
||
return sign * sqrt(1.0 - exp_term)
|
||
|
||
|
||
def erfc(x: double) -> double:
|
||
return 1.0 - erf(x)
|
||
|
||
def sqrtf(x: float) -> float:
|
||
if x < 0.0: return 0.0
|
||
if x == 0.0: return 0.0
|
||
guess: float = x / 2.0
|
||
for i in range(10):
|
||
guess = (guess + x / guess) / 2.0
|
||
return guess
|
||
|
||
def sinf(x: float) -> float:
|
||
while x > 3.14159265:
|
||
x -= 6.28318530
|
||
while x < -3.14159265:
|
||
x += 6.28318530
|
||
x2: float = x * x
|
||
result: float = x
|
||
term: float = x
|
||
sign: float = -1.0
|
||
coeffs: const [8]float = [6.0, 20.0, 42.0, 72.0, 110.0, 156.0, 210.0, 272.0]
|
||
for c in coeffs:
|
||
term = term * x2 / c
|
||
result += sign * term
|
||
sign = -sign
|
||
return result
|
||
|
||
|
||
def cosf(x: float) -> float:
|
||
while x > 3.14159265:
|
||
x -= 6.28318530
|
||
while x < -3.14159265:
|
||
x += 6.28318530
|
||
x2: float = x * x
|
||
result: float = 1.0
|
||
term: float = 1.0
|
||
sign: float = -1.0
|
||
coeffs: const [8]float = [2.0, 12.0, 30.0, 56.0, 90.0, 132.0, 182.0, 240.0]
|
||
for c in coeffs:
|
||
term = term * x2 / c
|
||
result += sign * term
|
||
sign = -sign
|
||
return result
|
||
|
||
def tanf(x: float) -> float:
|
||
c: float = cosf(x)
|
||
if c == 0.0: return 0.0
|
||
return sinf(x) / c
|
||
|
||
def fabsf(x: float) -> float: return x if x >= 0.0 else -x
|
||
|
||
def floorf(x: float) -> float:
|
||
if x >= 0.0: return float(int(x))
|
||
i: float = float(int(x))
|
||
if x < i: return i - 1
|
||
else: return i
|
||
|
||
def ceilf(x: float) -> float:
|
||
if x < 0.0: return float(int(x))
|
||
i: float = float(int(x))
|
||
if x > i: return i + 1
|
||
else: return i
|
||
|
||
|
||
union _U:
|
||
d: double
|
||
u: uint64
|
||
|
||
def isnan(x: double) -> int:
|
||
# IEEE-754 bit reinterpretation (can't use uint64(x) numerical conversion: that's fptosi, which loses the bit pattern completely
|
||
# → exp_part is always 0 → always returns 0). Bitwise copy
|
||
xv: double = x
|
||
u: uint64 = 0
|
||
*&u = *&xv
|
||
exp_part: uint64 = (u >> 52) & 0x7FF
|
||
mantissa_part: uint64 = u & 0xFFFFFFFFFFFFF
|
||
return 1 if exp_part == 0x7FF and mantissa_part != 0 else 0
|
||
|
||
def isinf(x: double) -> int:
|
||
xv2: double = x
|
||
u2: uint64 = 0
|
||
*&u2 = *&xv2
|
||
exp_part2: uint64 = (u2 >> 52) & 0x7FF
|
||
mantissa_part2: uint64 = u2 & 0xFFFFFFFFFFFFF
|
||
return 1 if exp_part2 == 0x7FF and mantissa_part2 == 0 else 0
|