fix(double): op_mod for subnormal double (moonbitlang#1303)

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.
 

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()
 }

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))
 }