1 /* mpfr_get_decimal128 -- convert a multiple precision floating-point number
2 to an IEEE 754-2008 decimal128 float
3
4 See https://gcc.gnu.org/legacy-ml/gcc/2006-06/msg00691.html,
5 https://gcc.gnu.org/onlinedocs/gcc/Decimal-Float.html,
6 and TR 24732 <https://www.open-std.org/jtc1/sc22/wg14/www/projects#24732>.
7
8 Copyright 2006-2023 Free Software Foundation, Inc.
9 Contributed by the AriC and Caramba projects, INRIA.
10
11 This file is part of the GNU MPFR Library.
12
13 The GNU MPFR Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU Lesser General Public License as published by
15 the Free Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 The GNU MPFR Library is distributed in the hope that it will be useful, but
19 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
20 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21 License for more details.
22
23 You should have received a copy of the GNU Lesser General Public License
24 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
25 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
26 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
27
28 /* Warning! Do not use any conversion between binary and decimal types,
29 * otherwise GCC will generate from 2 to 3 MB of code (depending on the
30 * GCC version) in the MPFR shared library when the _Decimal128 format
31 * is BID (e.g. on x86).
32 * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96173
33 * https://gforge.inria.fr/tracker/index.php?func=detail&aid=21849&group_id=136&atid=619
34 *
35 * FIXME: Try to save even more space in the MPFR library by avoiding
36 * _Decimal128 operations entirely. These operations now appear only in
37 * string_to_Decimal128(). In the case where the _Decimal128 format is
38 * recognized as BID, this function should be reimplemented directly by
39 * using the specification of the encoding of this format, as already
40 * done for _Decimal64 (see string_to_Decimal64 in get_d64.c).
41 * Or use strtod128 when available, making sure that the string is
42 * locale-independent? (Should one optionally use libdfp for that?)
43 */
44
45 #include "mpfr-impl.h"
46 #include "ieee_floats.h"
47
48 #define ISDIGIT(c) ('0' <= c && c <= '9')
49
50 #ifdef MPFR_WANT_DECIMAL_FLOATS
51
52 #ifndef DEC128_MAX
53 # define DEC128_MAX 9.999999999999999999999999999999999E6144dl
54 #endif
55
56 /* construct a decimal128 NaN */
57 static _Decimal128
get_decimal128_nan(void)58 get_decimal128_nan (void)
59 {
60 return 0.dl / 0.dl;
61 }
62
63 /* construct the decimal128 Inf with given sign */
64 static _Decimal128
get_decimal128_inf(int negative)65 get_decimal128_inf (int negative)
66 {
67 return negative ? - 1.dl / 0.dl : 1.dl / 0.dl;
68 }
69
70 /* construct the decimal128 zero with given sign */
71 static _Decimal128
get_decimal128_zero(int negative)72 get_decimal128_zero (int negative)
73 {
74 return negative ? - 0.dl : 0.dl;
75 }
76
77 /* construct the decimal128 smallest non-zero with given sign:
78 it is 10^emin * 10^(1-p). Since emax = 6144, emin = 1-emax = -6143,
79 and p = 34, we get 10^(-6176) */
80 static _Decimal128
get_decimal128_min(int negative)81 get_decimal128_min (int negative)
82 {
83 return negative ? - 1E-6176dl : 1E-6176dl;
84 }
85
86 /* construct the decimal128 largest finite number with given sign */
87 static _Decimal128
get_decimal128_max(int negative)88 get_decimal128_max (int negative)
89 {
90 return negative ? - DEC128_MAX : DEC128_MAX;
91 }
92
93 /* one-to-one conversion:
94 s is a decimal string representing a number x = m * 10^e which must be
95 exactly representable in the decimal128 format, i.e.
96 (a) the mantissa m has at most 34 decimal digits
97 (b1) -6143 <= e <= 6144 with m integer multiple of 10^(-33), |m| < 10
98 (b2) or -6176 <= e <= 6111 with m integer, |m| < 10^34.
99 Assumes s is neither NaN nor +Inf nor -Inf.
100 s = [-][0-9]+E[-][0-9]+
101
102 The decimal128 format (cf table 3.6 of IEEE 754-2008) has the following
103 parameters:
104 * k = 128 (number of bits of storage)
105 * p = 34 (precision in digits)
106 * emax = 6144
107 * bias = E-q = 6176
108 * sign bit has 1 bit
109 * w+5 = 17 bits (combination field width)
110 * t = 110 bits (trailing significand width)
111 We have k = 1 + 5 + w + t = 128.
112 */
113 static _Decimal128
string_to_Decimal128(char * s)114 string_to_Decimal128 (char *s) /* portable version */
115 {
116 long int exp = 0;
117 char m[35];
118 long n = 0; /* mantissa length */
119 char *endptr[1];
120 _Decimal128 x = 0;
121 int sign = 0;
122
123 /* read sign */
124 if (*s == '-')
125 {
126 sign = 1;
127 s ++;
128 }
129 /* read mantissa */
130 while (ISDIGIT (*s))
131 m[n++] = *s++;
132
133 /* as constructed in mpfr_get_decimal128, s cannot have any '.' separator */
134
135 /* we will consider an integer mantissa m*10^exp */
136 MPFR_ASSERTN(n <= 34);
137 /* s always has an exponent separator 'E' */
138 MPFR_ASSERTN(*s == 'E');
139 exp = strtol (s + 1, endptr, 10);
140 MPFR_ASSERTN(**endptr == '\0');
141 MPFR_ASSERTN(-6176 <= exp && exp <= (long) (6145 - n));
142 while (n < 34)
143 {
144 m[n++] = '0';
145 exp --;
146 }
147 /* now n=34 and -6176 <= exp <= 6111, cf (b2) */
148 m[n] = '\0';
149
150 /* the number to convert is m[] * 10^exp where the mantissa is a 34-digit
151 integer */
152
153 /* compute biased exponent */
154 exp += 6176;
155
156 MPFR_ASSERTN(exp >= -33);
157 if (exp < 0)
158 {
159 int i;
160 n = -exp;
161 /* check the last n digits of the mantissa are zero */
162 for (i = 1; i <= n; i++)
163 MPFR_ASSERTN(m[34 - n] == '0');
164 /* shift the first (34-n) digits to the right */
165 for (i = 34 - n - 1; i >= 0; i--)
166 m[i + n] = m[i];
167 /* zero the first n digits */
168 for (i = 0; i < n; i ++)
169 m[i] = '0';
170 exp = 0;
171 }
172
173 /* the number to convert is m[] * 10^(exp-6176) */
174 exp -= 6176;
175
176 for (n = 0; n < 34; n++)
177 x = (_Decimal128) 10 * x + (_Decimal128) (m[n] - '0');
178
179 /* multiply by 10^exp */
180 if (exp > 0)
181 {
182 _Decimal128 ten = 10;
183 _Decimal128 ten2 = ten * ten;
184 _Decimal128 ten4 = ten2 * ten2;
185 _Decimal128 ten8 = ten4 * ten4;
186 _Decimal128 ten16 = ten8 * ten8;
187 _Decimal128 ten32 = ten16 * ten16;
188 _Decimal128 ten64 = ten32 * ten32;
189 _Decimal128 ten128 = ten64 * ten64;
190 _Decimal128 ten256 = ten128 * ten128;
191 _Decimal128 ten512 = ten256 * ten256;
192 _Decimal128 ten1024 = ten512 * ten512;
193 _Decimal128 ten2048 = ten1024 * ten1024;
194 _Decimal128 ten4096 = ten2048 * ten2048;
195
196 if (exp >= 4096)
197 {
198 x *= ten4096;
199 exp -= 4096;
200 }
201 if (exp >= 2048)
202 {
203 x *= ten2048;
204 exp -= 2048;
205 }
206 if (exp >= 1024)
207 {
208 x *= ten1024;
209 exp -= 1024;
210 }
211 if (exp >= 512)
212 {
213 x *= ten512;
214 exp -= 512;
215 }
216 if (exp >= 256)
217 {
218 x *= ten256;
219 exp -= 256;
220 }
221 if (exp >= 128)
222 {
223 x *= ten128;
224 exp -= 128;
225 }
226 if (exp >= 64)
227 {
228 x *= ten64;
229 exp -= 64;
230 }
231 if (exp >= 32)
232 {
233 x *= ten32;
234 exp -= 32;
235 }
236 if (exp >= 16)
237 {
238 x *= ten16;
239 exp -= 16;
240 }
241 if (exp >= 8)
242 {
243 x *= ten8;
244 exp -= 8;
245 }
246 if (exp >= 4)
247 {
248 x *= ten4;
249 exp -= 4;
250 }
251 if (exp >= 2)
252 {
253 x *= ten2;
254 exp -= 2;
255 }
256 if (exp >= 1)
257 {
258 x *= ten;
259 exp -= 1;
260 }
261 }
262 else if (exp < 0)
263 {
264 _Decimal128 ten = 10;
265 _Decimal128 ten2 = ten * ten;
266 _Decimal128 ten4 = ten2 * ten2;
267 _Decimal128 ten8 = ten4 * ten4;
268 _Decimal128 ten16 = ten8 * ten8;
269 _Decimal128 ten32 = ten16 * ten16;
270 _Decimal128 ten64 = ten32 * ten32;
271 _Decimal128 ten128 = ten64 * ten64;
272 _Decimal128 ten256 = ten128 * ten128;
273 _Decimal128 ten512 = ten256 * ten256;
274 _Decimal128 ten1024 = ten512 * ten512;
275 _Decimal128 ten2048 = ten1024 * ten1024;
276 _Decimal128 ten4096 = ten2048 * ten2048;
277
278 if (exp <= -4096)
279 {
280 x /= ten4096;
281 exp += 4096;
282 }
283 if (exp <= -2048)
284 {
285 x /= ten2048;
286 exp += 2048;
287 }
288 if (exp <= -1024)
289 {
290 x /= ten1024;
291 exp += 1024;
292 }
293 if (exp <= -512)
294 {
295 x /= ten512;
296 exp += 512;
297 }
298 if (exp <= -256)
299 {
300 x /= ten256;
301 exp += 256;
302 }
303 if (exp <= -128)
304 {
305 x /= ten128;
306 exp += 128;
307 }
308 if (exp <= -64)
309 {
310 x /= ten64;
311 exp += 64;
312 }
313 if (exp <= -32)
314 {
315 x /= ten32;
316 exp += 32;
317 }
318 if (exp <= -16)
319 {
320 x /= ten16;
321 exp += 16;
322 }
323 if (exp <= -8)
324 {
325 x /= ten8;
326 exp += 8;
327 }
328 if (exp <= -4)
329 {
330 x /= ten4;
331 exp += 4;
332 }
333 if (exp <= -2)
334 {
335 x /= ten2;
336 exp += 2;
337 }
338 if (exp <= -1)
339 {
340 x /= ten;
341 exp += 1;
342 }
343 }
344
345 if (sign)
346 x = -x;
347
348 return x;
349 }
350
351 _Decimal128
mpfr_get_decimal128(mpfr_srcptr src,mpfr_rnd_t rnd_mode)352 mpfr_get_decimal128 (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
353 {
354 int negative;
355 mpfr_exp_t e;
356
357 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
358 {
359 if (MPFR_IS_NAN (src))
360 {
361 /* we don't propagate the sign bit */
362 return get_decimal128_nan ();
363 }
364
365 negative = MPFR_IS_NEG (src);
366
367 if (MPFR_IS_INF (src))
368 return get_decimal128_inf (negative);
369
370 MPFR_ASSERTD (MPFR_IS_ZERO(src));
371 return get_decimal128_zero (negative);
372 }
373
374 e = MPFR_GET_EXP (src);
375 negative = MPFR_IS_NEG (src);
376
377 MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (src));
378
379 /* now rnd_mode is RNDN, RNDF, RNDA or RNDZ */
380
381 /* the smallest decimal128 number is 10^(-6176),
382 with 2^(-20517) < 10^(-6176) < 2^(-20516) */
383 if (MPFR_UNLIKELY (e < -20517)) /* src <= 2^(-20518) < 1/2*10^(-6176) */
384 {
385 if (rnd_mode != MPFR_RNDA)
386 return get_decimal128_zero (negative);
387 else /* RNDA: return the smallest non-zero number */
388 return get_decimal128_min (negative);
389 }
390 /* the largest decimal128 number is just below 10^6145 < 2^20414 */
391 else if (MPFR_UNLIKELY (e > 20414)) /* then src >= 2^20414 */
392 {
393 if (rnd_mode == MPFR_RNDZ)
394 return get_decimal128_max (negative);
395 else /* RNDN, RNDA, RNDF: round away */
396 return get_decimal128_inf (negative);
397 }
398 else
399 {
400 /* we need to store the sign (1 character), the significand (at most 34
401 characters), the exponent part (at most 6 characters for "E-6176"),
402 and the terminating null character, thus we need at least 42
403 characters in s. */
404 char s[42];
405 mpfr_get_str (s, &e, 10, 34, src, rnd_mode);
406 /* the smallest normal number is 1.000...000E-6143,
407 which corresponds to s=[0.]1000...000 and e=-6142 */
408 if (e < -6142)
409 {
410 /* the smallest subnormal number is 0.000...001E-6143 = 1E-6176,
411 which corresponds to s=[0.]1000...000 and e=-6175 */
412 if (e < -6175)
413 {
414 if (rnd_mode == MPFR_RNDN && e == -6176)
415 {
416 /* If 0.5E-6176 < |src| < 1E-6176 (smallest subnormal),
417 src should round to +/- 1E-6176 in MPFR_RNDN. */
418 mpfr_get_str (s, &e, 10, 1, src, MPFR_RNDA);
419 return e == -6176 && s[negative] <= '5' ?
420 get_decimal128_zero (negative) :
421 get_decimal128_min (negative);
422 }
423 if (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDN)
424 return get_decimal128_zero (negative);
425 else /* RNDA or RNDF: return the smallest non-zero number */
426 return get_decimal128_min (negative);
427 }
428 else
429 {
430 mpfr_exp_t e2;
431 long digits = 34 - (-6142 - e);
432 /* if e = -6175 then 34 - (-6142 - e) = 1 */
433 mpfr_get_str (s, &e2, 10, digits, src, rnd_mode);
434 /* Warning: we can have e2 = e + 1 here, when rounding to
435 nearest or away from zero. */
436 s[negative + digits] = 'E';
437 sprintf (s + negative + digits + 1, "%ld",
438 (long int)e2 - digits);
439 return string_to_Decimal128 (s);
440 }
441 }
442 /* the largest number is 9.999...999E+6144,
443 which corresponds to s=[0.]9999...999 and e=6145 */
444 else if (e > 6145)
445 {
446 if (rnd_mode == MPFR_RNDZ)
447 return get_decimal128_max (negative);
448 else /* RNDN, RNDA, RNDF: round away */
449 return get_decimal128_inf (negative);
450 }
451 else /* -6142 <= e <= 6145 */
452 {
453 s[34 + negative] = 'E';
454 sprintf (s + 35 + negative, "%ld", (long int) e - 34);
455 return string_to_Decimal128 (s);
456 }
457 }
458 }
459
460 #endif /* MPFR_WANT_DECIMAL_FLOATS */
461