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 | |
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 | |