-
Notifications
You must be signed in to change notification settings - Fork 96
/
Copy pathexp_nonjs.mbt
149 lines (145 loc) · 3.95 KB
/
exp_nonjs.mbt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
// 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.
///|
/// Computes `e` raised to the power of a given number.
///
/// Parameters:
///
/// * `self` : The exponent value to compute `e^x`.
///
/// Returns the value of `e^x`, where `e` is Euler's number (approximately
/// 2.71828).
///
/// Special cases:
///
/// * Returns `x` if `x` is `infinity` (positive infinity)
/// * Returns `0` if `x` is negative infinity
/// * Returns `NaN` if `x` is `NaN`
/// * Returns `1` if `x` is `0`
///
/// Example:
///
/// ```moonbit
/// test "exp" {
/// inspect!(0.0.exp(), content="1")
/// inspect!(1.0.exp(), content="2.718281828459045")
/// inspect!((-1.0).exp(), content="0.36787944117144233")
/// inspect!(not_a_number.exp().is_nan(), content="true")
/// }
/// ```
pub fn exp(self : Double) -> 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 insert_words(ix0 : UInt64, ix1 : UInt64) -> Double {
let mut bits : UInt64 = 0
bits = bits | (ix0 << 32)
bits = bits | ix1
bits.reinterpret_as_double()
}
let mut x = self
let one = 1.0
let halF = [0.5, -0.5]
let o_threshold = 7.09782712893383973096e+02
let u_threshold = -7.45133219101941108420e+02
let ln2HI = [6.93147180369123816490e-01, -6.93147180369123816490e-01]
let ln2LO = [1.90821492927058770002e-10, -1.90821492927058770002e-10]
let invln2 = 1.44269504088896338700e+00
let p1 = 1.66666666666666019037e-01
let p2 = -2.77777777770155933842e-03
let p3 = 6.61375632143793436117e-05
let p4 = -1.65339022054652515390e-06
let p5 = 4.13813679705723846039e-08
let e = 2.718281828459045
let mut hi = 0.0
let mut lo = 0.0
let huge = 1.0e+300
let twom1000 = 9.33263618503218878990e-302
let two1023 = 8.988465674311579539e307
let mut k : Int = 0
let mut hx : UInt = get_high_word(self)
let xsb : Int = ((hx >> 31) & 1).reinterpret_as_int()
hx = hx & 0x7FFFFFFF
if hx >= 0x40862E42 {
if hx >= 0x7FF00000 {
let lx : UInt = get_low_word(self)
if ((hx & 0xFFFFF) | lx) != 0 {
return self + self
} else if xsb == 0 {
return self
} else {
return 0.0
}
}
if self > o_threshold {
return huge * huge
}
if self < u_threshold {
return twom1000 * twom1000
}
}
if hx > 0x3FD62E42 {
if hx < 0x3FF0A2B2 {
if self == 1.0 {
return e
}
hi = self - ln2HI[xsb]
lo = ln2LO[xsb]
k = 1 - xsb - xsb
} else {
k = (invln2 * self + halF[xsb]).to_int()
let t = k.to_double()
hi = self - t * ln2HI[0]
lo = t * ln2LO[0]
}
x = hi - lo
} else if hx < 0x3E300000 {
if huge + x > one {
return one + x
}
} else {
k = 0
}
let t = x * x
let twopk = if k >= -1021 {
insert_words(
(0x3FF00000 + (k.reinterpret_as_uint() << 20).reinterpret_as_int())
.to_int64()
.reinterpret_as_uint64(),
0,
)
} else {
insert_words(
0x3FF00000UL + ((k + 1000).reinterpret_as_uint() << 20).to_uint64(),
0,
)
}
let c = x - t * (p1 + t * (p2 + t * (p3 + t * (p4 + t * p5))))
if k == 0 {
return one - (x * c / (c - 2.0) - x)
}
let y = one - (lo - x * c / (2.0 - c) - hi)
if k >= -1021 {
if k == 1024 {
return y * 2.0 * two1023
} else {
return y * twopk
}
} else {
return y * twopk * twom1000
}
}