diff --git a/libmath/desk/lib/math-chebyshev.hoon b/libmath/desk/lib/math-chebyshev.hoon new file mode 100644 index 0000000..64ce4f1 --- /dev/null +++ b/libmath/desk/lib/math-chebyshev.hoon @@ -0,0 +1,2885 @@ +:: +:::: Chebyshev/Fixed-Polynomial Math Library +:: +:: Jet-compatible transcendental functions using fixed polynomial +:: evaluation. Every floating-point operation happens in the same +:: order with the same intermediate values in both Hoon and the C jet. +:: +:: Reference: musl libc (MIT license) +:: Coefficients expressed as hex literals for bit-exactness. +:: +:: For each precision level (@rh, @rs, @rd), we provide: +:: - Polynomial coefficients as hex constants +:: - Argument reduction helpers +:: - Kernel functions (sindf, cosdf) +:: - Public API (sin, cos, tan, exp) +:: - Chebyshev polynomial generator (+cheb) +:: +|% +:: +:: ================================================================ +:: HALF PRECISION (@rh) - 16-bit IEEE 754 +:: ================================================================ +:: +:: Half precision has limited range and precision, so we use +:: reduced-degree polynomials (fewer terms). +:: +++ rh + |% + :: + :: Constants (hex-exact) + :: + ++ pio2-hi `@rh`0x3e48 :: pi/2 high bits ~1.5703125 + ++ pio2-lo `@rh`0x0500 :: pi/2 low bits (correction) + ++ pio4 `@rh`0x3a48 :: pi/4 ~0.785 + ++ invpio2 `@rh`0x3518 :: 2/pi ~0.6367 + ++ ln2hi `@rh`0x398c :: ln(2) high ~0.693 + ++ ln2lo `@rh`0x0000 :: ln(2) low (negligible at half) + ++ invln2 `@rh`0x3dc5 :: 1/ln(2) ~1.4425 + :: + :: Sin coefficients for half precision (degree 5) + :: sin(x) ~ x - x^3/6 + x^5/120 + :: + ++ sin-s1 `@rh`0xb555 :: -1/6 = -0.1666... + ++ sin-s2 `@rh`0x2444 :: 1/120 = 0.00833... + :: + :: Cos coefficients for half precision (degree 4) + :: cos(x) ~ 1 - x^2/2 + x^4/24 + :: + ++ cos-c0 `@rh`0xb800 :: -1/2 = -0.5 + ++ cos-c1 `@rh`0x2955 :: 1/24 = 0.04166... + :: + :: Exp coefficients for half precision (degree 3) + :: + ++ exp-p1 `@rh`0x3800 :: 1/2 + ++ exp-p2 `@rh`0x3155 :: 1/6 ~0.1666 + :: + :: +cheb: Chebyshev polynomial T_n(x) + :: + :: Returns the Chebyshev polynomial of the first kind. + :: T_n(x) = 2x*T_{n-1}(x) - T_{n-2}(x) with T_0=1, T_1=x + :: + ++ cheb + |= n=@ud + |= x=@rh + ^- @rh + ?: =(0 n) `@rh`0x3c00 :: 1.0 + ?: =(1 n) x + %+ sub:^rh + (mul:^rh `@rh`0x4000 (mul:^rh x $(n (dec n)))) :: 2*x*T_{n-1} + $(n (sub n 2)) :: T_{n-2} + :: + :: +sindf: Kernel sin for |x| <= pi/4 + :: + ++ sindf + |= x=@rh + ^- @rh + =/ z (mul:^rh x x) + :: sin(x) = x + x^3*S1 + x^5*S2 + :: = x*(1 + z*(S1 + z*S2)) + =/ r (add:^rh sin-s1 (mul:^rh z sin-s2)) + (add:^rh x (mul:^rh (mul:^rh x z) r)) + :: + :: +cosdf: Kernel cos for |x| <= pi/4 + :: + ++ cosdf + |= x=@rh + ^- @rh + =/ z (mul:^rh x x) + :: cos(x) = 1 + z*C0 + z^2*C1 + :: = 1 + z*(C0 + z*C1) + =/ r (add:^rh cos-c0 (mul:^rh z cos-c1)) + (add:^rh `@rh`0x3c00 (mul:^rh z r)) + :: + :: +sin: Sine with argument reduction + :: + ++ sin + |= x=@rh + ^- @rh + =/ bits `@`x + =/ ax `@rh`(dis bits 0x7fff) + :: tiny: sin(x) ~ x + ?: (lth:^rh ax `@rh`0x1000) x + :: inf/nan + ?: (gte:^rh ax `@rh`0x7c00) (sub:^rh x x) + :: Argument reduction (simplified for half precision range) + =/ n (div:^rh ax pio2-hi) + =/ ni `@`(rsh [0 10] `@`n) + =/ y (sub:^rh x (mul:^rh (sun:^rh ni) pio2-hi)) + ?- (mod ni 4) + %0 (sindf y) + %1 (cosdf y) + %2 (sub:^rh `@rh`0x0 (sindf y)) + %3 (sub:^rh `@rh`0x0 (cosdf y)) + == + :: + :: +cos: Cosine with argument reduction + :: + ++ cos + |= x=@rh + ^- @rh + =/ bits `@`x + =/ ax `@rh`(dis bits 0x7fff) + :: tiny: cos(x) ~ 1 + ?: (lth:^rh ax `@rh`0x1000) `@rh`0x3c00 + :: inf/nan + ?: (gte:^rh ax `@rh`0x7c00) (sub:^rh x x) + :: Argument reduction + =/ n (div:^rh ax pio2-hi) + =/ ni `@`(rsh [0 10] `@`n) + =/ y (sub:^rh x (mul:^rh (sun:^rh ni) pio2-hi)) + ?- (mod ni 4) + %0 (cosdf y) + %1 (sub:^rh `@rh`0x0 (sindf y)) + %2 (sub:^rh `@rh`0x0 (cosdf y)) + %3 (sindf y) + == + :: + :: +tan: Tangent + :: + ++ tan + |= x=@rh + ^- @rh + (div:^rh (sin x) (cos x)) + :: + :: +exp: Exponential (simplified for half precision) + :: + ++ exp + |= x=@rh + ^- @rh + =/ bits `@`x + :: x = 0 -> 1 + ?: =(x `@rh`0x0) `@rh`0x3c00 + :: +inf -> +inf + ?: =(bits 0x7c00) `@rh`0x7c00 + :: -inf -> 0 + ?: =(bits 0xfc00) `@rh`0x0 + :: nan -> nan + ?: (gth:^rh `@rh`(dis bits 0x7fff) `@rh`0x7c00) x + :: overflow check (x > 11.09) + ?: (gth:^rh x `@rh`0x498a) `@rh`0x7c00 + :: underflow check (x < -17) + ?: (lth:^rh x `@rh`0xcc40) `@rh`0x0 + :: exp(r) ~ 1 + r + r^2/2 + r^3/6 + =/ r2 (mul:^rh x x) + =/ p (add:^rh exp-p1 (mul:^rh x exp-p2)) + (add:^rh `@rh`0x3c00 (mul:^rh x (add:^rh `@rh`0x3c00 (mul:^rh x p)))) + :: + :: +log: Natural logarithm (simplified for half precision) + :: + ++ log + |= x=@rh + ^- @rh + =/ bits `@`x + :: Special cases + ?: =(x `@rh`0x0) `@rh`0xfc00 :: log(0) = -inf + ?: =(bits 0x7c00) `@rh`0x7c00 :: log(+inf) = +inf + ?: (gth:^rh `@rh`(dis bits 0x7fff) `@rh`0x7c00) x :: nan + ?: !=(0 (rsh [0 15] bits)) (sub:^rh x x) :: log(negative) = nan + :: + :: Extract exponent and mantissa + :: Half: 1 sign, 5 exp (bias 15), 10 mantissa + :: + =/ exp-bits (dis (rsh [0 10] bits) 0x1f) + =/ mant-bits (dis bits 0x3ff) + ?: =(exp-bits 0) `@rh`0xfc00 :: subnormal + =/ k (sub exp-bits 15) + =/ m `@rh`(con 0x3c00 mant-bits) :: 1.0 + mantissa + :: + :: s = (m-1)/(m+1), log(m) ~ 2*s*(1 + s^2*Lg1) + :: + =/ one `@rh`0x3c00 + =/ two `@rh`0x4000 + =/ f (sub:^rh m one) + =/ s (div:^rh f (add:^rh two f)) + =/ s2 (mul:^rh s s) + =/ lg1 `@rh`0x3955 :: ~0.666 (2/3) + =/ log-m (mul:^rh (mul:^rh two s) (add:^rh one (mul:^rh s2 lg1))) + :: + =/ kf (sun:^rh k) + (add:^rh (mul:^rh kf ln2hi) log-m) + :: + ++ log-2 + |= x=@rh + ^- @rh + (mul:^rh (log x) invln2) + :: + ++ log-10 + |= x=@rh + ^- @rh + =/ invlog10 `@rh`0x36f3 :: ~0.4343 + (mul:^rh (log x) invlog10) + :: + :: +sqt: Square root (Newton-Raphson, 3 iterations for half) + :: + ++ sqt + |= x=@rh + ^- @rh + =/ bits `@`x + ?: =(x `@rh`0x0) `@rh`0x0 + ?: =(bits 0x7c00) `@rh`0x7c00 + ?: (gth:^rh `@rh`(dis bits 0x7fff) `@rh`0x7c00) x + ?: !=(0 (rsh [0 15] bits)) (sub:^rh x x) + :: Initial guess + =/ g0 `@rh`(add (rsh [0 1] bits) 0x1c00) + =/ half `@rh`0x3800 :: 0.5 + =/ g g0 + =. g (mul:^rh half (add:^rh g (div:^rh x g))) + =. g (mul:^rh half (add:^rh g (div:^rh x g))) + =. g (mul:^rh half (add:^rh g (div:^rh x g))) + g + :: + ++ sqrt sqt + :: + :: +atan: Inverse tangent (simplified for half precision) + :: + ++ atan + |= x=@rh + ^- @rh + =/ bits `@`x + =/ sign (rsh [0 15] bits) + =/ ax `@rh`(dis bits 0x7fff) + :: + ?: (gte:^rh ax `@rh`0x7c00) + ?: =(bits 0x7c00) `@rh`0x3e48 :: pi/2 + ?: =(bits 0xfc00) `@rh`0xbe48 :: -pi/2 + x + :: + :: Simple polynomial for small |x| + :: atan(x) ~ x - x^3/3 + x^5/5 + :: + =/ one `@rh`0x3c00 + =/ pio2 `@rh`0x3e48 :: pi/2 + =/ at1 `@rh`0xb555 :: -1/3 + =/ at2 `@rh`0x3266 :: ~0.2 (1/5) + :: + =/ result=@rh + ?: (lte:^rh ax one) + :: Small: polynomial directly + =/ x2 (mul:^rh ax ax) + =/ x3 (mul:^rh x2 ax) + =/ r (add:^rh at1 (mul:^rh x2 at2)) + (add:^rh ax (mul:^rh x3 r)) + :: Large: atan(x) = pi/2 - atan(1/x) + =/ t (div:^rh one ax) + =/ t2 (mul:^rh t t) + =/ t3 (mul:^rh t2 t) + =/ r (add:^rh at1 (mul:^rh t2 at2)) + (sub:^rh pio2 (add:^rh t (mul:^rh t3 r))) + :: + ?:(=(sign 0) result (sub:^rh `@rh`0x0 result)) + :: + ++ asin + |= x=@rh + ^- @rh + =/ bits `@`x + =/ ax `@rh`(dis bits 0x7fff) + =/ one `@rh`0x3c00 + ?: (gth:^rh ax one) (sub:^rh x x) + ?: =(ax one) + ?: =(0 (rsh [0 15] bits)) `@rh`0x3e48 :: pi/2 + `@rh`0xbe48 :: -pi/2 + (atan (div:^rh x (sqt (sub:^rh one (mul:^rh x x))))) + :: + ++ acos + |= x=@rh + ^- @rh + =/ bits `@`x + =/ ax `@rh`(dis bits 0x7fff) + =/ one `@rh`0x3c00 + =/ pi `@rh`0x4248 :: pi + ?: (gth:^rh ax one) (sub:^rh x x) + ?: =(x one) `@rh`0x0 + ?: =(x `@rh`0xbc00) pi :: acos(-1) = pi + =/ s (sqt (sub:^rh one (mul:^rh x x))) + ?: (gte:^rh x `@rh`0x0) + (atan (div:^rh s x)) + (add:^rh pi (atan (div:^rh s x))) + :: + ++ atan2 + |= [y=@rh x=@rh] + ^- @rh + =/ pi `@rh`0x4248 + =/ pio2 `@rh`0x3e48 + =/ zero `@rh`0x0 + ?: (gth:^rh x zero) + (atan (div:^rh y x)) + ?: &((lth:^rh x zero) (gte:^rh y zero)) + (add:^rh (atan (div:^rh y x)) pi) + ?: &((lth:^rh x zero) (lth:^rh y zero)) + (sub:^rh (atan (div:^rh y x)) pi) + ?: &(=(x zero) (gth:^rh y zero)) + pio2 + ?: &(=(x zero) (lth:^rh y zero)) + (sub:^rh zero pio2) + zero + :: + ++ pow-n + |= [x=@rh n=@rh] + ^- @rh + =/ one `@rh`0x3c00 + ?: =(n `@rh`0x0) one + =/ ni (abs:si (need (toi:^rh n))) + =/ neg (lth:^rh n `@rh`0x0) + =/ result one + =/ base x + =/ i 0 + |- + ?: |(=(ni 0) (gth i 15)) + ?:(neg (div:^rh one result) result) + =/ new-result + ?: =(1 (dis ni 1)) + (mul:^rh result base) + result + $(ni (rsh [0 1] ni), base (mul:^rh base base), result new-result, i +(i)) + :: + ++ pow + |= [x=@rh n=@rh] + ^- @rh + =/ ni (toi:^rh n) + ?: &(?=(^ ni) =(n (san:^rh (need ni)))) + (pow-n x n) + (exp (mul:^rh n (log x))) + :: + :: ============================================================ + :: HYPERBOLIC FUNCTIONS + :: ============================================================ + :: + :: +sinh: Hyperbolic sine + :: sinh(x) = (exp(x) - exp(-x)) / 2 + :: + ++ sinh + |= x=@rh + ^- @rh + =/ half `@rh`0x3800 :: 0.5 + =/ ep (exp x) + =/ em (exp (sub:^rh `@rh`0x0 x)) + (mul:^rh half (sub:^rh ep em)) + :: + :: +cosh: Hyperbolic cosine + :: cosh(x) = (exp(x) + exp(-x)) / 2 + :: + ++ cosh + |= x=@rh + ^- @rh + =/ half `@rh`0x3800 :: 0.5 + =/ ep (exp x) + =/ em (exp (sub:^rh `@rh`0x0 x)) + (mul:^rh half (add:^rh ep em)) + :: + :: +tanh: Hyperbolic tangent + :: tanh(x) = sinh(x) / cosh(x) + :: + ++ tanh + |= x=@rh + ^- @rh + =/ e2x (exp (mul:^rh `@rh`0x4000 x)) :: exp(2x) + =/ one `@rh`0x3c00 + :: tanh(x) = (e^2x - 1) / (e^2x + 1) + (div:^rh (sub:^rh e2x one) (add:^rh e2x one)) + :: + :: +asinh: Inverse hyperbolic sine + :: asinh(x) = ln(x + sqrt(x^2 + 1)) + :: + ++ asinh + |= x=@rh + ^- @rh + =/ one `@rh`0x3c00 + (log (add:^rh x (sqt (add:^rh (mul:^rh x x) one)))) + :: + :: +acosh: Inverse hyperbolic cosine + :: acosh(x) = ln(x + sqrt(x^2 - 1)), x >= 1 + :: + ++ acosh + |= x=@rh + ^- @rh + =/ one `@rh`0x3c00 + ?: (lth:^rh x one) (sub:^rh x x) :: NaN for x < 1 + (log (add:^rh x (sqt (sub:^rh (mul:^rh x x) one)))) + :: + :: +atanh: Inverse hyperbolic tangent + :: atanh(x) = 0.5 * ln((1+x)/(1-x)), |x| < 1 + :: + ++ atanh + |= x=@rh + ^- @rh + =/ one `@rh`0x3c00 + =/ half `@rh`0x3800 + =/ ax `@rh`(dis `@`x 0x7fff) + ?: (gte:^rh ax one) (sub:^rh x x) :: NaN for |x| >= 1 + (mul:^rh half (log (div:^rh (add:^rh one x) (sub:^rh one x)))) + :: + :: ============================================================ + :: SPECIAL FUNCTIONS + :: ============================================================ + :: + :: +erf: Error function (Horner approximation) + :: Uses Abramowitz & Stegun approximation 7.1.26 + :: + ++ erf + |= x=@rh + ^- @rh + =/ bits `@`x + =/ sign (rsh [0 15] bits) + =/ ax `@rh`(dis bits 0x7fff) + =/ one `@rh`0x3c00 + :: Coefficients (simplified for half precision) + =/ a1 `@rh`0x3528 :: ~0.254829592 + =/ a2 `@rh`0xaccc :: ~-0.284496736 + =/ a3 `@rh`0x3b85 :: ~1.421413741 + =/ a4 `@rh`0xbbc8 :: ~-1.453152027 + =/ a5 `@rh`0x3a51 :: ~1.061405429 + =/ p `@rh`0x352b :: ~0.3275911 + :: + =/ t (div:^rh one (add:^rh one (mul:^rh p ax))) + =/ t2 (mul:^rh t t) + =/ t3 (mul:^rh t2 t) + =/ t4 (mul:^rh t3 t) + =/ t5 (mul:^rh t4 t) + :: + =/ poly %+ add:^rh (mul:^rh a1 t) + %+ add:^rh (mul:^rh a2 t2) + %+ add:^rh (mul:^rh a3 t3) + %+ add:^rh (mul:^rh a4 t4) + (mul:^rh a5 t5) + :: + =/ result (sub:^rh one (mul:^rh poly (exp (sub:^rh `@rh`0x0 (mul:^rh ax ax))))) + ?:(=(sign 0) result (sub:^rh `@rh`0x0 result)) + :: + :: +gamma: Gamma function (simplified Lanczos for half precision) + :: Γ(n) = (n-1)! for positive integers + :: + ++ gamma + |= x=@rh + ^- @rh + =/ one `@rh`0x3c00 + =/ half `@rh`0x3800 + =/ two `@rh`0x4000 + :: Handle special cases + ?: (lte:^rh x `@rh`0x0) (sub:^rh x x) :: NaN for x <= 0 + ?: =(x one) one :: Γ(1) = 1 + ?: =(x two) one :: Γ(2) = 1 + :: For x > 2, use Γ(x) = (x-1) * Γ(x-1) + :: For 0 < x < 1, use reflection or direct + :: Simplified: use exp-log for small half precision + :: Γ(x) ≈ sqrt(2π/x) * (x/e + 1/(12ex))^x (Stirling) + =/ pi `@rh`0x4248 :: π ≈ 3.14159 + =/ e `@rh`0x4170 :: e ≈ 2.718 + =/ twopi (mul:^rh two pi) + :: Stirling: Γ(x) ≈ sqrt(2π/x) * (x/e)^x + =/ base (div:^rh x e) + =/ term1 (sqt (div:^rh twopi x)) + =/ term2 (pow base x) + (mul:^rh term1 term2) + :: + :: +lgamma: Log-gamma function + :: + ++ lgamma + |= x=@rh + ^- @rh + (log (gamma x)) + :: + :: +j0: Bessel function of first kind, order 0 + :: Polynomial approximation for |x| <= 8 + :: + ++ j0 + |= x=@rh + ^- @rh + =/ ax `@rh`(dis `@`x 0x7fff) + =/ one `@rh`0x3c00 + :: For small x: J0(x) ≈ 1 - x²/4 + x⁴/64 + =/ x2 (mul:^rh ax ax) + =/ x4 (mul:^rh x2 x2) + =/ c1 `@rh`0xb800 :: -0.25 + =/ c2 `@rh`0x2c00 :: 0.015625 (1/64) + (add:^rh one (add:^rh (mul:^rh c1 x2) (mul:^rh c2 x4))) + :: + :: +j1: Bessel function of first kind, order 1 + :: Polynomial approximation for |x| <= 8 + :: + ++ j1 + |= x=@rh + ^- @rh + =/ bits `@`x + =/ sign (rsh [0 15] bits) + =/ ax `@rh`(dis bits 0x7fff) + :: For small x: J1(x) ≈ x/2 - x³/16 + x⁵/384 + =/ half `@rh`0x3800 + =/ x2 (mul:^rh ax ax) + =/ x3 (mul:^rh x2 ax) + =/ c1 `@rh`0xac00 :: -0.0625 (1/16) + =/ term1 (mul:^rh half ax) + =/ term2 (mul:^rh c1 x3) + =/ result (add:^rh term1 term2) + ?:(=(sign 0) result (sub:^rh `@rh`0x0 result)) + :: + :: +y0: Bessel function of second kind, order 0 + :: + ++ y0 + |= x=@rh + ^- @rh + ?: (lte:^rh x `@rh`0x0) (sub:^rh x x) :: NaN for x <= 0 + =/ twoopi `@rh`0x3a22 :: 2/π ≈ 0.6366 + =/ euler `@rh`0x394a :: γ ≈ 0.5772 + :: Y0(x) ≈ (2/π) * (ln(x/2) + γ) * J0(x) for small x + =/ half `@rh`0x3800 + =/ lnterm (add:^rh (log (mul:^rh half x)) euler) + (mul:^rh twoopi (mul:^rh lnterm (j0 x))) + :: + :: +y1: Bessel function of second kind, order 1 + :: + ++ y1 + |= x=@rh + ^- @rh + ?: (lte:^rh x `@rh`0x0) (sub:^rh x x) :: NaN for x <= 0 + =/ twoopi `@rh`0x3a22 :: 2/π + =/ one `@rh`0x3c00 + :: Y1(x) ≈ (2/π) * (J1(x)*ln(x) - 1/x) for small x + =/ j1x (j1 x) + =/ lnx (log x) + =/ invx (div:^rh one x) + (mul:^rh twoopi (sub:^rh (mul:^rh j1x lnx) invx)) + -- +:: +:: ================================================================ +:: SINGLE PRECISION (@rs) - 32-bit IEEE 754 +:: ================================================================ +:: +:: Coefficients from musl __sindf.c, __cosdf.c, expf.c +:: All values as hex literals for bit-exact matching with C jets. +:: +++ rs + |% + :: + :: ============================================================ + :: CONSTANTS - Must match exactly between Hoon and C jet + :: ============================================================ + :: + :: Pi-related constants for argument reduction + :: From musl __rem_pio2f.c + :: + ++ pio2-1 `@rs`0x3fc9.0fda :: 1.5707963267341256e+00 + ++ pio2-1t `@rs`0x33a2.2168 :: 1.5893254773528196e-08 + ++ pio4 `@rs`0x3f49.0fda :: pi/4 = 0.7853981633974483 + ++ invpio2 `@rs`0x3f22.f983 :: 2/pi = 0.6366197723675814 + :: + :: Exp constants + :: + ++ ln2hi `@rs`0x3f31.7200 :: 0.6931381225585938 + ++ ln2lo `@rs`0x3717.f7d1 :: 9.058001e-06 + ++ invln2 `@rs`0x3fb8.aa3b :: 1.4426950216293335 + ++ exp-huge `@rs`0x42b1.7218 :: 88.72283935546875 (overflow threshold) + ++ exp-tiny `@rs`0xc2cf.f1b5 :: -103.97208404541016 (underflow threshold) + :: + :: ============================================================ + :: POLYNOMIAL COEFFICIENTS - from musl (hex-exact) + :: ============================================================ + :: + :: Sin coefficients: sin(x) = x + x^3*S1 + x^5*S2 + x^7*S3 + x^9*S4 + :: From musl __sindf.c + :: + ++ sin-s1 `@rs`0xbe2a.aaab :: -1.6666666641831398e-01 + ++ sin-s2 `@rs`0x3c08.8889 :: 8.3333293930888176e-03 + ++ sin-s3 `@rs`0xb950.0d01 :: -1.9839334836483002e-04 + ++ sin-s4 `@rs`0x3638.ef1b :: 2.7181216275692873e-06 + :: + :: Cos coefficients: cos(x) = 1 - x^2/2 + x^4*C0 + x^6*C1 + x^8*C2 + x^10*C3 + :: From musl __cosdf.c + :: + ++ cos-c0 `@rs`0x3d2a.aaab :: 4.1666668653488159e-02 + ++ cos-c1 `@rs`0xbab6.0b61 :: -1.3888889225125313e-03 + ++ cos-c2 `@rs`0x37d0.0d01 :: 2.4801587642822415e-05 + ++ cos-c3 `@rs`0xb493.f27c :: -2.7557314603984356e-07 + :: + :: Exp polynomial coefficients + :: exp(r) ~ 1 + r + r^2*(P1 + r*(P2 + r*(P3 + r*P4))) + :: + ++ exp-p1 `@rs`0x3f00.0000 :: 0.5 + ++ exp-p2 `@rs`0x3e2a.aaab :: 0.16666667163372040 + ++ exp-p3 `@rs`0x3d2a.aaab :: 0.04166668653488159 + ++ exp-p4 `@rs`0x3c08.8889 :: 0.008333338864892721 + :: + :: ============================================================ + :: CHEBYSHEV POLYNOMIAL + :: ============================================================ + :: + :: +cheb: Chebyshev polynomial T_n(x) of the first kind + :: + :: The Chebyshev polynomials are defined by the recurrence: + :: T_n(x) = 2x*T_{n-1}(x) - T_{n-2}(x) + :: with T_0(x) = 1 and T_1(x) = x. + :: + :: Examples + :: > ((cheb 0) .0.25) + :: .1 + :: > ((cheb 1) .0.25) + :: .0.25 + :: > ((cheb 2) .0.25) + :: .-0.875 + :: > ((cheb 3) .0.25) + :: .-0.6875 + :: + ++ cheb + |= n=@ud + |= x=@rs + ^- @rs + ?: =(0 n) .1 + ?: =(1 n) x + (sub:^rs (mul:^rs .2 (mul:^rs x $(n (dec n)))) $(n (sub n 2))) + :: + :: ============================================================ + :: HELPER FUNCTIONS + :: ============================================================ + :: + :: +scalbn: multiply x by 2^n (ldexp equivalent) + :: + ++ scalbn + |= [x=@rs n=@sd] + ^- @rs + =/ bits `@`x + =/ exp (dis (rsh [0 23] bits) 0xff) + =/ sign (rsh [0 31] bits) + :: Special cases + ?: =(exp 0xff) x :: inf or nan + ?: =(exp 0) x :: zero or subnormal + :: Compute new exponent + =/ new-exp (sum:si (sun:si exp) n) + :: Overflow + ?: (gth:si new-exp (sun:si 254)) + ?:(=(sign 0) `@rs`0x7f80.0000 `@rs`0xff80.0000) + :: Underflow + ?: (lth:si new-exp (sun:si 1)) + .0 + :: Reconstruct + =/ new-bits + %+ con (lsh [0 31] sign) + %+ con (lsh [0 23] (abs:si new-exp)) + (dis bits 0x7f.ffff) + `@rs`new-bits + :: + :: +floor-int: floor to signed integer + :: + ++ floor-int + |= x=@rs + ^- @sd + =/ bits `@`x + =/ sign (rsh [0 31] bits) + =/ bexp (dis (rsh [0 23] bits) 0xff) + ?: (^lth bexp 127) :: |x| < 1 + ?:(=(sign 0) --0 -1) + =/ exp (sub bexp 127) + =/ mant (con (dis bits 0x7f.ffff) 0x80.0000) + ?: (^gte exp 23) :: No fractional part + ?: =(sign 0) + (sun:si (lsh [0 (sub exp 23)] mant)) + (new:si %.n (lsh [0 (sub exp 23)] mant)) + =/ shift (sub 23 exp) + =/ int-mant (rsh [0 shift] mant) + =/ frac-mask (dec (lsh [0 shift] 1)) + =/ had-frac !=(0 (dis mant frac-mask)) + =/ result ?:(=(sign 0) (sun:si int-mant) (new:si %.n int-mant)) + ?: &(=(sign 1) had-frac) + (dif:si result --1) + result + :: + :: ============================================================ + :: ARGUMENT REDUCTION + :: ============================================================ + :: + :: +rem-pio2: Reduce x to [-pi/4, pi/4], return quadrant + :: + :: Returns [n y] where x = n*(pi/2) + y, |y| <= pi/4 + :: + ++ rem-pio2 + |= x=@rs + ^- [n=@ y=@rs] + =/ bits `@`x + =/ sign (rsh [0 31] bits) + =/ ax `@rs`(dis bits 0x7fff.ffff) + :: Small argument: no reduction needed + ?: (lte:^rs ax pio4) + [0 x] + :: Compute n = round(x * 2/pi) + =/ fn-raw (mul:^rs x invpio2) + =/ fn ?: (gte:^rs fn-raw .0) + (add:^rs fn-raw .0.5) + (sub:^rs fn-raw .0.5) + =/ n (floor-int fn) + =/ fn (san:^rs n) + :: y = x - n*pio2_1 - n*pio2_1t (extended precision) + =/ y (sub:^rs x (mul:^rs fn pio2-1)) + =. y (sub:^rs y (mul:^rs fn pio2-1t)) + [(mod (abs:si n) 4) y] + :: + :: ============================================================ + :: KERNEL FUNCTIONS - Fixed polynomial evaluation + :: ============================================================ + :: + :: +sindf: Kernel sin for reduced argument |x| <= pi/4 + :: + :: sin(x) = x + x^3*S1 + x^5*S2 + x^7*S3 + x^9*S4 + :: = x + x*z*(S1 + z*(S2 + z*(S3 + z*S4))) + :: where z = x^2 + :: + ++ sindf + |= x=@rs + ^- @rs + =/ z (mul:^rs x x) :: z = x^2 + :: r = S3 + z*S4 + =/ r (add:^rs sin-s3 (mul:^rs z sin-s4)) + :: r = S2 + z*r + =. r (add:^rs sin-s2 (mul:^rs z r)) + :: r = S1 + z*r + =. r (add:^rs sin-s1 (mul:^rs z r)) + :: result = x + x*z*r + (add:^rs x (mul:^rs (mul:^rs x z) r)) + :: + :: +cosdf: Kernel cos for reduced argument |x| <= pi/4 + :: + :: cos(x) = 1 - z/2 + z^2*(C0 + z*(C1 + z*(C2 + z*C3))) + :: where z = x^2 + :: + ++ cosdf + |= x=@rs + ^- @rs + =/ z (mul:^rs x x) :: z = x^2 + =/ w (mul:^rs z z) :: w = z^2 + :: r = C2 + z*C3 + =/ r (add:^rs cos-c2 (mul:^rs z cos-c3)) + :: r = C1 + z*r + =. r (add:^rs cos-c1 (mul:^rs z r)) + :: r = C0 + z*r + =. r (add:^rs cos-c0 (mul:^rs z r)) + :: result = 1 - z/2 + w*r + =/ hz (mul:^rs .0.5 z) + (add:^rs (sub:^rs .1 hz) (mul:^rs w r)) + :: + :: ============================================================ + :: PUBLIC API + :: ============================================================ + :: + :: +sin: Sine + :: + ++ sin + |= x=@rs + ^- @rs + =/ bits `@`x + =/ ax `@rs`(dis bits 0x7fff.ffff) + :: Tiny: sin(x) ~ x + ?: (lth:^rs ax `@rs`0x3980.0000) x :: |x| < 2^-12 + :: Inf/NaN + ?: (gte:^rs ax `@rs`0x7f80.0000) (sub:^rs x x) + :: Argument reduction + =/ [n=@ y=@rs] (rem-pio2 x) + :: Quadrant dispatch + ?- n + %0 (sindf y) + %1 (cosdf y) + %2 (sub:^rs .0 (sindf y)) + %3 (sub:^rs .0 (cosdf y)) + == + :: + :: +cos: Cosine + :: + ++ cos + |= x=@rs + ^- @rs + =/ bits `@`x + =/ ax `@rs`(dis bits 0x7fff.ffff) + :: Tiny: cos(x) ~ 1 + ?: (lth:^rs ax `@rs`0x3980.0000) .1 + :: Inf/NaN + ?: (gte:^rs ax `@rs`0x7f80.0000) (sub:^rs x x) + :: Argument reduction + =/ [n=@ y=@rs] (rem-pio2 x) + :: Quadrant dispatch + ?- n + %0 (cosdf y) + %1 (sub:^rs .0 (sindf y)) + %2 (sub:^rs .0 (cosdf y)) + %3 (sindf y) + == + :: + :: +tan: Tangent + :: + ++ tan + |= x=@rs + ^- @rs + (div:^rs (sin x) (cos x)) + :: + :: +exp: Exponential + :: + :: Algorithm: + :: 1. Handle special cases + :: 2. Reduce: x = k*ln2 + r, |r| <= ln2/2 + :: 3. Compute exp(r) via polynomial + :: 4. Return 2^k * exp(r) + :: + ++ exp + |= x=@rs + ^- @rs + =/ bits `@`x + =/ ax `@rs`(dis bits 0x7fff.ffff) + :: x = 0 -> 1 + ?: =(x .0) .1 + :: +inf -> +inf + ?: =(bits 0x7f80.0000) `@rs`0x7f80.0000 + :: -inf -> 0 + ?: =(bits 0xff80.0000) .0 + :: NaN -> NaN + ?: (gth:^rs ax `@rs`0x7f80.0000) x + :: Overflow + ?: (gth:^rs x exp-huge) `@rs`0x7f80.0000 + :: Underflow + ?: (lth:^rs x exp-tiny) .0 + :: + :: Argument reduction: x = k*ln2 + r + :: + =/ fn-raw (mul:^rs x invln2) + =/ fn ?: (gte:^rs fn-raw .0) + (add:^rs fn-raw .0.5) + (sub:^rs fn-raw .0.5) + =/ k (floor-int fn) + =/ kf (san:^rs k) + :: r = x - k*ln2hi - k*ln2lo (extended precision) + =/ r (sub:^rs x (mul:^rs kf ln2hi)) + =. r (sub:^rs r (mul:^rs kf ln2lo)) + :: + :: Polynomial: exp(r) ~ 1 + r + r^2*(P1 + r*(P2 + r*(P3 + r*P4))) + :: + =/ r2 (mul:^rs r r) + :: p = P3 + r*P4 + =/ p (add:^rs exp-p3 (mul:^rs r exp-p4)) + :: p = P2 + r*p + =. p (add:^rs exp-p2 (mul:^rs r p)) + :: p = P1 + r*p + =. p (add:^rs exp-p1 (mul:^rs r p)) + :: exp(r) = 1 + r + r^2*p + =/ expr (add:^rs .1 (add:^rs r (mul:^rs r2 p))) + :: + :: Scale by 2^k + :: + (scalbn expr k) + :: + :: +log: Natural logarithm + :: + :: Algorithm (from FreeBSD e_logf.c): + :: 1. Handle special cases (0, negative, inf, nan) + :: 2. Reduce: x = 2^k * m where 1 <= m < 2 + :: 3. Further reduce: s = (m-1)/(m+1), so m = (1+s)/(1-s) + :: 4. Compute log(m) = 2*s + 2*s^3*(Lg1 + s^2*(Lg2 + s^2*(Lg3 + s^2*Lg4))) + :: 5. Return k*ln2 + log(m) + :: + ++ log + |= x=@rs + ^- @rs + =/ bits `@`x + :: Special cases + ?: =(x .0) `@rs`0xff80.0000 :: log(0) = -inf + ?: =(bits 0x7f80.0000) `@rs`0x7f80.0000 :: log(+inf) = +inf + ?: (gth:^rs `@rs`(dis bits 0x7fff.ffff) `@rs`0x7f80.0000) x :: nan + ?: !=(0 (rsh [0 31] bits)) (sub:^rs x x) :: log(negative) = nan + :: + :: Extract exponent and mantissa + :: x = 2^k * m where 1 <= m < 2 + :: + =/ exp-bits (dis (rsh [0 23] bits) 0xff) + =/ mant-bits (dis bits 0x7f.ffff) + :: Handle subnormals (simplified: treat as very small) + ?: =(exp-bits 0) `@rs`0xff80.0000 + :: k = exponent - bias + =/ k (sub exp-bits 127) + :: m = 1.mantissa (reconstruct normalized mantissa) + =/ m `@rs`(con 0x3f80.0000 mant-bits) + :: + :: Compute s = (m-1)/(m+1) + :: + =/ f (sub:^rs m .1) + =/ s (div:^rs f (add:^rs .2 f)) + =/ s2 (mul:^rs s s) + =/ s4 (mul:^rs s2 s2) + :: + :: Polynomial: R = Lg2 + s2*(Lg4) [even terms] + :: + Lg1 + s2*(Lg3) [odd terms] + :: Coefficients from FreeBSD e_logf.c + :: + =/ lg1 `@rs`0x3f2a.aaaa :: 0.66666662693 + =/ lg2 `@rs`0x3ecc.ce13 :: 0.40000972152 + =/ lg3 `@rs`0x3e91.e9ee :: 0.28498786688 + =/ lg4 `@rs`0x3e78.9e26 :: 0.24279078841 + :: + =/ t1 (mul:^rs s4 (add:^rs lg2 (mul:^rs s4 lg4))) + =/ t2 (mul:^rs s2 (add:^rs lg1 (mul:^rs s4 lg3))) + =/ r (add:^rs t1 t2) + :: + :: log(m) = 2*s + 2*s^3*R = 2*s*(1 + s^2*R) + :: + =/ log-m (mul:^rs (mul:^rs .2 s) (add:^rs .1 (mul:^rs s2 r))) + :: + :: Result = k*ln2 + log(m) + :: + =/ kf (san:^rs (sun:si k)) + (add:^rs (mul:^rs kf ln2hi) log-m) + :: + :: +log-2: Base-2 logarithm + :: + ++ log-2 + |= x=@rs + ^- @rs + (mul:^rs (log x) invln2) + :: + :: +log-10: Base-10 logarithm + :: + ++ log-10 + |= x=@rs + ^- @rs + =/ invlog10 `@rs`0x3ede.5bd9 :: 1/log(10) = 0.43429448 + (mul:^rs (log x) invlog10) + :: + :: +sqt: Square root (Newton-Raphson with fixed iterations) + :: + :: Algorithm: + :: 1. Handle special cases + :: 2. Initial guess from bit manipulation + :: 3. Fixed 4 Newton-Raphson iterations: g = (g + x/g) / 2 + :: + ++ sqt + |= x=@rs + ^- @rs + =/ bits `@`x + :: Special cases + ?: =(x .0) .0 + ?: =(bits 0x7f80.0000) `@rs`0x7f80.0000 :: sqrt(+inf) = +inf + ?: (gth:^rs `@rs`(dis bits 0x7fff.ffff) `@rs`0x7f80.0000) x :: nan + ?: !=(0 (rsh [0 31] bits)) (sub:^rs x x) :: sqrt(negative) = nan + :: + :: Initial guess: use bit manipulation for good starting point + :: g0 = bits/2 + 0x1fc00000 (magic constant for float sqrt) + :: + =/ g0 `@rs`(add (rsh [0 1] bits) 0x1fc0.0000) + :: + :: Newton-Raphson: g = (g + x/g) / 2 + :: Fixed 4 iterations for single precision + :: + =/ g g0 + =. g (mul:^rs .0.5 (add:^rs g (div:^rs x g))) :: iter 1 + =. g (mul:^rs .0.5 (add:^rs g (div:^rs x g))) :: iter 2 + =. g (mul:^rs .0.5 (add:^rs g (div:^rs x g))) :: iter 3 + =. g (mul:^rs .0.5 (add:^rs g (div:^rs x g))) :: iter 4 + g + :: + ++ sqrt sqt + :: + :: +atan: Inverse tangent + :: + :: Algorithm (from FreeBSD s_atanf.c): + :: 1. Reduce to [0, 0.4375] using atan(x) = pi/2 - atan(1/x) for |x| > 1 + :: and atan(x) = atan(c) + atan((x-c)/(1+x*c)) for thresholds + :: 2. Polynomial approximation in reduced range + :: + ++ atan + |= x=@rs + ^- @rs + =/ bits `@`x + =/ sign (rsh [0 31] bits) + =/ ax `@rs`(dis bits 0x7fff.ffff) + :: + :: Special cases + ?: (gte:^rs ax `@rs`0x7f80.0000) :: inf or nan + ?: =(bits 0x7f80.0000) :: +inf + `@rs`0x3fc9.0fdb :: pi/2 + ?: =(bits 0xff80.0000) :: -inf + `@rs`0xbfc9.0fdb :: -pi/2 + x :: nan + :: + :: Polynomial coefficients from FreeBSD s_atanf.c + :: + =/ at0 `@rs`0x3eaa.aaaa :: 0.33333334 + =/ at1 `@rs`0xbe4c.cccc :: -0.20000000 + =/ at2 `@rs`0x3e12.4925 :: 0.14285715 + =/ at3 `@rs`0xbdc3.b8a4 :: -0.11111089 + :: + :: Reduction thresholds + :: + =/ lo-thresh `@rs`0x3ee0.0000 :: 0.4375 + =/ hi-thresh `@rs`0x401c.0000 :: 2.4375 + :: + :: atan constants for reduction + :: + =/ pio2 `@rs`0x3fc9.0fdb :: pi/2 + =/ pio4 `@rs`0x3f49.0fdb :: pi/4 + :: + :: Core polynomial evaluation for |x| <= 0.4375 + :: + =/ eval-poly + |= z=@rs + ^- @rs + =/ w (mul:^rs z z) + =/ s1 (add:^rs at2 (mul:^rs w at3)) + =/ s2 (add:^rs at0 (mul:^rs w at1)) + =/ r (add:^rs s2 (mul:^rs (mul:^rs w w) s1)) + (sub:^rs z (mul:^rs z (mul:^rs w r))) + :: + =/ result=@rs + ?: (lte:^rs ax lo-thresh) + :: Small: use polynomial directly + (eval-poly ax) + ?: (lte:^rs ax .1) + :: Medium-small: atan(x) = atan(0.5) + atan((x-0.5)/(1+x*0.5)) + =/ atan-half `@rs`0x3eed.6338 :: atan(0.5) = 0.46364760 + =/ t (div:^rs (sub:^rs ax .0.5) (add:^rs .1 (mul:^rs ax .0.5))) + (add:^rs atan-half (eval-poly t)) + ?: (lte:^rs ax hi-thresh) + :: Medium: atan(x) = pi/4 + atan((x-1)/(1+x)) + =/ t (div:^rs (sub:^rs ax .1) (add:^rs .1 ax)) + (add:^rs pio4 (eval-poly t)) + :: Large: atan(x) = pi/2 - atan(1/x) + =/ t (div:^rs .1 ax) + (sub:^rs pio2 (eval-poly t)) + :: + :: Apply sign + ?:(=(sign 0) result (sub:^rs .0 result)) + :: + :: +asin: Inverse sine + :: + ++ asin + |= x=@rs + ^- @rs + =/ bits `@`x + =/ ax `@rs`(dis bits 0x7fff.ffff) + ?: (gth:^rs ax .1) (sub:^rs x x) :: |x| > 1: nan + ?: =(ax .1) :: |x| = 1: +/- pi/2 + ?: =(0 (rsh [0 31] bits)) + `@rs`0x3fc9.0fdb :: pi/2 + `@rs`0xbfc9.0fdb :: -pi/2 + :: asin(x) = atan(x / sqrt(1 - x^2)) + (atan (div:^rs x (sqt (sub:^rs .1 (mul:^rs x x))))) + :: + :: +acos: Inverse cosine + :: + ++ acos + |= x=@rs + ^- @rs + =/ bits `@`x + =/ ax `@rs`(dis bits 0x7fff.ffff) + =/ pi `@rs`0x4049.0fdb :: pi + ?: (gth:^rs ax .1) (sub:^rs x x) :: |x| > 1: nan + ?: =(x .1) .0 :: acos(1) = 0 + ?: =(x .-1) pi :: acos(-1) = pi + :: acos(x) = atan(sqrt(1 - x^2) / x) with quadrant adjustment + =/ s (sqt (sub:^rs .1 (mul:^rs x x))) + ?: (gte:^rs x .0) + (atan (div:^rs s x)) + (add:^rs pi (atan (div:^rs s x))) + :: + :: +atan2: Two-argument inverse tangent + :: + ++ atan2 + |= [y=@rs x=@rs] + ^- @rs + =/ pi `@rs`0x4049.0fdb :: pi + =/ pio2 `@rs`0x3fc9.0fdb :: pi/2 + ?: (gth:^rs x .0) + (atan (div:^rs y x)) + ?: &((lth:^rs x .0) (gte:^rs y .0)) + (add:^rs (atan (div:^rs y x)) pi) + ?: &((lth:^rs x .0) (lth:^rs y .0)) + (sub:^rs (atan (div:^rs y x)) pi) + ?: &(=(x .0) (gth:^rs y .0)) + pio2 + ?: &(=(x .0) (lth:^rs y .0)) + (sub:^rs .0 pio2) + .0 :: undefined (0,0) + :: + :: +pow-n: Integer power (fixed iteration) + :: + ++ pow-n + |= [x=@rs n=@rs] + ^- @rs + ?: =(n .0) .1 + :: Convert n to integer (must be positive integer) + =/ ni (abs:si (need (toi:^rs n))) + =/ neg (lth:^rs n .0) + :: Binary exponentiation (fixed max 32 iterations for 32-bit range) + =/ result .1 + =/ base x + =/ i 0 + |- + ?: |(=(ni 0) (gth i 31)) + ?:(neg (div:^rs .1 result) result) + =/ new-result + ?: =(1 (dis ni 1)) + (mul:^rs result base) + result + $(ni (rsh [0 1] ni), base (mul:^rs base base), result new-result, i +(i)) + :: + :: +pow: General power function + :: + ++ pow + |= [x=@rs n=@rs] + ^- @rs + :: Special case: integer exponent + =/ ni (toi:^rs n) + ?: &(?=(^ ni) =(n (san:^rs (need ni)))) + (pow-n x n) + :: General case: x^n = exp(n * log(x)) + (exp (mul:^rs n (log x))) + :: + :: ============================================================ + :: HYPERBOLIC FUNCTIONS + :: ============================================================ + :: + :: +sinh: Hyperbolic sine + :: sinh(x) = (exp(x) - exp(-x)) / 2 + :: + ++ sinh + |= x=@rs + ^- @rs + =/ ep (exp x) + =/ em (exp (sub:^rs .0 x)) + (mul:^rs .0.5 (sub:^rs ep em)) + :: + :: +cosh: Hyperbolic cosine + :: cosh(x) = (exp(x) + exp(-x)) / 2 + :: + ++ cosh + |= x=@rs + ^- @rs + =/ ep (exp x) + =/ em (exp (sub:^rs .0 x)) + (mul:^rs .0.5 (add:^rs ep em)) + :: + :: +tanh: Hyperbolic tangent + :: tanh(x) = (e^2x - 1) / (e^2x + 1) + :: + ++ tanh + |= x=@rs + ^- @rs + =/ e2x (exp (mul:^rs .2 x)) + (div:^rs (sub:^rs e2x .1) (add:^rs e2x .1)) + :: + :: +asinh: Inverse hyperbolic sine + :: asinh(x) = ln(x + sqrt(x^2 + 1)) + :: + ++ asinh + |= x=@rs + ^- @rs + (log (add:^rs x (sqt (add:^rs (mul:^rs x x) .1)))) + :: + :: +acosh: Inverse hyperbolic cosine + :: acosh(x) = ln(x + sqrt(x^2 - 1)), x >= 1 + :: + ++ acosh + |= x=@rs + ^- @rs + ?: (lth:^rs x .1) (sub:^rs x x) :: NaN for x < 1 + (log (add:^rs x (sqt (sub:^rs (mul:^rs x x) .1)))) + :: + :: +atanh: Inverse hyperbolic tangent + :: atanh(x) = 0.5 * ln((1+x)/(1-x)), |x| < 1 + :: + ++ atanh + |= x=@rs + ^- @rs + =/ ax (abs x) + ?: (gte:^rs ax .1) (sub:^rs x x) :: NaN for |x| >= 1 + (mul:^rs .0.5 (log (div:^rs (add:^rs .1 x) (sub:^rs .1 x)))) + :: + :: ============================================================ + :: SPECIAL FUNCTIONS + :: ============================================================ + :: + :: +erf: Error function + :: Uses Abramowitz & Stegun approximation 7.1.26 (max error ~1.5e-7) + :: + ++ erf + |= x=@rs + ^- @rs + =/ bits `@`x + =/ sign (rsh [0 31] bits) + =/ ax (abs x) + :: Coefficients from A&S + =/ a1 .0.254829592 + =/ a2 .-0.284496736 + =/ a3 .1.421413741 + =/ a4 .-1.453152027 + =/ a5 .1.061405429 + =/ p .0.3275911 + :: + =/ t (div:^rs .1 (add:^rs .1 (mul:^rs p ax))) + =/ t2 (mul:^rs t t) + =/ t3 (mul:^rs t2 t) + =/ t4 (mul:^rs t3 t) + =/ t5 (mul:^rs t4 t) + :: + =/ poly %+ add:^rs (mul:^rs a1 t) + %+ add:^rs (mul:^rs a2 t2) + %+ add:^rs (mul:^rs a3 t3) + %+ add:^rs (mul:^rs a4 t4) + (mul:^rs a5 t5) + :: + =/ result (sub:^rs .1 (mul:^rs poly (exp (sub:^rs .0 (mul:^rs ax ax))))) + ?:(=(sign 0) result (sub:^rs .0 result)) + :: + :: +gamma: Gamma function using Lanczos approximation + :: Γ(z) = sqrt(2π) * (z + g + 0.5)^(z+0.5) * e^-(z+g+0.5) * A(z) + :: + ++ gamma + |= x=@rs + ^- @rs + :: Handle special cases + ?: (lte:^rs x .0) (sub:^rs x x) :: NaN for x <= 0 + ?: =(x .1) .1 :: Γ(1) = 1 + ?: =(x .2) .1 :: Γ(2) = 1 + :: Lanczos coefficients (g=5, n=7) + =/ g .5 + =/ c0 .1.000000000190015 + =/ c1 .76.18009172947146 + =/ c2 .-86.50532032941677 + =/ c3 .24.01409824083091 + =/ c4 .-1.231739572450155 + =/ c5 .0.001208650973866179 + =/ c6 .-0.000005395239384953 + :: + =/ xm1 (sub:^rs x .1) + =/ sum c0 + =. sum (add:^rs sum (div:^rs c1 (add:^rs xm1 .1))) + =. sum (add:^rs sum (div:^rs c2 (add:^rs xm1 .2))) + =. sum (add:^rs sum (div:^rs c3 (add:^rs xm1 .3))) + =. sum (add:^rs sum (div:^rs c4 (add:^rs xm1 .4))) + =. sum (add:^rs sum (div:^rs c5 (add:^rs xm1 .5))) + =. sum (add:^rs sum (div:^rs c6 (add:^rs xm1 .6))) + :: + =/ sqrt2pi .2.5066282746310002 + =/ t (add:^rs xm1 (add:^rs g .0.5)) + =/ term1 (mul:^rs sqrt2pi sum) + =/ term2 (pow t (add:^rs xm1 .0.5)) + =/ term3 (exp (sub:^rs .0 t)) + (mul:^rs term1 (mul:^rs term2 term3)) + :: + :: +lgamma: Log-gamma function + :: + ++ lgamma + |= x=@rs + ^- @rs + (log (gamma x)) + :: + :: +j0: Bessel function of first kind, order 0 + :: Uses polynomial approximation from Abramowitz & Stegun + :: + ++ j0 + |= x=@rs + ^- @rs + =/ ax (abs x) + ?: (lth:^rs ax .8) + :: Polynomial for |x| < 8 + =/ x2 (mul:^rs ax ax) + =/ x4 (mul:^rs x2 x2) + =/ x6 (mul:^rs x4 x2) + =/ x8 (mul:^rs x4 x4) + :: J0(x) ≈ 1 - x²/4 + x⁴/64 - x⁶/2304 + x⁸/147456 + =/ c1 .-0.25 + =/ c2 .0.015625 + =/ c3 .-0.000434028 + =/ c4 .0.0000067817 + %+ add:^rs .1 + %+ add:^rs (mul:^rs c1 x2) + %+ add:^rs (mul:^rs c2 x4) + %+ add:^rs (mul:^rs c3 x6) + (mul:^rs c4 x8) + :: Asymptotic for |x| >= 8 + =/ z (div:^rs .8 ax) + =/ z2 (mul:^rs z z) + =/ p0 .1 + =/ q0 .-0.125 + =/ p (add:^rs p0 (mul:^rs z2 .-0.00021)) + =/ q (add:^rs q0 (mul:^rs z2 .0.00066)) + =/ xx (sub:^rs ax .0.785398163) + =/ amp (div:^rs .0.79788456 (sqt ax)) + (mul:^rs amp (sub:^rs (mul:^rs p (cos xx)) (mul:^rs (mul:^rs z q) (sin xx)))) + :: + :: +j1: Bessel function of first kind, order 1 + :: + ++ j1 + |= x=@rs + ^- @rs + =/ bits `@`x + =/ sign (rsh [0 31] bits) + =/ ax (abs x) + =/ result=@rs + ?: (lth:^rs ax .8) + :: Polynomial for |x| < 8 + =/ x2 (mul:^rs ax ax) + =/ x3 (mul:^rs x2 ax) + =/ x5 (mul:^rs x3 x2) + =/ x7 (mul:^rs x5 x2) + :: J1(x) ≈ x/2 - x³/16 + x⁵/384 - x⁷/18432 + =/ c0 .0.5 + =/ c1 .-0.0625 + =/ c2 .0.00260417 + =/ c3 .-0.0000543 + %+ add:^rs (mul:^rs c0 ax) + %+ add:^rs (mul:^rs c1 x3) + %+ add:^rs (mul:^rs c2 x5) + (mul:^rs c3 x7) + :: Asymptotic for |x| >= 8 + =/ z (div:^rs .8 ax) + =/ z2 (mul:^rs z z) + =/ p1 .1 + =/ q1 .0.375 + =/ p (add:^rs p1 (mul:^rs z2 .0.00084)) + =/ q (add:^rs q1 (mul:^rs z2 .-0.00195)) + =/ xx (sub:^rs ax .2.35619449) + =/ amp (div:^rs .0.79788456 (sqt ax)) + (mul:^rs amp (sub:^rs (mul:^rs p (cos xx)) (mul:^rs (mul:^rs z q) (sin xx)))) + ?:(=(sign 0) result (sub:^rs .0 result)) + :: + :: +y0: Bessel function of second kind, order 0 + :: + ++ y0 + |= x=@rs + ^- @rs + ?: (lte:^rs x .0) (sub:^rs x x) :: NaN for x <= 0 + =/ twoopi .0.6366197723675814 + =/ euler .0.5772156649015329 + ?: (lth:^rs x .8) + :: Y0(x) ≈ (2/π) * (J0(x) * ln(x/2) + (x²/4 - ...)) + =/ j0x (j0 x) + =/ lnhx (log (mul:^rs .0.5 x)) + (add:^rs (mul:^rs twoopi (mul:^rs j0x lnhx)) (mul:^rs .0.36746691 (mul:^rs x x))) + :: Asymptotic for x >= 8 + =/ z (div:^rs .8 x) + =/ z2 (mul:^rs z z) + =/ p0 .1 + =/ q0 .-0.125 + =/ p (add:^rs p0 (mul:^rs z2 .-0.00021)) + =/ q (add:^rs q0 (mul:^rs z2 .0.00066)) + =/ xx (sub:^rs x .0.785398163) + =/ amp (div:^rs .0.79788456 (sqt x)) + (mul:^rs amp (add:^rs (mul:^rs p (sin xx)) (mul:^rs (mul:^rs z q) (cos xx)))) + :: + :: +y1: Bessel function of second kind, order 1 + :: + ++ y1 + |= x=@rs + ^- @rs + ?: (lte:^rs x .0) (sub:^rs x x) :: NaN for x <= 0 + =/ twoopi .0.6366197723675814 + ?: (lth:^rs x .8) + =/ j1x (j1 x) + =/ lnx (log x) + =/ invx (div:^rs .1 x) + (mul:^rs twoopi (sub:^rs (mul:^rs j1x lnx) invx)) + :: Asymptotic for x >= 8 + =/ z (div:^rs .8 x) + =/ z2 (mul:^rs z z) + =/ p1 .1 + =/ q1 .0.375 + =/ p (add:^rs p1 (mul:^rs z2 .0.00084)) + =/ q (add:^rs q1 (mul:^rs z2 .-0.00195)) + =/ xx (sub:^rs x .2.35619449) + =/ amp (div:^rs .0.79788456 (sqt x)) + (mul:^rs amp (add:^rs (mul:^rs p (sin xx)) (mul:^rs (mul:^rs z q) (cos xx)))) + -- +:: +:: ================================================================ +:: DOUBLE PRECISION (@rd) - 64-bit IEEE 754 +:: ================================================================ +:: +:: Coefficients from musl __sin.c, __cos.c, exp.c +:: Using higher-degree polynomials for double precision accuracy. +:: +++ rd + |% + :: + :: ============================================================ + :: CONSTANTS (hex-exact, from musl) + :: ============================================================ + :: + :: Pi-related constants for argument reduction + :: + ++ pio2-1 `@rd`0x3ff9.21fb.5440.0000 :: pi/2 high bits + ++ pio2-1t `@rd`0x3dd0.b461.1a62.6331 :: pi/2 - pio2_1 + ++ pio2-2 `@rd`0x3dd0.b461.1a60.0000 :: second 33 bits of pi/2 + ++ pio2-2t `@rd`0x3ba3.198a.2e03.7073 :: pi/2 - pio2_1 - pio2_2 + ++ pio4 `@rd`0x3fe9.21fb.5444.2d18 :: pi/4 + ++ invpio2 `@rd`0x3fe4.5f30.6dc9.c883 :: 2/pi + :: + :: Exp constants + :: + ++ ln2hi `@rd`0x3fe6.2e42.fee0.0000 :: ln(2) high + ++ ln2lo `@rd`0x3dea.39ef.3579.3c76 :: ln(2) low + ++ invln2 `@rd`0x3ff7.1547.652b.82fe :: 1/ln(2) + ++ exp-huge `@rd`0x4086.2e42.fefa.39ef :: 709.782... (overflow) + ++ exp-tiny `@rd`0xc087.4910.d52d.3051 :: -745.133... (underflow) + :: + :: ============================================================ + :: POLYNOMIAL COEFFICIENTS (hex-exact, from musl) + :: ============================================================ + :: + :: Sin coefficients: degree 13 polynomial + :: sin(x) = x + x^3*(S1 + x^2*(S2 + x^2*(S3 + x^2*(S4 + x^2*(S5 + x^2*S6))))) + :: + ++ sin-s1 `@rd`0xbfc5.5555.5555.5549 :: -1.66666666666666324348e-01 + ++ sin-s2 `@rd`0x3f81.1111.1110.f8a6 :: 8.33333333332248946124e-03 + ++ sin-s3 `@rd`0xbf2a.01a0.19c1.61d5 :: -1.98412698298579493134e-04 + ++ sin-s4 `@rd`0x3ec7.1de3.57b1.fe7d :: 2.75573137070700676789e-06 + ++ sin-s5 `@rd`0xbe5a.e5e6.8a2b.9ceb :: -2.50507602534068634195e-08 + ++ sin-s6 `@rd`0x3de5.d93a.5acf.d57c :: 1.58969099521155010221e-10 + :: + :: Cos coefficients: degree 14 polynomial + :: cos(x) = 1 - x^2/2 + x^4*(C1 + x^2*(C2 + x^2*(C3 + x^2*(C4 + x^2*(C5 + x^2*C6))))) + :: + ++ cos-c1 `@rd`0x3fa5.5555.5555.554c :: 4.16666666666666019037e-02 + ++ cos-c2 `@rd`0xbf56.c16c.16c1.5177 :: -1.38888888888741095749e-03 + ++ cos-c3 `@rd`0x3efa.01a0.19cb.1590 :: 2.48015872894767294178e-05 + ++ cos-c4 `@rd`0xbe92.7e4f.809c.52ad :: -2.75573143513906633035e-07 + ++ cos-c5 `@rd`0x3e21.ee9e.bdb4.b1c4 :: 2.08757232129817482790e-09 + ++ cos-c6 `@rd`0xbda8.fae9.be88.38d4 :: -1.13596475577881948265e-11 + :: + :: Exp polynomial coefficients + :: From musl exp_data.c + :: + ++ exp-p1 `@rd`0x3fe0.0000.0000.0000 :: 0.5 + ++ exp-p2 `@rd`0x3fc5.5555.5555.5555 :: 0.16666666666666666 + ++ exp-p3 `@rd`0x3fa5.5555.5555.5555 :: 0.041666666666666664 + ++ exp-p4 `@rd`0x3f81.1111.1111.1111 :: 0.008333333333333333 + ++ exp-p5 `@rd`0x3f56.c16c.16c1.6c17 :: 0.001388888888888889 + :: + :: ============================================================ + :: CHEBYSHEV POLYNOMIAL + :: ============================================================ + :: + :: +cheb: Chebyshev polynomial T_n(x) of the first kind + :: + ++ cheb + |= n=@ud + |= x=@rd + ^- @rd + ?: =(0 n) .~1 + ?: =(1 n) x + (sub:^rd (mul:^rd .~2 (mul:^rd x $(n (dec n)))) $(n (sub n 2))) + :: + :: ============================================================ + :: HELPER FUNCTIONS + :: ============================================================ + :: + :: +scalbn: multiply x by 2^n + :: + ++ scalbn + |= [x=@rd n=@sd] + ^- @rd + =/ bits `@`x + =/ exp (dis (rsh [0 52] bits) 0x7ff) + =/ sign (rsh [0 63] bits) + ?: =(exp 0x7ff) x :: inf or nan + ?: =(exp 0) x :: zero or subnormal + =/ new-exp (sum:si (sun:si exp) n) + ?: (gth:si new-exp (sun:si 2.046)) + ?:(=(sign 0) `@rd`0x7ff0.0000.0000.0000 `@rd`0xfff0.0000.0000.0000) + ?: (lth:si new-exp (sun:si 1)) + .~0 + =/ new-bits + %+ con (lsh [0 63] sign) + %+ con (lsh [0 52] (abs:si new-exp)) + (dis bits 0xf.ffff.ffff.ffff) + `@rd`new-bits + :: + :: +floor-int: floor to signed integer + :: + ++ floor-int + |= x=@rd + ^- @sd + =/ bits `@`x + =/ sign (rsh [0 63] bits) + =/ bexp (dis (rsh [0 52] bits) 0x7ff) + ?: (^lth bexp 1.023) + ?:(=(sign 0) --0 -1) + =/ exp (sub bexp 1.023) + =/ mant (con (dis bits 0xf.ffff.ffff.ffff) 0x10.0000.0000.0000) + ?: (^gte exp 52) + ?: =(sign 0) + (sun:si (lsh [0 (sub exp 52)] mant)) + (new:si %.n (lsh [0 (sub exp 52)] mant)) + =/ shift (sub 52 exp) + =/ int-mant (rsh [0 shift] mant) + =/ frac-mask (dec (lsh [0 shift] 1)) + =/ had-frac !=(0 (dis mant frac-mask)) + =/ result ?:(=(sign 0) (sun:si int-mant) (new:si %.n int-mant)) + ?: &(=(sign 1) had-frac) + (dif:si result --1) + result + :: + :: ============================================================ + :: ARGUMENT REDUCTION + :: ============================================================ + :: + :: +rem-pio2: Reduce x to [-pi/4, pi/4] + :: + ++ rem-pio2 + |= x=@rd + ^- [n=@ y=@rd] + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + ?: (lte:^rd ax pio4) + [0 x] + =/ fn-raw (mul:^rd x invpio2) + =/ fn ?: (gte:^rd fn-raw .~0) + (add:^rd fn-raw .~0.5) + (sub:^rd fn-raw .~0.5) + =/ n (floor-int fn) + =/ fn (san:^rd n) + :: Extended precision reduction + =/ y (sub:^rd x (mul:^rd fn pio2-1)) + =. y (sub:^rd y (mul:^rd fn pio2-1t)) + [(mod (abs:si n) 4) y] + :: + :: ============================================================ + :: KERNEL FUNCTIONS + :: ============================================================ + :: + :: +sindf: Kernel sin (degree 13) + :: + ++ sindf + |= x=@rd + ^- @rd + =/ z (mul:^rd x x) + :: Horner: S5 + z*S6 + =/ r (add:^rd sin-s5 (mul:^rd z sin-s6)) + :: S4 + z*r + =. r (add:^rd sin-s4 (mul:^rd z r)) + :: S3 + z*r + =. r (add:^rd sin-s3 (mul:^rd z r)) + :: S2 + z*r + =. r (add:^rd sin-s2 (mul:^rd z r)) + :: S1 + z*r + =. r (add:^rd sin-s1 (mul:^rd z r)) + :: x + x*z*r + (add:^rd x (mul:^rd (mul:^rd x z) r)) + :: + :: +cosdf: Kernel cos (degree 14) + :: + ++ cosdf + |= x=@rd + ^- @rd + =/ z (mul:^rd x x) + =/ w (mul:^rd z z) + :: Horner: C5 + z*C6 + =/ r (add:^rd cos-c5 (mul:^rd z cos-c6)) + :: C4 + z*r + =. r (add:^rd cos-c4 (mul:^rd z r)) + :: C3 + z*r + =. r (add:^rd cos-c3 (mul:^rd z r)) + :: C2 + z*r + =. r (add:^rd cos-c2 (mul:^rd z r)) + :: C1 + z*r + =. r (add:^rd cos-c1 (mul:^rd z r)) + :: 1 - z/2 + w*r + =/ hz (mul:^rd .~0.5 z) + (add:^rd (sub:^rd .~1 hz) (mul:^rd w r)) + :: + :: ============================================================ + :: PUBLIC API + :: ============================================================ + :: + ++ sin + |= x=@rd + ^- @rd + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + :: Tiny + ?: (lth:^rd ax `@rd`0x3e40.0000.0000.0000) x + :: Inf/NaN + ?: (gte:^rd ax `@rd`0x7ff0.0000.0000.0000) (sub:^rd x x) + =/ [n=@ y=@rd] (rem-pio2 x) + ?- n + %0 (sindf y) + %1 (cosdf y) + %2 (sub:^rd .~0 (sindf y)) + %3 (sub:^rd .~0 (cosdf y)) + == + :: + ++ cos + |= x=@rd + ^- @rd + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + ?: (lth:^rd ax `@rd`0x3e40.0000.0000.0000) .~1 + ?: (gte:^rd ax `@rd`0x7ff0.0000.0000.0000) (sub:^rd x x) + =/ [n=@ y=@rd] (rem-pio2 x) + ?- n + %0 (cosdf y) + %1 (sub:^rd .~0 (sindf y)) + %2 (sub:^rd .~0 (cosdf y)) + %3 (sindf y) + == + :: + ++ tan + |= x=@rd + ^- @rd + (div:^rd (sin x) (cos x)) + :: + ++ exp + |= x=@rd + ^- @rd + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + ?: =(x .~0) .~1 + ?: =(bits 0x7ff0.0000.0000.0000) `@rd`0x7ff0.0000.0000.0000 + ?: =(bits 0xfff0.0000.0000.0000) .~0 + ?: (gth:^rd ax `@rd`0x7ff0.0000.0000.0000) x + ?: (gth:^rd x exp-huge) `@rd`0x7ff0.0000.0000.0000 + ?: (lth:^rd x exp-tiny) .~0 + :: Argument reduction + =/ fn-raw (mul:^rd x invln2) + =/ fn ?: (gte:^rd fn-raw .~0) + (add:^rd fn-raw .~0.5) + (sub:^rd fn-raw .~0.5) + =/ k (floor-int fn) + =/ kf (san:^rd k) + =/ r (sub:^rd x (mul:^rd kf ln2hi)) + =. r (sub:^rd r (mul:^rd kf ln2lo)) + :: Polynomial (degree 5) + =/ r2 (mul:^rd r r) + =/ p (add:^rd exp-p4 (mul:^rd r exp-p5)) + =. p (add:^rd exp-p3 (mul:^rd r p)) + =. p (add:^rd exp-p2 (mul:^rd r p)) + =. p (add:^rd exp-p1 (mul:^rd r p)) + =/ expr (add:^rd .~1 (add:^rd r (mul:^rd r2 p))) + (scalbn expr k) + :: + :: +log: Natural logarithm (double precision) + :: + ++ log + |= x=@rd + ^- @rd + =/ bits `@`x + :: Special cases + ?: =(x .~0) `@rd`0xfff0.0000.0000.0000 :: log(0) = -inf + ?: =(bits 0x7ff0.0000.0000.0000) `@rd`0x7ff0.0000.0000.0000 + ?: (gth:^rd `@rd`(dis bits 0x7fff.ffff.ffff.ffff) `@rd`0x7ff0.0000.0000.0000) x + ?: !=(0 (rsh [0 63] bits)) (sub:^rd x x) + :: + :: Extract exponent and mantissa + :: + =/ exp-bits (dis (rsh [0 52] bits) 0x7ff) + =/ mant-bits (dis bits 0xf.ffff.ffff.ffff) + ?: =(exp-bits 0) `@rd`0xfff0.0000.0000.0000 + =/ k (sub exp-bits 1.023) + =/ m `@rd`(con 0x3ff0.0000.0000.0000 mant-bits) + :: + :: s = (m-1)/(m+1) + :: + =/ f (sub:^rd m .~1) + =/ s (div:^rd f (add:^rd .~2 f)) + =/ s2 (mul:^rd s s) + =/ s4 (mul:^rd s2 s2) + :: + :: Lg1-Lg7 from FreeBSD e_log.c (hex-exact) + :: + =/ lg1 `@rd`0x3fe5.5555.5555.5593 :: 0.6666666666666735130 + =/ lg2 `@rd`0x3fd9.9999.9997.fa04 :: 0.3999999999940941908 + =/ lg3 `@rd`0x3fd2.4924.9422.9359 :: 0.2857142874366239149 + =/ lg4 `@rd`0x3fcc.71c5.1d8e.78af :: 0.2222219843214978396 + =/ lg5 `@rd`0x3fc7.4664.96cb.03de :: 0.1818357216161805012 + =/ lg6 `@rd`0x3fc3.9a09.d078.c69f :: 0.1531383769920937332 + =/ lg7 `@rd`0x3fc2.f112.df3e.5244 :: 0.1479819860511658591 + :: + =/ t1 (mul:^rd s4 (add:^rd lg2 (mul:^rd s4 (add:^rd lg4 (mul:^rd s4 lg6))))) + =/ t2 (mul:^rd s2 (add:^rd lg1 (mul:^rd s4 (add:^rd lg3 (mul:^rd s4 (add:^rd lg5 (mul:^rd s4 lg7))))))) + =/ r (add:^rd t1 t2) + =/ log-m (mul:^rd (mul:^rd .~2 s) (add:^rd .~1 (mul:^rd s2 r))) + =/ kf (san:^rd (sun:si k)) + (add:^rd (mul:^rd kf ln2hi) log-m) + :: + ++ log-2 + |= x=@rd + ^- @rd + (mul:^rd (log x) invln2) + :: + ++ log-10 + |= x=@rd + ^- @rd + =/ invlog10 `@rd`0x3fdb.cb7b.1526.e50e :: 0.43429448190325176 + (mul:^rd (log x) invlog10) + :: + :: +sqt: Square root (Newton-Raphson, 6 iterations for double) + :: + ++ sqt + |= x=@rd + ^- @rd + =/ bits `@`x + ?: =(x .~0) .~0 + ?: =(bits 0x7ff0.0000.0000.0000) `@rd`0x7ff0.0000.0000.0000 + ?: (gth:^rd `@rd`(dis bits 0x7fff.ffff.ffff.ffff) `@rd`0x7ff0.0000.0000.0000) x + ?: !=(0 (rsh [0 63] bits)) (sub:^rd x x) + :: Initial guess + =/ g0 `@rd`(add (rsh [0 1] bits) 0x1ff0.0000.0000.0000) + =/ g g0 + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + =. g (mul:^rd .~0.5 (add:^rd g (div:^rd x g))) + g + :: + ++ sqrt sqt + :: + :: +atan: Inverse tangent (double precision) + :: + ++ atan + |= x=@rd + ^- @rd + =/ bits `@`x + =/ sign (rsh [0 63] bits) + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + :: + ?: (gte:^rd ax `@rd`0x7ff0.0000.0000.0000) + ?: =(bits 0x7ff0.0000.0000.0000) `@rd`0x3ff9.21fb.5444.2d18 :: pi/2 + ?: =(bits 0xfff0.0000.0000.0000) `@rd`0xbff9.21fb.5444.2d18 + x + :: + :: Coefficients from FreeBSD s_atan.c + :: + =/ at0 `@rd`0x3fd5.5555.5555.550d :: 0.33333333333329318027 + =/ at1 `@rd`0xbfc9.9999.9998.ebc4 :: -0.19999999998764832476 + =/ at2 `@rd`0x3fc2.4924.9200.83ff :: 0.14285714272503466371 + =/ at3 `@rd`0xbfbc.71c6.fe23.1671 :: -0.11111104054623557880 + =/ at4 `@rd`0x3fb7.45cd.c54c.206e :: 0.09090887133436506962 + =/ at5 `@rd`0xbfb3.b0f2.af74.9a6d :: -0.07691876205044829950 + =/ at6 `@rd`0x3fb1.0d66.a0d0.3d51 :: 0.06661073137387531207 + =/ at7 `@rd`0xbfad.de2d.52de.fd9a :: -0.05833570133790573485 + =/ at8 `@rd`0x3fa9.7b4b.2476.0deb :: 0.04976877994615936017 + =/ at9 `@rd`0xbfa2.b444.2c6a.6c2f :: -0.03653157274421691527 + =/ at10 `@rd`0x3f90.ad3a.e322.da11 :: 0.01628582015365782362 + :: + =/ pio2 `@rd`0x3ff9.21fb.5444.2d18 :: pi/2 + =/ pio4 `@rd`0x3fe9.21fb.5444.2d18 :: pi/4 + :: + =/ eval-poly + |= z=@rd + ^- @rd + =/ w (mul:^rd z z) + =/ w2 (mul:^rd w w) + =/ s1 (add:^rd at0 (mul:^rd w (add:^rd at1 (mul:^rd w (add:^rd at2 (mul:^rd w (add:^rd at3 (mul:^rd w at4)))))))) + =/ s2 (add:^rd at5 (mul:^rd w (add:^rd at6 (mul:^rd w (add:^rd at7 (mul:^rd w (add:^rd at8 (mul:^rd w (add:^rd at9 (mul:^rd w at10)))))))))) + =/ r (add:^rd s1 (mul:^rd w2 (mul:^rd w s2))) + (sub:^rd z (mul:^rd z (mul:^rd w r))) + :: + =/ lo-thresh `@rd`0x3fdc.0000.0000.0000 :: 0.4375 + =/ hi-thresh `@rd`0x4003.8000.0000.0000 :: 2.4375 + :: + =/ result=@rd + ?: (lte:^rd ax lo-thresh) + (eval-poly ax) + ?: (lte:^rd ax .~1) + =/ atan-half `@rd`0x3fdd.ac67.0561.bb4f + =/ t (div:^rd (sub:^rd ax .~0.5) (add:^rd .~1 (mul:^rd ax .~0.5))) + (add:^rd atan-half (eval-poly t)) + ?: (lte:^rd ax hi-thresh) + =/ t (div:^rd (sub:^rd ax .~1) (add:^rd .~1 ax)) + (add:^rd pio4 (eval-poly t)) + =/ t (div:^rd .~1 ax) + (sub:^rd pio2 (eval-poly t)) + :: + ?:(=(sign 0) result (sub:^rd .~0 result)) + :: + ++ asin + |= x=@rd + ^- @rd + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + ?: (gth:^rd ax .~1) (sub:^rd x x) + ?: =(ax .~1) + ?: =(0 (rsh [0 63] bits)) + `@rd`0x3ff9.21fb.5444.2d18 + `@rd`0xbff9.21fb.5444.2d18 + (atan (div:^rd x (sqt (sub:^rd .~1 (mul:^rd x x))))) + :: + ++ acos + |= x=@rd + ^- @rd + =/ bits `@`x + =/ ax `@rd`(dis bits 0x7fff.ffff.ffff.ffff) + =/ pi `@rd`0x4009.21fb.5444.2d18 + ?: (gth:^rd ax .~1) (sub:^rd x x) + ?: =(x .~1) .~0 + ?: =(x .~-1) pi + =/ s (sqt (sub:^rd .~1 (mul:^rd x x))) + ?: (gte:^rd x .~0) + (atan (div:^rd s x)) + (add:^rd pi (atan (div:^rd s x))) + :: + ++ atan2 + |= [y=@rd x=@rd] + ^- @rd + =/ pi `@rd`0x4009.21fb.5444.2d18 + =/ pio2 `@rd`0x3ff9.21fb.5444.2d18 + ?: (gth:^rd x .~0) + (atan (div:^rd y x)) + ?: &((lth:^rd x .~0) (gte:^rd y .~0)) + (add:^rd (atan (div:^rd y x)) pi) + ?: &((lth:^rd x .~0) (lth:^rd y .~0)) + (sub:^rd (atan (div:^rd y x)) pi) + ?: &(=(x .~0) (gth:^rd y .~0)) + pio2 + ?: &(=(x .~0) (lth:^rd y .~0)) + (sub:^rd .~0 pio2) + .~0 + :: + ++ pow-n + |= [x=@rd n=@rd] + ^- @rd + ?: =(n .~0) .~1 + =/ ni (abs:si (need (toi:^rd n))) + =/ neg (lth:^rd n .~0) + =/ result .~1 + =/ base x + =/ i 0 + |- + ?: |(=(ni 0) (gth i 63)) + ?:(neg (div:^rd .~1 result) result) + =/ new-result + ?: =(1 (dis ni 1)) + (mul:^rd result base) + result + $(ni (rsh [0 1] ni), base (mul:^rd base base), result new-result, i +(i)) + :: + ++ pow + |= [x=@rd n=@rd] + ^- @rd + =/ ni (toi:^rd n) + ?: &(?=(^ ni) =(n (san:^rd (need ni)))) + (pow-n x n) + (exp (mul:^rd n (log x))) + :: + :: ============================================================ + :: HYPERBOLIC FUNCTIONS + :: ============================================================ + :: + :: +sinh: Hyperbolic sine + :: sinh(x) = (exp(x) - exp(-x)) / 2 + :: + ++ sinh + |= x=@rd + ^- @rd + =/ ep (exp x) + =/ em (exp (sub:^rd .~0 x)) + (mul:^rd .~0.5 (sub:^rd ep em)) + :: + :: +cosh: Hyperbolic cosine + :: cosh(x) = (exp(x) + exp(-x)) / 2 + :: + ++ cosh + |= x=@rd + ^- @rd + =/ ep (exp x) + =/ em (exp (sub:^rd .~0 x)) + (mul:^rd .~0.5 (add:^rd ep em)) + :: + :: +tanh: Hyperbolic tangent + :: tanh(x) = (e^2x - 1) / (e^2x + 1) + :: + ++ tanh + |= x=@rd + ^- @rd + =/ e2x (exp (mul:^rd .~2 x)) + (div:^rd (sub:^rd e2x .~1) (add:^rd e2x .~1)) + :: + :: +asinh: Inverse hyperbolic sine + :: asinh(x) = ln(x + sqrt(x^2 + 1)) + :: + ++ asinh + |= x=@rd + ^- @rd + (log (add:^rd x (sqt (add:^rd (mul:^rd x x) .~1)))) + :: + :: +acosh: Inverse hyperbolic cosine + :: acosh(x) = ln(x + sqrt(x^2 - 1)), x >= 1 + :: + ++ acosh + |= x=@rd + ^- @rd + ?: (lth:^rd x .~1) (sub:^rd x x) :: NaN for x < 1 + (log (add:^rd x (sqt (sub:^rd (mul:^rd x x) .~1)))) + :: + :: +atanh: Inverse hyperbolic tangent + :: atanh(x) = 0.5 * ln((1+x)/(1-x)), |x| < 1 + :: + ++ atanh + |= x=@rd + ^- @rd + =/ ax (abs x) + ?: (gte:^rd ax .~1) (sub:^rd x x) :: NaN for |x| >= 1 + (mul:^rd .~0.5 (log (div:^rd (add:^rd .~1 x) (sub:^rd .~1 x)))) + :: + :: ============================================================ + :: SPECIAL FUNCTIONS + :: ============================================================ + :: + :: +erf: Error function + :: Uses Abramowitz & Stegun approximation 7.1.26 + :: + ++ erf + |= x=@rd + ^- @rd + =/ bits `@`x + =/ sign (rsh [0 63] bits) + =/ ax (abs x) + :: Coefficients from A&S (double precision) + =/ a1 .~0.254829592 + =/ a2 .~-0.284496736 + =/ a3 .~1.421413741 + =/ a4 .~-1.453152027 + =/ a5 .~1.061405429 + =/ p .~0.3275911 + :: + =/ t (div:^rd .~1 (add:^rd .~1 (mul:^rd p ax))) + =/ t2 (mul:^rd t t) + =/ t3 (mul:^rd t2 t) + =/ t4 (mul:^rd t3 t) + =/ t5 (mul:^rd t4 t) + :: + =/ poly %+ add:^rd (mul:^rd a1 t) + %+ add:^rd (mul:^rd a2 t2) + %+ add:^rd (mul:^rd a3 t3) + %+ add:^rd (mul:^rd a4 t4) + (mul:^rd a5 t5) + :: + =/ result (sub:^rd .~1 (mul:^rd poly (exp (sub:^rd .~0 (mul:^rd ax ax))))) + ?:(=(sign 0) result (sub:^rd .~0 result)) + :: + :: +gamma: Gamma function using Lanczos approximation + :: + ++ gamma + |= x=@rd + ^- @rd + ?: (lte:^rd x .~0) (sub:^rd x x) + ?: =(x .~1) .~1 + ?: =(x .~2) .~1 + :: Lanczos coefficients (g=7, n=9) for double precision + =/ g .~7 + =/ c0 .~0.99999999999980993 + =/ c1 .~676.5203681218851 + =/ c2 .~-1259.1392167224028 + =/ c3 .~771.32342877765313 + =/ c4 .~-176.61502916214059 + =/ c5 .~12.507343278686905 + =/ c6 .~-0.13857109526572012 + =/ c7 .~0.0000099843695780195716 + =/ c8 .~0.00000015056327351493116 + :: + =/ xm1 (sub:^rd x .~1) + =/ sum c0 + =. sum (add:^rd sum (div:^rd c1 (add:^rd xm1 .~1))) + =. sum (add:^rd sum (div:^rd c2 (add:^rd xm1 .~2))) + =. sum (add:^rd sum (div:^rd c3 (add:^rd xm1 .~3))) + =. sum (add:^rd sum (div:^rd c4 (add:^rd xm1 .~4))) + =. sum (add:^rd sum (div:^rd c5 (add:^rd xm1 .~5))) + =. sum (add:^rd sum (div:^rd c6 (add:^rd xm1 .~6))) + =. sum (add:^rd sum (div:^rd c7 (add:^rd xm1 .~7))) + =. sum (add:^rd sum (div:^rd c8 (add:^rd xm1 .~8))) + :: + =/ sqrt2pi .~2.5066282746310002 + =/ t (add:^rd xm1 (add:^rd g .~0.5)) + =/ term1 (mul:^rd sqrt2pi sum) + =/ term2 (pow t (add:^rd xm1 .~0.5)) + =/ term3 (exp (sub:^rd .~0 t)) + (mul:^rd term1 (mul:^rd term2 term3)) + :: + :: +lgamma: Log-gamma function + :: + ++ lgamma + |= x=@rd + ^- @rd + (log (gamma x)) + :: + :: +j0: Bessel function of first kind, order 0 + :: + ++ j0 + |= x=@rd + ^- @rd + =/ ax (abs x) + ?: (lth:^rd ax .~8) + =/ x2 (mul:^rd ax ax) + =/ x4 (mul:^rd x2 x2) + =/ x6 (mul:^rd x4 x2) + =/ x8 (mul:^rd x4 x4) + =/ x10 (mul:^rd x8 x2) + =/ c1 .~-0.25 + =/ c2 .~0.015625 + =/ c3 .~-0.00043402777777777775 + =/ c4 .~0.0000067816840277777775 + =/ c5 .~-0.00000006781684027777778 + %+ add:^rd .~1 + %+ add:^rd (mul:^rd c1 x2) + %+ add:^rd (mul:^rd c2 x4) + %+ add:^rd (mul:^rd c3 x6) + %+ add:^rd (mul:^rd c4 x8) + (mul:^rd c5 x10) + =/ z (div:^rd .~8 ax) + =/ z2 (mul:^rd z z) + =/ p0 .~1 + =/ p2 .~-0.00021 + =/ q0 .~-0.125 + =/ q2 .~0.00066 + =/ p (add:^rd p0 (mul:^rd z2 p2)) + =/ q (add:^rd q0 (mul:^rd z2 q2)) + =/ xx (sub:^rd ax .~0.7853981633974483) + =/ amp (div:^rd .~0.7978845608028654 (sqt ax)) + (mul:^rd amp (sub:^rd (mul:^rd p (cos xx)) (mul:^rd (mul:^rd z q) (sin xx)))) + :: + :: +j1: Bessel function of first kind, order 1 + :: + ++ j1 + |= x=@rd + ^- @rd + =/ bits `@`x + =/ sign (rsh [0 63] bits) + =/ ax (abs x) + =/ result=@rd + ?: (lth:^rd ax .~8) + =/ x2 (mul:^rd ax ax) + =/ x3 (mul:^rd x2 ax) + =/ x5 (mul:^rd x3 x2) + =/ x7 (mul:^rd x5 x2) + =/ x9 (mul:^rd x7 x2) + =/ c0 .~0.5 + =/ c1 .~-0.0625 + =/ c2 .~0.00260416666666666667 + =/ c3 .~-0.0000542534722222222 + =/ c4 .~0.000000677168402777778 + %+ add:^rd (mul:^rd c0 ax) + %+ add:^rd (mul:^rd c1 x3) + %+ add:^rd (mul:^rd c2 x5) + %+ add:^rd (mul:^rd c3 x7) + (mul:^rd c4 x9) + =/ z (div:^rd .~8 ax) + =/ z2 (mul:^rd z z) + =/ p1 .~1 + =/ p3 .~0.00084 + =/ q1 .~0.375 + =/ q3 .~-0.00195 + =/ p (add:^rd p1 (mul:^rd z2 p3)) + =/ q (add:^rd q1 (mul:^rd z2 q3)) + =/ xx (sub:^rd ax .~2.356194490192345) + =/ amp (div:^rd .~0.7978845608028654 (sqt ax)) + (mul:^rd amp (sub:^rd (mul:^rd p (cos xx)) (mul:^rd (mul:^rd z q) (sin xx)))) + ?:(=(sign 0) result (sub:^rd .~0 result)) + :: + :: +y0: Bessel function of second kind, order 0 + :: + ++ y0 + |= x=@rd + ^- @rd + ?: (lte:^rd x .~0) (sub:^rd x x) + =/ twoopi .~0.6366197723675814 + ?: (lth:^rd x .~8) + =/ j0x (j0 x) + =/ lnhx (log (mul:^rd .~0.5 x)) + (add:^rd (mul:^rd twoopi (mul:^rd j0x lnhx)) (mul:^rd .~0.36746691 (mul:^rd x x))) + =/ z (div:^rd .~8 x) + =/ z2 (mul:^rd z z) + =/ p0 .~1 + =/ q0 .~-0.125 + =/ p (add:^rd p0 (mul:^rd z2 .~-0.00021)) + =/ q (add:^rd q0 (mul:^rd z2 .~0.00066)) + =/ xx (sub:^rd x .~0.7853981633974483) + =/ amp (div:^rd .~0.7978845608028654 (sqt x)) + (mul:^rd amp (add:^rd (mul:^rd p (sin xx)) (mul:^rd (mul:^rd z q) (cos xx)))) + :: + :: +y1: Bessel function of second kind, order 1 + :: + ++ y1 + |= x=@rd + ^- @rd + ?: (lte:^rd x .~0) (sub:^rd x x) + =/ twoopi .~0.6366197723675814 + ?: (lth:^rd x .~8) + =/ j1x (j1 x) + =/ lnx (log x) + =/ invx (div:^rd .~1 x) + (mul:^rd twoopi (sub:^rd (mul:^rd j1x lnx) invx)) + =/ z (div:^rd .~8 x) + =/ z2 (mul:^rd z z) + =/ p1 .~1 + =/ q1 .~0.375 + =/ p (add:^rd p1 (mul:^rd z2 .~0.00084)) + =/ q (add:^rd q1 (mul:^rd z2 .~-0.00195)) + =/ xx (sub:^rd x .~2.356194490192345) + =/ amp (div:^rd .~0.7978845608028654 (sqt x)) + (mul:^rd amp (add:^rd (mul:^rd p (sin xx)) (mul:^rd (mul:^rd z q) (cos xx)))) + -- +:: +:: ================================================================ +:: QUAD PRECISION (@rq) - 128-bit IEEE 754 +:: ================================================================ +:: +:: Coefficients from FreeBSD msun ld128/k_sinl.c, k_cosl.c +:: Quad precision uses degree 25 sin, degree 22 cos for ~113 bits accuracy. +:: +:: Note: 128-bit hex literals are 32 hex digits each. +:: Format: 1 sign + 15 exponent (bias 16383) + 112 mantissa bits +:: +++ rq + |% + :: + :: ============================================================ + :: CONSTANTS (hex-exact for binary128) + :: ============================================================ + :: + :: Pi-related constants + :: + ++ pio2-1 `@rq`0x3fff.921f.b544.42d1.8469.898c.c517.01b8 :: pi/2 high + ++ pio2-1t `@rq`0x3f8c.c51c.e0ae.9ca4.e68e.6e36.30d1.5dc4 :: pi/2 - pio2_1 + ++ pio4 `@rq`0x3ffe.921f.b544.42d1.8469.898c.c517.01b8 :: pi/4 + ++ invpio2 `@rq`0x3ffe.45f3.06dc.9c88.2a53.f84e.afa3.ea6a :: 2/pi + :: + :: Exp constants + :: + ++ ln2hi `@rq`0x3ffe.62e4.2fef.a39e.f357.93c7.6730.07e6 :: ln(2) high + ++ ln2lo `@rq`0x3fbc.7abc.9e3b.3998.0f2f.707c.9ca4.4b10 :: ln(2) low + ++ invln2 `@rq`0x3fff.7154.7652.b82f.e177.7d0f.fda0.d23a :: 1/ln(2) + ++ exp-huge `@rq`0x400c.62e4.2fef.a39e.f357.93c7.6730.07e6 :: ~11356.5 (overflow) + ++ exp-tiny `@rq`0xc00c.6760.f0c7.f77a.8a37.8c9c.e5ca.10bc :: ~-11433.5 (underflow) + :: + :: ============================================================ + :: POLYNOMIAL COEFFICIENTS + :: ============================================================ + :: + :: Sin coefficients: degree 25 polynomial (S1-S12) + :: From FreeBSD ld128/k_sinl.c + :: sin(x) = x + x^3*(S1 + x^2*(S2 + ...)) + :: + ++ sin-s1 `@rq`0xbffc.5555.5555.5555.5555.5555.5555.5555 :: -1/6 + ++ sin-s2 `@rq`0x3ff8.1111.1111.1111.1111.1111.1111.1111 :: 1/120 + ++ sin-s3 `@rq`0xbff2.a01a.01a0.1a01.a01a.01a0.1a01.a01a :: -1/5040 + ++ sin-s4 `@rq`0x3fec.71de.3a55.6c73.38fa.ac1c.88e5.0017 :: 1/362880 + ++ sin-s5 `@rq`0xbfe5.ae64.567f.544e.38fe.747e.4b60.c9a6 :: -1/39916800 + ++ sin-s6 `@rq`0x3fde.6124.613a.86d0.9393.30de.0b58.0e66 :: 1/6227020800 + ++ sin-s7 `@rq`0xbfd6.ae7f.3e73.3b81.f7be.ced5.1769.e17e :: ~-1.15e-15 + ++ sin-s8 `@rq`0x3fce.952c.7703.0ad4.a71d.fc88.40d2.35fc :: ~8.22e-18 + :: + :: Cos coefficients: degree 22 polynomial (C1-C11) + :: From FreeBSD ld128/k_cosl.c + :: cos(x) = 1 - x^2/2 + x^4*(C1 + x^2*(C2 + ...)) + :: + ++ cos-c1 `@rq`0x3ffa.5555.5555.5555.5555.5555.5555.5555 :: 1/24 + ++ cos-c2 `@rq`0xbff5.6c16.c16c.16c1.6c16.c16c.16c1.6c17 :: -1/720 + ++ cos-c3 `@rq`0x3fef.a01a.01a0.1a01.a01a.01a0.1a01.a01a :: 1/40320 + ++ cos-c4 `@rq`0xbfe9.27e4.fb77.89f5.c72e.f016.d3ea.6679 :: -1/3628800 + ++ cos-c5 `@rq`0x3fe2.1eed.8eff.8d89.7b5e.33bf.1a9b.a5c9 :: 1/479001600 + ++ cos-c6 `@rq`0xbfda.93f2.7dbb.c4fb.7459.4a58.2e61.6a96 :: ~-8.9e-14 + ++ cos-c7 `@rq`0x3fd2.ae7f.3e73.3b94.0c61.7a76.bec4.4b09 :: ~4.78e-16 + ++ cos-c8 `@rq`0xbfca.6827.863b.97d9.7e4c.9bbd.c044.2faa :: ~-1.56e-18 + :: + :: Exp polynomial coefficients (degree 7 for quad) + :: + ++ exp-p1 `@rq`0x3ffe.0000.0000.0000.0000.0000.0000.0000 :: 0.5 + ++ exp-p2 `@rq`0x3ffc.5555.5555.5555.5555.5555.5555.5555 :: 1/6 + ++ exp-p3 `@rq`0x3ffa.5555.5555.5555.5555.5555.5555.5555 :: 1/24 + ++ exp-p4 `@rq`0x3ff8.1111.1111.1111.1111.1111.1111.1111 :: 1/120 + ++ exp-p5 `@rq`0x3ff5.6c16.c16c.16c1.6c16.c16c.16c1.6c17 :: 1/720 + ++ exp-p6 `@rq`0x3ff2.a01a.01a0.1a01.a01a.01a0.1a01.a01a :: 1/5040 + ++ exp-p7 `@rq`0x3fef.a01a.01a0.1a01.a01a.01a0.1a01.a01a :: 1/40320 + :: + :: ============================================================ + :: CHEBYSHEV POLYNOMIAL + :: ============================================================ + :: + ++ cheb + |= n=@ud + |= x=@rq + ^- @rq + ?: =(0 n) .~~~1 + ?: =(1 n) x + (sub:^rq (mul:^rq .~~~2 (mul:^rq x $(n (dec n)))) $(n (sub n 2))) + :: + :: ============================================================ + :: HELPER FUNCTIONS + :: ============================================================ + :: + ++ scalbn + |= [x=@rq n=@sd] + ^- @rq + =/ bits `@`x + =/ exp (dis (rsh [0 112] bits) 0x7fff) + =/ sign (rsh [0 127] bits) + ?: =(exp 0x7fff) x :: inf or nan + ?: =(exp 0) x :: zero or subnormal + =/ new-exp (sum:si (sun:si exp) n) + ?: (gth:si new-exp (sun:si 32.766)) + ?: =(sign 0) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000 + `@rq`0xffff.0000.0000.0000.0000.0000.0000.0000 + ?: (lth:si new-exp (sun:si 1)) + .~~~0 + =/ new-bits + %+ con (lsh [0 127] sign) + %+ con (lsh [0 112] (abs:si new-exp)) + (dis bits 0xffff.ffff.ffff.ffff.ffff.ffff.ffff) + `@rq`new-bits + :: + ++ floor-int + |= x=@rq + ^- @sd + =/ bits `@`x + =/ sign (rsh [0 127] bits) + =/ bexp (dis (rsh [0 112] bits) 0x7fff) + ?: (^lth bexp 16.383) + ?:(=(sign 0) --0 -1) + =/ exp (sub bexp 16.383) + =/ mant (con (dis bits 0xffff.ffff.ffff.ffff.ffff.ffff.ffff) 0x1.0000.0000.0000.0000.0000.0000.0000) + ?: (^gte exp 112) + ?: =(sign 0) + (sun:si (lsh [0 (sub exp 112)] mant)) + (new:si %.n (lsh [0 (sub exp 112)] mant)) + =/ shift (sub 112 exp) + =/ int-mant (rsh [0 shift] mant) + =/ frac-mask (dec (lsh [0 shift] 1)) + =/ had-frac !=(0 (dis mant frac-mask)) + =/ result ?:(=(sign 0) (sun:si int-mant) (new:si %.n int-mant)) + ?: &(=(sign 1) had-frac) + (dif:si result --1) + result + :: + :: ============================================================ + :: ARGUMENT REDUCTION + :: ============================================================ + :: + ++ rem-pio2 + |= x=@rq + ^- [n=@ y=@rq] + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: (lte:^rq ax pio4) + [0 x] + =/ fn-raw (mul:^rq x invpio2) + =/ fn ?: (gte:^rq fn-raw .~~~0) + (add:^rq fn-raw .~~~0.5) + (sub:^rq fn-raw .~~~0.5) + =/ n (floor-int fn) + =/ fn (san:^rq n) + =/ y (sub:^rq x (mul:^rq fn pio2-1)) + =. y (sub:^rq y (mul:^rq fn pio2-1t)) + [(mod (abs:si n) 4) y] + :: + :: ============================================================ + :: KERNEL FUNCTIONS + :: ============================================================ + :: + :: +sindf: Kernel sin (degree 17 for practical quad precision) + :: + ++ sindf + |= x=@rq + ^- @rq + =/ z (mul:^rq x x) + :: Horner form with 8 coefficients + =/ r (add:^rq sin-s7 (mul:^rq z sin-s8)) + =. r (add:^rq sin-s6 (mul:^rq z r)) + =. r (add:^rq sin-s5 (mul:^rq z r)) + =. r (add:^rq sin-s4 (mul:^rq z r)) + =. r (add:^rq sin-s3 (mul:^rq z r)) + =. r (add:^rq sin-s2 (mul:^rq z r)) + =. r (add:^rq sin-s1 (mul:^rq z r)) + (add:^rq x (mul:^rq (mul:^rq x z) r)) + :: + :: +cosdf: Kernel cos (degree 16 for practical quad precision) + :: + ++ cosdf + |= x=@rq + ^- @rq + =/ z (mul:^rq x x) + =/ w (mul:^rq z z) + :: Horner form with 8 coefficients + =/ r (add:^rq cos-c7 (mul:^rq z cos-c8)) + =. r (add:^rq cos-c6 (mul:^rq z r)) + =. r (add:^rq cos-c5 (mul:^rq z r)) + =. r (add:^rq cos-c4 (mul:^rq z r)) + =. r (add:^rq cos-c3 (mul:^rq z r)) + =. r (add:^rq cos-c2 (mul:^rq z r)) + =. r (add:^rq cos-c1 (mul:^rq z r)) + =/ hz (mul:^rq .~~~0.5 z) + (add:^rq (sub:^rq .~~~1 hz) (mul:^rq w r)) + :: + :: ============================================================ + :: PUBLIC API + :: ============================================================ + :: + ++ sin + |= x=@rq + ^- @rq + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + :: Tiny + ?: (lth:^rq ax `@rq`0x3f8f.0000.0000.0000.0000.0000.0000.0000) x + :: Inf/NaN + ?: (gte:^rq ax `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) + (sub:^rq x x) + =/ [n=@ y=@rq] (rem-pio2 x) + ?- n + %0 (sindf y) + %1 (cosdf y) + %2 (sub:^rq .~~~0 (sindf y)) + %3 (sub:^rq .~~~0 (cosdf y)) + == + :: + ++ cos + |= x=@rq + ^- @rq + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: (lth:^rq ax `@rq`0x3f8f.0000.0000.0000.0000.0000.0000.0000) .~~~1 + ?: (gte:^rq ax `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) + (sub:^rq x x) + =/ [n=@ y=@rq] (rem-pio2 x) + ?- n + %0 (cosdf y) + %1 (sub:^rq .~~~0 (sindf y)) + %2 (sub:^rq .~~~0 (cosdf y)) + %3 (sindf y) + == + :: + ++ tan + |= x=@rq + ^- @rq + (div:^rq (sin x) (cos x)) + :: + ++ exp + |= x=@rq + ^- @rq + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: =(x .~~~0) .~~~1 + ?: =(bits 0x7fff.0000.0000.0000.0000.0000.0000.0000) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000 + ?: =(bits 0xffff.0000.0000.0000.0000.0000.0000.0000) .~~~0 + ?: (gth:^rq ax `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) x + ?: (gth:^rq x exp-huge) `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000 + ?: (lth:^rq x exp-tiny) .~~~0 + :: Argument reduction + =/ fn-raw (mul:^rq x invln2) + =/ fn ?: (gte:^rq fn-raw .~~~0) + (add:^rq fn-raw .~~~0.5) + (sub:^rq fn-raw .~~~0.5) + =/ k (floor-int fn) + =/ kf (san:^rq k) + =/ r (sub:^rq x (mul:^rq kf ln2hi)) + =. r (sub:^rq r (mul:^rq kf ln2lo)) + :: Polynomial (degree 7) + =/ r2 (mul:^rq r r) + =/ p (add:^rq exp-p6 (mul:^rq r exp-p7)) + =. p (add:^rq exp-p5 (mul:^rq r p)) + =. p (add:^rq exp-p4 (mul:^rq r p)) + =. p (add:^rq exp-p3 (mul:^rq r p)) + =. p (add:^rq exp-p2 (mul:^rq r p)) + =. p (add:^rq exp-p1 (mul:^rq r p)) + =/ expr (add:^rq .~~~1 (add:^rq r (mul:^rq r2 p))) + (scalbn expr k) + :: + :: +log: Natural logarithm (quad precision) + :: + :: Uses higher-degree polynomial for 113-bit mantissa accuracy + :: + ++ log + |= x=@rq + ^- @rq + =/ bits `@`x + :: Special cases + ?: =(x .~~~0) `@rq`0xffff.0000.0000.0000.0000.0000.0000.0000 :: -inf + ?: =(bits 0x7fff.0000.0000.0000.0000.0000.0000.0000) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000 :: +inf + ?: (gth:^rq `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) + x :: nan + ?: !=(0 (rsh [0 127] bits)) (sub:^rq x x) :: negative + :: + :: Extract exponent and mantissa + :: Quad: 1 sign, 15 exp (bias 16383), 112 mantissa + :: + =/ exp-bits (dis (rsh [0 112] bits) 0x7fff) + =/ mant-bits (dis bits 0xffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: =(exp-bits 0) `@rq`0xffff.0000.0000.0000.0000.0000.0000.0000 + =/ k (sub exp-bits 16.383) + =/ m `@rq`(con 0x3fff.0000.0000.0000.0000.0000.0000.0000 mant-bits) + :: + :: s = (m-1)/(m+1) + :: + =/ f (sub:^rq m .~~~1) + =/ s (div:^rq f (add:^rq .~~~2 f)) + =/ s2 (mul:^rq s s) + =/ s4 (mul:^rq s2 s2) + :: + :: Lg1-Lg8 for quad precision (higher accuracy) + :: + =/ lg1 `@rq`0x3ffc.5555.5555.5555.5555.5555.5555.5555 :: 2/3 + =/ lg2 `@rq`0x3ffb.9999.9999.9999.9999.9999.9999.999a :: 2/5 + =/ lg3 `@rq`0x3ffb.2492.4924.9249.2492.4924.9249.2492 :: 2/7 + =/ lg4 `@rq`0x3ffa.c71c.71c7.1c71.c71c.71c7.1c71.c71c :: 2/9 + =/ lg5 `@rq`0x3ffa.745d.1745.d174.5d17.45d1.745d.1746 :: 2/11 + =/ lg6 `@rq`0x3ffa.3a83.a83a.83a8.3a83.a83a.83a8.3a84 :: 2/13 + =/ lg7 `@rq`0x3ffa.0f0f.0f0f.0f0f.0f0f.0f0f.0f0f.0f0f :: 2/15 + =/ lg8 `@rq`0x3ff9.e1e1.e1e1.e1e1.e1e1.e1e1.e1e1.e1e2 :: 2/17 + :: + =/ t1 (mul:^rq s4 (add:^rq lg2 (mul:^rq s4 (add:^rq lg4 (mul:^rq s4 (add:^rq lg6 (mul:^rq s4 lg8))))))) + =/ t2 (mul:^rq s2 (add:^rq lg1 (mul:^rq s4 (add:^rq lg3 (mul:^rq s4 (add:^rq lg5 (mul:^rq s4 lg7))))))) + =/ r (add:^rq t1 t2) + =/ log-m (mul:^rq (mul:^rq .~~~2 s) (add:^rq .~~~1 (mul:^rq s2 r))) + :: + =/ kf (san:^rq (sun:si k)) + (add:^rq (mul:^rq kf ln2hi) log-m) + :: + ++ log-2 + |= x=@rq + ^- @rq + (mul:^rq (log x) invln2) + :: + ++ log-10 + |= x=@rq + ^- @rq + =/ invlog10 `@rq`0x3ffd.bcb7.b152.6e50.e32a.6ab7.555f.5a68 :: 1/ln(10) + (mul:^rq (log x) invlog10) + :: + :: +sqt: Square root (Newton-Raphson, 8 iterations for quad) + :: + ++ sqt + |= x=@rq + ^- @rq + =/ bits `@`x + ?: =(x .~~~0) .~~~0 + ?: =(bits 0x7fff.0000.0000.0000.0000.0000.0000.0000) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000 + ?: (gth:^rq `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) + x + ?: !=(0 (rsh [0 127] bits)) (sub:^rq x x) + :: Initial guess + =/ g0 `@rq`(add (rsh [0 1] bits) 0x1fff.8000.0000.0000.0000.0000.0000.0000) + =/ g g0 + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + =. g (mul:^rq .~~~0.5 (add:^rq g (div:^rq x g))) + g + :: + ++ sqrt sqt + :: + :: +atan: Inverse tangent (quad precision) + :: + ++ atan + |= x=@rq + ^- @rq + =/ bits `@`x + =/ sign (rsh [0 127] bits) + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + :: + ?: (gte:^rq ax `@rq`0x7fff.0000.0000.0000.0000.0000.0000.0000) + ?: =(bits 0x7fff.0000.0000.0000.0000.0000.0000.0000) + `@rq`0x3fff.921f.b544.42d1.8469.898c.c517.01b8 :: pi/2 + ?: =(bits 0xffff.0000.0000.0000.0000.0000.0000.0000) + `@rq`0xbfff.921f.b544.42d1.8469.898c.c517.01b8 :: -pi/2 + x + :: + :: Coefficients for atan polynomial (higher degree for quad) + :: + =/ at0 `@rq`0x3ffd.5555.5555.5555.5555.5555.5555.5555 :: 1/3 + =/ at1 `@rq`0xbffc.cccc.cccc.cccc.cccc.cccc.cccc.cccd :: -1/5 + =/ at2 `@rq`0x3ffc.4924.9249.2492.4924.9249.2492.4924 :: 1/7 + =/ at3 `@rq`0xbffc.0000.0000.0000.0000.0000.0000.0000 :: -1/9 + =/ at4 `@rq`0x3ffb.a2e8.ba2e.8ba2.e8ba.2e8b.a2e8.ba2f :: 1/11 + =/ at5 `@rq`0xbffb.5555.5555.5555.5555.5555.5555.5555 :: -1/13 + :: + =/ pio2 `@rq`0x3fff.921f.b544.42d1.8469.898c.c517.01b8 :: pi/2 + =/ pio4 `@rq`0x3ffe.921f.b544.42d1.8469.898c.c517.01b8 :: pi/4 + :: + =/ eval-poly + |= z=@rq + ^- @rq + =/ w (mul:^rq z z) + =/ w2 (mul:^rq w w) + =/ s1 (add:^rq at0 (mul:^rq w (add:^rq at1 (mul:^rq w at2)))) + =/ s2 (add:^rq at3 (mul:^rq w (add:^rq at4 (mul:^rq w at5)))) + =/ r (add:^rq s1 (mul:^rq w2 (mul:^rq w s2))) + (sub:^rq z (mul:^rq z (mul:^rq w r))) + :: + =/ lo-thresh `@rq`0x3ffe.0000.0000.0000.0000.0000.0000.0000 :: 0.5 + =/ hi-thresh `@rq`0x4000.8000.0000.0000.0000.0000.0000.0000 :: 2.0 + :: + =/ result=@rq + ?: (lte:^rq ax lo-thresh) + (eval-poly ax) + ?: (lte:^rq ax .~~~1) + =/ atan-half `@rq`0x3ffd.dac6.7056.1bb4.f68a.dec4.befd.cd8b + =/ t (div:^rq (sub:^rq ax .~~~0.5) (add:^rq .~~~1 (mul:^rq ax .~~~0.5))) + (add:^rq atan-half (eval-poly t)) + ?: (lte:^rq ax hi-thresh) + =/ t (div:^rq (sub:^rq ax .~~~1) (add:^rq .~~~1 ax)) + (add:^rq pio4 (eval-poly t)) + =/ t (div:^rq .~~~1 ax) + (sub:^rq pio2 (eval-poly t)) + :: + ?:(=(sign 0) result (sub:^rq .~~~0 result)) + :: + ++ asin + |= x=@rq + ^- @rq + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: (gth:^rq ax .~~~1) (sub:^rq x x) + ?: =(ax .~~~1) + ?: =(0 (rsh [0 127] bits)) + `@rq`0x3fff.921f.b544.42d1.8469.898c.c517.01b8 + `@rq`0xbfff.921f.b544.42d1.8469.898c.c517.01b8 + (atan (div:^rq x (sqt (sub:^rq .~~~1 (mul:^rq x x))))) + :: + ++ acos + |= x=@rq + ^- @rq + =/ bits `@`x + =/ ax `@rq`(dis bits 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + =/ pi `@rq`0x4000.921f.b544.42d1.8469.898c.c517.01b8 + ?: (gth:^rq ax .~~~1) (sub:^rq x x) + ?: =(x .~~~1) .~~~0 + ?: =(x .~~~-1) pi + =/ s (sqt (sub:^rq .~~~1 (mul:^rq x x))) + ?: (gte:^rq x .~~~0) + (atan (div:^rq s x)) + (add:^rq pi (atan (div:^rq s x))) + :: + ++ atan2 + |= [y=@rq x=@rq] + ^- @rq + =/ pi `@rq`0x4000.921f.b544.42d1.8469.898c.c517.01b8 + =/ pio2 `@rq`0x3fff.921f.b544.42d1.8469.898c.c517.01b8 + ?: (gth:^rq x .~~~0) + (atan (div:^rq y x)) + ?: &((lth:^rq x .~~~0) (gte:^rq y .~~~0)) + (add:^rq (atan (div:^rq y x)) pi) + ?: &((lth:^rq x .~~~0) (lth:^rq y .~~~0)) + (sub:^rq (atan (div:^rq y x)) pi) + ?: &(=(x .~~~0) (gth:^rq y .~~~0)) + pio2 + ?: &(=(x .~~~0) (lth:^rq y .~~~0)) + (sub:^rq .~~~0 pio2) + .~~~0 + :: + ++ pow-n + |= [x=@rq n=@rq] + ^- @rq + ?: =(n .~~~0) .~~~1 + =/ ni (abs:si (need (toi:^rq n))) + =/ neg (lth:^rq n .~~~0) + =/ result .~~~1 + =/ base x + =/ i 0 + |- + ?: |(=(ni 0) (gth i 127)) + ?:(neg (div:^rq .~~~1 result) result) + =/ new-result + ?: =(1 (dis ni 1)) + (mul:^rq result base) + result + $(ni (rsh [0 1] ni), base (mul:^rq base base), result new-result, i +(i)) + :: + ++ pow + |= [x=@rq n=@rq] + ^- @rq + =/ ni (toi:^rq n) + ?: &(?=(^ ni) =(n (san:^rq (need ni)))) + (pow-n x n) + (exp (mul:^rq n (log x))) + :: + :: ============================================================ + :: HYPERBOLIC FUNCTIONS + :: ============================================================ + :: + :: +sinh: Hyperbolic sine + :: sinh(x) = (exp(x) - exp(-x)) / 2 + :: + ++ sinh + |= x=@rq + ^- @rq + =/ ep (exp x) + =/ em (exp (sub:^rq .~~~0 x)) + (mul:^rq .~~~0.5 (sub:^rq ep em)) + :: + :: +cosh: Hyperbolic cosine + :: cosh(x) = (exp(x) + exp(-x)) / 2 + :: + ++ cosh + |= x=@rq + ^- @rq + =/ ep (exp x) + =/ em (exp (sub:^rq .~~~0 x)) + (mul:^rq .~~~0.5 (add:^rq ep em)) + :: + :: +tanh: Hyperbolic tangent + :: tanh(x) = (e^2x - 1) / (e^2x + 1) + :: + ++ tanh + |= x=@rq + ^- @rq + =/ e2x (exp (mul:^rq .~~~2 x)) + (div:^rq (sub:^rq e2x .~~~1) (add:^rq e2x .~~~1)) + :: + :: +asinh: Inverse hyperbolic sine + :: asinh(x) = ln(x + sqrt(x^2 + 1)) + :: + ++ asinh + |= x=@rq + ^- @rq + (log (add:^rq x (sqt (add:^rq (mul:^rq x x) .~~~1)))) + :: + :: +acosh: Inverse hyperbolic cosine + :: acosh(x) = ln(x + sqrt(x^2 - 1)), x >= 1 + :: + ++ acosh + |= x=@rq + ^- @rq + ?: (lth:^rq x .~~~1) (sub:^rq x x) :: NaN for x < 1 + (log (add:^rq x (sqt (sub:^rq (mul:^rq x x) .~~~1)))) + :: + :: +atanh: Inverse hyperbolic tangent + :: atanh(x) = 0.5 * ln((1+x)/(1-x)), |x| < 1 + :: + ++ atanh + |= x=@rq + ^- @rq + =/ ax (abs x) + ?: (gte:^rq ax .~~~1) (sub:^rq x x) :: NaN for |x| >= 1 + (mul:^rq .~~~0.5 (log (div:^rq (add:^rq .~~~1 x) (sub:^rq .~~~1 x)))) + :: + :: ============================================================ + :: SPECIAL FUNCTIONS + :: ============================================================ + :: + :: +erf: Error function + :: Uses Abramowitz & Stegun approximation 7.1.26 + :: + ++ erf + |= x=@rq + ^- @rq + =/ bits `@`x + =/ sign (rsh [0 127] bits) + =/ ax (abs x) + :: Coefficients from A&S (quad precision) + =/ a1 .~~~0.254829592 + =/ a2 .~~~-0.284496736 + =/ a3 .~~~1.421413741 + =/ a4 .~~~-1.453152027 + =/ a5 .~~~1.061405429 + =/ p .~~~0.3275911 + :: + =/ t (div:^rq .~~~1 (add:^rq .~~~1 (mul:^rq p ax))) + =/ t2 (mul:^rq t t) + =/ t3 (mul:^rq t2 t) + =/ t4 (mul:^rq t3 t) + =/ t5 (mul:^rq t4 t) + :: + =/ poly %+ add:^rq (mul:^rq a1 t) + %+ add:^rq (mul:^rq a2 t2) + %+ add:^rq (mul:^rq a3 t3) + %+ add:^rq (mul:^rq a4 t4) + (mul:^rq a5 t5) + :: + =/ result (sub:^rq .~~~1 (mul:^rq poly (exp (sub:^rq .~~~0 (mul:^rq ax ax))))) + ?:(=(sign 0) result (sub:^rq .~~~0 result)) + :: + :: +gamma: Gamma function using Lanczos approximation + :: + ++ gamma + |= x=@rq + ^- @rq + ?: (lte:^rq x .~~~0) (sub:^rq x x) + ?: =(x .~~~1) .~~~1 + ?: =(x .~~~2) .~~~1 + :: Lanczos coefficients (g=7, n=9) + =/ g .~~~7 + =/ c0 .~~~0.99999999999980993 + =/ c1 .~~~676.5203681218851 + =/ c2 .~~~-1259.1392167224028 + =/ c3 .~~~771.32342877765313 + =/ c4 .~~~-176.61502916214059 + =/ c5 .~~~12.507343278686905 + =/ c6 .~~~-0.13857109526572012 + =/ c7 .~~~0.0000099843695780195716 + =/ c8 .~~~0.00000015056327351493116 + :: + =/ xm1 (sub:^rq x .~~~1) + =/ sum c0 + =. sum (add:^rq sum (div:^rq c1 (add:^rq xm1 .~~~1))) + =. sum (add:^rq sum (div:^rq c2 (add:^rq xm1 .~~~2))) + =. sum (add:^rq sum (div:^rq c3 (add:^rq xm1 .~~~3))) + =. sum (add:^rq sum (div:^rq c4 (add:^rq xm1 .~~~4))) + =. sum (add:^rq sum (div:^rq c5 (add:^rq xm1 .~~~5))) + =. sum (add:^rq sum (div:^rq c6 (add:^rq xm1 .~~~6))) + =. sum (add:^rq sum (div:^rq c7 (add:^rq xm1 .~~~7))) + =. sum (add:^rq sum (div:^rq c8 (add:^rq xm1 .~~~8))) + :: + =/ sqrt2pi .~~~2.5066282746310002 + =/ t (add:^rq xm1 (add:^rq g .~~~0.5)) + =/ term1 (mul:^rq sqrt2pi sum) + =/ term2 (pow t (add:^rq xm1 .~~~0.5)) + =/ term3 (exp (sub:^rq .~~~0 t)) + (mul:^rq term1 (mul:^rq term2 term3)) + :: + :: +lgamma: Log-gamma function + :: + ++ lgamma + |= x=@rq + ^- @rq + (log (gamma x)) + :: + :: +j0: Bessel function of first kind, order 0 + :: + ++ j0 + |= x=@rq + ^- @rq + =/ ax (abs x) + ?: (lth:^rq ax .~~~8) + =/ x2 (mul:^rq ax ax) + =/ x4 (mul:^rq x2 x2) + =/ x6 (mul:^rq x4 x2) + =/ x8 (mul:^rq x4 x4) + =/ x10 (mul:^rq x8 x2) + =/ x12 (mul:^rq x6 x6) + =/ c1 .~~~-0.25 + =/ c2 .~~~0.015625 + =/ c3 .~~~-0.00043402777777777775 + =/ c4 .~~~0.0000067816840277777775 + =/ c5 .~~~-0.00000006781684027777778 + =/ c6 .~~~0.0000000004709502797067901 + %+ add:^rq .~~~1 + %+ add:^rq (mul:^rq c1 x2) + %+ add:^rq (mul:^rq c2 x4) + %+ add:^rq (mul:^rq c3 x6) + %+ add:^rq (mul:^rq c4 x8) + %+ add:^rq (mul:^rq c5 x10) + (mul:^rq c6 x12) + =/ z (div:^rq .~~~8 ax) + =/ z2 (mul:^rq z z) + =/ p0 .~~~1 + =/ p2 .~~~-0.00021 + =/ q0 .~~~-0.125 + =/ q2 .~~~0.00066 + =/ p (add:^rq p0 (mul:^rq z2 p2)) + =/ q (add:^rq q0 (mul:^rq z2 q2)) + =/ xx (sub:^rq ax .~~~0.7853981633974483) + =/ amp (div:^rq .~~~0.7978845608028654 (sqt ax)) + (mul:^rq amp (sub:^rq (mul:^rq p (cos xx)) (mul:^rq (mul:^rq z q) (sin xx)))) + :: + :: +j1: Bessel function of first kind, order 1 + :: + ++ j1 + |= x=@rq + ^- @rq + =/ bits `@`x + =/ sign (rsh [0 127] bits) + =/ ax (abs x) + =/ result=@rq + ?: (lth:^rq ax .~~~8) + =/ x2 (mul:^rq ax ax) + =/ x3 (mul:^rq x2 ax) + =/ x5 (mul:^rq x3 x2) + =/ x7 (mul:^rq x5 x2) + =/ x9 (mul:^rq x7 x2) + =/ x11 (mul:^rq x9 x2) + =/ c0 .~~~0.5 + =/ c1 .~~~-0.0625 + =/ c2 .~~~0.00260416666666666667 + =/ c3 .~~~-0.0000542534722222222 + =/ c4 .~~~0.000000677168402777778 + =/ c5 .~~~-0.00000000564307002314815 + %+ add:^rq (mul:^rq c0 ax) + %+ add:^rq (mul:^rq c1 x3) + %+ add:^rq (mul:^rq c2 x5) + %+ add:^rq (mul:^rq c3 x7) + %+ add:^rq (mul:^rq c4 x9) + (mul:^rq c5 x11) + =/ z (div:^rq .~~~8 ax) + =/ z2 (mul:^rq z z) + =/ p1 .~~~1 + =/ p3 .~~~0.00084 + =/ q1 .~~~0.375 + =/ q3 .~~~-0.00195 + =/ p (add:^rq p1 (mul:^rq z2 p3)) + =/ q (add:^rq q1 (mul:^rq z2 q3)) + =/ xx (sub:^rq ax .~~~2.356194490192345) + =/ amp (div:^rq .~~~0.7978845608028654 (sqt ax)) + (mul:^rq amp (sub:^rq (mul:^rq p (cos xx)) (mul:^rq (mul:^rq z q) (sin xx)))) + ?:(=(sign 0) result (sub:^rq .~~~0 result)) + :: + :: +y0: Bessel function of second kind, order 0 + :: + ++ y0 + |= x=@rq + ^- @rq + ?: (lte:^rq x .~~~0) (sub:^rq x x) + =/ twoopi .~~~0.6366197723675814 + ?: (lth:^rq x .~~~8) + =/ j0x (j0 x) + =/ lnhx (log (mul:^rq .~~~0.5 x)) + (add:^rq (mul:^rq twoopi (mul:^rq j0x lnhx)) (mul:^rq .~~~0.36746691 (mul:^rq x x))) + =/ z (div:^rq .~~~8 x) + =/ z2 (mul:^rq z z) + =/ p0 .~~~1 + =/ q0 .~~~-0.125 + =/ p (add:^rq p0 (mul:^rq z2 .~~~-0.00021)) + =/ q (add:^rq q0 (mul:^rq z2 .~~~0.00066)) + =/ xx (sub:^rq x .~~~0.7853981633974483) + =/ amp (div:^rq .~~~0.7978845608028654 (sqt x)) + (mul:^rq amp (add:^rq (mul:^rq p (sin xx)) (mul:^rq (mul:^rq z q) (cos xx)))) + :: + :: +y1: Bessel function of second kind, order 1 + :: + ++ y1 + |= x=@rq + ^- @rq + ?: (lte:^rq x .~~~0) (sub:^rq x x) + =/ twoopi .~~~0.6366197723675814 + ?: (lth:^rq x .~~~8) + =/ j1x (j1 x) + =/ lnx (log x) + =/ invx (div:^rq .~~~1 x) + (mul:^rq twoopi (sub:^rq (mul:^rq j1x lnx) invx)) + =/ z (div:^rq .~~~8 x) + =/ z2 (mul:^rq z z) + =/ p1 .~~~1 + =/ q1 .~~~0.375 + =/ p (add:^rq p1 (mul:^rq z2 .~~~0.00084)) + =/ q (add:^rq q1 (mul:^rq z2 .~~~-0.00195)) + =/ xx (sub:^rq x .~~~2.356194490192345) + =/ amp (div:^rq .~~~0.7978845608028654 (sqt x)) + (mul:^rq amp (add:^rq (mul:^rq p (sin xx)) (mul:^rq (mul:^rq z q) (cos xx)))) + -- +:: +:: ================================================================ +:: JET IMPLEMENTATION GUIDE +:: ================================================================ +:: +:: For C jets to produce identical results: +:: +:: 1. Use SoftFloat for all FP operations (IEEE 754 compliance) +:: +:: 2. Define coefficients as hex literals: +:: static const double S1 = 0xbfc5555555555549; // reinterpreted +:: Or use union-based bit casting. +:: +:: 3. Match operation order exactly: +:: // sindf for @rd +:: double z = x * x; +:: double r = S5 + z * S6; +:: r = S4 + z * r; +:: r = S3 + z * r; +:: r = S2 + z * r; +:: r = S1 + z * r; +:: return x + x * z * r; +:: +:: 4. Use volatile or explicit stores if compiler reorders. +:: +-- diff --git a/libmath/desk/tests/lib/math-chebyshev.hoon b/libmath/desk/tests/lib/math-chebyshev.hoon new file mode 100644 index 0000000..59588f5 --- /dev/null +++ b/libmath/desk/tests/lib/math-chebyshev.hoon @@ -0,0 +1,733 @@ +/+ *test +/+ *math-chebyshev +:: +:::: Tests for Chebyshev/fixed-polynomial math library +:: +^| +|_ $: atol=_.1e-3 :: absolute tolerance + rtol=_.1e-5 :: relative tolerance + == +:: +:: Helper functions +:: +++ expect-near-rs + |= [expected=@rs actual=@rs] ^- tang + =/ diff (abs:rs:^math-chebyshev (sub:^rs expected actual)) + ?: (lth:^rs diff atol) + ~ + :~ [%palm [": " ~ ~ ~] [leaf+"expected" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"actual" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"diff" "{}"]] + == +:: +++ expect-near-rd + |= [expected=@rd actual=@rd] ^- tang + =/ diff `@rd`(dis (sub:^rd expected actual) 0x7fff.ffff.ffff.ffff) + ?: (lth:^rd diff .~1e-10) + ~ + :~ [%palm [": " ~ ~ ~] [leaf+"expected" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"actual" "{}"]] + == +:: +++ expect-near-rh + |= [expected=@rh actual=@rh] ^- tang + :: Half precision: compare with ~1e-2 tolerance + =/ diff `@rh`(dis (sub:^rh expected actual) 0x7fff) + ?: (lth:^rh diff `@rh`0x2e66) :: ~0.01 + ~ + :~ [%palm [": " ~ ~ ~] [leaf+"expected" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"actual" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"diff" "{}"]] + == +:: +++ expect-near-rq + |= [expected=@rq actual=@rq] ^- tang + =/ diff `@rq`(dis (sub:^rq expected actual) 0x7fff.ffff.ffff.ffff.ffff.ffff.ffff.ffff) + ?: (lth:^rq diff .~~~1e-20) + ~ + :~ [%palm [": " ~ ~ ~] [leaf+"expected" "{}"]] + [%palm [": " ~ ~ ~] [leaf+"actual" "{}"]] + == +:: +:: ============================================================ +:: SINGLE PRECISION (@rs) TESTS +:: ============================================================ +:: +++ test-rs-sin ^- tang + ;: weld + %+ expect-near-rs .0 (sin:rs .0) + %+ expect-near-rs .0.84147096 (sin:rs .1) + %+ expect-near-rs .0.9092974 (sin:rs .2) + %+ expect-near-rs .-0.7568025 (sin:rs .4) + %+ expect-near-rs .0 (sin:rs .3.1415927) :: sin(pi) ~ 0 + == +:: +++ test-rs-cos ^- tang + ;: weld + %+ expect-near-rs .1 (cos:rs .0) + %+ expect-near-rs .0.5403023 (cos:rs .1) + %+ expect-near-rs .-0.4161468 (cos:rs .2) + %+ expect-near-rs .-0.6536436 (cos:rs .4) + %+ expect-near-rs .-1 (cos:rs .3.1415927) :: cos(pi) ~ -1 + == +:: +++ test-rs-tan ^- tang + ;: weld + %+ expect-near-rs .0 (tan:rs .0) + %+ expect-near-rs .1.5574077 (tan:rs .1) + == +:: +++ test-rs-exp ^- tang + ;: weld + %+ expect-near-rs .1 (exp:rs .0) + %+ expect-near-rs .2.7182817 (exp:rs .1) + %+ expect-near-rs .7.389056 (exp:rs .2) + %+ expect-near-rs .0.36787945 (exp:rs .-1) + == +:: +++ test-rs-log ^- tang + ;: weld + %+ expect-near-rs .0 (log:rs .1) + %+ expect-near-rs .0.6931472 (log:rs .2) + %+ expect-near-rs .2.302585 (log:rs .10) + %+ expect-near-rs .-0.6931472 (log:rs .0.5) + == +:: +++ test-rs-sqrt ^- tang + ;: weld + %+ expect-near-rs .0 (sqt:rs .0) + %+ expect-near-rs .1 (sqt:rs .1) + %+ expect-near-rs .1.4142135 (sqt:rs .2) + %+ expect-near-rs .2 (sqt:rs .4) + %+ expect-near-rs .3 (sqt:rs .9) + %+ expect-near-rs .10 (sqt:rs .100) + == +:: +++ test-rs-atan ^- tang + ;: weld + %+ expect-near-rs .0 (atan:rs .0) + %+ expect-near-rs .0.7853982 (atan:rs .1) :: atan(1) = pi/4 + %+ expect-near-rs .1.1071488 (atan:rs .2) + %+ expect-near-rs .-0.7853982 (atan:rs .-1) + == +:: +++ test-rs-asin ^- tang + ;: weld + %+ expect-near-rs .0 (asin:rs .0) + %+ expect-near-rs .0.5235988 (asin:rs .0.5) :: pi/6 + %+ expect-near-rs .1.5707964 (asin:rs .1) :: pi/2 + %+ expect-near-rs .-1.5707964 (asin:rs .-1) + == +:: +++ test-rs-acos ^- tang + ;: weld + %+ expect-near-rs .1.5707964 (acos:rs .0) :: pi/2 + %+ expect-near-rs .0 (acos:rs .1) + %+ expect-near-rs .3.1415927 (acos:rs .-1) :: pi + == +:: +++ test-rs-pow-n ^- tang + ;: weld + %+ expect-near-rs .1 (pow-n:rs .2 .0) + %+ expect-near-rs .2 (pow-n:rs .2 .1) + %+ expect-near-rs .4 (pow-n:rs .2 .2) + %+ expect-near-rs .8 (pow-n:rs .2 .3) + %+ expect-near-rs .27 (pow-n:rs .3 .3) + == +:: +++ test-rs-pow ^- tang + ;: weld + %+ expect-near-rs .1 (pow:rs .2 .0) + %+ expect-near-rs .8 (pow:rs .2 .3) + %+ expect-near-rs .2.828427 (pow:rs .2 .1.5) :: 2^1.5 = 2*sqrt(2) + %+ expect-near-rs .0.5 (pow:rs .2 .-1) + == +:: +++ test-rs-cheb ^- tang + :: Chebyshev polynomials: T_0=1, T_1=x, T_2=2x^2-1, T_3=4x^3-3x + ;: weld + %+ expect-near-rs .1 ((cheb:rs 0) .0.5) + %+ expect-near-rs .0.5 ((cheb:rs 1) .0.5) + %+ expect-near-rs .-0.5 ((cheb:rs 2) .0.5) :: 2*(0.5)^2 - 1 = -0.5 + %+ expect-near-rs .-1 ((cheb:rs 3) .0.5) :: 4*(0.5)^3 - 3*0.5 = -1 + == +:: +:: ============================================================ +:: DOUBLE PRECISION (@rd) TESTS +:: ============================================================ +:: +++ test-rd-sin ^- tang + ;: weld + %+ expect-near-rd .~0 (sin:rd .~0) + %+ expect-near-rd .~0.8414709848078965 (sin:rd .~1) + %+ expect-near-rd .~0.9092974268256817 (sin:rd .~2) + == +:: +++ test-rd-cos ^- tang + ;: weld + %+ expect-near-rd .~1 (cos:rd .~0) + %+ expect-near-rd .~0.5403023058681398 (cos:rd .~1) + == +:: +++ test-rd-exp ^- tang + ;: weld + %+ expect-near-rd .~1 (exp:rd .~0) + %+ expect-near-rd .~2.718281828459045 (exp:rd .~1) + == +:: +++ test-rd-log ^- tang + ;: weld + %+ expect-near-rd .~0 (log:rd .~1) + %+ expect-near-rd .~0.6931471805599453 (log:rd .~2) + %+ expect-near-rd .~2.302585092994046 (log:rd .~10) + == +:: +++ test-rd-sqrt ^- tang + ;: weld + %+ expect-near-rd .~1 (sqt:rd .~1) + %+ expect-near-rd .~1.4142135623730951 (sqt:rd .~2) + %+ expect-near-rd .~2 (sqt:rd .~4) + == +:: +++ test-rd-atan ^- tang + ;: weld + %+ expect-near-rd .~0 (atan:rd .~0) + %+ expect-near-rd .~0.7853981633974483 (atan:rd .~1) :: pi/4 + == +:: +++ test-rd-pow-n ^- tang + ;: weld + %+ expect-near-rd .~8 (pow-n:rd .~2 .~3) + %+ expect-near-rd .~27 (pow-n:rd .~3 .~3) + == +:: +++ test-rd-cheb ^- tang + ;: weld + %+ expect-near-rd .~1 ((cheb:rd 0) .~0.5) + %+ expect-near-rd .~0.5 ((cheb:rd 1) .~0.5) + %+ expect-near-rd .~-0.5 ((cheb:rd 2) .~0.5) + %+ expect-near-rd .~-1 ((cheb:rd 3) .~0.5) + == +:: +:: ============================================================ +:: HALF PRECISION (@rh) TESTS +:: ============================================================ +:: +:: Note: @rh uses hex literals. Common values: +:: 0x0000 = 0.0, 0x3c00 = 1.0, 0x4000 = 2.0 +:: 0x4200 = 3.0, 0x4400 = 4.0, 0x4900 = 10.0 +:: 0x3800 = 0.5, 0xbc00 = -1.0 +:: +++ test-rh-sin ^- tang + ;: weld + :: sin(0) = 0 + %+ expect-near-rh `@rh`0x0 (sin:rh `@rh`0x0) + :: sin(1) ~ 0.8415 + %+ expect-near-rh `@rh`0x3abb (sin:rh `@rh`0x3c00) + :: sin(2) ~ 0.9093 + %+ expect-near-rh `@rh`0x3b46 (sin:rh `@rh`0x4000) + == +:: +++ test-rh-cos ^- tang + ;: weld + :: cos(0) = 1 + %+ expect-near-rh `@rh`0x3c00 (cos:rh `@rh`0x0) + :: cos(1) ~ 0.5403 + %+ expect-near-rh `@rh`0x3854 (cos:rh `@rh`0x3c00) + == +:: +++ test-rh-exp ^- tang + ;: weld + :: exp(0) = 1 + %+ expect-near-rh `@rh`0x3c00 (exp:rh `@rh`0x0) + :: exp(1) ~ 2.718 + %+ expect-near-rh `@rh`0x4170 (exp:rh `@rh`0x3c00) + :: exp(2) ~ 7.389 + %+ expect-near-rh `@rh`0x4763 (exp:rh `@rh`0x4000) + == +:: +++ test-rh-log ^- tang + ;: weld + :: log(1) = 0 + %+ expect-near-rh `@rh`0x0 (log:rh `@rh`0x3c00) + :: log(2) ~ 0.6931 + %+ expect-near-rh `@rh`0x398c (log:rh `@rh`0x4000) + :: log(10) ~ 2.303 + %+ expect-near-rh `@rh`0x409b (log:rh `@rh`0x4900) + == +:: +++ test-rh-sqrt ^- tang + ;: weld + :: sqrt(1) = 1 + %+ expect-near-rh `@rh`0x3c00 (sqt:rh `@rh`0x3c00) + :: sqrt(2) ~ 1.4142 + %+ expect-near-rh `@rh`0x3da8 (sqt:rh `@rh`0x4000) + :: sqrt(4) = 2 + %+ expect-near-rh `@rh`0x4000 (sqt:rh `@rh`0x4400) + :: sqrt(9) = 3 + %+ expect-near-rh `@rh`0x4200 (sqt:rh `@rh`0x4880) + == +:: +++ test-rh-atan ^- tang + ;: weld + :: atan(0) = 0 + %+ expect-near-rh `@rh`0x0 (atan:rh `@rh`0x0) + :: atan(1) ~ 0.7854 (pi/4) + %+ expect-near-rh `@rh`0x3a48 (atan:rh `@rh`0x3c00) + == +:: +++ test-rh-pow-n ^- tang + ;: weld + :: 2^0 = 1 + %+ expect-near-rh `@rh`0x3c00 (pow-n:rh `@rh`0x4000 `@rh`0x0) + :: 2^3 = 8 + %+ expect-near-rh `@rh`0x4800 (pow-n:rh `@rh`0x4000 `@rh`0x4200) + :: 3^3 = 27 + %+ expect-near-rh `@rh`0x4ec0 (pow-n:rh `@rh`0x4200 `@rh`0x4200) + == +:: +++ test-rh-cheb ^- tang + :: Chebyshev polynomials: T_0=1, T_1=x, T_2=2x^2-1, T_3=4x^3-3x + :: At x=0.5: T_0=1, T_1=0.5, T_2=-0.5, T_3=-1 + ;: weld + %+ expect-near-rh `@rh`0x3c00 ((cheb:rh 0) `@rh`0x3800) :: T_0(0.5) = 1 + %+ expect-near-rh `@rh`0x3800 ((cheb:rh 1) `@rh`0x3800) :: T_1(0.5) = 0.5 + %+ expect-near-rh `@rh`0xb800 ((cheb:rh 2) `@rh`0x3800) :: T_2(0.5) = -0.5 + %+ expect-near-rh `@rh`0xbc00 ((cheb:rh 3) `@rh`0x3800) :: T_3(0.5) = -1 + == +:: +:: ============================================================ +:: QUAD PRECISION (@rq) TESTS +:: ============================================================ +:: +++ test-rq-sin ^- tang + ;: weld + %+ expect-near-rq .~~~0 (sin:rq .~~~0) + %+ expect-near-rq .~~~0.8414709848078965 (sin:rq .~~~1) + %+ expect-near-rq .~~~0.9092974268256817 (sin:rq .~~~2) + == +:: +++ test-rq-cos ^- tang + ;: weld + %+ expect-near-rq .~~~1 (cos:rq .~~~0) + %+ expect-near-rq .~~~0.5403023058681398 (cos:rq .~~~1) + == +:: +++ test-rq-exp ^- tang + ;: weld + %+ expect-near-rq .~~~1 (exp:rq .~~~0) + %+ expect-near-rq .~~~2.718281828459045 (exp:rq .~~~1) + == +:: +++ test-rq-log ^- tang + ;: weld + %+ expect-near-rq .~~~0 (log:rq .~~~1) + %+ expect-near-rq .~~~0.6931471805599453 (log:rq .~~~2) + %+ expect-near-rq .~~~2.302585092994046 (log:rq .~~~10) + == +:: +++ test-rq-sqrt ^- tang + ;: weld + %+ expect-near-rq .~~~1 (sqt:rq .~~~1) + %+ expect-near-rq .~~~1.4142135623730951 (sqt:rq .~~~2) + %+ expect-near-rq .~~~2 (sqt:rq .~~~4) + == +:: +++ test-rq-atan ^- tang + ;: weld + %+ expect-near-rq .~~~0 (atan:rq .~~~0) + %+ expect-near-rq .~~~0.7853981633974483 (atan:rq .~~~1) :: pi/4 + == +:: +++ test-rq-pow-n ^- tang + ;: weld + %+ expect-near-rq .~~~8 (pow-n:rq .~~~2 .~~~3) + %+ expect-near-rq .~~~27 (pow-n:rq .~~~3 .~~~3) + == +:: +++ test-rq-cheb ^- tang + ;: weld + %+ expect-near-rq .~~~1 ((cheb:rq 0) .~~~0.5) + %+ expect-near-rq .~~~0.5 ((cheb:rq 1) .~~~0.5) + %+ expect-near-rq .~~~-0.5 ((cheb:rq 2) .~~~0.5) + %+ expect-near-rq .~~~-1 ((cheb:rq 3) .~~~0.5) + == +:: +:: ============================================================ +:: IDENTITY TESTS - verify mathematical relationships +:: ============================================================ +:: +++ test-rs-sin-cos-identity ^- tang + :: sin^2(x) + cos^2(x) = 1 + =/ x .1.234 + =/ s (sin:rs x) + =/ c (cos:rs x) + =/ sum (add:^rs (mul:^rs s s) (mul:^rs c c)) + (expect-near-rs .1 sum) +:: +++ test-rs-exp-log-identity ^- tang + :: exp(log(x)) = x for x > 0 + =/ x .2.5 + =/ result (exp:rs (log:rs x)) + (expect-near-rs x result) +:: +++ test-rs-log-exp-identity ^- tang + :: log(exp(x)) = x + =/ x .1.5 + =/ result (log:rs (exp:rs x)) + (expect-near-rs x result) +:: +++ test-rs-sqrt-square-identity ^- tang + :: sqrt(x)^2 = x for x >= 0 + =/ x .7.5 + =/ s (sqt:rs x) + =/ result (mul:^rs s s) + (expect-near-rs x result) +:: +++ test-rs-atan-tan-identity ^- tang + :: atan(tan(x)) = x for |x| < pi/2 + =/ x .0.5 + =/ result (atan:rs (tan:rs x)) + (expect-near-rs x result) +:: +:: @rd identity tests +:: +++ test-rd-sin-cos-identity ^- tang + :: sin^2(x) + cos^2(x) = 1 + =/ x .~1.234 + =/ s (sin:rd x) + =/ c (cos:rd x) + =/ sum (add:^rd (mul:^rd s s) (mul:^rd c c)) + (expect-near-rd .~1 sum) +:: +++ test-rd-exp-log-identity ^- tang + :: exp(log(x)) = x for x > 0 + =/ x .~2.5 + =/ result (exp:rd (log:rd x)) + (expect-near-rd x result) +:: +++ test-rd-sqrt-square-identity ^- tang + :: sqrt(x)^2 = x for x >= 0 + =/ x .~7.5 + =/ s (sqt:rd x) + =/ result (mul:^rd s s) + (expect-near-rd x result) +:: +:: @rh identity tests (wider tolerance due to half precision) +:: +++ test-rh-sin-cos-identity ^- tang + :: sin^2(x) + cos^2(x) = 1 + :: x = 1.25 (0x3d00) + =/ x `@rh`0x3d00 + =/ one `@rh`0x3c00 + =/ s (sin:rh x) + =/ c (cos:rh x) + =/ sum (add:^rh (mul:^rh s s) (mul:^rh c c)) + (expect-near-rh one sum) +:: +++ test-rh-exp-log-identity ^- tang + :: exp(log(x)) = x for x > 0 + :: x = 2.5 (0x4100) + =/ x `@rh`0x4100 + =/ result (exp:rh (log:rh x)) + (expect-near-rh x result) +:: +++ test-rh-sqrt-square-identity ^- tang + :: sqrt(x)^2 = x for x >= 0 + :: x = 4.0 (0x4400) + =/ x `@rh`0x4400 + =/ s (sqt:rh x) + =/ result (mul:^rh s s) + (expect-near-rh x result) +:: +:: @rq identity tests +:: +++ test-rq-sin-cos-identity ^- tang + :: sin^2(x) + cos^2(x) = 1 + =/ x .~~~1.234 + =/ s (sin:rq x) + =/ c (cos:rq x) + =/ sum (add:^rq (mul:^rq s s) (mul:^rq c c)) + (expect-near-rq .~~~1 sum) +:: +++ test-rq-exp-log-identity ^- tang + :: exp(log(x)) = x for x > 0 + =/ x .~~~2.5 + =/ result (exp:rq (log:rq x)) + (expect-near-rq x result) +:: +++ test-rq-sqrt-square-identity ^- tang + :: sqrt(x)^2 = x for x >= 0 + =/ x .~~~7.5 + =/ s (sqt:rq x) + =/ result (mul:^rq s s) + (expect-near-rq x result) +:: +++ test-rq-atan-tan-identity ^- tang + :: atan(tan(x)) = x for |x| < pi/2 + =/ x .~~~0.5 + =/ result (atan:rq (tan:rq x)) + (expect-near-rq x result) +:: +:: ============================================================ +:: HYPERBOLIC FUNCTION TESTS +:: ============================================================ +:: +++ test-rs-sinh ^- tang + ;: weld + %+ expect-near-rs .0 (sinh:rs .0) + %+ expect-near-rs .1.1752012 (sinh:rs .1) :: sinh(1) + %+ expect-near-rs .3.6268604 (sinh:rs .2) :: sinh(2) + == +:: +++ test-rs-cosh ^- tang + ;: weld + %+ expect-near-rs .1 (cosh:rs .0) + %+ expect-near-rs .1.5430806 (cosh:rs .1) :: cosh(1) + %+ expect-near-rs .3.7621956 (cosh:rs .2) :: cosh(2) + == +:: +++ test-rs-tanh ^- tang + ;: weld + %+ expect-near-rs .0 (tanh:rs .0) + %+ expect-near-rs .0.7615942 (tanh:rs .1) :: tanh(1) + %+ expect-near-rs .0.9640276 (tanh:rs .2) :: tanh(2) + == +:: +++ test-rs-asinh ^- tang + ;: weld + %+ expect-near-rs .0 (asinh:rs .0) + %+ expect-near-rs .0.88137358 (asinh:rs .1) :: asinh(1) + %+ expect-near-rs .1.4436355 (asinh:rs .2) :: asinh(2) + == +:: +++ test-rs-acosh ^- tang + ;: weld + %+ expect-near-rs .0 (acosh:rs .1) + %+ expect-near-rs .1.3169579 (acosh:rs .2) :: acosh(2) + %+ expect-near-rs .1.7627472 (acosh:rs .3) :: acosh(3) + == +:: +++ test-rs-atanh ^- tang + ;: weld + %+ expect-near-rs .0 (atanh:rs .0) + %+ expect-near-rs .0.5493061 (atanh:rs .0.5) :: atanh(0.5) + == +:: +++ test-rd-sinh ^- tang + ;: weld + %+ expect-near-rd .~0 (sinh:rd .~0) + %+ expect-near-rd .~1.1752011936438014 (sinh:rd .~1) + == +:: +++ test-rd-cosh ^- tang + ;: weld + %+ expect-near-rd .~1 (cosh:rd .~0) + %+ expect-near-rd .~1.5430806348152437 (cosh:rd .~1) + == +:: +++ test-rd-tanh ^- tang + ;: weld + %+ expect-near-rd .~0 (tanh:rd .~0) + %+ expect-near-rd .~0.7615941559557649 (tanh:rd .~1) + == +:: +++ test-rq-sinh ^- tang + ;: weld + %+ expect-near-rq .~~~0 (sinh:rq .~~~0) + %+ expect-near-rq .~~~1.1752011936438014 (sinh:rq .~~~1) + == +:: +++ test-rq-cosh ^- tang + ;: weld + %+ expect-near-rq .~~~1 (cosh:rq .~~~0) + %+ expect-near-rq .~~~1.5430806348152437 (cosh:rq .~~~1) + == +:: +++ test-rq-tanh ^- tang + ;: weld + %+ expect-near-rq .~~~0 (tanh:rq .~~~0) + %+ expect-near-rq .~~~0.7615941559557649 (tanh:rq .~~~1) + == +:: +:: Hyperbolic identities +:: +++ test-rs-cosh-sinh-identity ^- tang + :: cosh^2(x) - sinh^2(x) = 1 + =/ x .1.5 + =/ sh (sinh:rs x) + =/ ch (cosh:rs x) + =/ result (sub:^rs (mul:^rs ch ch) (mul:^rs sh sh)) + (expect-near-rs .1 result) +:: +++ test-rd-cosh-sinh-identity ^- tang + :: cosh^2(x) - sinh^2(x) = 1 + =/ x .~1.5 + =/ sh (sinh:rd x) + =/ ch (cosh:rd x) + =/ result (sub:^rd (mul:^rd ch ch) (mul:^rd sh sh)) + (expect-near-rd .~1 result) +:: +++ test-rs-asinh-sinh-identity ^- tang + :: asinh(sinh(x)) = x + =/ x .1.5 + =/ result (asinh:rs (sinh:rs x)) + (expect-near-rs x result) +:: +++ test-rs-acosh-cosh-identity ^- tang + :: acosh(cosh(x)) = |x| for x > 0 + =/ x .1.5 + =/ result (acosh:rs (cosh:rs x)) + (expect-near-rs x result) +:: +++ test-rs-atanh-tanh-identity ^- tang + :: atanh(tanh(x)) = x for small x + =/ x .0.5 + =/ result (atanh:rs (tanh:rs x)) + (expect-near-rs x result) +:: +:: ============================================================ +:: SPECIAL FUNCTION TESTS +:: ============================================================ +:: +++ test-rs-erf ^- tang + ;: weld + %+ expect-near-rs .0 (erf:rs .0) + %+ expect-near-rs .0.8427008 (erf:rs .1) :: erf(1) + %+ expect-near-rs .0.9953223 (erf:rs .2) :: erf(2) + %+ expect-near-rs .-0.8427008 (erf:rs .-1) :: erf(-1) + == +:: +++ test-rd-erf ^- tang + ;: weld + %+ expect-near-rd .~0 (erf:rd .~0) + %+ expect-near-rd .~0.8427007929497149 (erf:rd .~1) + == +:: +++ test-rq-erf ^- tang + ;: weld + %+ expect-near-rq .~~~0 (erf:rq .~~~0) + %+ expect-near-rq .~~~0.8427007929497149 (erf:rq .~~~1) + == +:: +++ test-rs-erf-symmetry ^- tang + :: erf(-x) = -erf(x) + =/ x .1.5 + =/ ep (erf:rs x) + =/ em (erf:rs (sub:^rs .0 x)) + (expect-near-rs (sub:^rs .0 ep) em) +:: +:: ============================================================ +:: GAMMA FUNCTION TESTS +:: ============================================================ +:: +++ test-rs-gamma ^- tang + ;: weld + %+ expect-near-rs .1 (gamma:rs .1) :: Γ(1) = 1 + %+ expect-near-rs .1 (gamma:rs .2) :: Γ(2) = 1! = 1 + %+ expect-near-rs .2 (gamma:rs .3) :: Γ(3) = 2! = 2 + %+ expect-near-rs .6 (gamma:rs .4) :: Γ(4) = 3! = 6 + %+ expect-near-rs .24 (gamma:rs .5) :: Γ(5) = 4! = 24 + %+ expect-near-rs .1.7724539 (gamma:rs .1.5) :: Γ(1.5) = √π/2 + == +:: +++ test-rd-gamma ^- tang + ;: weld + %+ expect-near-rd .~1 (gamma:rd .~1) + %+ expect-near-rd .~1 (gamma:rd .~2) + %+ expect-near-rd .~2 (gamma:rd .~3) + %+ expect-near-rd .~6 (gamma:rd .~4) + %+ expect-near-rd .~24 (gamma:rd .~5) + %+ expect-near-rd .~1.7724538509055159 (gamma:rd .~1.5) + == +:: +++ test-rq-gamma ^- tang + ;: weld + %+ expect-near-rq .~~~1 (gamma:rq .~~~1) + %+ expect-near-rq .~~~1 (gamma:rq .~~~2) + %+ expect-near-rq .~~~2 (gamma:rq .~~~3) + %+ expect-near-rq .~~~6 (gamma:rq .~~~4) + %+ expect-near-rq .~~~1.7724538509055159 (gamma:rq .~~~1.5) + == +:: +++ test-rs-gamma-factorial-identity ^- tang + :: Γ(n+1) = n * Γ(n) + =/ n .3.5 + =/ gn (gamma:rs n) + =/ gnp1 (gamma:rs (add:^rs n .1)) + (expect-near-rs gnp1 (mul:^rs n gn)) +:: +:: ============================================================ +:: BESSEL FUNCTION TESTS +:: ============================================================ +:: +++ test-rs-j0 ^- tang + ;: weld + %+ expect-near-rs .1 (j0:rs .0) :: J0(0) = 1 + %+ expect-near-rs .0.7651977 (j0:rs .1) :: J0(1) + %+ expect-near-rs .0.22389078 (j0:rs .2) :: J0(2) + %+ expect-near-rs .-0.26005195 (j0:rs .3) :: J0(3) + == +:: +++ test-rs-j1 ^- tang + ;: weld + %+ expect-near-rs .0 (j1:rs .0) :: J1(0) = 0 + %+ expect-near-rs .0.44005059 (j1:rs .1) :: J1(1) + %+ expect-near-rs .0.5767248 (j1:rs .2) :: J1(2) + %+ expect-near-rs .0.33905896 (j1:rs .3) :: J1(3) + == +:: +++ test-rd-j0 ^- tang + ;: weld + %+ expect-near-rd .~1 (j0:rd .~0) + %+ expect-near-rd .~0.7651976865579666 (j0:rd .~1) + %+ expect-near-rd .~0.22389077914123567 (j0:rd .~2) + == +:: +++ test-rd-j1 ^- tang + ;: weld + %+ expect-near-rd .~0 (j1:rd .~0) + %+ expect-near-rd .~0.44005058574493355 (j1:rd .~1) + %+ expect-near-rd .~0.5767248077568734 (j1:rd .~2) + == +:: +++ test-rq-j0 ^- tang + ;: weld + %+ expect-near-rq .~~~1 (j0:rq .~~~0) + %+ expect-near-rq .~~~0.7651976865579666 (j0:rq .~~~1) + == +:: +++ test-rq-j1 ^- tang + ;: weld + %+ expect-near-rq .~~~0 (j1:rq .~~~0) + %+ expect-near-rq .~~~0.44005058574493355 (j1:rq .~~~1) + == +:: +++ test-rs-j0-j1-identity ^- tang + :: J1(x) = -J0'(x), approximated by (J0(x-h) - J0(x+h))/(2h) + :: Using small x for stability + =/ x .1 + =/ h .0.001 + =/ j0m (j0:rs (sub:^rs x h)) + =/ j0p (j0:rs (add:^rs x h)) + =/ deriv (div:^rs (sub:^rs j0m j0p) (mul:^rs .2 h)) + =/ j1x (j1:rs x) + :: Derivative should approximately equal J1 + (expect-near-rs deriv j1x) +:: +++ test-rs-y0 ^- tang + ;: weld + :: Y0(1) ~ 0.0882569642 + %+ expect-near-rs .0.088256965 (y0:rs .1) + :: Y0(2) ~ 0.5103756726 + %+ expect-near-rs .0.5103757 (y0:rs .2) + == +:: +++ test-rs-y1 ^- tang + ;: weld + :: Y1(1) ~ -0.7812128213 + %+ expect-near-rs .-0.78121282 (y1:rs .1) + :: Y1(2) ~ -0.1070324315 + %+ expect-near-rs .-0.10703243 (y1:rs .2) + == +--