From 128116254ac4601d14daa4e2c0e6b280e6f5beb8 Mon Sep 17 00:00:00 2001 From: liuly Date: Wed, 25 Dec 2024 15:37:13 +0800 Subject: [PATCH 1/3] fix(double): op_mod for subnormal double (#1303) --- NOTICE | 5 +- double/mod_nonjs.mbt | 121 +++++++++++++++++++++++-------------------- double/mod_test.mbt | 24 +++++++-- 3 files changed, 90 insertions(+), 60 deletions(-) diff --git a/NOTICE b/NOTICE index d522bde9c..3e28daab8 100644 --- a/NOTICE +++ b/NOTICE @@ -71,13 +71,14 @@ THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -File `float/exp.mbt`, `float/log.mbt`, and `double/scalbn.mbt` is adapted from +File `float/exp.mbt`, `float/log.mbt`, `double/mod_nonjs.mbt` and `double/scalbn.mbt` is adapted from [musl](https://www.musl-libc.org/). Specifically: - `float/log.mbt` is adapted from `src/math/logf.c`, `src/math/logf_data.h`, and `src/math/logf_data.c`. - `float/exp.mbt` is adapted from `src/math/expf.c`, `src/math/exp2f_data.h`, and `src/math/exp2f_data.c`. +- `double/mod_nonjs.mbt` is adapted from `src/math/fmod.c`. - `double/scalbn.mbt` is adapted from `src/math/scalbn.c`. `float/log.mbt` and `float/exp.mbt` files are Copyright (c) 2017-2018 Arm Limited and licensed under the MIT @@ -104,7 +105,7 @@ COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -`double/scalbn.mbt` is licensed under the following license: +`double/mod_nonjs.mbt` and `double/scalbn.mbt` is licensed under the following license: Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. diff --git a/double/mod_nonjs.mbt b/double/mod_nonjs.mbt index cbea93441..3f6fa1964 100644 --- a/double/mod_nonjs.mbt +++ b/double/mod_nonjs.mbt @@ -13,69 +13,80 @@ // limitations under the License. ///| -fn copy_sign(x : Double, y : Double) -> Double { - if x < 0 { - -y.abs() - } else { - y.abs() +pub fn op_mod(self : Double, other : Double) -> Double { + let x = self + let y = other + let mut uint64_x = x.reinterpret_as_uint64() + let mut uint64_y = y.reinterpret_as_uint64() + let mut ex = ((uint64_x >> 52) & 0x7FF).to_int() + let mut ey = ((uint64_y >> 52) & 0x7FF).to_int() + let sign_x = uint64_x >> 63 + let mut i : UInt64 = 0 + if (uint64_y << 1) == 0 || y.is_nan() || ex == 0x7ff { + return x * y / (x * y) } -} - -///| -fn ldexp(frac : Double, exp_ : Double) -> Double { - if frac == 0 { - 0 - } else if frac.is_inf() || frac.is_nan() { - frac - } else { - let (frac, e) = normalize(frac) - let mut exp = exp_ - exp += e.to_double() - let mut x = frac.reinterpret_as_uint64() - exp += (((x >> 52).to_int() & 0x7FF) - 1023).to_double() - if exp < -1075 { - copy_sign(0, frac) - } else if exp > 1023 { - if frac < 0 { - neg_infinity - } else { - infinity - } - } else { - let mut m : Double = 1 - if exp < -1022 { - exp += 53 - m = 1.0 / (1 << 53).to_double() - } - x = x & (0x7FFUL << 52).lnot() - x = x | ((exp + 1023.0).to_int64().reinterpret_as_uint64() << 52) - m * x.reinterpret_as_double() + if (uint64_x << 1) <= (uint64_y << 1) { + if (uint64_x << 1) == (uint64_y << 1) { + return 0.0 * x } + return x } -} -///| -pub fn op_mod(self : Double, other : Double) -> Double { - if other == 0 || self.is_inf() || self.is_nan() || other.is_nan() { - not_a_number + // normalize x and y + if ex == 0 { + i = uint64_x << 12 + while (i >> 63) == 0 { + ex -= 1 + i = i << 1 + } + uint64_x = uint64_x << (-ex + 1) } else { - let y = other.abs() - let (yfr, yexp) = frexp(y) - let mut r = self - if self < 0 { - r = -self + uint64_x = uint64_x & (18446744073709551615UL >> 12) + uint64_x = uint64_x | (1UL << 52) + } + if ey == 0 { + i = uint64_y << 12 + while (i >> 63) == 0 { + ey -= 1 + i = i << 1 } - while r >= y { - let (rfr, rexp_) = frexp(r) - let mut rexp = rexp_ - if rfr < yfr { - rexp -= 1 + uint64_y = uint64_y << (-ey + 1) + } else { + uint64_y = uint64_y & (18446744073709551615UL >> 12) + uint64_y = uint64_y | (1UL << 52) + } + + // x mod y + while ex > ey { + i = uint64_x - uint64_y + if (i >> 63) == 0 { + if i == 0 { + return 0.0 * x } - r = r - ldexp(y, (rexp - yexp).to_double()) + uint64_x = i } - if self < 0 { - r = -r + uint64_x = uint64_x << 1 + ex -= 1 + } + i = uint64_x - uint64_y + if (i >> 63) == 0 { + if i == 0 { + return 0.0 * x } - r + uint64_x = i + } + while (uint64_x >> 52) == 0 { + uint64_x = uint64_x << 1 + ex -= 1 + } + + // scale result + if ex > 0 { + uint64_x = uint64_x - (1UL << 52) + uint64_x = uint64_x | (ex.to_uint64() << 52) + } else { + uint64_x = uint64_x >> (-ex + 1) } + uint64_x = uint64_x | (sign_x << 63) + uint64_x.reinterpret_as_double() } diff --git a/double/mod_test.mbt b/double/mod_test.mbt index d167aa51f..a5aa1c865 100644 --- a/double/mod_test.mbt +++ b/double/mod_test.mbt @@ -13,15 +13,33 @@ // limitations under the License. ///| +/// Check if two floating point numbers are equal within a small epsilon. +/// +/// Only used for testing when the result can't be precisely represented +/// by IEEE 754 double. fn check(x : Double, y : Double) -> Bool { (x - y).abs() < 1.0e-15 } test "mod" { + // normal numbers inspect!(check(7.5 % 2.3, 0.6), content="true") - inspect!(check(5.75 % 5.75, 0), content="true") inspect!(check(-3.6 % 1.4, -0.8), content="true") - inspect!(check(15.25 % 4.5, 1.75), content="true") inspect!(check(0.7 % 0.2, 0.1), content="true") - inspect!(check(-8.4 % 3.2, -2), content="true") + inspect!(5.75 % 5.75, content="0") + inspect!(15.25 % 4.5, content="1.75") + inspect!(-8.4 % 3.2, content="-2") + // subnormal numbers + inspect!(5.0e-324 % 5.0e-324, content="0") + inspect!(0.5 % 1.5e-323, content="1e-323") + // inf + inspect!(1.0 % @double.infinity, content="1") + assert_true!(@double.is_nan(@double.infinity % @double.infinity)) + assert_true!(@double.is_nan(@double.infinity % 1)) + // division by zero + assert_true!(@double.is_nan(1.0 % 0)) + assert_true!(@double.is_nan(1.0 % -0)) + // nan + assert_true!(@double.is_nan(@double.not_a_number % 1.0)) + assert_true!(@double.is_nan(1.0 % @double.not_a_number)) } From 11f86f552c65fa341b8e5314e48e6aa1dca57625 Mon Sep 17 00:00:00 2001 From: liuly Date: Wed, 25 Dec 2024 16:04:57 +0800 Subject: [PATCH 2/3] fix typo --- NOTICE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NOTICE b/NOTICE index 3e28daab8..89c0b1075 100644 --- a/NOTICE +++ b/NOTICE @@ -105,7 +105,7 @@ COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -`double/mod_nonjs.mbt` and `double/scalbn.mbt` is licensed under the following license: +`double/mod_nonjs.mbt` and `double/scalbn.mbt` are licensed under the following license: Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. From a77d0125163dad4ff73ff362e397c10a43cdf0ad Mon Sep 17 00:00:00 2001 From: liuly Date: Wed, 25 Dec 2024 16:13:55 +0800 Subject: [PATCH 3/3] improve coverage --- double/mod_test.mbt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/double/mod_test.mbt b/double/mod_test.mbt index a5aa1c865..ed5d72c48 100644 --- a/double/mod_test.mbt +++ b/double/mod_test.mbt @@ -29,9 +29,12 @@ test "mod" { inspect!(5.75 % 5.75, content="0") inspect!(15.25 % 4.5, content="1.75") inspect!(-8.4 % 3.2, content="-2") + inspect!(8.0 % 4.0, content="0") + inspect!(15.0 % 3.0, content="0") // subnormal numbers inspect!(5.0e-324 % 5.0e-324, content="0") inspect!(0.5 % 1.5e-323, content="1e-323") + inspect!(4.0e-310 % 3.0e-310, content="1e-310") // inf inspect!(1.0 % @double.infinity, content="1") assert_true!(@double.is_nan(@double.infinity % @double.infinity))