1#include <float.h>
2#include <math.h>
3#include <stdint.h>
4
5double 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