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