Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_PrintDouble.cpp
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 
42 #include "Teuchos_PrintDouble.hpp"
43 #include "Teuchos_BigUInt.hpp"
44 
45 #include <cstring>
46 
47 namespace Teuchos {
48 
49 namespace {
50 
51 int 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 
97 void 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 }
void print_double(std::ostream &os, double v)
Prints a double-precision floating-point number exactly using minimal characters. ...
Arbitrary-precision unsigned integer class.
Arbitrary-precision unsigned integer definition.
Declares Teuchos::print_double.