From 22e45b1edad2664e5cd8cc109324e3296b6ce106 Mon Sep 17 00:00:00 2001 From: Kaida-Amethyst Date: Thu, 23 Jan 2025 16:54:28 +0800 Subject: [PATCH 1/3] better trigonometric functions' implementation --- double/double.mbti | 3 + double/moon.pkg.json | 2 + double/trig_js.mbt | 91 +++++ double/trig_nonjs.mbt | 770 ++++++++++++++++++++++++++++++++++++ double/trig_test.mbt | 94 +++++ math/trigonometric.mbt | 167 +------- math/trigonometric_test.mbt | 105 ----- 7 files changed, 963 insertions(+), 269 deletions(-) create mode 100644 double/trig_js.mbt create mode 100644 double/trig_nonjs.mbt create mode 100644 double/trig_test.mbt diff --git a/double/double.mbti b/double/double.mbti index 5bbda7491..bfd8c8688 100644 --- a/double/double.mbti +++ b/double/double.mbti @@ -19,6 +19,7 @@ let not_a_number : Double impl Double { abs(Double) -> Double ceil(Double) -> Double + cos(Double) -> Double exp(Double) -> Double floor(Double) -> Double from_int(Int) -> Double @@ -37,6 +38,8 @@ impl Double { pow(Double, Double) -> Double round(Double) -> Double signum(Double) -> Double + sin(Double) -> Double + tan(Double) -> Double to_be_bytes(Double) -> Bytes to_le_bytes(Double) -> Bytes to_string(Double) -> String diff --git a/double/moon.pkg.json b/double/moon.pkg.json index 34b0e672d..653cc2454 100644 --- a/double/moon.pkg.json +++ b/double/moon.pkg.json @@ -14,6 +14,8 @@ "mod_nonjs.mbt": ["not", "js"], "pow_js.mbt": ["js"], "pow_nonjs.mbt": ["not", "js"], + "trig_js.mbt" : ["js"], + "trig_nonjs.mbt" : ["not", "js"], "round_js.mbt": ["js"], "round_wasm.mbt": ["wasm", "wasm-gc"], "round.mbt": ["not", "js", "wasm", "wasm-gc"], diff --git a/double/trig_js.mbt b/double/trig_js.mbt new file mode 100644 index 000000000..601530942 --- /dev/null +++ b/double/trig_js.mbt @@ -0,0 +1,91 @@ +// Copyright 2025 International Digital Economy Academy +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +///| +/// Calculates the sine of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the sine. +/// +/// Returns the sine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "sin" { +/// inspect!(0.0.sin(), content="0") +/// inspect!(1.570796326794897.sin(), content="1") // pi / 2 +/// inspect!(2.0.sin(), content="0.9092974268256817") +/// inspect!(-5.0.sin(), content="0.9589242746631385") +/// inspect!(31415926535897.9323846.sin(), content="0.0012091232715481885") +/// inspect!(@double.not_a_number.sin(), content="NaN") +/// inspect!(@double.infinity.sin(), content="NaN") +/// inspect!(@double.neg_infinity.sin(), content="NaN") +/// } +/// ``` +pub fn sin(self : Double) -> Double = "Math" "sin" + +///| +/// Calculates the cosine of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the cosine. +/// +/// Returns the cosine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "cos" { +/// inspect!(0.0.cos(), content="1") +/// inspect!(2.5.cos(), content="-0.8011436155469337") +/// inspect!((-3.141592653589793).cos(), content="-1") // -pi +/// inspect!((-5.0).cos(), content="0.28366218546322625") +/// inspect!(31415926535897.9323846.cos(), content="0.9999992690101899") +/// inspect!(@double.not_a_number.cos(), content="NaN") +/// inspect!(@double.infinity.cos(), content="NaN") +/// inspect!(@double.neg_infinity.cos(), content="NaN") +/// } +/// ``` +pub fn cos(self : Double) -> Double = "Math" "cos" + +///| +/// Calculates the tangent of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the tangent. +/// +/// Returns the tangent of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "tan" { +/// inspect!(0.0.tan(), content="0") +/// inspect!(0.7853981633974483.tan(), content="0.9999999999999999") +/// inspect!(4.0.tan(), content="1.1578212823495777") +/// inspect!(5.0.tan(), content="-3.380515006246586") +/// inspect!(31415926535897.9323846.tan(), content="0.0012091241554056254") +/// inspect!(@double.not_a_number.tan(), content="NaN") +/// inspect!(@double.infinity.tan(), content="NaN") +/// inspect!(@double.neg_infinity.tan(), content="NaN") +/// } +/// ``` +pub fn tan(self : Double) -> Double = "Math" "tan" diff --git a/double/trig_nonjs.mbt b/double/trig_nonjs.mbt new file mode 100644 index 000000000..23329b849 --- /dev/null +++ b/double/trig_nonjs.mbt @@ -0,0 +1,770 @@ +// Copyright 2025 International Digital Economy Academy +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +///| +let two_over_pi : Array[Int] = [ + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, + 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, + 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, + 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, + 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, + 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, + 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, + 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, + 0x60E27B, 0xC08C6B, +] + +///| +let pi_over_2 : Array[Double] = [ + 1.57079625129699707031e+00, // 0x3FF921FB, 0x40000000 */ + 7.54978941586159635335e-08, // 0x3E74442D, 0x00000000 */ + 5.39030252995776476554e-15, // 0x3CF84698, 0x80000000 */ + 3.28200341580791294123e-22, // 0x3B78CC51, 0x60000000 */ + 1.27065575308067607349e-29, // 0x39F01B83, 0x80000000 */ + 1.22933308981111328932e-36, // 0x387A2520, 0x40000000 */ + 2.73370053816464559624e-44, // 0x36E38222, 0x80000000 */ + 2.16741683877804819444e-51, // 0x3569F31D, 0x00000000 */ +] + +///| +let npio2_hw : Array[Int] = [ + 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, + 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, + 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 0x403DD85A, 0x403F6A7A, 0x40407E4C, + 0x4041475C, 0x4042106C, 0x4042D97C, 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, + 0x4046C6CB, 0x40478FDB, 0x404858EB, 0x404921FB, +] + +///| +const PI_OVER_4 : Double = 0.785398163397448309616 + +///| +const PIO2_1 : Double = 1.5707963267341256e+00 // 0x3FF921FB, 0x54400000 */ + +///| +const PIO2_1T : Double = 6.0771005065061922e-11 // 0x3DD0B461, 0x1A600000 */ + +///| +const PIO2_2 : Double = 6.0771005063039660e-11 // 0x3DD0B461, 0x1A600000 */ + +///| +const PIO2_2T : Double = 2.0222662487959506e-21 // 0x3BA3198A, 0x2E037073 */ + +///| +const PIO2_3 : Double = 2.0222662487111665e-21 // 0x3BA3198A, 0x2E037073 */ + +///| +const PIO2_3T : Double = 8.4784276603688996e-32 // 0x397B839A, 0x252049C1 */ + +///| +const INV_PIO2 : Double = 6.3661977236758138e-01 // 0x3FE45F30, 0x6DC9C883 */ + +///| +const HALF : Double = 0.5 + +///| +const TWO24 : Double = 16777216.0 // 0x41700000, 0x00000000 */ + +///| +fn set_low_word(d : Double, v : UInt) -> Double { + let bits : UInt64 = d.reinterpret_as_uint64() + let bits = bits & 0xFFFF_FFFF_0000_0000 + let bits = bits | v.to_uint64() + bits.reinterpret_as_double() +} + +///| +fn set_high_word(d : Double, v : UInt) -> Double { + let bits : UInt64 = d.reinterpret_as_uint64() + let bits = bits & 0x0000_0000_FFFF_FFFF + let bits = bits | (v.to_uint64() << 32) + bits.reinterpret_as_double() +} + +///| +fn get_high_word(x : Double) -> UInt { + (x.reinterpret_as_uint64() >> 32).to_uint() +} + +///| +fn get_low_word(x : Double) -> UInt { + x.reinterpret_as_uint64().to_uint() +} + +///| +fn rem_pio2(x : Double, y : Array[Double]) -> Int { + let hx = get_high_word(x).reinterpret_as_int() + let ix : Int = hx & 0x7fffffff + let mut z = 0.0 + if ix <= 0x3fe921fb { + // |x| <= pi/4, no reduction needed + y[0] = x + y[1] = 0.0 + return 0 + } + if ix < 0x4002d97c { + // |x| < 3pi/4, special case with n = +-1 + if hx > 0 { + z = x - PIO2_1 + if ix != 0x3ff921fb { + // 33+53 bit pi is good enough + y[0] = z - PIO2_1T + y[1] = z - y[0] - PIO2_1T + } else { + // Near pi/2, use 33+33+53 bit pi + z = z - PIO2_2 + y[0] = z - PIO2_2T + y[1] = z - y[0] - PIO2_2T + } + return 1 + } else { + // Negative x + z = x + PIO2_1 + if ix != 0x3ff921fb { + // 33+53 bit pi is good enough + y[0] = z + PIO2_1T + y[1] = z - y[0] + PIO2_1T + } else { + // Near pi/2, use 33+33+53 bit pi + let z = z + PIO2_2 + y[0] = z + PIO2_2T + y[1] = z - y[0] + PIO2_2T + } + return -1 + } + } + if ix <= 0x413921fb { + // |x| <= 2^19 * (pi/2), medium size + let t = x.abs() + let n = (t * INV_PIO2 + HALF).to_int() + let fn_ = n.to_double() + let mut r = t - fn_ * PIO2_1 + let mut w = fn_ * PIO2_1T + if n < 32 && ix != npio2_hw[n - 1] { + y[0] = r - w + } else { + let j = ix >> 20 + y[0] = r - w + let i = j - ((get_high_word(y[0]) >> 20).reinterpret_as_int() & 0x7ff) + if i > 16 { + // 2nd iteration needed, good to 118 bits + let t = r + w = fn_ * PIO2_2 + r = t - w + w = fn_ * PIO2_2T - (t - r - w) + y[0] = r - w + let i = j - ((get_high_word(y[0]) >> 20).reinterpret_as_int() & 0x7ff) + if i > 49 { + // 3rd iteration needed, 151 bits accuracy + let t = r + w = fn_ * PIO2_3 + r = t - w + w = fn_ * PIO2_3T - (t - r - w) + y[0] = r - w + } + } + } + y[1] = r - y[0] - w + if hx > 0 { + return n + } else { + y[0] = -y[0] + y[1] = -y[1] + return -n + } + } + + // All other (large) arguments + if ix >= 0x7ff00000 { + // x is inf or NaN + y[0] = x - x + y[1] = y[0] + return 0 + } + + // Set z = scalbn(|x|, ilogb(x) - 23) + z = set_low_word(z, get_low_word(x)) + let e0 = (ix >> 20) - 1046 // e0 = ilogb(z) - 23 + z = set_high_word(z, (ix - (e0 << 20)).reinterpret_as_uint()) + let tx = [0.0, 0.0, 0.0] + for i = 0; i < 2; i = i + 1 { + tx[i] = z.to_int().to_double() + z = (z - tx[i]) * TWO24 + } + tx[2] = z + let mut nx = 3 + while tx[nx - 1] == 0.0 { + nx -= 1 + } + let n = __kernel_rem_pio2(tx, y, e0, nx, 2) + if hx > 0 { + n + } else { + y[0] = -y[0] + y[1] = -y[1] + -n + } +} + +///| +fn __kernel_rem_pio2( + x : Array[Double], + y : Array[Double], + e0 : Int, + nx : Int, + prec : Int +) -> Int { + let init_jk = [2, 3, 4, 6] + let two24 : Double = 16777216.0 // 0x41700000, 0x00000000 + let twon24 : Double = 5.96046447753906250000e-08 // 0x3E700000, 0x00000000 + + // Declare variables at the beginning, as in the original C code + let mut jz : Int = 0 + let mut jx : Int = 0 + let mut jv : Int = 0 + let mut jp : Int = 0 + let mut jk : Int = 0 + let mut carry : Int = 0 + let mut n : Int = 0 + let iq : Array[Int] = Array::make(20, 0) + let mut i : Int = 0 + let mut j : Int = 0 + let mut k : Int = 0 + let mut m : Int = 0 + let mut q0 : Int = 0 + let mut ih : Int = 0 + let mut z : Double = 0 + let mut fw : Double = 0 + let f : Array[Double] = Array::make(20, 0.0) + let fq : Array[Double] = Array::make(20, 0.0) + let q : Array[Double] = Array::make(20, 0.0) + + // initialize jk + jk = init_jk[prec] + jp = jk + + // determine jx, jv, q0, note that 3 > q0 + jx = nx - 1 + jv = (e0 - 3) / 24 + if jv < 0 { + jv = 0 + } + q0 = e0 - 24 * (jv + 1) + + // set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] + j = jv - jx + m = jx + jk + i = 0 + while i <= m { + f[i] = if j < 0 { 0.0 } else { two_over_pi[j].to_double() } + i += 1 + j += 1 + } + + // compute q[0], q[1], ... q[jk] + i = 0 + while i <= jk { + j = 0 + fw = 0.0 + while j <= jx { + fw += x[j] * f[jx + i - j] + j += 1 + } + q[i] = fw + i += 1 + } + jz = jk + + // Use a loop to replace the `goto` logic + let mut recompute = true + while recompute { + recompute = false + + // distill q[] into iq[] reversingly + i = 0 + j = jz + z = q[jz] + while j > 0 { + fw = (twon24 * z).floor() + iq[i] = (z - two24 * fw).to_int() + z = q[j - 1] + fw + i += 1 + j -= 1 + } + + // compute n + z = scalbn(z, q0) // actual value of z + z -= 8.0 * (z * 0.125).floor() // trim off integer >= 8 + n = z.to_int() + z -= n.to_double() + ih = 0 + if q0 > 0 { + // need iq[jz-1] to determine n + i = iq[jz - 1] >> (24 - q0) + n += i + iq[jz - 1] -= i << (24 - q0) + ih = iq[jz - 1] >> (23 - q0) + } else if q0 == 0 { + ih = iq[jz - 1] >> 23 + } else if z >= 0.5 { + ih = 2 + } + if ih > 0 { + // q > 0.5 + n += 1 + carry = 0 + i = 0 + while i < jz { + j = iq[i] + if carry == 0 { + if j != 0 { + carry = 1 + iq[i] = 0x1000000 - j + } + } else { + iq[i] = 0xffffff - j + } + i += 1 + } + if q0 > 0 { + // rare case: chance is 1 in 12 + match q0 { + 1 => iq[jz - 1] = iq[jz - 1] & 0x7fffff + 2 => iq[jz - 1] = iq[jz - 1] & 0x3fffff + _ => () + } + } + if ih == 2 { + z = 1.0 - z + if carry != 0 { + z -= scalbn(1.0, q0) + } + } + } + + // check if recomputation is needed + if z == 0.0 { + j = 0 + i = jz - 1 + while i >= jk { + j = j | iq[i] + i -= 1 + } + if j == 0 { + // need recomputation + k = 1 + while iq[jk - k] == 0 { + k += 1 + } + i = jz + 1 + while i <= jz + k { + f[jx + i] = two_over_pi[jv + i].to_double() + j = 0 + fw = 0.0 + while j <= jx { + fw += x[j] * f[jx + i - j] + j += 1 + } + q[i] = fw + i += 1 + } + jz += k + recompute = true // Continue to recompute + continue + } + } // Skip the rest of the loop and recompute + + // chop off zero terms + if z == 0.0 { + jz -= 1 + q0 -= 24 + while iq[jz] == 0 { + jz -= 1 + q0 -= 24 + } + } else { + // break z into 24-bit if necessary + z = scalbn(z, -q0) + if z >= two24 { + fw = (twon24 * z).floor() + iq[jz] = (z - two24 * fw).to_int() + jz += 1 + q0 += 24 + iq[jz] = fw.to_int() + } else { + iq[jz] = z.to_int() + } + } + + // convert integer "bit" chunk to floating-point value + fw = scalbn(1.0, q0) + i = jz + while i >= 0 { + q[i] = fw * iq[i].to_double() + fw *= twon24 + i -= 1 + } + + // compute PIo2[0,...,jp] * q[jz,...,0] + i = jz + while i >= 0 { + fw = 0.0 + k = 0 + while k <= jp && k <= jz - i { + fw += pi_over_2[k] * q[i + k] + k += 1 + } + fq[jz - i] = fw + i -= 1 + } + + // compress fq[] into y[] + match prec { + 0 => { + fw = 0.0 + i = jz + while i >= 0 { + fw += fq[i] + i -= 1 + } + y[0] = if ih == 0 { fw } else { -fw } + } + 1 | 2 => { + fw = 0.0 + i = jz + while i >= 0 { + fw += fq[i] + i -= 1 + } + y[0] = if ih == 0 { fw } else { -fw } + fw = fq[0] - fw + i = 1 + while i <= jz { + fw += fq[i] + i += 1 + } + y[1] = if ih == 0 { fw } else { -fw } + } + 3 => { + i = jz + while i > 0 { + fw = fq[i - 1] + fq[i] + fq[i] += fq[i - 1] - fw + fq[i - 1] = fw + i -= 1 + } + i = jz + while i > 1 { + fw = fq[i - 1] + fq[i] + fq[i] += fq[i - 1] - fw + fq[i - 1] = fw + i -= 1 + } + fw = 0.0 + i = jz + while i >= 2 { + fw += fq[i] + i -= 1 + } + if ih == 0 { + y[0] = fq[0] + y[1] = fq[1] + y[2] = fw + } else { + y[0] = -fq[0] + y[1] = -fq[1] + y[2] = -fw + } + } + _ => () + } + } + n & 7 +} + +///| +fn __kernel_sin(x : Double, y : Double, iy : Int) -> Double { + let s1 = -1.66666666666666324348e-01 + let s2 = 8.33333333332248946124e-03 + let s3 = -1.98412698298579493134e-04 + let s4 = 2.75573137070700676789e-06 + let s5 = -2.50507602534068634195e-08 + let s6 = 1.58969099521155010221e-10 + let mut z = 0.0 + let mut r = 0.0 + let mut v = 0.0 + let ix = get_high_word(x) & 0x7fffffff + if ix < 0x3e400000 { + if x.to_int() == 0 { + return x + } + } + z = x * x + v = z * x + r = s2 + z * (s3 + z * (s4 + z * (s5 + z * s6))) + if iy == 0 { + x + v * (s1 + z * r) + } else { + x - (z * (0.5 * y - v * r) - y - v * s1) + } +} + +///| +fn __kernel_cos(x : Double, y : Double) -> Double { + let one = 1.00000000000000000000e+00 + let c1 = 4.16666666666666019037e-02 + let c2 = -1.38888888888741095749e-03 + let c3 = 2.48015872894767294178e-05 + let c4 = -2.75573143513906633035e-07 + let c5 = 2.08757232129817482790e-09 + let c6 = -1.13596475577881948265e-11 + let mut a = 0.0 + let mut hz = 0.0 + let mut z = 0.0 + let mut r = 0.0 + let mut qx = 0.0 + let ix = get_high_word(x) & 0x7fffffff + if ix < 0x3e400000 { + if x.to_int() == 0 { + return one + } + } + z = x * x + r = z * (c1 + z * (c2 + z * (c3 + z * (c4 + z * (c5 + z * c6))))) + if ix < 0x3fd33333 { + return one - (0.5 * z - (z * r - x * y)) + } else { + if ix > 0x3fe90000 { + qx = 0.28125 + } else { + qx = ((ix - 0x00200000).to_uint64() << 32).reinterpret_as_double() + } + hz = 0.5 * z - qx + a = one - qx + return a - (hz - (z * r - x * y)) + } +} + +///| +fn __kernal_tan(x : Double, y : Double, iy : Int) -> Double { + let one = 1.0 + let pio4 = 7.85398163397448278999e-01 + let pio4lo = 3.06161699786838301793e-17 + let mut x = x + let mut y = y + let mut z = 0.0 + let mut r = 0.0 + let mut v = 0.0 + let mut w = 0.0 + let mut s = 0.0 + let t = [ + 3.33333333333334091986e-01, // 3FD55555, 55555563 */ + 1.33333333333201242699e-01, // 3FC11111, 1110FE7A */ + 5.39682539762260521377e-02, // 3FABA1BA, 1BB341FE */ + 2.18694882948595424599e-02, // 3F9664F4, 8406D637 */ + 8.86323982359930005737e-03, // 3F8226E3, E96E8493 */ + 3.59207910759131235356e-03, // 3F6D6D22, C9560328 */ + 1.45620945432529025516e-03, // 3F57DBC8, FEE08315 */ + 5.88041240820264096874e-04, // 3F4344D8, F2F26501 */ + 2.46463134818469906812e-04, // 3F3026F7, 1A8D1068 */ + 7.81794442939557092300e-05, // 3F147E88, A03792A6 */ + 7.14072491382608190305e-05, // 3F12B80F, 32F0A7E9 */ + -1.85586374855275456654e-05, // BEF375CB, DB605373 */ + 2.59073051863633712884e-05, // 3EFB2A70, 74BF7AD4 */ + 1.00000000000000000000e+00, // 3FF00000, 00000000 (one) */ + 7.85398163397448278999e-01, // 3FE921FB, 54442D18 (pio4) */ + 3.06161699786838301793e-17, // 3C81A626, 33145C07 (pio4lo) */ + ] + let hx = get_high_word(x).reinterpret_as_int() + let ix = hx & 0x7fffffff + if ix < 0x3e300000 { + if x.to_int() == 0 { + if (ix | get_low_word(x).reinterpret_as_int() | (iy + 1)) == 0 { + return one / x.abs() + } else if iy == 1 { + return x + } else { + w = x + y + z = w + z = set_low_word(z, 0) + v = y - (z - x) + let a = -one / w + let mut t = a + t = set_low_word(t, 0) + s = one + t * z + return t + a * (s + t * v) + } + } + } + if ix >= 0x3fe59428 { + if hx < 0 { + x = -x + y = -y + } + z = pio4 - x + w = pio4lo - y + x = z + w + y = 0.0 + } + z = x * x + w = z * z + r = t[1] + w * (t[3] + w * (t[5] + w * (t[7] + w * (t[9] + w * t[11])))) + v = z * + (t[2] + w * (t[4] + w * (t[6] + w * (t[8] + w * (t[10] + w * t[12]))))) + s = z * x + r = y + z * (s * (r + v) + y) + r += t[0] * s + w = x + r + if ix >= 0x3fe59428 { + v = iy.to_double() + return (1 - ((hx >> 30) & 2)).to_double() * + (v - 2.0 * (x - (w * w / (w + v) - r))) + } + if iy == 1 { + w + } else { + z = w + z = set_low_word(z, 0) + v = r - (z - x) + let a = -1.0 / w + let mut t = a + t = set_low_word(t, 0) + s = 1.0 + t * z + t + a * (s + t * v) + } +} + +///| +/// Calculates the sine of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the sine. +/// +/// Returns the sine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "sin" { +/// inspect!(0.0.sin(), content="0") +/// inspect!(1.570796326794897.sin(), content="1") // pi / 2 +/// inspect!(2.0.sin(), content="0.9092974268256817") +/// inspect!(-5.0.sin(), content="0.9589242746631385") +/// inspect!(31415926535897.9323846.sin(), content="0.0012091232715481885") +/// inspect!(@double.not_a_number.sin(), content="NaN") +/// inspect!(@double.infinity.sin(), content="NaN") +/// inspect!(@double.neg_infinity.sin(), content="NaN") +/// } +/// ``` +pub fn sin(self : Double) -> Double { + if self.is_inf() || self.is_nan() { + return not_a_number + } + let y = [0.0, 0.0] + let z = 0.0 + if self.abs() <= PI_OVER_4 { + return __kernel_sin(self, z, 0) + } else { + let n = rem_pio2(self, y) + match n & 3 { + 0 => __kernel_sin(y[0], y[1], 1) + 1 => __kernel_cos(y[0], y[1]) + 2 => -__kernel_sin(y[0], y[1], 1) + _ => -__kernel_cos(y[0], y[1]) + } + } +} + +///| +/// Calculates the cosine of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the cosine. +/// +/// Returns the cosine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "cos" { +/// inspect!(0.0.cos(), content="1") +/// inspect!(2.5.cos(), content="-0.8011436155469337") +/// inspect!((-3.141592653589793).cos(), content="-1") // -pi +/// inspect!((-5.0).cos(), content="0.28366218546322625") +/// inspect!(31415926535897.9323846.cos(), content="0.9999992690101899") +/// inspect!(@double.not_a_number.cos(), content="NaN") +/// inspect!(@double.infinity.cos(), content="NaN") +/// inspect!(@double.neg_infinity.cos(), content="NaN") +/// } +/// ``` +pub fn cos(self : Double) -> Double { + if self.is_inf() || self.is_nan() { + return not_a_number + } + let y = [0.0, 0.0] + let z = 0.0 + if self.abs() <= PI_OVER_4 { + return __kernel_cos(self, z) + } else { + let n = rem_pio2(self, y) + match n & 3 { + 0 => __kernel_cos(y[0], y[1]) + 1 => -__kernel_sin(y[0], y[1], 1) + 2 => -__kernel_cos(y[0], y[1]) + _ => __kernel_sin(y[0], y[1], 1) + } + } +} + +///| +/// Calculates the tangent of a number in radians. Handles special cases and edge +/// conditions according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the tangent. +/// +/// Returns the tangent of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "tan" { +/// inspect!(0.0.tan(), content="0") +/// inspect!(0.7853981633974483.tan(), content="0.9999999999999999") +/// inspect!(4.0.tan(), content="1.1578212823495777") +/// inspect!(5.0.tan(), content="-3.380515006246586") +/// inspect!(31415926535897.9323846.tan(), content="0.0012091241554056254") +/// inspect!(@double.not_a_number.tan(), content="NaN") +/// inspect!(@double.infinity.tan(), content="NaN") +/// inspect!(@double.neg_infinity.tan(), content="NaN") +/// } +/// ``` +pub fn tan(self : Double) -> Double { + if self.is_inf() || self.is_nan() { + return not_a_number + } + let y = Array::make(2, 0.0) + let z = 0.0 + if self.abs() <= PI_OVER_4 { + __kernal_tan(self, z, 1) + } else { + let n = rem_pio2(self, y) + __kernal_tan(y[0], y[1], 1 - ((n & 1) << 1)) + } +} diff --git a/double/trig_test.mbt b/double/trig_test.mbt new file mode 100644 index 000000000..b4c727314 --- /dev/null +++ b/double/trig_test.mbt @@ -0,0 +1,94 @@ +// Copyright 2025 International Digital Economy Academy +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +test "trigonometric functions" { + // sin + inspect!(0.0.sin(), content="0") + inspect!(1.570796326794897.sin(), content="1") // pi / 2 + inspect!(2.0.sin(), content="0.9092974268256817") + inspect!(-5.0.sin(), content="0.9589242746631385") + inspect!(31415926535897.9323846.sin(), content="0.0012091232715481885") + inspect!(@double.not_a_number.sin(), content="NaN") + inspect!(@double.infinity.sin(), content="NaN") + inspect!(@double.neg_infinity.sin(), content="NaN") + + // cos + inspect!(0.0.cos(), content="1") + inspect!(2.5.cos(), content="-0.8011436155469337") + inspect!((-3.141592653589793).cos(), content="-1") // -pi + inspect!((-5.0).cos(), content="0.28366218546322625") + inspect!(31415926535897.9323846.cos(), content="0.9999992690101899") + inspect!(@double.not_a_number.cos(), content="NaN") + inspect!(@double.infinity.cos(), content="NaN") + inspect!(@double.neg_infinity.cos(), content="NaN") + + // tan + inspect!(0.0.tan(), content="0") + inspect!(0.7853981633974483.tan(), content="0.9999999999999999") + inspect!(4.0.tan(), content="1.1578212823495777") + inspect!(5.0.tan(), content="-3.380515006246586") + inspect!(31415926535897.9323846.tan(), content="0.0012091241554056254") + inspect!(@double.not_a_number.tan(), content="NaN") + inspect!(@double.infinity.tan(), content="NaN") + inspect!(@double.neg_infinity.tan(), content="NaN") + + // Consistency Test + let vals = [ + 1.1, 3.3, 4.4, 5.5, 11.22, 33.44, 55.66, 77.88, 99.99, -1.0, -2.1, -3.21, -4.321, + -5.4321, -6.54321, -7.654321, -8.7654321, -9.87654321, -10.987654321, 123456789.1, + 234567891.2, 345678912.3, 456789123.4, 567891234.5, 678912345.6, 789123456.7, + 891234567.8, 912345678.9, 123456789.01, + ] + + // Results get from JavaScript V8 + let sin_res = [ + 0.8912073600614354, -0.1577456941432482, -0.9516020738895161, -0.7055403255703919, + -0.9749220735246146, 0.8990168278954119, -0.7762325332088996, 0.612971751458885, + -0.5149633680424761, -0.8414709848078965, -0.8632093666488737, 0.06835400612104778, + 0.9243800935781268, 0.7519962469828865, -0.2571044147079568, -0.9801339430243416, + -0.6125999559240065, 0.43655436914142864, 0.999968637228749, 0.9991709051289052, + 0.7065650094118109, -0.6621352221182213, -0.36973782493802004, 0.6337066920505413, + -0.990764758712146, -0.8321035463409452, -0.4831904888511182, -0.6152185206743059, + 0.9914678206150951, + ] + let cos_res = [ + 0.4535961214255773, -0.9874797699088649, -0.30733286997841935, 0.70866977429126, + 0.22254651323816765, -0.43791408193945003, 0.6304467101889696, -0.7901048233705619, + 0.8572121846861191, 0.5403023058681398, -0.5048461045998576, -0.9976611297666176, + -0.38147272850951436, 0.6591673873331824, 0.966383629795993, 0.19833671806187747, + -0.7903931262364985, -0.8996778772335848, -0.007919883766735825, 0.040712434757508714, + -0.7076481381837216, -0.7493843790942362, 0.9291361260924591, 0.7735734150359359, + -0.13559200896093865, 0.5546203099119456, 0.8755152491440783, -0.7883566273072857, + 0.13035168078990628, + ] + let tan_res = [ + 1.9647596572486523, 0.15974574766003222, 3.0963237806497457, -0.995584052213885, + -4.380756451040238, -2.052952542457217, -1.2312421028832592, -0.7758106688223559, + -0.6007420067541824, -1.5574077246549023, 1.7098465429045073, -0.06851425206576686, + -2.4231878834166034, 1.1408274460077057, -0.2660479821685643, -4.941767478065043, + 0.775057291857959, -0.48523408231820436, -126.26051930568741, 24.542155512932695, + -0.9984693964225064, 0.883572223534375, -0.3979371962351479, 0.8191939895208301, + 7.306955375206259, -1.500312072006621, -0.5518927161160188, 0.7803809841437486, + 7.606099243270125, + ] + fn high_accuracy_test(expect : Double, actual : Double) -> Unit!Error { + assert_eq!((expect - actual).abs() < 1.0e-14, true) + } + + for i in 0.. Double { - if Double::is_inf(x) || Double::is_nan(x) { - @double.not_a_number - } else { - let mut sign = false - let mut x_ = x - if x_ < 0.0 { - x_ = -x_ - sign = true - } - let mut j = (x_ * pi4).to_int64() - if (j & 1L) == 1L { - j += 1L - } - let y = j.to_double() - j = j & 7L - if j > 3L { - sign = not(sign) - j -= 4L - } - let z = x_ - y * pi4a - y * pi4b - y * pi4c - let zz = z * z - (if j == 1L || j == 2L { - 1.0 - - 0.5 * zz + - zz * - zz * - ( - ( - (((cos_[0] * zz + cos_[1]) * zz + cos_[2]) * zz + cos_[3]) * zz + - cos_[4] - ) * - zz + - cos_[5] - ) - } else { - z + - z * - zz * - ( - ( - (((sin_[0] * zz + sin_[1]) * zz + sin_[2]) * zz + sin_[3]) * zz + - sin_[4] - ) * - zz + - sin_[5] - ) - }) * - (if sign { -1.0 } else { 1.0 }) - } + x.sin() } -///| -let cos_ = [ - -0x0.000000000C7D24D03C8, 0x0.00000008F74EBDA72B71, -0x0.0000049F93DFAB12FDC6, - 0x0.0001A01A019C844F3E07, -0x0.005B05B05B053E453C34, 0x0.0AAAAAAAAAAAA55A47E1, -] - ///| pub fn cos(x : Double) -> Double { - if Double::is_inf(x) || Double::is_nan(x) { - @double.not_a_number - } else { - let mut sign = false - let x_ = x.abs() - let mut j = (x_ * pi4).to_int64() - if (j & 1L) == 1L { - j += 1L - } - let y = j.to_double() - j = j & 7L - if j > 3L { - sign = not(sign) - j -= 4L - } - let z = x_ - y * pi4a - y * pi4b - y * pi4c - if j > 1L { - sign = not(sign) - } - let zz = z * z - (if j == 1L || j == 2L { - z + - z * - zz * - ( - ( - (((sin_[0] * zz + sin_[1]) * zz + sin_[2]) * zz + sin_[3]) * zz + - sin_[4] - ) * - zz + - sin_[5] - ) - } else { - 1.0 - - 0.5 * zz + - zz * - zz * - ( - ( - (((cos_[0] * zz + cos_[1]) * zz + cos_[2]) * zz + cos_[3]) * zz + - cos_[4] - ) * - zz + - cos_[5] - ) - }) * - (if sign { -1.0 } else { 1.0 }) - } + x.cos() } -///| -let tan_p = [ - -0x3325.B1A49E7E6F070EEA7F33, 0x1199EC.A5FC9DDCB62A659C2F5C, -0x111FEAD.3299175F331F9B486E14, -] - -///| -let tan_q = [ - 1.0, 0x3571.4BDD66CAE3B4F57353EC, -0x1427BC.582ABC95B9AA3C9963BF, 0x17D98FC.2EAD8EF5BE3D475ECC36, - -0x335FC07.97CB461BD8709633FE8F, -] - ///| pub fn tan(x : Double) -> Double { - if Double::is_inf(x) || Double::is_nan(x) { - @double.not_a_number - } else { - let x_ = x.abs() - let mut j = (x_ * pi4).to_int64() - if (j & 1L) == 1L { - j += 1L - } - let y = j.to_double() - let z = x_ - y * pi4a - y * pi4b - y * pi4c - let mut zz = z * z - zz = if zz > 1.0e-14 { - z + - z * - ( - zz * - ((tan_p[0] * zz + tan_p[1]) * zz + tan_p[2]) / - ((((zz + tan_q[1]) * zz + tan_q[2]) * zz + tan_q[3]) * zz + tan_q[4]) - ) - } else { - z - } - if (j & 2L) == 2L { - zz = -1.0 / zz - } - if x.signum() < 0.0 { - -zz - } else { - zz - } - } + x.tan() } ///| diff --git a/math/trigonometric_test.mbt b/math/trigonometric_test.mbt index b8358fefc..9e91f462a 100644 --- a/math/trigonometric_test.mbt +++ b/math/trigonometric_test.mbt @@ -30,111 +30,6 @@ let vf = [ 1.8253080916808550e+00, -8.6859247685756013e+00, ] -test "sin" { - let sin_res : FixedArray[_] = [ - -9.6466616586009283766724726e-01, 9.9338225271646545763467022e-01, -2.7335587039794393342449301e-01, - 9.5586257685042792878173752e-01, -2.099421066779969164496634e-01, 2.135578780799860532750616e-01, - -8.694568971167362743327708e-01, 4.019566681155577786649878e-01, 9.6778633541687993721617774e-01, - -6.734405869050344734943028e-01, - ] - let sin_res_large : FixedArray[_] = [ - -9.646661658548936063912e-01, 9.933822527198506903752e-01, -2.7335587036246899796e-01, - 9.55862576853689321268e-01, -2.099421066862688873691e-01, 2.13557878070308981163e-01, - -8.694568970959221300497e-01, 4.01956668098863248917e-01, 9.67786335404528727927e-01, - -6.7344058693131973066e-01, - ] - inspect!(@math.sin(@double.not_a_number), content="NaN") - inspect!(@math.sin(@double.infinity), content="NaN") - inspect!(@math.sin(@double.neg_infinity), content="NaN") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!(imprecise_equal(@math.sin(vf[i]), sin_res[i]), content="true") - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.sin(vf[i]), sin_res[i]), - ) - inspect!( - imprecise_equal(@math.sin(100000.0 * PI + vf[i]), sin_res_large[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.sin(100000.0 * PI + vf[i]), sin_res_large[i]), - ) - } - inspect!(max_inaccuracy, content="5.6910670620524684e-11") -} - -test "cos" { - let cos_res : FixedArray[_] = [ - 2.634752140995199110787593e-01, 1.148551260848219865642039e-01, 9.6191297325640768154550453e-01, - 2.938141150061714816890637e-01, -9.777138189897924126294461e-01, -9.7693041344303219127199518e-01, - 4.940088096948647263961162e-01, -9.1565869021018925545016502e-01, -2.517729313893103197176091e-01, - -7.39241351595676573201918e-01, - ] - let cos_res_large : FixedArray[_] = [ - 2.634752141185559426744e-01, 1.14855126055543100712e-01, 9.61912973266488928113e-01, - 2.9381411499556122552e-01, -9.777138189880161924641e-01, -9.76930413445147608049e-01, - 4.940088097314976789841e-01, -9.15658690217517835002e-01, -2.51772931436786954751e-01, - -7.3924135157173099849e-01, - ] - inspect!(@math.cos(@double.not_a_number), content="NaN") - inspect!(@math.cos(@double.infinity), content="NaN") - inspect!(@math.cos(@double.neg_infinity), content="NaN") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!(imprecise_equal(@math.cos(vf[i]), cos_res[i]), content="true") - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.cos(vf[i]), cos_res[i]), - ) - inspect!( - imprecise_equal(@math.cos(100000.0 * PI + vf[i]), cos_res_large[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.cos(100000.0 * PI + vf[i]), cos_res_large[i]), - ) - } - inspect!(max_inaccuracy, content="5.7822691079678634e-11") -} - -test "tan" { - let tan_res : FixedArray[_] = [ - -3.661316565040227801781974e+00, 8.64900232648597589369854e+00, -2.8417941955033612725238097e-01, - 3.253290185974728640827156e+00, 2.147275640380293804770778e-01, -2.18600910711067004921551e-01, - -1.760002817872367935518928e+00, -4.389808914752818126249079e-01, -3.843885560201130679995041e+00, - 9.10988793377685105753416e-01, - ] - let tan_res_large : FixedArray[_] = [ - -3.66131656475596512705e+00, 8.6490023287202547927e+00, -2.841794195104782406e-01, - 3.2532901861033120983e+00, 2.14727564046880001365e-01, -2.18600910700688062874e-01, - -1.760002817699722747043e+00, -4.38980891453536115952e-01, -3.84388555942723509071e+00, - 9.1098879344275101051e-01, - ] - inspect!(@math.tan(@double.not_a_number), content="NaN") - inspect!(@math.tan(@double.infinity), content="NaN") - inspect!(@math.tan(@double.neg_infinity), content="NaN") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!(imprecise_equal(@math.tan(vf[i]), tan_res[i]), content="true") - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.tan(vf[i]), tan_res[i]), - ) - inspect!( - imprecise_equal(@math.tan(100000.0 * PI + vf[i]), tan_res_large[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.tan(100000.0 * PI + vf[i]), tan_res_large[i]), - ) - } - inspect!(max_inaccuracy, content="4.4124632836428646e-9") -} - test "atan" { let atan_res : FixedArray[_] = [ 1.372590262129621651920085e+00, 1.442290609645298083020664e+00, -2.7011324359471758245192595e-01, From f32edd50b9830c1eb02dc9eaa75749ff78ec8190 Mon Sep 17 00:00:00 2001 From: Kaida-Amethyst Date: Fri, 31 Jan 2025 15:22:22 +0800 Subject: [PATCH 2/3] better asin, acos, atan, atan2 implementation --- double/double.mbti | 4 + double/trig_js.mbt | 112 +++++++++++ double/trig_nonjs.mbt | 380 ++++++++++++++++++++++++++++++++++++ double/trig_test.mbt | 139 +++++++++++++ math/trigonometric.mbt | 124 +----------- math/trigonometric_test.mbt | 182 ----------------- 6 files changed, 639 insertions(+), 302 deletions(-) delete mode 100644 math/trigonometric_test.mbt diff --git a/double/double.mbti b/double/double.mbti index bfd8c8688..1812e7dd9 100644 --- a/double/double.mbti +++ b/double/double.mbti @@ -18,6 +18,10 @@ let not_a_number : Double // Types and methods impl Double { abs(Double) -> Double + acos(Double) -> Double + asin(Double) -> Double + atan(Double) -> Double + atan2(Double, Double) -> Double ceil(Double) -> Double cos(Double) -> Double exp(Double) -> Double diff --git a/double/trig_js.mbt b/double/trig_js.mbt index 601530942..4b8dbdc22 100644 --- a/double/trig_js.mbt +++ b/double/trig_js.mbt @@ -89,3 +89,115 @@ pub fn cos(self : Double) -> Double = "Math" "cos" /// } /// ``` pub fn tan(self : Double) -> Double = "Math" "tan" + +///| +/// Calculates the arcsine of a number. Handles special cases and edge conditions +/// according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arcsine. +/// +/// Returns the arcsine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "asin" { +/// inspect!(0.0.asin(), content="0") +/// inspect!(1.0.asin(), content="1.5707963267948966") +/// inspect!((-1.0).asin(), content="-1.5707963267948966") +/// inspect!(@double.not_a_number.asin(), content="NaN") +/// inspect!(@double.infinity.asin(), content="NaN") +/// inspect!(@double.neg_infinity.asin(), content="NaN") +/// } +/// ``` +pub fn asin(self : Double) -> Double = "Math" "asin" + +///| +/// Calculates the arccosine of a number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arccosine. +/// +/// Returns the arccosine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "acos" { +/// inspect!(0.0.acos(), content="1.5707963267948966") +/// inspect!(1.0.acos(), content="0") +/// inspect!((-1.0).acos(), content="3.141592653589793") +/// inspect!(@double.not_a_number.acos(), content="NaN") +/// inspect!(@double.infinity.acos(), content="NaN") +/// inspect!(@double.neg_infinity.acos(), content="NaN") +/// inspect!(0.0.acos(), content="1.5707963267948966") +/// } +/// ``` +pub fn acos(self : Double) -> Double = "Math" "acos" + +///| +/// Calculates the arctangent of a number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arctangent. +/// +/// Returns the arctangent of the number `x`. +/// +/// Example: +/// +/// * Returns NaN if the input is NaN. +/// +/// ```moonbit +/// test "atan" { +/// inspect!(0.0.atan(), content="0") +/// inspect!(1.0.atan(), content="0.7853981633974483") +/// inspect!((-1.0).atan(), content="-0.7853981633974483") +/// inspect!(@double.not_a_number.atan(), content="NaN") +/// inspect!(@double.infinity.atan(), content="1.5707963267948966") +/// inspect!(@double.neg_infinity.atan(), content="-1.5707963267948966") +/// } +/// ``` +pub fn atan(self : Double) -> Double = "Math" "atan" + +///| +/// Calculates the arctangent of the quotient of two numbers. +/// +/// Parameters: +/// +/// * `self` : The numerator of the quotient. +/// * `x` : The denominator of the quotient. +/// +/// Returns the arctangent of the quotient `self / x`. +/// +/// * Returns NaN if self or x is NaN. +/// +/// Example: +/// +/// ```moonbit +/// test "atan2" { +/// inspect!(0.0.atan2(-1.0), content="3.141592653589793") +/// inspect!(1.0.atan2(0.0), content="1.5707963267948966") +/// inspect!(1.0.atan2(1.0), content="0.7853981633974483") +/// inspect!(not_a_number.atan2(1.0), content="NaN") +/// inspect!(1.0.atan2(not_a_number), content="NaN") +/// inspect!(infinity.atan2(1.0), content="1.5707963267948966") +/// inspect!(1.0.atan2(infinity), content="0") +/// inspect!(neg_infinity.atan2(1.0), content="-1.5707963267948966") +/// inspect!(1.0.atan2(neg_infinity), content="3.141592653589793") +/// inspect!(infinity.atan2(infinity), content="0.7853981633974483") +/// inspect!(neg_infinity.atan2(neg_infinity), content="-2.356194490192345") +/// inspect!(infinity.atan2(neg_infinity), content="2.356194490192345") +/// inspect!(neg_infinity.atan2(infinity), content="-0.7853981633974483") +/// } +/// ``` +pub fn atan2(self : Double, x : Double) -> Double = "Math" "atan2" diff --git a/double/trig_nonjs.mbt b/double/trig_nonjs.mbt index 23329b849..b1ba30ed1 100644 --- a/double/trig_nonjs.mbt +++ b/double/trig_nonjs.mbt @@ -12,6 +12,24 @@ // See the License for the specific language governing permissions and // limitations under the License. +// +// origin: FreeBSD /usr/src/lib/msun/src/e_sin.c +// origin: FreeBSD /usr/src/lib/msun/src/e_cos.c +// origin: FreeBSD /usr/src/lib/msun/src/e_tan.c +// origin: FreeBSD /usr/src/lib/msun/src/e_asin.c +// origin: FreeBSD /usr/src/lib/msun/src/e_acos.c +// origin: FreeBSD /usr/src/lib/msun/src/s_atan.c +// origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunSoft, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// + ///| let two_over_pi : Array[Int] = [ 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, @@ -768,3 +786,365 @@ pub fn tan(self : Double) -> Double { __kernal_tan(y[0], y[1], 1 - ((n & 1) << 1)) } } + +///| +/// Calculates the arcsine of a number. Handles special cases and edge conditions +/// according to IEEE 754 standards. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arcsine. +/// +/// Returns the arcsine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "asin" { +/// inspect!(0.0.asin(), content="0") +/// inspect!(1.0.asin(), content="1.5707963267948966") +/// inspect!((-1.0).asin(), content="-1.5707963267948966") +/// inspect!(@double.not_a_number.asin(), content="NaN") +/// inspect!(@double.infinity.asin(), content="NaN") +/// inspect!(@double.neg_infinity.asin(), content="NaN") +/// } +/// ``` +pub fn asin(self : Double) -> Double { + let huge = 1.0e+300 + let pio4_hi = 7.85398163397448278999e-01 + let pio2_hi = 1.57079632679489655800 + let pio2_lo = 6.12323399573676603587e-17 + let ps0 = 1.66666666666666657415e-01 + let ps1 = -3.25565818622400915405e-01 + let ps2 = 2.01212532134862925881e-01 + let ps3 = -4.00555345006794114027e-02 + let ps4 = 7.91534994289814532176e-04 + let ps5 = 3.47933107596021167570e-05 + let qs1 = -2.40339491173441421878e+00 + let qs2 = 2.02094576023350569471e+00 + let qs3 = -6.88283971605453293030e-01 + let qs4 = 7.70381505559019352791e-02 + let ix = get_high_word(self).reinterpret_as_int() & 0x7fffffff + let absx = self.abs() + if absx >= 1.0 { + if absx == 1.0 { + return self * pio2_hi + self * pio2_lo + } else { + return not_a_number + } + } else if absx < 0.5 { + if ix < 0x3e400000 { + if huge + self > 1.0 { + return self + } + } else { + let t = self * self + let p = t * + (ps0 + t * (ps1 + t * (ps2 + t * (ps3 + t * (ps4 + t * ps5))))) + let q = 1.0 + t * (qs1 + t * (qs2 + t * (qs3 + t * qs4))) + let w = p / q + return self + self * w + } + } + let w = 1.0 - absx + let t = w * 0.5 + let p = t * (ps0 + t * (ps1 + t * (ps2 + t * (ps3 + t * (ps4 + t * ps5))))) + let q = 1.0 + t * (qs1 + t * (qs2 + t * (qs3 + t * qs4))) + let s = t.sqrt() + if ix >= 0x3FEF3333 { + let w = p / q + let t = pio2_hi - (2.0 * (s + s * w) - pio2_lo) + return if self > 0.0 { t } else { -t } + } else { + let mut w = s + w = set_low_word(w, 0) + let c = (t - w * w) / (s + w) + let r = p / q + let p = 2.0 * s * r - (pio2_lo - 2.0 * c) + let q = pio4_hi - 2.0 * w + let t = pio4_hi - (p - q) + return if self > 0.0 { t } else { -t } + } +} + +///| +/// Calculates the arccosine of a number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arccosine. +/// +/// Returns the arccosine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "acos" { +/// inspect!(0.0.acos(), content="1.5707963267948966") +/// inspect!(1.0.acos(), content="0") +/// inspect!((-1.0).acos(), content="3.141592653589793") +/// inspect!(@double.not_a_number.acos(), content="NaN") +/// inspect!(@double.infinity.acos(), content="NaN") +/// inspect!(@double.neg_infinity.acos(), content="NaN") +/// inspect!(0.0.acos(), content="1.5707963267948966") +/// } +/// ``` +pub fn acos(self : Double) -> Double { + let one : Double = 1.0 + let pi : Double = 3.14159265358979311600 + let pio2_hi : Double = 1.57079632679489655800 + let pio2_lo : Double = 6.12323399573676603587e-17 + let ps0 : Double = 1.66666666666666657415e-01 + let ps1 : Double = -3.25565818622400915405e-01 + let ps2 : Double = 2.01212532134862925881e-01 + let ps3 : Double = -4.00555345006794114027e-02 + let ps4 : Double = 7.91534994289814532176e-04 + let ps5 : Double = 3.47933107596021167570e-05 + let qs1 : Double = -2.40339491173441421878e+00 + let qs2 : Double = 2.02094576023350569471e+00 + let qs3 : Double = -6.88283971605453293030e-01 + let qs4 : Double = 7.70381505559019352791e-02 + let ix = get_high_word(self).reinterpret_as_int() & 0x7fffffff + let absx = self.abs() + if absx >= 1.0 { + if absx == 1.0 { + if self > 0 { + return 0.0 + } else { + return pi + 2.0 * pio2_lo + } + } + return not_a_number + } + if absx < 0.5 { + if ix <= 0x3c600000 { + return pio2_hi + pio2_lo + } + let z = self * self + let p = z * (ps0 + z * (ps1 + z * (ps2 + z * (ps3 + z * (ps4 + z * ps5))))) + let q = one + z * (qs1 + z * (qs2 + z * (qs3 + z * qs4))) + let r = p / q + pio2_hi - (self - (pio2_lo - self * r)) + } else if self < 0 { + let z = (one + self) * 0.5 + let p = z * (ps0 + z * (ps1 + z * (ps2 + z * (ps3 + z * (ps4 + z * ps5))))) + let q = one + z * (qs1 + z * (qs2 + z * (qs3 + z * qs4))) + let s = z.sqrt() + let r = p / q + let w = r * s - pio2_lo + pi - 2.0 * (s + w) + } else { + let z = (one - self) * 0.5 + let s = z.sqrt() + let df = s + let c = (z - df * df) / (s + df) + let p = z * (ps0 + z * (ps1 + z * (ps2 + z * (ps3 + z * (ps4 + z * ps5))))) + let q = one + z * (qs1 + z * (qs2 + z * (qs3 + z * qs4))) + let r = p / q + let w = r * s + c + 2.0 * (df + w) + } +} + +///| +/// Calculates the arctangent of a number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arctangent. +/// +/// Returns the arctangent of the number `x`. +/// +/// Example: +/// +/// * Returns NaN if the input is NaN. +/// +/// ```moonbit +/// test "atan" { +/// inspect!(0.0.atan(), content="0") +/// inspect!(1.0.atan(), content="0.7853981633974483") +/// inspect!((-1.0).atan(), content="-0.7853981633974483") +/// inspect!(@double.not_a_number.atan(), content="NaN") +/// inspect!(@double.infinity.atan(), content="1.5707963267948966") +/// inspect!(@double.neg_infinity.atan(), content="-1.5707963267948966") +/// } +/// ``` +pub fn atan(self : Double) -> Double { + if self.is_nan() || self == 0.0 { + return self + } + let atan_hi = [ + 4.63647609000806093515e-01, 7.85398163397448278999e-01, 9.82793723247329054082e-01, + 1.57079632679489655800e+00, + ] + let atan_lo = [ + 2.26987774529616870924e-17, 3.06161699786838301793e-17, 1.39033110312309984516e-17, + 6.12323399573676603587e-17, + ] + let a_t = [ + 3.33333333333329318027e-01, -1.99999999998764832476e-01, 1.42857142725034663711e-01, + -1.11111104054623557880e-01, 9.09088713343650656196e-02, -7.69187620504482999495e-02, + 6.66107313738753120669e-02, -5.83357013379057348645e-02, 4.97687799461593236017e-02, + -3.65315727442169155270e-02, 1.62858201153657823623e-02, + ] + let one = 1.0 + let huge = 1.0e300 + let ix = get_high_word(self).reinterpret_as_int() & 0x7fffffff + let mut id = 0 + let mut z = 0.0 + let mut w = 0.0 + let mut self = self + let x_is_neg = self < 0.0 + if ix >= 0x44100000 { + if self > 0 { + return atan_hi[3] + atan_lo[3] + } else { + return -atan_hi[3] - atan_lo[3] + } + } + if ix < 0x3fdc0000 { + if ix < 0x3e200000 { + if huge + self > one { + return self + } + } + id = -1 + } else { + self = self.abs() + if ix < 0x3ff30000 { + if ix < 0x3fe60000 { + id = 0 + self = (2.0 * self - one) / (2.0 + self) + } else { + id = 1 + self = (self - one) / (self + one) + } + } else if ix < 0x40038000 { + id = 2 + self = (self - 1.5) / (one + 1.5 * self) + } else { + id = 3 + self = -1.0 / self + } + } + z = self * self + w = z * z + let s1 = z * + ( + a_t[0] + + w * (a_t[2] + w * (a_t[4] + w * (a_t[6] + w * (a_t[8] + w * a_t[10])))) + ) + let s2 = w * + (a_t[1] + w * (a_t[3] + w * (a_t[5] + w * (a_t[7] + w * a_t[9])))) + if id < 0 { + self - self * (s1 + s2) + } else { + z = atan_hi[id] - (self * (s1 + s2) - atan_lo[id] - self) + if x_is_neg { + -z + } else { + z + } + } +} + +///| +/// Calculates the arctangent of the quotient of two numbers. +/// +/// Parameters: +/// +/// * `self` : The numerator of the quotient. +/// * `x` : The denominator of the quotient. +/// +/// Returns the arctangent of the quotient `self / x`. +/// +/// * Returns NaN if self or x is NaN. +/// +/// Example: +/// +/// ```moonbit +/// test "atan2" { +/// inspect!(0.0.atan2(-1.0), content="3.141592653589793") +/// inspect!(1.0.atan2(0.0), content="1.5707963267948966") +/// inspect!(1.0.atan2(1.0), content="0.7853981633974483") +/// inspect!(not_a_number.atan2(1.0), content="NaN") +/// inspect!(1.0.atan2(not_a_number), content="NaN") +/// inspect!(infinity.atan2(1.0), content="1.5707963267948966") +/// inspect!(1.0.atan2(infinity), content="0") +/// inspect!(neg_infinity.atan2(1.0), content="-1.5707963267948966") +/// inspect!(1.0.atan2(neg_infinity), content="3.141592653589793") +/// inspect!(infinity.atan2(infinity), content="0.7853981633974483") +/// inspect!(neg_infinity.atan2(neg_infinity), content="-2.356194490192345") +/// inspect!(infinity.atan2(neg_infinity), content="2.356194490192345") +/// inspect!(neg_infinity.atan2(infinity), content="-0.7853981633974483") +/// } +/// ``` +pub fn atan2(self : Double, x : Double) -> Double { + if x.is_nan() || self.is_nan() { + return not_a_number + } + let tiny = 1.0e-300 + let zero = 0.0 + let pi_o_4 = 7.8539816339744827900E-01 + let pi_o_2 = 1.5707963267948965580E+00 + let pi = 3.1415926535897931160E+00 + let pi_lo = 1.2246467991473531772E-16 + let hx = get_high_word(x).reinterpret_as_int() + let hy = get_high_word(self).reinterpret_as_int() + let ix = hx & 0x7fffffff + let iy = hy & 0x7fffffff + if x == 1.0 { + return self.atan() + } + let m = ((hy >> 31) & 1) | ((hx >> 30) & 2) + if self == 0 { + match m { + 0 | 1 => return self + 2 => return pi + tiny + _ => return -pi - tiny + } + } + if x == 0 { + return if hy < 0 { -pi_o_2 - tiny } else { pi_o_2 + tiny } + } + if x.is_inf() { + if self.is_inf() { + match m { + 0 => return pi_o_4 + tiny + 1 => return -pi_o_4 - tiny + 2 => return 3.0 * pi_o_4 + tiny + _ => return -3.0 * pi_o_4 - tiny + } + } else { + match m { + 0 => return zero + 1 => return -zero + 2 => return pi + tiny + _ => return -pi - tiny + } + } + } + if self.is_inf() { + return if hy < 0 { -pi_o_2 - tiny } else { pi_o_2 + tiny } + } + let k = (iy - ix) >> 20 + let z = if k > 60 { + pi_o_2 + 0.5 * pi_lo + } else if hx < 0 && k < -60 { + 0.0 + } else { + (self / x).abs().atan() + } + match m { + 0 => z + 1 => -z + 2 => pi - (z - pi_lo) + _ => z - pi_lo - pi + } +} diff --git a/double/trig_test.mbt b/double/trig_test.mbt index b4c727314..36db8cde3 100644 --- a/double/trig_test.mbt +++ b/double/trig_test.mbt @@ -92,3 +92,142 @@ test "trigonometric functions" { high_accuracy_test!(vals[i].tan(), tan_res[i]) } } + +test "Inverse trigonometric functions" { + // asin + inspect!(0.0.asin(), content="0") + inspect!(1.0.asin(), content="1.5707963267948966") + inspect!((-1.0).asin(), content="-1.5707963267948966") + inspect!(0.5.asin(), content="0.5235987755982989") + inspect!(0.45.asin(), content="0.4667653390472964") + inspect!(@double.not_a_number.asin(), content="NaN") + inspect!(@double.infinity.asin(), content="NaN") + inspect!(@double.neg_infinity.asin(), content="NaN") + + // acos + inspect!(0.0.acos(), content="1.5707963267948966") + inspect!(1.0.acos(), content="0") + inspect!((-1.0).acos(), content="3.141592653589793") + inspect!(0.5.acos(), content="1.0471975511965979") + inspect!(0.45.acos(), content="1.1040309877476002") + inspect!(@double.not_a_number.acos(), content="NaN") + inspect!(@double.infinity.acos(), content="NaN") + inspect!(@double.neg_infinity.acos(), content="NaN") + + // High accuracy test + let vals = [ + 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.123, 0.234, 0.345, 0.456, 0.567, + 0.678, 0.789, 0.891, 0.912, 0.1234, -0.135, -0.246, -0.357, -0.468, -0.579, -0.681, + -0.792, -0.893, -0.994, -0.12345, + ] + let asin_res = [ + 0.1001674211615598, 0.2013579207903308, 0.3046926540153975, 0.41151684606748806, + 0.5235987755982989, 0.6435011087932844, 0.775397496610753, 0.9272952180016123, + 1.1197695149986342, 0.12331227519187199, 0.23618988434810217, 0.35223878509706474, + 0.47349551215005636, 0.6028592463891203, 0.745038349213386, 0.9091796699171293, + 1.0995430581994015, 1.148133721990592, 0.12371534584255098, -0.13541346246655556, + -0.24855126213339576, -0.36505428193897027, -0.4870262796643988, -0.61750165481717, + -0.7491273605521073, -0.9140779620461925, -1.1039675362433812, -1.4611969689632571, + -0.12376573109305462, + ] + let acos_res = [ + 1.4706289056333368, 1.369438406004566, 1.2661036727794992, 1.1592794807274085, + 1.0471975511965979, 0.9272952180016123, 0.7953988301841436, 0.6435011087932843, + 0.45102681179626236, 1.4474840516030247, 1.3346064424467945, 1.2185575416978318, + 1.0973008146448402, 0.9679370804057764, 0.8257579775815106, 0.6616166568777674, + 0.47125326859549527, 0.4226626048043045, 1.4470809809523457, 1.706209789261452, + 1.8193475889282924, 1.935850608733867, 2.0578226064592955, 2.1882979816120667, + 2.319923687347004, 2.4848742888410893, 2.6747638630382777, 3.031993295758154, + 1.6945620578879512, + ] + fn high_accuracy_test(expect : Double, actual : Double) -> Unit!Error { + let ulp = (expect.reinterpret_as_int64() - actual.reinterpret_as_int64()).abs() + assert_eq!(ulp <= 1, true) + } + + for i in 0.. Unit!Error { + let ulp = (expect.reinterpret_as_int64() - actual.reinterpret_as_int64()).abs() + assert_eq!(ulp <= 1, true) + } + + for i in 0.. Unit!Error { + let ulp = (expect.reinterpret_as_int64() - actual.reinterpret_as_int64()).abs() + assert_eq!(ulp <= 1, true) + } + + let mut idx = 0 + for y in vals { + for x in vals { + high_accuracy_test!(y.atan2(x), atan2_res[idx]) + idx = idx + 1 + } + } +} diff --git a/math/trigonometric.mbt b/math/trigonometric.mbt index b271660c2..eebdfe318 100644 --- a/math/trigonometric.mbt +++ b/math/trigonometric.mbt @@ -34,138 +34,22 @@ pub fn tan(x : Double) -> Double { x.tan() } -///| -let atan_p = [ - -8.750608600031904122785e-01, -1.615753718733365076637e+01, -7.500855792314704667340e+01, - -1.228866684490136173410e+02, -6.485021904942025371773e+01, -] - -///| -let atan_q = [ - 2.485846490142306297962e+01, 1.650270098316988542046e+02, 4.328810604912902668951e+02, - 4.853903996359136964868e+02, 1.945506571482613964425e+02, -] - -///| -let morebits = 6.123233995736765886130e-17 - -///| -let tan3pio8 = 2.41421356237309504880 - -///| -fn xatan(x : Double) -> Double { - let z = x * x - x * - ( - z * - ( - (((atan_p[0] * z + atan_p[1]) * z + atan_p[2]) * z + atan_p[3]) * z + - atan_p[4] - ) / - ( - ((((z + atan_q[0]) * z + atan_q[1]) * z + atan_q[2]) * z + atan_q[3]) * z + - atan_q[4] - ) - ) + - x -} - -///| -fn satan(x : Double) -> Double { - if x <= 0.66 { - xatan(x) - } else if x > tan3pio8 { - PI / 2.0 - xatan(1.0 / x) + morebits - } else { - PI / 4.0 + xatan((x - 1.0) / (x + 1.0)) + 0.5 * morebits - } -} - ///| pub fn atan(x : Double) -> Double { - if Double::is_nan(x) || x == 0.0 { - x - } else if x > 0.0 { - satan(x) - } else { - -satan(-x) - } + x.atan() } ///| pub fn asin(x : Double) -> Double { - if Double::is_nan(x) || x == 0.0 { - x - } else { - let x_ = x.abs() - if x_ > 1.0 { - @double.not_a_number - } else { - let temp = (1.0 - x_ * x_).sqrt() - (if x > 0.7 { PI / 2.0 - satan(temp / x_) } else { satan(x_ / temp) }) * - x.signum() - } - } + x.asin() } ///| pub fn acos(x : Double) -> Double { - if Double::is_nan(x) { - x - } else { - PI / 2.0 - asin(x) - } -} - -///| -const SIGNBIT : UInt64 = 0x8000000000000000UL - -///| -fn copysign(x : Double, y : Double) -> Double { - // should check the sign bit to distinguish +0/-0 - if (y.reinterpret_as_uint64() & SIGNBIT) > 0 { - -x.abs() - } else { - x.abs() - } + x.acos() } ///| pub fn atan2(y : Double, x : Double) -> Double { - if Double::is_nan(x) || Double::is_nan(y) { - @double.not_a_number - } else if y == 0 { - if x >= 0 { - copysign(0, y) - } else { - copysign(PI, y) - } - } else if x == 0 { - copysign(PI / 2, y) - } else if x.is_inf() { - if x.is_pos_inf() { - if y.is_inf() { - copysign(PI / 4, y) - } else { - copysign(0, y) - } - } else if y.is_inf() { - copysign(3.0 * PI / 4, y) - } else { - copysign(PI, y) - } - } else if y.is_inf() { - copysign(PI / 2, y) - } else { - let q = atan(y / x) - if x < 0 { - if q <= 0 { - q + PI - } else { - q - PI - } - } else { - q - } - } + y.atan2(x) } diff --git a/math/trigonometric_test.mbt b/math/trigonometric_test.mbt deleted file mode 100644 index 1e4e4f23f..000000000 --- a/math/trigonometric_test.mbt +++ /dev/null @@ -1,182 +0,0 @@ -// Copyright 2025 International Digital Economy Academy -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -///| -//Compare two Double values for equality within the allowed error range of 4.5e-9. -fn imprecise_equal(x : Double, y : Double) -> Bool { - (x - y).abs() < 4.5e-9 -} - -///| -fn inaccuracy(x : Double, y : Double) -> Double { - (x - y).abs() -} - -///| -let vf = [ - 4.9790119248836735e+00, 7.7388724745781045e+00, -2.7688005719200159e-01, -5.0106036182710749e+00, - 9.6362937071984173e+00, 2.9263772392439646e+00, 5.2290834314593066e+00, 2.7279399104360102e+00, - 1.8253080916808550e+00, -8.6859247685756013e+00, -] - -test "atan" { - let atan_res : FixedArray[_] = [ - 1.372590262129621651920085e+00, 1.442290609645298083020664e+00, -2.7011324359471758245192595e-01, - -1.3738077684543379452781531e+00, 1.4673921193587666049154681e+00, 1.2415173565870168649117764e+00, - 1.3818396865615168979966498e+00, 1.2194305844639670701091426e+00, 1.0696031952318783760193244e+00, - -1.4561721938838084990898679e+00, - ] - inspect!(@math.atan(@double.not_a_number), content="NaN") - inspect!(@math.atan(@double.infinity), content="1.5707963267948966") - inspect!(@math.atan(@double.neg_infinity), content="-1.5707963267948966") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!(imprecise_equal(@math.atan(vf[i]), atan_res[i]), content="true") - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.atan(vf[i]), atan_res[i]), - ) - } - inspect!(max_inaccuracy, content="2.220446049250313e-16") -} - -test "asin" { - let asin_res : FixedArray[_] = [ - 5.2117697218417440497416805e-01, 8.8495619865825236751471477e-01, -02.769154466281941332086016e-02, - -5.2482360935268931351485822e-01, 1.3002662421166552333051524e+00, 2.9698415875871901741575922e-01, - 5.5025938468083370060258102e-01, 2.7629597861677201301553823e-01, 1.83559892257451475846656e-01, - -1.0523547536021497774980928e+00, - ] - inspect!(@math.asin(@double.not_a_number), content="NaN") - inspect!(@math.asin(@double.infinity), content="NaN") - inspect!(@math.asin(@double.neg_infinity), content="NaN") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!( - imprecise_equal(@math.asin(vf[i] / 10.0), asin_res[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.asin(vf[i] / 10.0), asin_res[i]), - ) - } - inspect!(max_inaccuracy, content="2.220446049250313e-16") -} - -test "acos" { - let acos_res : FixedArray[_] = [ - 1.0496193546107222142571536e+00, 6.8584012813664425171660692e-01, 1.5984878714577160325521819e+00, - 2.0956199361475859327461799e+00, 2.7053008467824138592616927e-01, 1.2738121680361776018155625e+00, - 1.0205369421140629186287407e+00, 1.2945003481781246062157835e+00, 1.3872364345374451433846657e+00, - 2.6231510803970463967294145e+00, - ] - inspect!(@math.acos(@double.not_a_number), content="NaN") - inspect!(@math.acos(@double.infinity), content="NaN") - inspect!(@math.acos(@double.neg_infinity), content="NaN") - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!( - imprecise_equal(@math.acos(vf[i] / 10.0), acos_res[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.acos(vf[i] / 10.0), acos_res[i]), - ) - } - inspect!(max_inaccuracy, content="2.220446049250313e-16") -} - -test "atan2" { - let atan2_res : FixedArray[Double] = [ - 1.1088291730037004444527075e+00, 9.1218183188715804018797795e-01, 1.5984772603216203736068915e+00, - 2.0352918654092086637227327e+00, 8.0391819139044720267356014e-01, 1.2861075249894661588866752e+00, - 1.0889904479131695712182587e+00, 1.3044821793397925293797357e+00, 1.3902530903455392306872261e+00, - 2.2859857424479142655411058e+00, - ] - let mut max_inaccuracy = 0.0 - for i = 0; i < vf.length(); i = i + 1 { - inspect!( - imprecise_equal(@math.atan2(10, vf[i]), atan2_res[i]), - content="true", - ) - max_inaccuracy = @math.maximum( - max_inaccuracy, - inaccuracy(@math.atan2(10, vf[i]), atan2_res[i]), - ) - } - inspect!(max_inaccuracy, content="4.440892098500626e-16") -} - -test "tan with very small value" { - let x = 1.0e-8 // x * x = 1.0e-16 < 1.0e-14 - inspect!(@math.tan(x), content="1e-8") -} - -test "asin with value greater than 1" { - inspect!(@math.asin(1.5), content="NaN") -} - -test "atan2 with NaN" { - let y = @double.not_a_number - let x = 1.0 - inspect!(atan2(y, x), content="NaN") - inspect!(atan2(1.0, @double.not_a_number), content="NaN") -} - -test "atan2 with zero x" { - inspect!(atan2(1.0, 0.0), content="1.5707963267948966") - inspect!(atan2(-1.0, 0.0), content="-1.5707963267948966") -} - -test "atan2 with zero y" { - inspect!(atan2(0.0, 1.0), content="0") - inspect!(atan2(0.0, -1.0), content="3.141592653589793") - inspect!(atan2(-0.0, -1.0), content="-3.141592653589793") -} - -test "atan2 with infinite y and finite x" { - inspect!(atan2(@double.infinity, 1.0), content="1.5707963267948966") - inspect!(atan2(-@double.infinity, 1.0), content="-1.5707963267948966") -} - -test "atan2 with infinite x and finite y" { - inspect!(@math.atan2(1.0, @double.infinity), content="0") - inspect!(@math.atan2(1.0, @double.neg_infinity), content="3.141592653589793") - inspect!(@math.atan2(-1.0, @double.infinity), content="0") - inspect!( - @math.atan2(-1.0, @double.neg_infinity), - content="-3.141592653589793", - ) -} - -test "atan2 with infinite x and infinite y" { - inspect!( - @math.atan2(@double.infinity, @double.infinity), - content="0.7853981633974483", - ) - inspect!( - @math.atan2(@double.neg_infinity, @double.infinity), - content="-0.7853981633974483", - ) - inspect!( - @math.atan2(@double.infinity, @double.neg_infinity), - content="2.356194490192345", - ) - inspect!( - @math.atan2(@double.neg_infinity, @double.neg_infinity), - content="-2.356194490192345", - ) -} From 5f1ea099a5e2f9270315a0757310635bfe987620 Mon Sep 17 00:00:00 2001 From: Kaida-Amethyst Date: Sat, 1 Feb 2025 11:33:50 +0800 Subject: [PATCH 3/3] add sin, cos, tan, asin, acos, atan, atan2 for floating-point --- double/trig_nonjs.mbt | 18 ++-- double/trig_test.mbt | 3 +- float/float.mbti | 7 ++ float/trig.mbt | 206 ++++++++++++++++++++++++++++++++++++++++++ float/trig_test.mbt | 88 ++++++++++++++++++ 5 files changed, 312 insertions(+), 10 deletions(-) create mode 100644 float/trig.mbt create mode 100644 float/trig_test.mbt diff --git a/double/trig_nonjs.mbt b/double/trig_nonjs.mbt index b1ba30ed1..0d46d93d1 100644 --- a/double/trig_nonjs.mbt +++ b/double/trig_nonjs.mbt @@ -13,6 +13,7 @@ // limitations under the License. // +// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c // origin: FreeBSD /usr/src/lib/msun/src/e_sin.c // origin: FreeBSD /usr/src/lib/msun/src/e_cos.c // origin: FreeBSD /usr/src/lib/msun/src/e_tan.c @@ -892,7 +893,6 @@ pub fn asin(self : Double) -> Double { /// inspect!(@double.not_a_number.acos(), content="NaN") /// inspect!(@double.infinity.acos(), content="NaN") /// inspect!(@double.neg_infinity.acos(), content="NaN") -/// inspect!(0.0.acos(), content="1.5707963267948966") /// } /// ``` pub fn acos(self : Double) -> Double { @@ -1073,16 +1073,16 @@ pub fn atan(self : Double) -> Double { /// inspect!(0.0.atan2(-1.0), content="3.141592653589793") /// inspect!(1.0.atan2(0.0), content="1.5707963267948966") /// inspect!(1.0.atan2(1.0), content="0.7853981633974483") -/// inspect!(not_a_number.atan2(1.0), content="NaN") +/// inspect!(@double.not_a_number.atan2(1.0), content="NaN") /// inspect!(1.0.atan2(not_a_number), content="NaN") -/// inspect!(infinity.atan2(1.0), content="1.5707963267948966") +/// inspect!(@double.infinity.atan2(1.0), content="1.5707963267948966") /// inspect!(1.0.atan2(infinity), content="0") -/// inspect!(neg_infinity.atan2(1.0), content="-1.5707963267948966") -/// inspect!(1.0.atan2(neg_infinity), content="3.141592653589793") -/// inspect!(infinity.atan2(infinity), content="0.7853981633974483") -/// inspect!(neg_infinity.atan2(neg_infinity), content="-2.356194490192345") -/// inspect!(infinity.atan2(neg_infinity), content="2.356194490192345") -/// inspect!(neg_infinity.atan2(infinity), content="-0.7853981633974483") +/// inspect!(@double.neg_infinity.atan2(1.0), content="-1.5707963267948966") +/// inspect!(1.0.atan2(@double.neg_infinity), content="3.141592653589793") +/// inspect!(@double.infinity.atan2(@double.infinity), content="0.7853981633974483") +/// inspect!(@double.neg_infinity.atan2(@double.neg_infinity), content="-2.356194490192345") +/// inspect!(@double.infinity.atan2(@double.neg_infinity), content="2.356194490192345") +/// inspect!(@double.neg_infinity.atan2(@double.infinity), content="-0.7853981633974483") /// } /// ``` pub fn atan2(self : Double, x : Double) -> Double { diff --git a/double/trig_test.mbt b/double/trig_test.mbt index 36db8cde3..3c96c2983 100644 --- a/double/trig_test.mbt +++ b/double/trig_test.mbt @@ -83,7 +83,8 @@ test "trigonometric functions" { 7.606099243270125, ] fn high_accuracy_test(expect : Double, actual : Double) -> Unit!Error { - assert_eq!((expect - actual).abs() < 1.0e-14, true) + let ulp = (expect.reinterpret_as_int64() - actual.reinterpret_as_int64()).abs() + assert_eq!(ulp <= 1, true) } for i in 0.. Float + acos(Float) -> Float + asin(Float) -> Float + atan(Float) -> Float + atan2(Float, Float) -> Float ceil(Float) -> Float + cos(Float) -> Float default() -> Float exp(Float) -> Float floor(Float) -> Float @@ -30,6 +35,8 @@ impl Float { op_mod(Float, Float) -> Float pow(Float, Float) -> Float round(Float) -> Float + sin(Float) -> Float + tan(Float) -> Float to_be_bytes(Float) -> Bytes to_int(Float) -> Int to_le_bytes(Float) -> Bytes diff --git a/float/trig.mbt b/float/trig.mbt new file mode 100644 index 000000000..a9065db83 --- /dev/null +++ b/float/trig.mbt @@ -0,0 +1,206 @@ +// Copyright 2025 International Digital Economy Academy +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +///| +/// Calculates the sine of a floating-point number in radians. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the sine. +/// +/// Returns the sine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "sin" { +/// inspect!((0.0 : Float).sin(), content="0") +/// inspect!((2.0 : Float).sin(), content="0.9092974066734314") +/// inspect!((-5.0 : Float).sin(), content="0.9589242935180664") +/// inspect!(@float.not_a_number.sin(), content="NaN") +/// inspect!(@float.infinity.sin(), content="NaN") +/// inspect!(@float.neg_infinity.sin(), content="NaN") +/// } +/// ``` +pub fn sin(self : Float) -> Float { + self.to_double().sin().to_float() +} + +///| +/// Calculates the cosine of a floating-point number in radians. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the cosine. +/// +/// Returns the cosine of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "cos" { +/// inspect!((0.0 : Float).cos(), content="1") +/// inspect!((2.0 : Float).cos(), content="-0.416146844625473") +/// inspect!((-5.0 : Float).cos(), content="0.28366219997406006") +/// inspect!(@float.not_a_number.cos(), content="NaN") +/// inspect!(@float.infinity.cos(), content="NaN") +/// inspect!(@float.neg_infinity.cos(), content="NaN") +/// } +/// ``` +pub fn cos(self : Float) -> Float { + self.to_double().cos().to_float() +} + +///| +/// Calculates the tangent of a floating-point number in radians. +/// +/// Parameters: +/// +/// * `x` : The angle in radians for which to calculate the tangent. +/// +/// Returns the tangent of the angle `x`. +/// +/// Example: +/// +/// ```moonbit +/// test "tan" { +/// inspect!((0.0 : Float).tan(), content="0") +/// inspect!((2.0 : Float).tan(), content="-2.185039758682251") +/// inspect!((-5.0 : Float).tan(), content="3.3805150985717773") +/// inspect!(@float.not_a_number.tan(), content="NaN") +/// inspect!(@float.infinity.tan(), content="NaN") +/// inspect!(@float.neg_infinity.tan(), content="NaN") +/// } +/// ``` +pub fn tan(self : Float) -> Float { + self.to_double().tan().to_float() +} + +///| +/// Calculates the arcsine of a floating-point number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arcsine. +/// +/// Returns the arcsine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "asin" { +/// inspect!((0.0 : Float).asin(), content="0") +/// inspect!((1.0 : Float).asin(), content="1.5707963705062866") +/// inspect!((-1.0 : Float).asin(), content="-1.5707963705062866") +/// inspect!(@float.not_a_number.asin(), content="NaN") +/// inspect!(@float.infinity.asin(), content="NaN") +/// inspect!(@float.neg_infinity.asin(), content="NaN") +/// } +/// ``` +pub fn asin(self : Float) -> Float { + self.to_double().asin().to_float() +} + +///| +/// Calculates the arccosine of a floating-point number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arccosine. +/// +/// Returns the arccosine of the number `x`. +/// +/// * Returns NaN if the input is NaN. +/// * Returns NaN if the input is less than -1 or greater than 1. +/// +/// Example: +/// +/// ```moonbit +/// test "acos" { +/// inspect!((0.0 : Float).acos(), content="1.5707963705062866") +/// inspect!((1.0 : Float).acos(), content="0") +/// inspect!((-1.0 : Float).acos(), content="3.1415927410125732") +/// inspect!(@float.not_a_number.acos(), content="NaN") +/// inspect!(@float.infinity.acos(), content="NaN") +/// inspect!(@float.neg_infinity.acos(), content="NaN") +/// } +/// ``` +pub fn acos(self : Float) -> Float { + self.to_double().acos().to_float() +} + +///| +/// Calculates the arctangent of a floating-point number. +/// +/// Parameters: +/// +/// * `x` : The number for which to calculate the arctangent. +/// +/// Returns the arctangent of the number `x`. +/// +/// Example: +/// +/// * Returns NaN if the input is NaN. +/// +/// ```moonbit +/// test "atan" { +/// inspect!((0.0 : Float).atan(), content="0") +/// inspect!((1.0 : Float).atan(), content="0.7853981852531433") +/// inspect!((-1.0 : Float).atan(), content="-0.7853981852531433") +/// inspect!(@float.not_a_number.atan(), content="NaN") +/// inspect!(@float.infinity.atan(), content="1.5707963705062866") +/// inspect!(@float.neg_infinity.atan(), content="-1.5707963705062866") +/// } +/// ``` +pub fn atan(self : Float) -> Float { + self.to_double().atan().to_float() +} + +///| +/// Calculates the arctangent of the quotient of two floating-point numbers. +/// +/// Parameters: +/// +/// * `self` : The numerator of the quotient. +/// * `x` : The denominator of the quotient. +/// +/// Returns the arctangent of the quotient `self / x`. +/// +/// * Returns NaN if self or x is NaN. +/// +/// Example: +/// +/// ```moonbit +/// test "atan2" { +/// inspect!((0.0 : Float).atan2(-1.0), content="3.1415927410125732") +/// inspect!((1.0 : Float).atan2(0.0), content="1.5707963705062866") +/// inspect!((1.0 : Float).atan2(1.0), content="0.7853981852531433") +/// inspect!(@float.not_a_number.atan2(1.0), content="NaN") +/// inspect!((1.0 : Float).atan2(not_a_number), content="NaN") +/// inspect!(@float.infinity.atan2(1.0), content="1.5707963705062866") +/// inspect!((1.0 : Float).atan2(infinity), content="0") +/// inspect!(@float.neg_infinity.atan2(1.0), content="-1.5707963705062866") +/// inspect!((1.0 : Float).atan2(@float.neg_infinity), content="3.1415927410125732") +/// inspect!(@float.infinity.atan2(@float.infinity), content="0.7853981852531433") +/// inspect!(@float.neg_infinity.atan2(@float.neg_infinity), content="-2.356194496154785") +/// inspect!(@float.infinity.atan2(@float.neg_infinity), content="2.356194496154785") +/// inspect!(@float.neg_infinity.atan2(@float.infinity), content="-0.7853981852531433") +/// } +/// ``` +pub fn atan2(self : Float, other : Float) -> Float { + self.to_double().atan2(other.to_double()).to_float() +} diff --git a/float/trig_test.mbt b/float/trig_test.mbt new file mode 100644 index 000000000..cc5564917 --- /dev/null +++ b/float/trig_test.mbt @@ -0,0 +1,88 @@ +// Copyright 2025 International Digital Economy Academy +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +test "floating point trignometry" { + // sin + inspect!((0.0 : Float).sin(), content="0") + inspect!((2.0 : Float).sin(), content="0.9092974066734314") + inspect!((-5.0 : Float).sin(), content="0.9589242935180664") + inspect!(@float.not_a_number.sin(), content="NaN") + inspect!(@float.infinity.sin(), content="NaN") + inspect!(@float.neg_infinity.sin(), content="NaN") + + // cos + inspect!((0.0 : Float).cos(), content="1") + inspect!((2.0 : Float).cos(), content="-0.416146844625473") + inspect!((-5.0 : Float).cos(), content="0.28366219997406006") //< + inspect!(@float.not_a_number.cos(), content="NaN") + inspect!(@float.infinity.cos(), content="NaN") + inspect!(@float.neg_infinity.cos(), content="NaN") + + // tan + inspect!((0.0 : Float).tan(), content="0") + inspect!((2.0 : Float).tan(), content="-2.185039758682251") + inspect!((-5.0 : Float).tan(), content="3.3805150985717773") + inspect!(@float.not_a_number.tan(), content="NaN") + inspect!(@float.infinity.tan(), content="NaN") + inspect!(@float.neg_infinity.tan(), content="NaN") + + // asin + inspect!((0.0 : Float).asin(), content="0") + inspect!((1.0 : Float).asin(), content="1.5707963705062866") + inspect!((-1.0 : Float).asin(), content="-1.5707963705062866") + inspect!(@float.not_a_number.asin(), content="NaN") + inspect!(@float.infinity.asin(), content="NaN") + inspect!(@float.neg_infinity.asin(), content="NaN") + inspect!((0.0 : Float).acos(), content="1.5707963705062866") + inspect!((1.0 : Float).acos(), content="0") + inspect!((-1.0 : Float).acos(), content="3.1415927410125732") //< + inspect!(@float.not_a_number.acos(), content="NaN") + inspect!(@float.infinity.acos(), content="NaN") + inspect!(@float.neg_infinity.acos(), content="NaN") + + // atan + inspect!((0.0 : Float).atan(), content="0") + inspect!((1.0 : Float).atan(), content="0.7853981852531433") + inspect!((-1.0 : Float).atan(), content="-0.7853981852531433") + inspect!(@float.not_a_number.atan(), content="NaN") + inspect!(@float.infinity.atan(), content="1.5707963705062866") + inspect!(@float.neg_infinity.atan(), content="-1.5707963705062866") + + // atan2 + inspect!((0.0 : Float).atan2(-1.0), content="3.1415927410125732") + inspect!((1.0 : Float).atan2(0.0), content="1.5707963705062866") + inspect!((1.0 : Float).atan2(1.0), content="0.7853981852531433") + inspect!(@float.not_a_number.atan2(1.0), content="NaN") + inspect!((1.0 : Float).atan2(not_a_number), content="NaN") + inspect!(@float.infinity.atan2(1.0), content="1.5707963705062866") + inspect!((1.0 : Float).atan2(infinity), content="0") + inspect!(@float.neg_infinity.atan2(1.0), content="-1.5707963705062866") + inspect!( + (1.0 : Float).atan2(@float.neg_infinity), + content="3.1415927410125732", + ) + inspect!(@float.infinity.atan2(@float.infinity), content="0.7853981852531433") + inspect!( + @float.neg_infinity.atan2(@float.neg_infinity), + content="-2.356194496154785", + ) + inspect!( + @float.infinity.atan2(@float.neg_infinity), + content="2.356194496154785", + ) + inspect!( + @float.neg_infinity.atan2(@float.infinity), + content="-0.7853981852531433", + ) +}