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