]>
Commit | Line | Data |
---|---|---|
f964490f | 1 | /* s_frexpl.c -- long double version of s_frexp.c. |
f964490f RM |
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 | ||
1ed0291c RH |
29 | #include <math.h> |
30 | #include <math_private.h> | |
f964490f RM |
31 | #include <math_ldbl_opt.h> |
32 | ||
8db21882 | 33 | long double __frexpl(long double x, int *eptr) |
f964490f | 34 | { |
aa5f0ff1 AM |
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; | |
88283451 | 49 | return x + x; |
aa5f0ff1 AM |
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; | |
765714ca | 72 | |
aa5f0ff1 AM |
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 | { | |
706688aa | 81 | hx += 1LL << 52; |
aa5f0ff1 | 82 | expon -= 1; |
f964490f | 83 | } |
f964490f | 84 | |
aa5f0ff1 AM |
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; | |
f964490f | 91 | |
aa5f0ff1 AM |
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. */ | |
706688aa | 128 | hx -= 1LL << 52; |
aa5f0ff1 AM |
129 | expon += 1; |
130 | } | |
131 | } | |
132 | explo = 0; | |
133 | } | |
134 | lx |= (explo << 52) | ixl; | |
135 | } | |
f964490f | 136 | |
aa5f0ff1 AM |
137 | INSERT_WORDS64 (xhi, hx); |
138 | INSERT_WORDS64 (xlo, lx); | |
139 | x = ldbl_pack (xhi, xlo); | |
140 | *eptr = expon; | |
141 | return x; | |
f964490f | 142 | } |
a109996e | 143 | #if IS_IN (libm) |
f964490f RM |
144 | long_double_symbol (libm, __frexpl, frexpl); |
145 | #else | |
146 | long_double_symbol (libc, __frexpl, frexpl); | |
147 | #endif |