@@ -53,71 +53,38 @@ object Math {
5353
5454// Wasm intrinsic
5555def rint (a : scala.Double ): scala.Double = {
56- /* Is the integer-valued `x` odd? Fused by hand of `(x.toLong & 1L) != 0L`.
57- * Corner cases: returns false for Infinities and NaN.
58- */
59- @ inlinedef isOdd (x : scala.Double ): scala.Boolean =
60- (x.asInstanceOf [js.Dynamic ]& 1 .asInstanceOf [js.Dynamic ]).asInstanceOf [Int ]!= 0
61-
62- /* js.Math.round(a) does *almost* what we want. It rounds to nearest,
63- * breaking ties *up*. We need to break ties to *even*. So we need to
64- * detect ties, and for them, detect if we rounded to odd instead of even.
65- *
66- * The reasons why the apparently simple algorithm below works are subtle,
67- * and vary a lot depending on the range of `a`:
68- *
69- * - a is NaN
70- * r is NaN, then the == is false
71- * -> return r
72- *
73- * - a is +-Infinity
74- * r == a, then == is true! but isOdd(r) is false
75- * -> return r
76- *
77- * - 2**53 <= abs(a) < Infinity
78- * r == a, r - 0.5 rounds back to a so == is true!
79- * fortunately, isOdd(r) is false because all a >= 2**53 are even
80- * -> return r
56+ /* We apply the technique described in Section II of
57+ * Claude-Pierre Jeannerod, Jean-Michel Muller, Paul Zimmermann.
58+ * On various ways to split a floating-point number.
59+ * ARITH 2018 - 25th IEEE Symposium on Computer Arithmetic,
60+ * Jun 2018, Amherst (MA), United States.
61+ * pp.53-60, 10.1109/ARITH.2018.8464793. hal-01774587v2
62+ * available at
63+ * https://hal.inria.fr/hal-01774587v2/document
64+ * with β = 2, p = 53, and C = 2^(p-1) = 2^52.
8165 *
82- * - 2**52 <= abs(a) < 2**53
83- * r == a (because all a's are integers in that range)
84- * - a is odd
85- * r - 0.5 rounds down (towards even) to r - 1.0
86- * so a == r - 0.5 is false
87- * -> return r
88- * - a is even
89- * r - 0.5 rounds back up! (towards even) to r
90- * so a == r - 0.5 is true!
91- * but, isOdd(r) is false
92- * -> return r
66+ * That is only valid for values x <= 2^52. Fortunately, all values that
67+ * are >= 2^52 are already integers, so we can return them as is.
9368 *
94- * - 0.5 < abs(a) < 2**52
95- * then -2**52 + 0.5 <= a <= 2**52 - 0.5 (because values in-between are not representable)
96- * since Math.round rounds *up* on ties, r is an integer in the range (-2**52, 2**52]
97- * r - 0.5 is therefore lossless
98- * so a == r - 0.5 accurately detects ties, and isOdd(r) breaks ties
99- * -> return `r`` or `r - 1.0`
100- *
101- * - a == +0.5
102- * r == 1.0
103- * a == r - 0.5 is true and isOdd(r) is true
104- * -> return `r - 1.0`, which is +0.0
105- *
106- * - a == -0.5
107- * r == -0.0
108- * a == r - 0.5 is true and isOdd(r) is false
109- * -> return `r`, which is -0.0
110- *
111- * - 0.0 <= abs(a) < 0.5
112- * r == 0.0 with the same sign as a
113- * a == r - 0.5 is false
114- * -> return r
69+ * We cannot use "the 1.5 trick" with C = 2^(p-1) + 2^(p-2) to handle
70+ * negative numbers, because that would further reduce the range of valid
71+ * `x` to maximum 2^51, although we actually need it up to 2^52. Therefore,
72+ * we have a separate branch for negative numbers. This also allows to
73+ * gracefully deal with the fact that we need to return -0.0 for values in
74+ * the range [-0.5,-0.0).
11575*/
116- val r = js.Math .round(a)
117- if ((a== r- 0.5 )&& isOdd(r))
118- r- 1.0
119- else
120- r
76+ val C = 4503599627370496.0 // 2^52
77+ if (a> 0 ) {
78+ if (a>= C ) a
79+ else (C + a)- C
80+ }else if (a< 0 ) {
81+ // do not "optimize" as `C - (C - a)`, as it would return +0.0 where it should return -0.0
82+ if (a<= - C ) a
83+ else - ((C - a)- C )
84+ }else {
85+ // Handling zeroes here avoids the need to distinguish +0.0 from -0.0
86+ a// 0.0, -0.0 and NaN
87+ }
12188 }
12289
12390@ inlinedef round (a : scala.Float ): scala.Int = js.Math .round(a).toInt