1 | /* Split a double into fraction and mantissa, for hexadecimal printf. |
2 | Copyright (C) 2007, 2009-2019 Free Software Foundation, Inc. |
3 | |
4 | This program is free software: you can redistribute it and/or modify |
5 | it under the terms of the GNU General Public License as published by |
6 | the Free Software Foundation; either version 3 of the License, or |
7 | (at your option) any later version. |
8 | |
9 | This program is distributed in the hope that it will be useful, |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
12 | GNU General Public License for more details. |
13 | |
14 | You should have received a copy of the GNU General Public License |
15 | along with this program. If not, see <https://www.gnu.org/licenses/>. */ |
16 | |
17 | #if ! defined USE_LONG_DOUBLE |
18 | # include <config.h> |
19 | #endif |
20 | |
21 | /* Specification. */ |
22 | #ifdef USE_LONG_DOUBLE |
23 | # include "printf-frexpl.h" |
24 | #else |
25 | # include "printf-frexp.h" |
26 | #endif |
27 | |
28 | #include <float.h> |
29 | #include <math.h> |
30 | #ifdef USE_LONG_DOUBLE |
31 | # include "fpucw.h" |
32 | #endif |
33 | |
34 | /* This file assumes FLT_RADIX = 2. If FLT_RADIX is a power of 2 greater |
35 | than 2, or not even a power of 2, some rounding errors can occur, so that |
36 | then the returned mantissa is only guaranteed to be <= 2.0, not < 2.0. */ |
37 | |
38 | #ifdef USE_LONG_DOUBLE |
39 | # define FUNC printf_frexpl |
40 | # define DOUBLE long double |
41 | # define MIN_EXP LDBL_MIN_EXP |
42 | # if HAVE_FREXPL_IN_LIBC && HAVE_LDEXPL_IN_LIBC |
43 | # define USE_FREXP_LDEXP |
44 | # define FREXP frexpl |
45 | # define LDEXP ldexpl |
46 | # endif |
47 | # define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING |
48 | # define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING () |
49 | # define END_ROUNDING() END_LONG_DOUBLE_ROUNDING () |
50 | # define L_(literal) literal##L |
51 | #else |
52 | # define FUNC printf_frexp |
53 | # define DOUBLE double |
54 | # define MIN_EXP DBL_MIN_EXP |
55 | # if HAVE_FREXP_IN_LIBC && HAVE_LDEXP_IN_LIBC |
56 | # define USE_FREXP_LDEXP |
57 | # define FREXP frexp |
58 | # define LDEXP ldexp |
59 | # endif |
60 | # define DECL_ROUNDING |
61 | # define BEGIN_ROUNDING() |
62 | # define END_ROUNDING() |
63 | # define L_(literal) literal |
64 | #endif |
65 | |
66 | DOUBLE |
67 | FUNC (DOUBLE x, int *expptr) |
68 | { |
69 | int exponent; |
70 | DECL_ROUNDING |
71 | |
72 | BEGIN_ROUNDING (); |
73 | |
74 | #ifdef USE_FREXP_LDEXP |
75 | /* frexp and ldexp are usually faster than the loop below. */ |
76 | x = FREXP (x, &exponent); |
77 | |
78 | x = x + x; |
79 | exponent -= 1; |
80 | |
81 | if (exponent < MIN_EXP - 1) |
82 | { |
83 | x = LDEXP (x, exponent - (MIN_EXP - 1)); |
84 | exponent = MIN_EXP - 1; |
85 | } |
86 | #else |
87 | { |
88 | /* Since the exponent is an 'int', it fits in 64 bits. Therefore the |
89 | loops are executed no more than 64 times. */ |
90 | DOUBLE pow2[64]; /* pow2[i] = 2^2^i */ |
91 | DOUBLE powh[64]; /* powh[i] = 2^-2^i */ |
92 | int i; |
93 | |
94 | exponent = 0; |
95 | if (x >= L_(1.0)) |
96 | { |
97 | /* A nonnegative exponent. */ |
98 | { |
99 | DOUBLE pow2_i; /* = pow2[i] */ |
100 | DOUBLE powh_i; /* = powh[i] */ |
101 | |
102 | /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, |
103 | x * 2^exponent = argument, x >= 1.0. */ |
104 | for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); |
105 | ; |
106 | i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) |
107 | { |
108 | if (x >= pow2_i) |
109 | { |
110 | exponent += (1 << i); |
111 | x *= powh_i; |
112 | } |
113 | else |
114 | break; |
115 | |
116 | pow2[i] = pow2_i; |
117 | powh[i] = powh_i; |
118 | } |
119 | } |
120 | /* Here 1.0 <= x < 2^2^i. */ |
121 | } |
122 | else |
123 | { |
124 | /* A negative exponent. */ |
125 | { |
126 | DOUBLE pow2_i; /* = pow2[i] */ |
127 | DOUBLE powh_i; /* = powh[i] */ |
128 | |
129 | /* Invariants: pow2_i = 2^2^i, powh_i = 2^-2^i, |
130 | x * 2^exponent = argument, x < 1.0, exponent >= MIN_EXP - 1. */ |
131 | for (i = 0, pow2_i = L_(2.0), powh_i = L_(0.5); |
132 | ; |
133 | i++, pow2_i = pow2_i * pow2_i, powh_i = powh_i * powh_i) |
134 | { |
135 | if (exponent - (1 << i) < MIN_EXP - 1) |
136 | break; |
137 | |
138 | exponent -= (1 << i); |
139 | x *= pow2_i; |
140 | if (x >= L_(1.0)) |
141 | break; |
142 | |
143 | pow2[i] = pow2_i; |
144 | powh[i] = powh_i; |
145 | } |
146 | } |
147 | /* Here either x < 1.0 and exponent - 2^i < MIN_EXP - 1 <= exponent, |
148 | or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1. */ |
149 | |
150 | if (x < L_(1.0)) |
151 | /* Invariants: x * 2^exponent = argument, x < 1.0 and |
152 | exponent - 2^i < MIN_EXP - 1 <= exponent. */ |
153 | while (i > 0) |
154 | { |
155 | i--; |
156 | if (exponent - (1 << i) >= MIN_EXP - 1) |
157 | { |
158 | exponent -= (1 << i); |
159 | x *= pow2[i]; |
160 | if (x >= L_(1.0)) |
161 | break; |
162 | } |
163 | } |
164 | |
165 | /* Here either x < 1.0 and exponent = MIN_EXP - 1, |
166 | or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1. */ |
167 | } |
168 | |
169 | /* Invariants: x * 2^exponent = argument, and |
170 | either x < 1.0 and exponent = MIN_EXP - 1, |
171 | or 1.0 <= x < 2^2^i and exponent >= MIN_EXP - 1. */ |
172 | while (i > 0) |
173 | { |
174 | i--; |
175 | if (x >= pow2[i]) |
176 | { |
177 | exponent += (1 << i); |
178 | x *= powh[i]; |
179 | } |
180 | } |
181 | /* Here either x < 1.0 and exponent = MIN_EXP - 1, |
182 | or 1.0 <= x < 2.0 and exponent >= MIN_EXP - 1. */ |
183 | } |
184 | #endif |
185 | |
186 | END_ROUNDING (); |
187 | |
188 | *expptr = exponent; |
189 | return x; |
190 | } |
191 | |