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