1#include <math.h>
2#include <stdint.h>
3
4double scalbn(double x, int n)
5{
6 union {double f; uint64_t i;} u;
7 double_t y = x;
8
9 if (n > 1023) {
10 y *= 0x1p1023;
11 n -= 1023;
12 if (n > 1023) {
13 y *= 0x1p1023;
14 n -= 1023;
15 if (n > 1023)
16 n = 1023;
17 }
18 } else if (n < -1022) {
19 /* make sure final n < -53 to avoid double
20 rounding in the subnormal range */
21 y *= 0x1p-1022 * 0x1p53;
22 n += 1022 - 53;
23 if (n < -1022) {
24 y *= 0x1p-1022 * 0x1p53;
25 n += 1022 - 53;
26 if (n < -1022)
27 n = -1022;
28 }
29 }
30 u.i = (uint64_t)(0x3ff+n)<<52;
31 x = y * u.f;
32 return x;
33}