1 /* 2 * Copyright (c) 2005-2020 Rich Felker, et al. 3 * 4 * SPDX-License-Identifier: MIT 5 * 6 * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT 7 * for all contributors to musl. 8 */ 9 #include <float.h> 10 #include <math.h> 11 #include <stdint.h> 12 13 double scalbn(double x, int n) 14 { 15 union {double f; uint64_t i;} u; 16 double_t y = x; 17 18 if (n > 1023) { 19 y *= 0x1p1023; 20 n -= 1023; 21 if (n > 1023) { 22 y *= 0x1p1023; 23 n -= 1023; 24 if (n > 1023) 25 n = 1023; 26 } 27 } else if (n < -1022) { 28 /* make sure final n < -53 to avoid double 29 rounding in the subnormal range */ 30 y *= 0x1p-1022 * 0x1p53; 31 n += 1022 - 53; 32 if (n < -1022) { 33 y *= 0x1p-1022 * 0x1p53; 34 n += 1022 - 53; 35 if (n < -1022) 36 n = -1022; 37 } 38 } 39 u.i = (uint64_t)(0x3ff+n)<<52; 40 x = y * u.f; 41 return x; 42 } 43 44 #if (LDBL_MANT_DIG == 53) && !defined(scalbn) 45 __weak_reference(scalbn, ldexpl); 46 __weak_reference(scalbn, scalbnl); 47 #endif 48