Skip to content

Commit 448b7d7

Browse files
Kaida-AmethystKaida-Amethyst
committed
implement fp32 math functions
1 parent d955146 commit 448b7d7

26 files changed

+2067
-205
lines changed

‎float/cbrt.mbt

+58-1
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,21 @@
1212
// See the License for the specific language governing permissions and
1313
// limitations under the License.
1414

15+
// origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c
16+
//
17+
// Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
18+
// Debugged and optimized by Bruce D. Evans.
19+
//
20+
//
21+
// ====================================================
22+
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
23+
//
24+
// Developed at SunPro, a Sun Microsystems, Inc. business.
25+
// Permission to use, copy, modify, and distribute this
26+
// software is freely granted, provided that this notice
27+
// is preserved.
28+
// ====================================================
29+
1530
///|
1631
/// Calculates the cube root of a number.
1732
///
@@ -42,5 +57,47 @@
4257
/// }
4358
/// ```
4459
pub fn cbrt(self : Float) -> Float {
45-
self.to_double().cbrt().to_float()
60+
let b1 : UInt = 709958130 // B1 = (127-127.0/3-0.03306235651)*2**23 */
61+
let b2 : UInt = 642849266 // B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
62+
let mut ui : UInt = self.reinterpret_as_uint()
63+
let mut hx : UInt = ui & 0x7fffffff
64+
if hx >= 0x7f800000 {
65+
// cbrt(NaN,INF) is itself
66+
return self + self
67+
}
68+
69+
// rough cbrt to 5 bits
70+
if hx < 0x00800000 {
71+
// zero or subnormal?
72+
if hx == 0 {
73+
return self
74+
} // cbrt(+-0) is itself
75+
ui = (self * (0x1.0p24 : Float)).reinterpret_as_uint()
76+
hx = ui & 0x7fffffff
77+
hx = hx / 3 + b2
78+
} else {
79+
hx = hx / 3 + b1
80+
}
81+
ui = ui & 0x80000000
82+
ui = ui | hx
83+
84+
//
85+
// First step Newton iteration (solving t*t-self/t == 0) to 16 bits. In
86+
// double precision so that its terms can be arranged for efficiency
87+
// without causing overflow or underflow.
88+
//
89+
let dx = self.to_double()
90+
let t = ui.reinterpret_as_float().to_double()
91+
let r = t * t * t
92+
let t = t * (dx + dx + r) / (dx + r + r)
93+
94+
//
95+
// Second step Newton iteration to 47 bits. In double precision for
96+
// efficiency and accuracy.
97+
//
98+
let r = t * t * t
99+
let t = t * (dx + dx + r) / (dx + r + r)
100+
101+
// rounding to 24 bits is perfect in round-to-nearest mode
102+
t.to_float()
46103
}

‎float/exp.mbt

+117-2
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,22 @@
1212
// See the License for the specific language governing permissions and
1313
// limitations under the License.
1414

15+
// origin: FreeBSD /usr/src/lib/msun/src/e_expf.c
16+
// origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c
17+
//
18+
// Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
19+
//
20+
//
21+
// ====================================================
22+
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
23+
//
24+
// Developed at SunPro, a Sun Microsystems, Inc. business.
25+
// Permission to use, copy, modify, and distribute this
26+
// software is freely granted, provided that this notice
27+
// is preserved.
28+
// ====================================================
29+
//
30+
1531
///|
1632
fn top12(x : Float) -> UInt {
1733
x.reinterpret_as_int().reinterpret_as_uint() >> 20
@@ -142,13 +158,112 @@ pub fn exp(self : Float) -> Float {
142158
/// ```moonbit
143159
/// test "expm1" {
144160
/// inspect!((0 : Float).expm1(), content="0")
145-
/// inspect!((1 : Float).expm1(), content="1.718281865119934")
161+
/// inspect!((1 : Float).expm1(), content="1.7182817459106445")
146162
/// inspect!((-1 : Float).expm1(), content="-0.6321205496788025")
147163
/// inspect!(@float.not_a_number.expm1(), content="NaN")
148164
/// inspect!(@float.infinity.expm1(), content="Infinity")
149165
/// inspect!(@float.neg_infinity.expm1(), content="-1")
150166
/// }
151167
/// ```
152168
pub fn expm1(self : Float) -> Float {
153-
self.to_double().expm1().to_float()
169+
let float_ln2_hi : Float = 6.9314575195e-01 // 0x3f317200
170+
let float_ln2_lo : Float = 1.4286067653e-06 // 0x35bfbe8e
171+
let inv_ln2 : Float = 1.4426950216e+00 // 0x3fb8aa3b
172+
let mut x = self
173+
let q1 : Float = -3.3333212137e-2 // -0x888868.0p-28
174+
let q2 : Float = 1.5807170421e-3 // 0xcf3010.0p-33
175+
let mut hx = x.reinterpret_as_uint()
176+
let sign = (hx >> 31) != 0
177+
hx = hx & 0x7fffffff
178+
179+
// filter out huge and non-finite argument
180+
if hx >= 0x4195b844 {
181+
// if |x|>=27*ln2
182+
if hx > 0x7f800000 {
183+
// NaN
184+
return x
185+
}
186+
if sign {
187+
return -1.0
188+
}
189+
if hx > 0x42b17217 {
190+
x *= (0x1.0p127 : Float)
191+
return x
192+
}
193+
}
194+
let mut k : Int = 0
195+
let mut hi : Float = 0
196+
let mut lo : Float = 0
197+
let mut c : Float = 0
198+
// argument reduction
199+
if hx > 0x3eb17218 {
200+
// if |x| > 0.5 ln2
201+
if hx < 0x3F851592 {
202+
// and |x| < 1.5 ln2
203+
if not(sign) {
204+
hi = x - float_ln2_hi
205+
lo = float_ln2_lo
206+
k = 1
207+
} else {
208+
hi = x + float_ln2_hi
209+
lo = -float_ln2_lo
210+
k = -1
211+
}
212+
} else {
213+
k = (inv_ln2 * x + (if sign { -0.5 } else { 0.5 })).to_int()
214+
let t = k.to_float()
215+
hi = x - t * float_ln2_hi // t*ln2_hi is exact here
216+
lo = t * float_ln2_lo
217+
}
218+
x = hi - lo
219+
c = hi - x - lo
220+
} else if hx < 0x33000000 {
221+
// when |x|<2**-25, return x
222+
//if hx < 0x00800000 {
223+
// force_eval!(x * x);
224+
//}
225+
return x
226+
} else {
227+
k = 0
228+
}
229+
230+
// x is now in primary range
231+
let hfx = (0.5 : Float) * x
232+
let hxs = x * hfx
233+
let r1 = (1.0 : Float) + hxs * (q1 + hxs * q2)
234+
let t = (3.0 : Float) - r1 * hfx
235+
let mut e = hxs * ((r1 - t) / ((6.0 : Float) - x * t))
236+
if k == 0 {
237+
// c is 0
238+
return x - (x * e - hxs)
239+
}
240+
e = x * (e - c) - c
241+
e -= hxs
242+
// exp(x) ~ 2^k (x_reduced - e + 1)
243+
if k == -1 {
244+
return (0.5 : Float) * (x - e) - 0.5
245+
}
246+
if k == 1 {
247+
if x < -0.25 {
248+
return -(2.0 : Float) * (e - (x + 0.5))
249+
}
250+
return (1.0 : Float) + (2.0 : Float) * (x - e)
251+
}
252+
let twopk = ((0x7f + k) << 23).reinterpret_as_float() // 2^k
253+
if not(k is (0..=56)) {
254+
// suffice to return exp(x)-1
255+
let mut y = x - e + 1.0
256+
if k == 128 {
257+
y = y * 2.0 * (0x1.0p127 : Float)
258+
} else {
259+
y = y * twopk
260+
}
261+
return y - 1.0
262+
}
263+
let uf = ((0x7f - k) << 23).reinterpret_as_float() // 2^-k
264+
if k < 23 {
265+
(x - e + ((1.0 : Float) - uf)) * twopk
266+
} else {
267+
(x - (e + uf) + 1.0) * twopk
268+
}
154269
}

‎float/exp_test.mbt

+42-45
Original file line numberDiff line numberDiff line change
@@ -13,52 +13,49 @@
1313
// limitations under the License.
1414

1515
///|
16-
test "exp" {
17-
inspect!((0.5 : Float).exp(), content="1.6487212181091309")
18-
inspect!((1.0 : Float).exp(), content="2.7182817459106445")
19-
inspect!((0.0 : Float).exp(), content="1")
20-
inspect!((2.0 : Float).exp(), content="7.389056205749512")
21-
inspect!((-1.0 : Float).exp(), content="0.3678794503211975")
22-
inspect!(@float.infinity.exp(), content="Infinity")
23-
inspect!(@float.neg_infinity.exp(), content="0")
24-
inspect!((89 : Float).exp(), content="Infinity")
25-
inspect!((-104 : Float).exp(), content="0")
16+
test "expf Test" {
17+
let data : Array[Float] = [
18+
1, 2, 3, 4, 5, 6.5, 7.5, 8.5, 9.5, 10.5, 0.25, 1.25, 2.25, 3.25, 4.25, -0.125,
19+
-1.125, -2.125, -3.125, -4.125, -0.0625, -1.0625, -2.0625, -3.0625, -4.0625,
20+
]
21+
let libm_results : Array[Float] = [
22+
2.7182817459106445, 7.389056205749512, 20.08553695678711, 54.598148345947266,
23+
148.4131622314453, 665.1416625976562, 1808.0423583984375, 4914.76904296875, 13359.7265625,
24+
36315.50390625, 1.2840254306793213, 3.490342855453491, 9.487735748291016, 25.790340423583984,
25+
70.10541534423828, 0.8824968934059143, 0.32465246319770813, 0.1194329708814621,
26+
0.04393693432211876, 0.016163494437932968, 0.9394130706787109, 0.345590740442276,
27+
0.1271357387304306, 0.046770621091127396, 0.017205949872732162,
28+
]
29+
for i in 0..<data.length() {
30+
let x = data[i]
31+
let res = x.exp()
32+
let libm_res = libm_results[i]
33+
assert_true!(
34+
iabs(res.reinterpret_as_int() - libm_res.reinterpret_as_int()) <= 2,
35+
)
36+
}
2637
}
2738

2839
///|
29-
test "exp - quickcheck" {
30-
let double : Array[Double] = @quickcheck.samples(4)
31-
inspect!(
32-
double,
33-
content="[0.12125190562607557, 0.2656190379663609, 0.7248184591724188, 0.8261703932281662]",
34-
)
35-
inspect!(
36-
double.map(fn(x) { x.exp() }),
37-
content="[1.1289092551447675, 1.3042380988330389, 2.064356300995553, 2.2845530266177425]",
38-
)
39-
let float = double.map(Double::to_float)
40-
inspect!(
41-
float,
42-
content="[0.1212519034743309, 0.26561903953552246, 0.7248184680938721, 0.8261703848838806]",
43-
)
44-
inspect!(
45-
float.map(fn(x) { x.exp() }),
46-
content="[1.1289092302322388, 1.3042380809783936, 2.0643563270568848, 2.284553050994873]",
47-
)
48-
}
49-
50-
///|
51-
test "expm1" {
52-
inspect!((0 : Float).expm1(), content="0")
53-
inspect!((1 : Float).expm1(), content="1.718281865119934")
54-
inspect!((2 : Float).expm1(), content="6.389056205749512")
55-
inspect!((3 : Float).expm1(), content="19.08553695678711")
56-
inspect!((4 : Float).expm1(), content="53.598148345947266")
57-
inspect!((5 : Float).expm1(), content="147.4131622314453")
58-
inspect!((-0.5 : Float).expm1(), content="-0.39346933364868164")
59-
inspect!((-1 : Float).expm1(), content="-0.6321205496788025")
60-
inspect!((-2 : Float).expm1(), content="-0.8646647334098816")
61-
inspect!(@float.not_a_number.expm1(), content="NaN")
62-
inspect!(@float.infinity.expm1(), content="Infinity")
63-
inspect!(@float.neg_infinity.expm1(), content="-1")
40+
test "expm1f Test" {
41+
let data : Array[Float] = [
42+
1, 2, 3, 4, 5, 6.5, 7.5, 8.5, 9.5, 10.5, 0.25, 1.25, 2.25, 3.25, 4.25, -0.125,
43+
-1.125, -2.125, -3.125, -4.125, -0.0625, -1.0625, -2.0625, -3.0625, -4.0625,
44+
]
45+
let libm_results : Array[Float] = [
46+
1.718281865119934, 6.389056205749512, 19.08553695678711, 53.598148345947266,
47+
147.4131622314453, 664.1416625976562, 1807.0423583984375, 4913.76904296875, 13358.7265625,
48+
36314.50390625, 0.2840254306793213, 2.490342855453491, 8.487735748291016, 24.790340423583984,
49+
69.10541534423828, -0.1175030991435051, -0.6753475069999695, -0.8805670142173767,
50+
-0.9560630917549133, -0.9838365316390991, -0.06058693677186966, -0.6544092297554016,
51+
-0.8728642463684082, -0.9532293677330017, -0.9827940464019775,
52+
]
53+
for i in 0..<data.length() {
54+
let x = data[i]
55+
let res = x.expm1()
56+
let libm_res = libm_results[i]
57+
assert_true!(
58+
iabs(res.reinterpret_as_int() - libm_res.reinterpret_as_int()) <= 2,
59+
)
60+
}
6461
}

‎float/float.mbti

+3
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,8 @@ fn pow(Float, Float) -> Float
6363

6464
fn round(Float) -> Float
6565

66+
fn scalbn(Float, Int) -> Float
67+
6668
fn sin(Float) -> Float
6769

6870
fn sinh(Float) -> Float
@@ -105,6 +107,7 @@ impl Float {
105107
ln_1p(Float) -> Float
106108
pow(Float, Float) -> Float
107109
round(Float) -> Float
110+
scalbn(Float, Int) -> Float
108111
sin(Float) -> Float
109112
sinh(Float) -> Float
110113
tan(Float) -> Float

0 commit comments

Comments
 (0)