TermiNexus
db3be6968f
memhub: alloc 失败打印归属 site/ctx 与存活块 order 分布（新增 dump_alive_orders）；free 被 _is_valid_ptr 拒绝时首次打印现场；接通 big_alloc_count gzip: ub: uchar*= buf 补空格（*= 被词法合并为单 OP，注解止于 uchar） 新增 traceback 模块；删除 dbgstack.vp；其余模块随新语法与尺寸口径同步
563 lines
20 KiB
Plaintext
563 lines
20 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 O(1): a three-constant
|
||
# * Cody-Waite step for |x| < 2^20, and a Payne-Hanek step (20-word table of 2/π) above. Both yield r ∈ [-π/4, π/4]
|
||
# * as a double-double fed to the truncated Taylor series.
|
||
# * - 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
|
||
|
||
# Argument reduction constants for sin/cos/tan: pi/2 split into HI+MID+LO and 2/pi (all exact doubles),
|
||
# plus the 2^20 limit below which Cody-Waite applies (so each fn*HK product stays exact in 53 bits).
|
||
define U_M_PIO2_HI: double = 1.57079632673412561417e+00
|
||
define U_M_PIO2_MID: double = 6.07710050650619224932e-11
|
||
define U_M_PIO2_LO: double = 3.52155986518320018165e-27
|
||
define U_M_INV_PIO2: double = 6.36619772367581382433e-01
|
||
define U_M_REDUCE_LIMIT: double = 1048576.0
|
||
|
||
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_core(x: double) -> double:
|
||
"""sin(r) Taylor series, 8 terms, valid on |r| <= pi/4."""
|
||
x2: double = x * x
|
||
result: double = x
|
||
term: double = x
|
||
sign: double = -1.0
|
||
coeffs: const [8]double = [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 _cos_core(x: double) -> double:
|
||
"""cos(r) Taylor series, 8 terms, valid on |r| <= pi/4."""
|
||
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 _mul_add(a: uint64, b: uint64, c: uint64, lo: rawptr[uint64], hi: rawptr[uint64]):
|
||
"""Write (hi:lo) = a*b + c using 32-bit limbs for the full 128-bit product."""
|
||
a0: uint64 = a & 0xFFFFFFFF
|
||
a1: uint64 = a >> 32
|
||
b0: uint64 = b & 0xFFFFFFFF
|
||
b1: uint64 = b >> 32
|
||
w0: uint64 = a0 * b0
|
||
t: uint64 = a1 * b0 + (w0 >> 32)
|
||
w1: uint64 = t & 0xFFFFFFFF
|
||
w2: uint64 = t >> 32
|
||
w1 = a0 * b1 + w1
|
||
h: uint64 = a1 * b1 + w2 + (w1 >> 32)
|
||
l: uint64 = (w1 << 32) | (w0 & 0xFFFFFFFF)
|
||
l2: uint64 = l + c
|
||
if l2 < l: h = h + 1
|
||
*lo = l2
|
||
*hi = h
|
||
|
||
|
||
def _two_prod(a: double, b: double, ph: rawptr[double], pl: rawptr[double]):
|
||
"""Write a*b as the exact double-double (ph + pl) via Dekker splitting (no FMA required)."""
|
||
split: double = 134217729.0
|
||
p: double = a * b
|
||
ca: double = split * a
|
||
ahi: double = ca - (ca - a)
|
||
alo: double = a - ahi
|
||
cb: double = split * b
|
||
bhi: double = cb - (cb - b)
|
||
blo: double = b - bhi
|
||
*ph = p
|
||
*pl = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo
|
||
|
||
|
||
def _rem_pio2(x: double, rh: rawptr[double], rl: rawptr[double]) -> int:
|
||
"""O(1) reduction of finite x >= 0: x = n*(pi/2) + r with |r| <= pi/4; returns n mod 4, writes r = rh + rl."""
|
||
if x < U_M_REDUCE_LIMIT:
|
||
# Cody-Waite three-step with pi/2 = HI + MID + LO; each fn*HK product is exact for fn < 2^20
|
||
q: double = x * U_M_INV_PIO2
|
||
k: int = int(q)
|
||
if q - double(k) > 0.5: k = k + 1
|
||
fn: double = double(k)
|
||
r: double = x - fn * U_M_PIO2_HI
|
||
w1: double = fn * U_M_PIO2_MID
|
||
y0: double = r - w1
|
||
w2: double = fn * U_M_PIO2_LO
|
||
rhv: double = y0 - w2
|
||
*rh = rhv
|
||
*rl = ((r - y0) - w1) + ((y0 - rhv) - w2)
|
||
return k & 3
|
||
# Payne-Hanek: a 20-word (1280-bit) 2/pi table; multiply the 53-bit mantissa and read off the fraction
|
||
twopi: const [20]uint64 = [
|
||
0xA2F9836E4E441529, 0xFC2757D1F534DDC0, 0xDB6295993C439041, 0xFE5163ABDEBBC561,
|
||
0xB7246E3A424DD2E0, 0x06492EEA09D1921C, 0xFE1DEB1CB129A73E, 0xE88235F52EBB4484,
|
||
0xE99C7026B45F7E41, 0x3991D639835339F4, 0x9C845F8BBDF9283B, 0x1FF897FFDE05980F,
|
||
0xEF2F118B5A0A6D1F, 0x6D367ECF27CB09B7, 0x4F463F669E5FEA2D, 0x7527BAC7EBE5F17B,
|
||
0x3D0739F78A5292EA, 0x6BFB5FB11F8D5D08, 0x56033046FC7B6BAB, 0xF0CFBC209AF4361D
|
||
]
|
||
xv: double = x
|
||
xb: uint64 = 0
|
||
*&xb = *&xv
|
||
ef: uint64 = (xb >> 52) & 0x7FF
|
||
m: uint64 = (xb & 0xFFFFFFFFFFFFF) | 0x10000000000000
|
||
e: int = int(ef) - 1075
|
||
shift: int = 1280 - e
|
||
qw: int = shift >> 6
|
||
s: int = shift & 63
|
||
b: [21]uint64
|
||
carry: uint64 = 0
|
||
for j in range(20):
|
||
lo: uint64 = 0
|
||
hi: uint64 = 0
|
||
_mul_add(m, twopi[19 - j], carry, &lo, &hi)
|
||
b[j] = lo
|
||
carry = hi
|
||
b[20] = carry
|
||
lo2: uint64 = b[qw] >> s
|
||
if s == 63: lo2 = (b[qw] >> 63) | (b[qw + 1] << 1)
|
||
k: int = int(lo2 & 3)
|
||
a1: uint64 = b[qw] & ((uint64(1) << s) - 1)
|
||
a2: uint64 = b[qw - 1]
|
||
a3: uint64 = b[qw - 2]
|
||
fhi: uint64 = a2
|
||
flo: uint64 = a3
|
||
if s != 0:
|
||
sh: int = 64 - s
|
||
fhi = (a1 << sh) | (a2 >> s)
|
||
flo = (a2 << sh) | (a3 >> s)
|
||
# fractional part f = fhi*2^-64 + flo*2^-128; split off the top 53 bits so df is exact
|
||
f53: uint64 = fhi >> 11
|
||
df: double = double(f53)
|
||
de: double = double(fhi - (f53 << 11))
|
||
fh: double = df * 1.1102230246251565e-16
|
||
fl: double = de * 5.421010862427522e-20 + double(flo) * 2.9387358770557188e-39
|
||
r0h: double = fh
|
||
r0l: double = fl
|
||
n: int = k
|
||
if fh > 0.5 or (fh == 0.5 and fl > 0.0):
|
||
n = (k + 1) & 3
|
||
r0h = fh - 1.0
|
||
ph: double = 0.0
|
||
pl: double = 0.0
|
||
_two_prod(r0h, U_M_PIO2_HI, &ph, &pl)
|
||
*rh = ph
|
||
*rl = pl + r0h * U_M_PIO2_MID + r0h * U_M_PIO2_LO + r0l * U_M_PIO2_HI
|
||
return n
|
||
|
||
|
||
def sin(x: double) -> double:
|
||
u: uint64 = 0
|
||
xv: double = x
|
||
*&u = *&xv
|
||
ef: uint64 = (u >> 52) & 0x7FF
|
||
if ef == 0x7FF: return x - x # inf/nan -> NaN
|
||
neg: int = int(u >> 63)
|
||
ax: double = x if neg == 0 else -x
|
||
rh: double = 0.0
|
||
rl: double = 0.0
|
||
n: int = _rem_pio2(ax, &rh, &rl)
|
||
sr: double = _sin_core(rh) + _cos_core(rh) * rl
|
||
cr: double = _cos_core(rh) - _sin_core(rh) * rl
|
||
res: double = 0.0
|
||
if n == 0: res = sr
|
||
elif n == 1: res = cr
|
||
elif n == 2: res = -sr
|
||
else: res = -cr
|
||
return -res if neg != 0 else res
|
||
|
||
|
||
def cos(x: double) -> double:
|
||
u: uint64 = 0
|
||
xv: double = x
|
||
*&u = *&xv
|
||
ef: uint64 = (u >> 52) & 0x7FF
|
||
if ef == 0x7FF: return x - x # inf/nan -> NaN
|
||
neg: int = int(u >> 63)
|
||
ax: double = x if neg == 0 else -x
|
||
rh: double = 0.0
|
||
rl: double = 0.0
|
||
n: int = _rem_pio2(ax, &rh, &rl)
|
||
sr: double = _sin_core(rh) + _cos_core(rh) * rl
|
||
cr: double = _cos_core(rh) - _sin_core(rh) * rl
|
||
if n == 0: return cr
|
||
elif n == 1: return -sr
|
||
elif n == 2: return -cr
|
||
else: return sr
|
||
|
||
|
||
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: return float(sin(double(x)))
|
||
def cosf(x: float) -> float: return float(cos(double(x)))
|
||
|
||
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
|