1 #include <math.h> 2 #include <stdint.h> 3 scalbn(double x,int n)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 } 34