#include <math.h> | |
#include <stdint.h> | |
double scalbn(double x, int n) { | |
union { | |
double f; | |
uint64_t i; | |
} u; | |
double_t y = x; | |
if (n > 1023) { | |
y *= 0x1p1023; | |
n -= 1023; | |
if (n > 1023) { | |
y *= 0x1p1023; | |
n -= 1023; | |
if (n > 1023) | |
n = 1023; | |
} | |
} else if (n < -1022) { | |
y *= 0x1p-1022; | |
n += 1022; | |
if (n < -1022) { | |
y *= 0x1p-1022; | |
n += 1022; | |
if (n < -1022) | |
n = -1022; | |
} | |
} | |
u.i = (uint64_t)(0x3ff + n) << 52; | |
x = y * u.f; | |
return x; | |
} |