1 /* s_frexpl.c -- long double version of s_frexp.c.
2 */
3
4 /*
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
7 *
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
11 * is preserved.
12 * ====================================================
13 */
14
15 #if defined(LIBM_SCCS) && !defined(lint)
16 static char rcsid[] = "$NetBSD: $";
17 #endif
18
19 /*
20 * for non-zero x
21 * x = frexpl(arg,&exp);
22 * return a long double fp quantity x such that 0.5 <= |x| <1.0
23 * and the corresponding binary exponent "exp". That is
24 * arg = x*2^exp.
25 * If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg
26 * with *exp=0.
27 */
28
29 #include <math.h>
30 #include <math_private.h>
31 #include <math_ldbl_opt.h>
32
__frexpl(long double x,int * eptr)33 long double __frexpl(long double x, int *eptr)
34 {
35 uint64_t hx, lx, ix, ixl;
36 int64_t explo, expon;
37 double xhi, xlo;
38
39 ldbl_unpack (x, &xhi, &xlo);
40 EXTRACT_WORDS64 (hx, xhi);
41 EXTRACT_WORDS64 (lx, xlo);
42 ixl = 0x7fffffffffffffffULL & lx;
43 ix = 0x7fffffffffffffffULL & hx;
44 expon = 0;
45 if (ix >= 0x7ff0000000000000ULL || ix == 0)
46 {
47 /* 0,inf,nan. */
48 *eptr = expon;
49 return x + x;
50 }
51 expon = ix >> 52;
52 if (expon == 0)
53 {
54 /* Denormal high double, the low double must be 0.0. */
55 int cnt;
56
57 /* Normalize. */
58 if (sizeof (ix) == sizeof (long))
59 cnt = __builtin_clzl (ix);
60 else if ((ix >> 32) != 0)
61 cnt = __builtin_clzl ((long) (ix >> 32));
62 else
63 cnt = __builtin_clzl ((long) ix) + 32;
64 cnt = cnt - 12;
65 expon -= cnt;
66 ix <<= cnt + 1;
67 }
68 expon -= 1022;
69 ix &= 0x000fffffffffffffULL;
70 hx &= 0x8000000000000000ULL;
71 hx |= (1022LL << 52) | ix;
72
73 if (ixl != 0)
74 {
75 /* If the high double is an exact power of two and the low
76 double has the opposite sign, then the exponent calculated
77 from the high double is one too big. */
78 if (ix == 0
79 && (int64_t) (hx ^ lx) < 0)
80 {
81 hx += 1LL << 52;
82 expon -= 1;
83 }
84
85 explo = ixl >> 52;
86 if (explo == 0)
87 {
88 /* The low double started out as a denormal. Normalize its
89 mantissa and adjust the exponent. */
90 int cnt;
91
92 if (sizeof (ixl) == sizeof (long))
93 cnt = __builtin_clzl (ixl);
94 else if ((ixl >> 32) != 0)
95 cnt = __builtin_clzl ((long) (ixl >> 32));
96 else
97 cnt = __builtin_clzl ((long) ixl) + 32;
98 cnt = cnt - 12;
99 explo -= cnt;
100 ixl <<= cnt + 1;
101 }
102
103 /* With variable precision we can't assume much about the
104 magnitude of the returned low double. It may even be a
105 denormal. */
106 explo -= expon;
107 ixl &= 0x000fffffffffffffULL;
108 lx &= 0x8000000000000000ULL;
109 if (explo <= 0)
110 {
111 /* Handle denormal low double. */
112 if (explo > -52)
113 {
114 ixl |= 1LL << 52;
115 ixl >>= 1 - explo;
116 }
117 else
118 {
119 ixl = 0;
120 lx = 0;
121 if ((hx & 0x7ff0000000000000ULL) == (1023LL << 52))
122 {
123 /* Oops, the adjustment we made above for values a
124 little smaller than powers of two turned out to
125 be wrong since the returned low double will be
126 zero. This can happen if the input was
127 something weird like 0x1p1000 - 0x1p-1000. */
128 hx -= 1LL << 52;
129 expon += 1;
130 }
131 }
132 explo = 0;
133 }
134 lx |= (explo << 52) | ixl;
135 }
136
137 INSERT_WORDS64 (xhi, hx);
138 INSERT_WORDS64 (xlo, lx);
139 x = ldbl_pack (xhi, xlo);
140 *eptr = expon;
141 return x;
142 }
143 #if IS_IN (libm)
144 long_double_symbol (libm, __frexpl, frexpl);
145 #else
146 long_double_symbol (libc, __frexpl, frexpl);
147 #endif
148