1 | module math |
2 | |
3 | // scalbn scales x by FLT_RADIX raised to the power of n, returning the same as: |
4 | // scalbn(x,n) = x * FLT_RADIX ** n |
5 | pub fn scalbn(x f64, n_ int) f64 { |
6 | mut n := n_ |
7 | x1p1023 := f64_from_bits(u64(0x7fe0000000000000)) |
8 | x1p53 := f64_from_bits(u64(0x4340000000000000)) |
9 | x1p_1022 := f64_from_bits(u64(0x0010000000000000)) |
10 | |
11 | mut y := x |
12 | if n > 1023 { |
13 | y *= x1p1023 |
14 | n -= 1023 |
15 | if n > 1023 { |
16 | y *= x1p1023 |
17 | n -= 1023 |
18 | if n > 1023 { |
19 | n = 1023 |
20 | } |
21 | } |
22 | } else if n < -1022 { |
23 | /* |
24 | make sure final n < -53 to avoid double |
25 | rounding in the subnormal range |
26 | */ |
27 | y *= x1p_1022 * x1p53 |
28 | n += 1022 - 53 |
29 | if n < -1022 { |
30 | y *= x1p_1022 * x1p53 |
31 | n += 1022 - 53 |
32 | if n < -1022 { |
33 | n = -1022 |
34 | } |
35 | } |
36 | } |
37 | return y * f64_from_bits(u64((0x3ff + n)) << 52) |
38 | } |