1#include <math.h>
2#include <stdint.h>
3
4double scalbn(double x, int n) {
5    union {
6        double f;
7        uint64_t i;
8    } u;
9    double_t y = x;
10
11    if (n > 1023) {
12        y *= 0x1p1023;
13        n -= 1023;
14        if (n > 1023) {
15            y *= 0x1p1023;
16            n -= 1023;
17            if (n > 1023)
18                n = 1023;
19        }
20    } else if (n < -1022) {
21        y *= 0x1p-1022;
22        n += 1022;
23        if (n < -1022) {
24            y *= 0x1p-1022;
25            n += 1022;
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}
34