Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_PrintDouble.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Teuchos: Common Tools Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
43#include "Teuchos_BigUInt.hpp"
44
45#include <cstring>
46
47namespace Teuchos {
48
49namespace {
50
51int ndigits_for(std::uint32_t x) {
52 int n = 0;
53 while (x) {
54 ++n;
55 x /= 10;
56 }
57 return n;
58}
59
60}
61
97void print_double(std::ostream& os, double v) {
98 char buffer[64];
99 constexpr std::uint64_t one = 1;
100 constexpr std::uint64_t p = 53;
101 std::uint64_t pun;
102 std::memcpy(&pun, &v, sizeof(v));
103 auto sign = pun >> 63;
104 pun -= sign << 63;
105 auto be = pun >> (p - 1);
106 pun -= be << (p - 1);
107 auto m = pun;
108 int bp = 0;
109 if (be == 2047) {
110 if (m == 0) {
111 if (sign) buffer[bp++] = '-';
112 buffer[bp++] = 'i';
113 buffer[bp++] = 'n';
114 buffer[bp++] = 'f';
115 } else {
116 buffer[bp++] = 'n';
117 buffer[bp++] = 'a';
118 buffer[bp++] = 'n';
119 }
120 } else {
121 if (sign) buffer[bp++] = '-';
122 std::uint64_t f;
123 if (be == 0) {
124 f = m;
125 } else {
126 f = m + (one << (p - 1));
127 }
128 auto e = std::int64_t(be) - 1075;
129 BigUInt<34> r, s, mp, mm;
130 if (e >= 0) {
131 if (f != (one << (p - 1))) {
132 r = BigUInt<34>(f);
133 r <<= (e + 1);
134 s = BigUInt<34>(2);
135 mp = BigUInt<34>(1);
136 mp <<= e;
137 mm = BigUInt<34>(1);
138 mm <<= e;
139 } else {
140 r = BigUInt<34>(f);
141 r <<= (e + 2);
142 s = BigUInt<34>(2 * 2);
143 mp = BigUInt<34>(1);
144 mp <<= (e + 1);
145 mm = BigUInt<34>(1);
146 mm <<= e;
147 }
148 } else {
149 if ((be == 0) || (f != (one << (p - 1)))) {
150 r = BigUInt<34>(f);
151 r <<= 1;
152 s = BigUInt<34>(1);
153 s <<= (1 - e);
154 mp = BigUInt<34>(1);
155 mm = BigUInt<34>(1);
156 } else {
157 r = BigUInt<34>(f);
158 r <<= 2;
159 s = BigUInt<34>(1);
160 s <<= (2 - e);
161 mp = BigUInt<34>(2);
162 mm = BigUInt<34>(1);
163 }
164 }
165 std::int32_t k = 0;
166 BigUInt<34> B_k{1};
167 auto r_p_mp = r + mp;
168 auto r_p_mp_comp = comp(r_p_mp, s);
169 if (r_p_mp_comp == 0) {
170 } else if (r_p_mp_comp == 1) {
171 while (r_p_mp > (s * B_k)) {
172 ++k, B_k *= 10;
173 }
174 } else {
175 while ((r_p_mp * B_k) < s) {
176 --k, B_k *= 10;
177 }
178 ++k;
179 B_k = B_k / 10;
180 }
181 if (k >= 0) {
182 s = s * B_k;
183 } else {
184 r = r * B_k;
185 mp = mp * B_k;
186 mm = mm * B_k;
187 }
188 char last_d = '0';
189 int n;
190 for (n = 0; true; ++n) {
191 auto r_x_10 = r;
192 r_x_10 *= 10;
193 auto d_np1 = r_x_10 / s;
194 auto cond1 = r < mm;
195 auto cond2 = (r + mp) > s;
196 if (cond1 && cond2) {
197 r <<= 1;
198 if (r < s) {
199 buffer[bp++] = last_d;
200 } else {
201 buffer[bp++] = last_d + 1;
202 }
203 break;
204 } else if (cond1) {
205 buffer[bp++] = last_d;
206 break;
207 } else if (cond2) {
208 buffer[bp++] = last_d + 1;
209 break;
210 } else {
211 if (n) buffer[bp++] = last_d;
212 r = r_x_10;
213 r -= (s * d_np1);
214 mp *= 10;
215 mm *= 10;
216 last_d = char(d_np1[0]) + '0';
217 }
218 }
219 if (v == 0.0) {
220 k = 1;
221 ++n;
222 }
223 int dot_pos = -1;
224 bool do_scientific = false;
225 if (0 <= k && k <= n) {
226 // dot is touching significant digits
227 dot_pos = k;
228 } else if (k < 0) {
229 auto nchars_sci = ndigits_for(-k + 1) + 2;
230 if (n > 1) nchars_sci += 1; // add a dot to scientific notation if more than one digit
231 if (nchars_sci < (-k + 1)) {
232 // scientific notation requires fewer chars than trailing zeros
233 if (n > 1) dot_pos = 1;
234 do_scientific = true;
235 } else {
236 // trailing zeros are no more chars than scientific
237 for (int i = 0; i < n; ++i) {
238 buffer[bp + (-k) - i - 1] = buffer[bp - i - 1];
239 }
240 for (int i = 0; i < -k; ++i) {
241 buffer[bp - n + i] = '0';
242 }
243 dot_pos = bp - n;
244 bp += -k;
245 n += -k;
246 }
247 } else if (k > n) {
248 auto nchars_sci = ndigits_for(k - 1) + 1;
249 if (n > 1) nchars_sci += 1; // add a dot to scientific notation if more than one digit
250 if (nchars_sci < ((k-n)+1)) {
251 // scientific notation requires fewer chars than trailing zeros
252 if (n > 1) dot_pos = 1;
253 do_scientific = true;
254 } else {
255 // trailing zeros are no more chars than scientific
256 for (; n < k; ++n) buffer[bp++] = '0';
257 dot_pos = n;
258 }
259 }
260 if (dot_pos != -1) {
261 for (int i = 0; i < (n - dot_pos) && i < bp; ++i) buffer[bp - i] = buffer[bp - i - 1];
262 buffer[bp - n + dot_pos] = '.';
263 ++bp;
264 }
265 if (do_scientific) {
266 buffer[bp++] = 'e';
267 auto decimal_exponent = (k - 1);
268 if (decimal_exponent < 0) {
269 buffer[bp++] = '-';
270 decimal_exponent = -decimal_exponent;
271 }
272 int ne;
273 for (ne = 0; decimal_exponent; ++ne) {
274 buffer[bp++] = char(decimal_exponent % 10) + '0';
275 decimal_exponent /= 10;
276 }
277 for (int i = 0; i < ne / 2; ++i) {
278 auto tmp = buffer[bp - ne + i];
279 buffer[bp - ne + i] = buffer[bp - i - 1];
280 buffer[bp - i - 1] = tmp;
281 }
282 }
283 }
284 buffer[bp] = '\0';
285 os << buffer;
286}
287
288}
Arbitrary-precision unsigned integer definition.
Declares Teuchos::print_double.
Arbitrary-precision unsigned integer class.
int comp(BigUInt< n > const &a, BigUInt< n > const &b)
void print_double(std::ostream &os, double v)
Prints a double-precision floating-point number exactly using minimal characters.