1 /* Software floating-point emulation. Common operations.
2 Copyright (C) 1997-2014 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Richard Henderson (rth@cygnus.com),
5 Jakub Jelinek (jj@ultra.linux.cz),
6 David S. Miller (davem@redhat.com) and
7 Peter Maydell (pmaydell@chiark.greenend.org.uk).
9 The GNU C Library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
14 In addition to the permissions in the GNU Lesser General Public
15 License, the Free Software Foundation gives you unlimited
16 permission to link the compiled version of this file into
17 combinations with other programs, and to distribute those
18 combinations without any restriction coming from the use of this
19 file. (The Lesser General Public License restrictions do apply in
20 other respects; for example, they cover modification of the file,
21 and distribution when not linked into a combine executable.)
23 The GNU C Library is distributed in the hope that it will be useful,
24 but WITHOUT ANY WARRANTY; without even the implied warranty of
25 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
26 Lesser General Public License for more details.
28 You should have received a copy of the GNU Lesser General Public
29 License along with the GNU C Library; if not, see
30 <http://www.gnu.org/licenses/>. */
32 #define _FP_DECL(wc, X) \
33 _FP_I_TYPE X##_c __attribute__ ((unused)); \
34 _FP_I_TYPE X##_s __attribute__ ((unused)); \
35 _FP_I_TYPE X##_e __attribute__ ((unused)); \
36 _FP_FRAC_DECL_##wc (X)
38 /* Test whether the qNaN bit denotes a signaling NaN. */
39 #define _FP_FRAC_SNANP(fs, X) \
41 ? (_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs) \
42 : !(_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs))
43 #define _FP_FRAC_SNANP_SEMIRAW(fs, X) \
45 ? (_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs) \
46 : !(_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs))
48 /* Finish truly unpacking a native fp value by classifying the kind
49 of fp value and normalizing both the exponent and the fraction. */
51 #define _FP_UNPACK_CANONICAL(fs, wc, X) \
57 _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs; \
58 _FP_FRAC_SLL_##wc (X, _FP_WORKBITS); \
59 X##_e -= _FP_EXPBIAS_##fs; \
60 X##_c = FP_CLS_NORMAL; \
64 if (_FP_FRAC_ZEROP_##wc (X)) \
65 X##_c = FP_CLS_ZERO; \
68 /* A denormalized number. */ \
69 _FP_I_TYPE _FP_UNPACK_CANONICAL_shift; \
70 _FP_FRAC_CLZ_##wc (_FP_UNPACK_CANONICAL_shift, \
72 _FP_UNPACK_CANONICAL_shift -= _FP_FRACXBITS_##fs; \
73 _FP_FRAC_SLL_##wc (X, (_FP_UNPACK_CANONICAL_shift \
75 X##_e -= (_FP_EXPBIAS_##fs - 1 \
76 + _FP_UNPACK_CANONICAL_shift); \
77 X##_c = FP_CLS_NORMAL; \
78 FP_SET_EXCEPTION (FP_EX_DENORM); \
82 case _FP_EXPMAX_##fs: \
83 if (_FP_FRAC_ZEROP_##wc (X)) \
88 /* Check for signaling NaN. */ \
89 if (_FP_FRAC_SNANP (fs, X)) \
90 FP_SET_EXCEPTION (FP_EX_INVALID); \
97 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
98 shifted by _FP_WORKBITS but the implicit MSB is not inserted and
99 other classification is not done. */
100 #define _FP_UNPACK_SEMIRAW(fs, wc, X) _FP_FRAC_SLL_##wc (X, _FP_WORKBITS)
102 /* A semi-raw value has overflowed to infinity. Adjust the mantissa
103 and exponent appropriately. */
104 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X) \
107 if (FP_ROUNDMODE == FP_RND_NEAREST \
108 || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \
109 || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \
111 X##_e = _FP_EXPMAX_##fs; \
112 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
116 X##_e = _FP_EXPMAX_##fs - 1; \
117 _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc); \
119 FP_SET_EXCEPTION (FP_EX_INEXACT); \
120 FP_SET_EXCEPTION (FP_EX_OVERFLOW); \
124 /* Check for a semi-raw value being a signaling NaN and raise the
125 invalid exception if so. */
126 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X) \
129 if (X##_e == _FP_EXPMAX_##fs \
130 && !_FP_FRAC_ZEROP_##wc (X) \
131 && _FP_FRAC_SNANP_SEMIRAW (fs, X)) \
132 FP_SET_EXCEPTION (FP_EX_INVALID); \
136 /* Choose a NaN result from an operation on two semi-raw NaN
138 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP) \
141 /* _FP_CHOOSENAN expects raw values, so shift as required. */ \
142 _FP_FRAC_SRL_##wc (X, _FP_WORKBITS); \
143 _FP_FRAC_SRL_##wc (Y, _FP_WORKBITS); \
144 _FP_CHOOSENAN (fs, wc, R, X, Y, OP); \
145 _FP_FRAC_SLL_##wc (R, _FP_WORKBITS); \
149 /* Make the fractional part a quiet NaN, preserving the payload
150 if possible, otherwise make it the canonical quiet NaN and set
151 the sign bit accordingly. */
152 #define _FP_SETQNAN(fs, wc, X) \
155 if (_FP_QNANNEGATEDP) \
157 _FP_FRAC_HIGH_RAW_##fs (X) &= _FP_QNANBIT_##fs - 1; \
158 if (_FP_FRAC_ZEROP_##wc (X)) \
160 X##_s = _FP_NANSIGN_##fs; \
161 _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs); \
165 _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_QNANBIT_##fs; \
168 #define _FP_SETQNAN_SEMIRAW(fs, wc, X) \
171 if (_FP_QNANNEGATEDP) \
173 _FP_FRAC_HIGH_##fs (X) &= _FP_QNANBIT_SH_##fs - 1; \
174 if (_FP_FRAC_ZEROP_##wc (X)) \
176 X##_s = _FP_NANSIGN_##fs; \
177 _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs); \
178 _FP_FRAC_SLL_##wc (X, _FP_WORKBITS); \
182 _FP_FRAC_HIGH_##fs (X) |= _FP_QNANBIT_SH_##fs; \
186 /* Test whether a biased exponent is normal (not zero or maximum). */
187 #define _FP_EXP_NORMAL(fs, wc, X) (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
189 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
190 rounded and shifted right, with the rounding possibly increasing
191 the exponent (including changing a finite value to infinity). */
192 #define _FP_PACK_SEMIRAW(fs, wc, X) \
195 int _FP_PACK_SEMIRAW_is_tiny \
196 = X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X); \
197 if (_FP_TININESS_AFTER_ROUNDING \
198 && _FP_PACK_SEMIRAW_is_tiny) \
200 FP_DECL_##fs (_FP_PACK_SEMIRAW_T); \
201 _FP_FRAC_COPY_##wc (_FP_PACK_SEMIRAW_T, X); \
202 _FP_PACK_SEMIRAW_T##_s = X##_s; \
203 _FP_PACK_SEMIRAW_T##_e = X##_e; \
204 _FP_FRAC_SLL_##wc (_FP_PACK_SEMIRAW_T, 1); \
205 _FP_ROUND (wc, _FP_PACK_SEMIRAW_T); \
206 if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_SEMIRAW_T)) \
207 _FP_PACK_SEMIRAW_is_tiny = 0; \
210 if (_FP_PACK_SEMIRAW_is_tiny) \
212 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
213 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
214 FP_SET_EXCEPTION (FP_EX_UNDERFLOW); \
216 if (_FP_FRAC_HIGH_##fs (X) \
217 & (_FP_OVERFLOW_##fs >> 1)) \
219 _FP_FRAC_HIGH_##fs (X) &= ~(_FP_OVERFLOW_##fs >> 1); \
221 if (X##_e == _FP_EXPMAX_##fs) \
222 _FP_OVERFLOW_SEMIRAW (fs, wc, X); \
224 _FP_FRAC_SRL_##wc (X, _FP_WORKBITS); \
225 if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X)) \
227 if (!_FP_KEEPNANFRACP) \
229 _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs); \
230 X##_s = _FP_NANSIGN_##fs; \
233 _FP_SETQNAN (fs, wc, X); \
238 /* Before packing the bits back into the native fp result, take care
239 of such mundane things as rounding and overflow. Also, for some
240 kinds of fp values, the original parts may not have been fully
241 extracted -- but that is ok, we can regenerate them now. */
243 #define _FP_PACK_CANONICAL(fs, wc, X) \
248 case FP_CLS_NORMAL: \
249 X##_e += _FP_EXPBIAS_##fs; \
253 if (_FP_FRAC_OVERP_##wc (fs, X)) \
255 _FP_FRAC_CLEAR_OVERP_##wc (fs, X); \
258 _FP_FRAC_SRL_##wc (X, _FP_WORKBITS); \
259 if (X##_e >= _FP_EXPMAX_##fs) \
262 switch (FP_ROUNDMODE) \
264 case FP_RND_NEAREST: \
265 X##_c = FP_CLS_INF; \
269 X##_c = FP_CLS_INF; \
273 X##_c = FP_CLS_INF; \
276 if (X##_c == FP_CLS_INF) \
278 /* Overflow to infinity. */ \
279 X##_e = _FP_EXPMAX_##fs; \
280 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
284 /* Overflow to maximum normal. */ \
285 X##_e = _FP_EXPMAX_##fs - 1; \
286 _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc); \
288 FP_SET_EXCEPTION (FP_EX_OVERFLOW); \
289 FP_SET_EXCEPTION (FP_EX_INEXACT); \
294 /* We've got a denormalized number. */ \
295 int _FP_PACK_CANONICAL_is_tiny = 1; \
296 if (_FP_TININESS_AFTER_ROUNDING && X##_e == 0) \
298 FP_DECL_##fs (_FP_PACK_CANONICAL_T); \
299 _FP_FRAC_COPY_##wc (_FP_PACK_CANONICAL_T, X); \
300 _FP_PACK_CANONICAL_T##_s = X##_s; \
301 _FP_PACK_CANONICAL_T##_e = X##_e; \
302 _FP_ROUND (wc, _FP_PACK_CANONICAL_T); \
303 if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_CANONICAL_T)) \
304 _FP_PACK_CANONICAL_is_tiny = 0; \
306 X##_e = -X##_e + 1; \
307 if (X##_e <= _FP_WFRACBITS_##fs) \
309 _FP_FRAC_SRS_##wc (X, X##_e, _FP_WFRACBITS_##fs); \
311 if (_FP_FRAC_HIGH_##fs (X) \
312 & (_FP_OVERFLOW_##fs >> 1)) \
315 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
316 FP_SET_EXCEPTION (FP_EX_INEXACT); \
321 _FP_FRAC_SRL_##wc (X, _FP_WORKBITS); \
323 if (_FP_PACK_CANONICAL_is_tiny \
324 && ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
325 || (FP_TRAPPING_EXCEPTIONS \
326 & FP_EX_UNDERFLOW))) \
327 FP_SET_EXCEPTION (FP_EX_UNDERFLOW); \
331 /* Underflow to zero. */ \
333 if (!_FP_FRAC_ZEROP_##wc (X)) \
335 _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc); \
337 _FP_FRAC_LOW_##wc (X) >>= (_FP_WORKBITS); \
339 FP_SET_EXCEPTION (FP_EX_UNDERFLOW); \
346 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
350 X##_e = _FP_EXPMAX_##fs; \
351 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
355 X##_e = _FP_EXPMAX_##fs; \
356 if (!_FP_KEEPNANFRACP) \
358 _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs); \
359 X##_s = _FP_NANSIGN_##fs; \
362 _FP_SETQNAN (fs, wc, X); \
368 /* This one accepts raw argument and not cooked, returns
369 1 if X is a signaling NaN. */
370 #define _FP_ISSIGNAN(fs, wc, X) \
372 int _FP_ISSIGNAN_ret = 0; \
373 if (X##_e == _FP_EXPMAX_##fs) \
375 if (!_FP_FRAC_ZEROP_##wc (X) \
376 && _FP_FRAC_SNANP (fs, X)) \
377 _FP_ISSIGNAN_ret = 1; \
386 /* Addition on semi-raw values. */
387 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
390 if (X##_s == Y##_s) \
394 int _FP_ADD_INTERNAL_ediff = X##_e - Y##_e; \
395 if (_FP_ADD_INTERNAL_ediff > 0) \
400 /* Y is zero or denormalized. */ \
401 if (_FP_FRAC_ZEROP_##wc (Y)) \
403 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
404 _FP_FRAC_COPY_##wc (R, X); \
409 FP_SET_EXCEPTION (FP_EX_DENORM); \
410 _FP_ADD_INTERNAL_ediff--; \
411 if (_FP_ADD_INTERNAL_ediff == 0) \
413 _FP_FRAC_ADD_##wc (R, X, Y); \
416 if (X##_e == _FP_EXPMAX_##fs) \
418 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
419 _FP_FRAC_COPY_##wc (R, X); \
425 else if (X##_e == _FP_EXPMAX_##fs) \
427 /* X is NaN or Inf, Y is normal. */ \
428 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
429 _FP_FRAC_COPY_##wc (R, X); \
433 /* Insert implicit MSB of Y. */ \
434 _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs; \
437 /* Shift the mantissa of Y to the right \
438 _FP_ADD_INTERNAL_EDIFF steps; remember to account \
439 later for the implicit MSB of X. */ \
440 if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs) \
441 _FP_FRAC_SRS_##wc (Y, _FP_ADD_INTERNAL_ediff, \
442 _FP_WFRACBITS_##fs); \
443 else if (!_FP_FRAC_ZEROP_##wc (Y)) \
444 _FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc); \
445 _FP_FRAC_ADD_##wc (R, X, Y); \
447 else if (_FP_ADD_INTERNAL_ediff < 0) \
449 _FP_ADD_INTERNAL_ediff = -_FP_ADD_INTERNAL_ediff; \
453 /* X is zero or denormalized. */ \
454 if (_FP_FRAC_ZEROP_##wc (X)) \
456 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
457 _FP_FRAC_COPY_##wc (R, Y); \
462 FP_SET_EXCEPTION (FP_EX_DENORM); \
463 _FP_ADD_INTERNAL_ediff--; \
464 if (_FP_ADD_INTERNAL_ediff == 0) \
466 _FP_FRAC_ADD_##wc (R, Y, X); \
469 if (Y##_e == _FP_EXPMAX_##fs) \
471 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
472 _FP_FRAC_COPY_##wc (R, Y); \
478 else if (Y##_e == _FP_EXPMAX_##fs) \
480 /* Y is NaN or Inf, X is normal. */ \
481 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
482 _FP_FRAC_COPY_##wc (R, Y); \
486 /* Insert implicit MSB of X. */ \
487 _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs; \
490 /* Shift the mantissa of X to the right \
491 _FP_ADD_INTERNAL_EDIFF steps; remember to account \
492 later for the implicit MSB of Y. */ \
493 if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs) \
494 _FP_FRAC_SRS_##wc (X, _FP_ADD_INTERNAL_ediff, \
495 _FP_WFRACBITS_##fs); \
496 else if (!_FP_FRAC_ZEROP_##wc (X)) \
497 _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc); \
498 _FP_FRAC_ADD_##wc (R, Y, X); \
502 /* _FP_ADD_INTERNAL_ediff == 0. */ \
503 if (!_FP_EXP_NORMAL (fs, wc, X)) \
507 /* X and Y are zero or denormalized. */ \
509 if (_FP_FRAC_ZEROP_##wc (X)) \
511 if (!_FP_FRAC_ZEROP_##wc (Y)) \
512 FP_SET_EXCEPTION (FP_EX_DENORM); \
513 _FP_FRAC_COPY_##wc (R, Y); \
516 else if (_FP_FRAC_ZEROP_##wc (Y)) \
518 FP_SET_EXCEPTION (FP_EX_DENORM); \
519 _FP_FRAC_COPY_##wc (R, X); \
524 FP_SET_EXCEPTION (FP_EX_DENORM); \
525 _FP_FRAC_ADD_##wc (R, X, Y); \
526 if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
528 /* Normalized result. */ \
529 _FP_FRAC_HIGH_##fs (R) \
530 &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
538 /* X and Y are NaN or Inf. */ \
539 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
540 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
541 R##_e = _FP_EXPMAX_##fs; \
542 if (_FP_FRAC_ZEROP_##wc (X)) \
543 _FP_FRAC_COPY_##wc (R, Y); \
544 else if (_FP_FRAC_ZEROP_##wc (Y)) \
545 _FP_FRAC_COPY_##wc (R, X); \
547 _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP); \
551 /* The exponents of X and Y, both normal, are equal. The \
552 implicit MSBs will always add to increase the \
554 _FP_FRAC_ADD_##wc (R, X, Y); \
556 _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs); \
557 if (R##_e == _FP_EXPMAX_##fs) \
558 /* Overflow to infinity (depending on rounding mode). */ \
559 _FP_OVERFLOW_SEMIRAW (fs, wc, R); \
563 if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
566 _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
568 _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs); \
569 if (R##_e == _FP_EXPMAX_##fs) \
570 /* Overflow to infinity (depending on rounding mode). */ \
571 _FP_OVERFLOW_SEMIRAW (fs, wc, R); \
578 int _FP_ADD_INTERNAL_ediff = X##_e - Y##_e; \
579 if (_FP_ADD_INTERNAL_ediff > 0) \
585 /* Y is zero or denormalized. */ \
586 if (_FP_FRAC_ZEROP_##wc (Y)) \
588 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
589 _FP_FRAC_COPY_##wc (R, X); \
594 FP_SET_EXCEPTION (FP_EX_DENORM); \
595 _FP_ADD_INTERNAL_ediff--; \
596 if (_FP_ADD_INTERNAL_ediff == 0) \
598 _FP_FRAC_SUB_##wc (R, X, Y); \
601 if (X##_e == _FP_EXPMAX_##fs) \
603 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
604 _FP_FRAC_COPY_##wc (R, X); \
610 else if (X##_e == _FP_EXPMAX_##fs) \
612 /* X is NaN or Inf, Y is normal. */ \
613 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
614 _FP_FRAC_COPY_##wc (R, X); \
618 /* Insert implicit MSB of Y. */ \
619 _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs; \
622 /* Shift the mantissa of Y to the right \
623 _FP_ADD_INTERNAL_EDIFF steps; remember to account \
624 later for the implicit MSB of X. */ \
625 if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs) \
626 _FP_FRAC_SRS_##wc (Y, _FP_ADD_INTERNAL_ediff, \
627 _FP_WFRACBITS_##fs); \
628 else if (!_FP_FRAC_ZEROP_##wc (Y)) \
629 _FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc); \
630 _FP_FRAC_SUB_##wc (R, X, Y); \
632 else if (_FP_ADD_INTERNAL_ediff < 0) \
634 _FP_ADD_INTERNAL_ediff = -_FP_ADD_INTERNAL_ediff; \
639 /* X is zero or denormalized. */ \
640 if (_FP_FRAC_ZEROP_##wc (X)) \
642 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
643 _FP_FRAC_COPY_##wc (R, Y); \
648 FP_SET_EXCEPTION (FP_EX_DENORM); \
649 _FP_ADD_INTERNAL_ediff--; \
650 if (_FP_ADD_INTERNAL_ediff == 0) \
652 _FP_FRAC_SUB_##wc (R, Y, X); \
655 if (Y##_e == _FP_EXPMAX_##fs) \
657 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
658 _FP_FRAC_COPY_##wc (R, Y); \
664 else if (Y##_e == _FP_EXPMAX_##fs) \
666 /* Y is NaN or Inf, X is normal. */ \
667 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
668 _FP_FRAC_COPY_##wc (R, Y); \
672 /* Insert implicit MSB of X. */ \
673 _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs; \
676 /* Shift the mantissa of X to the right \
677 _FP_ADD_INTERNAL_EDIFF steps; remember to account \
678 later for the implicit MSB of Y. */ \
679 if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs) \
680 _FP_FRAC_SRS_##wc (X, _FP_ADD_INTERNAL_ediff, \
681 _FP_WFRACBITS_##fs); \
682 else if (!_FP_FRAC_ZEROP_##wc (X)) \
683 _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc); \
684 _FP_FRAC_SUB_##wc (R, Y, X); \
689 if (!_FP_EXP_NORMAL (fs, wc, X)) \
693 /* X and Y are zero or denormalized. */ \
695 if (_FP_FRAC_ZEROP_##wc (X)) \
697 _FP_FRAC_COPY_##wc (R, Y); \
698 if (_FP_FRAC_ZEROP_##wc (Y)) \
699 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
702 FP_SET_EXCEPTION (FP_EX_DENORM); \
707 else if (_FP_FRAC_ZEROP_##wc (Y)) \
709 FP_SET_EXCEPTION (FP_EX_DENORM); \
710 _FP_FRAC_COPY_##wc (R, X); \
716 FP_SET_EXCEPTION (FP_EX_DENORM); \
717 _FP_FRAC_SUB_##wc (R, X, Y); \
719 if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
721 /* |X| < |Y|, negate result. */ \
722 _FP_FRAC_SUB_##wc (R, Y, X); \
725 else if (_FP_FRAC_ZEROP_##wc (R)) \
726 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
732 /* X and Y are NaN or Inf, of opposite signs. */ \
733 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X); \
734 _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y); \
735 R##_e = _FP_EXPMAX_##fs; \
736 if (_FP_FRAC_ZEROP_##wc (X)) \
738 if (_FP_FRAC_ZEROP_##wc (Y)) \
741 R##_s = _FP_NANSIGN_##fs; \
742 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
743 _FP_FRAC_SLL_##wc (R, _FP_WORKBITS); \
744 FP_SET_EXCEPTION (FP_EX_INVALID); \
750 _FP_FRAC_COPY_##wc (R, Y); \
755 if (_FP_FRAC_ZEROP_##wc (Y)) \
759 _FP_FRAC_COPY_##wc (R, X); \
764 _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP); \
770 /* The exponents of X and Y, both normal, are equal. The \
771 implicit MSBs cancel. */ \
773 _FP_FRAC_SUB_##wc (R, X, Y); \
775 if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
777 /* |X| < |Y|, negate result. */ \
778 _FP_FRAC_SUB_##wc (R, Y, X); \
781 else if (_FP_FRAC_ZEROP_##wc (R)) \
784 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
790 if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
792 int _FP_ADD_INTERNAL_diff; \
793 /* Carry into most significant bit of larger one of X and Y, \
794 canceling it; renormalize. */ \
795 _FP_FRAC_HIGH_##fs (R) &= _FP_IMPLBIT_SH_##fs - 1; \
797 _FP_FRAC_CLZ_##wc (_FP_ADD_INTERNAL_diff, R); \
798 _FP_ADD_INTERNAL_diff -= _FP_WFRACXBITS_##fs; \
799 _FP_FRAC_SLL_##wc (R, _FP_ADD_INTERNAL_diff); \
800 if (R##_e <= _FP_ADD_INTERNAL_diff) \
802 /* R is denormalized. */ \
803 _FP_ADD_INTERNAL_diff \
804 = _FP_ADD_INTERNAL_diff - R##_e + 1; \
805 _FP_FRAC_SRS_##wc (R, _FP_ADD_INTERNAL_diff, \
806 _FP_WFRACBITS_##fs); \
811 R##_e -= _FP_ADD_INTERNAL_diff; \
812 _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
820 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL (fs, wc, R, X, Y, '+')
821 #define _FP_SUB(fs, wc, R, X, Y) \
824 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))) \
826 _FP_ADD_INTERNAL (fs, wc, R, X, Y, '-'); \
831 /* Main negation routine. The input value is raw. */
833 #define _FP_NEG(fs, wc, R, X) \
836 _FP_FRAC_COPY_##wc (R, X); \
843 /* Main multiplication routine. The input values should be cooked. */
845 #define _FP_MUL(fs, wc, R, X, Y) \
848 R##_s = X##_s ^ Y##_s; \
849 R##_e = X##_e + Y##_e + 1; \
850 switch (_FP_CLS_COMBINE (X##_c, Y##_c)) \
852 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL): \
853 R##_c = FP_CLS_NORMAL; \
855 _FP_MUL_MEAT_##fs (R, X, Y); \
857 if (_FP_FRAC_OVERP_##wc (fs, R)) \
858 _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs); \
863 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN): \
864 _FP_CHOOSENAN (fs, wc, R, X, Y, '*'); \
867 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL): \
868 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF): \
869 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO): \
872 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF): \
873 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL): \
874 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL): \
875 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO): \
876 _FP_FRAC_COPY_##wc (R, X); \
880 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN): \
881 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN): \
882 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN): \
885 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF): \
886 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO): \
887 _FP_FRAC_COPY_##wc (R, Y); \
891 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO): \
892 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF): \
893 R##_s = _FP_NANSIGN_##fs; \
894 R##_c = FP_CLS_NAN; \
895 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
896 FP_SET_EXCEPTION (FP_EX_INVALID); \
906 /* Fused multiply-add. The input values should be cooked. */
908 #define _FP_FMA(fs, wc, dwc, R, X, Y, Z) \
911 FP_DECL_##fs (_FP_FMA_T); \
912 _FP_FMA_T##_s = X##_s ^ Y##_s; \
913 _FP_FMA_T##_e = X##_e + Y##_e + 1; \
914 switch (_FP_CLS_COMBINE (X##_c, Y##_c)) \
916 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL): \
922 _FP_FRAC_COPY_##wc (R, Z); \
927 R##_c = FP_CLS_NORMAL; \
928 R##_s = _FP_FMA_T##_s; \
929 R##_e = _FP_FMA_T##_e; \
931 _FP_MUL_MEAT_##fs (R, X, Y); \
933 if (_FP_FRAC_OVERP_##wc (fs, R)) \
934 _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs); \
939 case FP_CLS_NORMAL:; \
940 _FP_FRAC_DECL_##dwc (_FP_FMA_TD); \
941 _FP_FRAC_DECL_##dwc (_FP_FMA_ZD); \
942 _FP_FRAC_DECL_##dwc (_FP_FMA_RD); \
943 _FP_MUL_MEAT_DW_##fs (_FP_FMA_TD, X, Y); \
944 R##_e = _FP_FMA_T##_e; \
946 = _FP_FRAC_HIGHBIT_DW_##dwc (fs, _FP_FMA_TD) == 0; \
947 _FP_FMA_T##_e -= _FP_FMA_tsh; \
948 int _FP_FMA_ediff = _FP_FMA_T##_e - Z##_e; \
949 if (_FP_FMA_ediff >= 0) \
952 = _FP_WFRACBITS_##fs - _FP_FMA_tsh - _FP_FMA_ediff; \
953 if (_FP_FMA_shift <= -_FP_WFRACBITS_##fs) \
954 _FP_FRAC_SET_##dwc (_FP_FMA_ZD, _FP_MINFRAC_##dwc); \
957 _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z); \
958 if (_FP_FMA_shift < 0) \
959 _FP_FRAC_SRS_##dwc (_FP_FMA_ZD, -_FP_FMA_shift, \
960 _FP_WFRACBITS_DW_##fs); \
961 else if (_FP_FMA_shift > 0) \
962 _FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_FMA_shift); \
964 R##_s = _FP_FMA_T##_s; \
965 if (_FP_FMA_T##_s == Z##_s) \
966 _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_TD, \
970 _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_TD, \
972 if (_FP_FRAC_NEGP_##dwc (_FP_FMA_RD)) \
975 _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD, \
984 _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z); \
985 _FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_WFRACBITS_##fs); \
986 int _FP_FMA_shift = -_FP_FMA_ediff - _FP_FMA_tsh; \
987 if (_FP_FMA_shift >= _FP_WFRACBITS_DW_##fs) \
988 _FP_FRAC_SET_##dwc (_FP_FMA_TD, _FP_MINFRAC_##dwc); \
989 else if (_FP_FMA_shift > 0) \
990 _FP_FRAC_SRS_##dwc (_FP_FMA_TD, _FP_FMA_shift, \
991 _FP_WFRACBITS_DW_##fs); \
992 if (Z##_s == _FP_FMA_T##_s) \
993 _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_ZD, \
996 _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD, \
999 if (_FP_FRAC_ZEROP_##dwc (_FP_FMA_RD)) \
1001 if (_FP_FMA_T##_s == Z##_s) \
1004 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
1005 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1006 R##_c = FP_CLS_ZERO; \
1011 _FP_FRAC_CLZ_##dwc (_FP_FMA_rlz, _FP_FMA_RD); \
1012 _FP_FMA_rlz -= _FP_WFRACXBITS_DW_##fs; \
1013 R##_e -= _FP_FMA_rlz; \
1014 int _FP_FMA_shift = _FP_WFRACBITS_##fs - _FP_FMA_rlz; \
1015 if (_FP_FMA_shift > 0) \
1016 _FP_FRAC_SRS_##dwc (_FP_FMA_RD, _FP_FMA_shift, \
1017 _FP_WFRACBITS_DW_##fs); \
1018 else if (_FP_FMA_shift < 0) \
1019 _FP_FRAC_SLL_##dwc (_FP_FMA_RD, -_FP_FMA_shift); \
1020 _FP_FRAC_COPY_##wc##_##dwc (R, _FP_FMA_RD); \
1021 R##_c = FP_CLS_NORMAL; \
1027 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN): \
1028 _FP_CHOOSENAN (fs, wc, _FP_FMA_T, X, Y, '*'); \
1031 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL): \
1032 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF): \
1033 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO): \
1034 _FP_FMA_T##_s = X##_s; \
1036 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF): \
1037 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL): \
1038 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL): \
1039 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO): \
1040 _FP_FRAC_COPY_##wc (_FP_FMA_T, X); \
1041 _FP_FMA_T##_c = X##_c; \
1044 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN): \
1045 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN): \
1046 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN): \
1047 _FP_FMA_T##_s = Y##_s; \
1049 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF): \
1050 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO): \
1051 _FP_FRAC_COPY_##wc (_FP_FMA_T, Y); \
1052 _FP_FMA_T##_c = Y##_c; \
1055 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO): \
1056 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF): \
1057 _FP_FMA_T##_s = _FP_NANSIGN_##fs; \
1058 _FP_FMA_T##_c = FP_CLS_NAN; \
1059 _FP_FRAC_SET_##wc (_FP_FMA_T, _FP_NANFRAC_##fs); \
1060 FP_SET_EXCEPTION (FP_EX_INVALID); \
1067 /* T = X * Y is zero, infinity or NaN. */ \
1068 switch (_FP_CLS_COMBINE (_FP_FMA_T##_c, Z##_c)) \
1070 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN): \
1071 _FP_CHOOSENAN (fs, wc, R, _FP_FMA_T, Z, '+'); \
1074 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL): \
1075 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF): \
1076 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO): \
1077 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL): \
1078 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO): \
1079 R##_s = _FP_FMA_T##_s; \
1080 _FP_FRAC_COPY_##wc (R, _FP_FMA_T); \
1081 R##_c = _FP_FMA_T##_c; \
1084 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN): \
1085 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN): \
1086 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL): \
1087 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF): \
1089 _FP_FRAC_COPY_##wc (R, Z); \
1093 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF): \
1094 if (_FP_FMA_T##_s == Z##_s) \
1097 _FP_FRAC_COPY_##wc (R, Z); \
1102 R##_s = _FP_NANSIGN_##fs; \
1103 R##_c = FP_CLS_NAN; \
1104 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
1105 FP_SET_EXCEPTION (FP_EX_INVALID); \
1109 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO): \
1110 if (_FP_FMA_T##_s == Z##_s) \
1113 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
1114 _FP_FRAC_COPY_##wc (R, Z); \
1126 /* Main division routine. The input values should be cooked. */
1128 #define _FP_DIV(fs, wc, R, X, Y) \
1131 R##_s = X##_s ^ Y##_s; \
1132 R##_e = X##_e - Y##_e; \
1133 switch (_FP_CLS_COMBINE (X##_c, Y##_c)) \
1135 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL): \
1136 R##_c = FP_CLS_NORMAL; \
1138 _FP_DIV_MEAT_##fs (R, X, Y); \
1141 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN): \
1142 _FP_CHOOSENAN (fs, wc, R, X, Y, '/'); \
1145 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL): \
1146 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF): \
1147 case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO): \
1149 _FP_FRAC_COPY_##wc (R, X); \
1153 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN): \
1154 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN): \
1155 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN): \
1157 _FP_FRAC_COPY_##wc (R, Y); \
1161 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF): \
1162 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF): \
1163 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL): \
1164 R##_c = FP_CLS_ZERO; \
1167 case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO): \
1168 FP_SET_EXCEPTION (FP_EX_DIVZERO); \
1169 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO): \
1170 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL): \
1171 R##_c = FP_CLS_INF; \
1174 case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF): \
1175 case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO): \
1176 R##_s = _FP_NANSIGN_##fs; \
1177 R##_c = FP_CLS_NAN; \
1178 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
1179 FP_SET_EXCEPTION (FP_EX_INVALID); \
1189 /* Main differential comparison routine. The inputs should be raw not
1190 cooked. The return is -1,0,1 for normal values, 2 otherwise. */
1192 #define _FP_CMP(fs, wc, ret, X, Y, un) \
1195 /* NANs are unordered. */ \
1196 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X)) \
1197 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))) \
1203 int _FP_CMP_is_zero_x; \
1204 int _FP_CMP_is_zero_y; \
1207 = (!X##_e && _FP_FRAC_ZEROP_##wc (X)) ? 1 : 0; \
1209 = (!Y##_e && _FP_FRAC_ZEROP_##wc (Y)) ? 1 : 0; \
1211 if (_FP_CMP_is_zero_x && _FP_CMP_is_zero_y) \
1213 else if (_FP_CMP_is_zero_x) \
1214 ret = Y##_s ? 1 : -1; \
1215 else if (_FP_CMP_is_zero_y) \
1216 ret = X##_s ? -1 : 1; \
1217 else if (X##_s != Y##_s) \
1218 ret = X##_s ? -1 : 1; \
1219 else if (X##_e > Y##_e) \
1220 ret = X##_s ? -1 : 1; \
1221 else if (X##_e < Y##_e) \
1222 ret = X##_s ? 1 : -1; \
1223 else if (_FP_FRAC_GT_##wc (X, Y)) \
1224 ret = X##_s ? -1 : 1; \
1225 else if (_FP_FRAC_GT_##wc (Y, X)) \
1226 ret = X##_s ? 1 : -1; \
1234 /* Simplification for strict equality. */
1236 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
1239 /* NANs are unordered. */ \
1240 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X)) \
1241 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))) \
1247 ret = !(X##_e == Y##_e \
1248 && _FP_FRAC_EQ_##wc (X, Y) \
1249 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc (X)))); \
1254 /* Version to test unordered. */
1256 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
1259 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X)) \
1260 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))); \
1264 /* Main square root routine. The input value should be cooked. */
1266 #define _FP_SQRT(fs, wc, R, X) \
1269 _FP_FRAC_DECL_##wc (_FP_SQRT_T); \
1270 _FP_FRAC_DECL_##wc (_FP_SQRT_S); \
1271 _FP_W_TYPE _FP_SQRT_q; \
1275 _FP_FRAC_COPY_##wc (R, X); \
1277 R##_c = FP_CLS_NAN; \
1282 R##_s = _FP_NANSIGN_##fs; \
1283 R##_c = FP_CLS_NAN; /* NAN */ \
1284 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
1285 FP_SET_EXCEPTION (FP_EX_INVALID); \
1290 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
1295 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
1297 case FP_CLS_NORMAL: \
1301 R##_c = FP_CLS_NAN; /* NAN */ \
1302 R##_s = _FP_NANSIGN_##fs; \
1303 _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs); \
1304 FP_SET_EXCEPTION (FP_EX_INVALID); \
1307 R##_c = FP_CLS_NORMAL; \
1309 _FP_FRAC_SLL_##wc (X, 1); \
1310 R##_e = X##_e >> 1; \
1311 _FP_FRAC_SET_##wc (_FP_SQRT_S, _FP_ZEROFRAC_##wc); \
1312 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1313 _FP_SQRT_q = _FP_OVERFLOW_##fs >> 1; \
1314 _FP_SQRT_MEAT_##wc (R, _FP_SQRT_S, _FP_SQRT_T, X, \
1320 /* Convert from FP to integer. Input is raw. */
1322 /* RSIGNED can have following values:
1323 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1324 the result is either 0 or (2^rsize)-1 depending on the sign in such
1326 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1327 NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1328 depending on the sign in such case.
1329 -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1330 set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1331 depending on the sign in such case. */
1332 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1335 if (X##_e < _FP_EXPBIAS_##fs) \
1340 if (!_FP_FRAC_ZEROP_##wc (X)) \
1342 FP_SET_EXCEPTION (FP_EX_INEXACT); \
1343 FP_SET_EXCEPTION (FP_EX_DENORM); \
1347 FP_SET_EXCEPTION (FP_EX_INEXACT); \
1349 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1350 || (!rsigned && X##_s)) \
1352 /* Overflow or converting to the most negative integer. */ \
1366 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1368 /* Possibly converting to most negative integer; check the \
1370 int _FP_TO_INT_inexact = 0; \
1371 (void) ((_FP_FRACBITS_##fs > rsize) \
1373 _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact, \
1374 _FP_FRACBITS_##fs - rsize, \
1375 _FP_FRACBITS_##fs); \
1379 if (!_FP_FRAC_ZEROP_##wc (X)) \
1380 FP_SET_EXCEPTION (FP_EX_INVALID); \
1381 else if (_FP_TO_INT_inexact) \
1382 FP_SET_EXCEPTION (FP_EX_INEXACT); \
1385 FP_SET_EXCEPTION (FP_EX_INVALID); \
1389 _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs; \
1390 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1392 _FP_FRAC_ASSEMBLE_##wc (r, X, rsize); \
1393 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1397 int _FP_TO_INT_inexact; \
1398 _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact, \
1399 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1401 _FP_FRACBITS_##fs); \
1402 if (_FP_TO_INT_inexact) \
1403 FP_SET_EXCEPTION (FP_EX_INEXACT); \
1404 _FP_FRAC_ASSEMBLE_##wc (r, X, rsize); \
1406 if (rsigned && X##_s) \
1412 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1414 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1419 rtype _FP_FROM_INT_ur; \
1421 if ((X##_s = (r < 0))) \
1424 _FP_FROM_INT_ur = (rtype) r; \
1425 (void) ((rsize <= _FP_W_TYPE_SIZE) \
1427 int _FP_FROM_INT_lz; \
1428 __FP_CLZ (_FP_FROM_INT_lz, \
1429 (_FP_W_TYPE) _FP_FROM_INT_ur); \
1430 X##_e = (_FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 \
1431 - _FP_FROM_INT_lz); \
1433 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1435 int _FP_FROM_INT_lz; \
1436 __FP_CLZ_2 (_FP_FROM_INT_lz, \
1437 (_FP_W_TYPE) (_FP_FROM_INT_ur \
1438 >> _FP_W_TYPE_SIZE), \
1439 (_FP_W_TYPE) _FP_FROM_INT_ur); \
1440 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1441 - _FP_FROM_INT_lz); \
1443 : (abort (), 0))); \
1445 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1446 && X##_e >= _FP_EXPMAX_##fs) \
1448 /* Exponent too big; overflow to infinity. (May also \
1449 happen after rounding below.) */ \
1450 _FP_OVERFLOW_SEMIRAW (fs, wc, X); \
1451 goto pack_semiraw; \
1454 if (rsize <= _FP_FRACBITS_##fs \
1455 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1457 /* Exactly representable; shift left. */ \
1458 _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, rsize); \
1459 if (_FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1 - X##_e > 0) \
1460 _FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs \
1461 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1465 /* More bits in integer than in floating type; need to \
1467 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1469 = ((_FP_FROM_INT_ur >> (X##_e - _FP_EXPBIAS_##fs \
1470 - _FP_WFRACBITS_##fs + 1)) \
1471 | ((_FP_FROM_INT_ur \
1472 << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1473 - _FP_WFRACBITS_##fs + 1))) \
1475 _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, rsize); \
1476 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1477 _FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs \
1478 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1479 _FP_FRAC_HIGH_##fs (X) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
1481 _FP_PACK_SEMIRAW (fs, wc, X); \
1488 _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc); \
1494 /* Extend from a narrower floating-point format to a wider one. Input
1495 and output are raw. */
1496 #define FP_EXTEND(dfs, sfs, dwc, swc, D, S) \
1499 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1500 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1501 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1502 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1503 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
1506 _FP_FRAC_COPY_##dwc##_##swc (D, S); \
1507 if (_FP_EXP_NORMAL (sfs, swc, S)) \
1509 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1510 _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1516 if (_FP_FRAC_ZEROP_##swc (S)) \
1518 else if (_FP_EXPBIAS_##dfs \
1519 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1521 FP_SET_EXCEPTION (FP_EX_DENORM); \
1522 _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs \
1523 - _FP_FRACBITS_##sfs)); \
1525 if (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW) \
1526 FP_SET_EXCEPTION (FP_EX_UNDERFLOW); \
1531 FP_SET_EXCEPTION (FP_EX_DENORM); \
1532 _FP_FRAC_CLZ_##swc (FP_EXTEND_lz, S); \
1533 _FP_FRAC_SLL_##dwc (D, \
1534 FP_EXTEND_lz + _FP_FRACBITS_##dfs \
1535 - _FP_FRACTBITS_##sfs); \
1536 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1537 + _FP_FRACXBITS_##sfs - FP_EXTEND_lz); \
1542 D##_e = _FP_EXPMAX_##dfs; \
1543 if (!_FP_FRAC_ZEROP_##swc (S)) \
1545 if (_FP_FRAC_SNANP (sfs, S)) \
1546 FP_SET_EXCEPTION (FP_EX_INVALID); \
1547 _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs \
1548 - _FP_FRACBITS_##sfs)); \
1549 _FP_SETQNAN (dfs, dwc, D); \
1556 /* Truncate from a wider floating-point format to a narrower one.
1557 Input and output are semi-raw. */
1558 #define FP_TRUNC(dfs, sfs, dwc, swc, D, S) \
1561 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1562 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1563 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1566 if (_FP_EXP_NORMAL (sfs, swc, S)) \
1568 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1569 if (D##_e >= _FP_EXPMAX_##dfs) \
1570 _FP_OVERFLOW_SEMIRAW (dfs, dwc, D); \
1575 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1577 _FP_FRAC_SET_##swc (S, _FP_ZEROFRAC_##swc); \
1578 _FP_FRAC_LOW_##swc (S) |= 1; \
1582 _FP_FRAC_HIGH_##sfs (S) |= _FP_IMPLBIT_SH_##sfs; \
1583 _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs \
1584 - _FP_WFRACBITS_##dfs \
1586 _FP_WFRACBITS_##sfs); \
1591 _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs \
1592 - _FP_WFRACBITS_##dfs), \
1593 _FP_WFRACBITS_##sfs); \
1594 _FP_FRAC_COPY_##dwc##_##swc (D, S); \
1602 if (_FP_FRAC_ZEROP_##swc (S)) \
1603 _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc); \
1606 FP_SET_EXCEPTION (FP_EX_DENORM); \
1607 if (_FP_EXPBIAS_##sfs \
1608 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1610 _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs \
1611 - _FP_WFRACBITS_##dfs), \
1612 _FP_WFRACBITS_##sfs); \
1613 _FP_FRAC_COPY_##dwc##_##swc (D, S); \
1617 _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc); \
1618 _FP_FRAC_LOW_##dwc (D) |= 1; \
1624 D##_e = _FP_EXPMAX_##dfs; \
1625 if (_FP_FRAC_ZEROP_##swc (S)) \
1626 _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc); \
1629 _FP_CHECK_SIGNAN_SEMIRAW (sfs, swc, S); \
1630 _FP_FRAC_SRL_##swc (S, (_FP_WFRACBITS_##sfs \
1631 - _FP_WFRACBITS_##dfs)); \
1632 _FP_FRAC_COPY_##dwc##_##swc (D, S); \
1633 /* Semi-raw NaN must have all workbits cleared. */ \
1634 _FP_FRAC_LOW_##dwc (D) \
1635 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1636 _FP_SETQNAN_SEMIRAW (dfs, dwc, D); \
1643 /* Helper primitives. */
1645 /* Count leading zeros in a word. */
1648 /* GCC 3.4 and later provide the builtins for us. */
1649 # define __FP_CLZ(r, x) \
1652 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1653 r = __builtin_clz (x); \
1654 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1655 r = __builtin_clzl (x); \
1656 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1657 r = __builtin_clzll (x); \
1662 #endif /* ndef __FP_CLZ */
1664 #define _FP_DIV_HELP_imm(q, r, n, d) \
1667 q = n / d, r = n % d; \
1672 /* A restoring bit-by-bit division primitive. */
1674 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1677 int _FP_DIV_MEAT_N_loop_count = _FP_WFRACBITS_##fs; \
1678 _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_u); \
1679 _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_v); \
1680 _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_u, X); \
1681 _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_v, Y); \
1682 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1683 /* Normalize _FP_DIV_MEAT_N_LOOP_U and _FP_DIV_MEAT_N_LOOP_V. */ \
1684 _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, _FP_WFRACXBITS_##fs); \
1685 _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_v, _FP_WFRACXBITS_##fs); \
1686 /* First round. Since the operands are normalized, either the \
1687 first or second bit will be set in the fraction. Produce a \
1688 normalized result by checking which and adjusting the loop \
1689 count and exponent accordingly. */ \
1690 if (_FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u, _FP_DIV_MEAT_N_loop_v)) \
1692 _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u, \
1693 _FP_DIV_MEAT_N_loop_u, \
1694 _FP_DIV_MEAT_N_loop_v); \
1695 _FP_FRAC_LOW_##wc (R) |= 1; \
1696 _FP_DIV_MEAT_N_loop_count--; \
1700 /* Subsequent rounds. */ \
1703 int _FP_DIV_MEAT_N_loop_msb \
1704 = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (_FP_DIV_MEAT_N_loop_u) < 0; \
1705 _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, 1); \
1706 _FP_FRAC_SLL_##wc (R, 1); \
1707 if (_FP_DIV_MEAT_N_loop_msb \
1708 || _FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u, \
1709 _FP_DIV_MEAT_N_loop_v)) \
1711 _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u, \
1712 _FP_DIV_MEAT_N_loop_u, \
1713 _FP_DIV_MEAT_N_loop_v); \
1714 _FP_FRAC_LOW_##wc (R) |= 1; \
1717 while (--_FP_DIV_MEAT_N_loop_count > 0); \
1718 /* If there's anything left in _FP_DIV_MEAT_N_LOOP_U, the result \
1720 _FP_FRAC_LOW_##wc (R) \
1721 |= !_FP_FRAC_ZEROP_##wc (_FP_DIV_MEAT_N_loop_u); \
1725 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1726 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1727 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)