1 /* Software floating-point emulation. Common operations.
2 Copyright (C) 1997-2013 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)); \
39 * Finish truely unpacking a native fp value by classifying the kind
40 * of fp value and normalizing both the exponent and the fraction.
43 #define _FP_UNPACK_CANONICAL(fs, wc, X) \
48 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
49 _FP_FRAC_SLL_##wc(X, _FP_WORKBITS); \
50 X##_e -= _FP_EXPBIAS_##fs; \
51 X##_c = FP_CLS_NORMAL; \
55 if (_FP_FRAC_ZEROP_##wc(X)) \
56 X##_c = FP_CLS_ZERO; \
59 /* a denormalized number */ \
61 _FP_FRAC_CLZ_##wc(_shift, X); \
62 _shift -= _FP_FRACXBITS_##fs; \
63 _FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS)); \
64 X##_e -= _FP_EXPBIAS_##fs - 1 + _shift; \
65 X##_c = FP_CLS_NORMAL; \
66 FP_SET_EXCEPTION(FP_EX_DENORM); \
70 case _FP_EXPMAX_##fs: \
71 if (_FP_FRAC_ZEROP_##wc(X)) \
76 /* Check for signaling NaN */ \
77 if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
78 FP_SET_EXCEPTION(FP_EX_INVALID); \
84 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
85 shifted by _FP_WORKBITS but the implicit MSB is not inserted and
86 other classification is not done. */
87 #define _FP_UNPACK_SEMIRAW(fs, wc, X) _FP_FRAC_SLL_##wc(X, _FP_WORKBITS)
89 /* A semi-raw value has overflowed to infinity. Adjust the mantissa
90 and exponent appropriately. */
91 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X) \
93 if (FP_ROUNDMODE == FP_RND_NEAREST \
94 || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \
95 || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \
97 X##_e = _FP_EXPMAX_##fs; \
98 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
102 X##_e = _FP_EXPMAX_##fs - 1; \
103 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
105 FP_SET_EXCEPTION(FP_EX_INEXACT); \
106 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
109 /* Check for a semi-raw value being a signaling NaN and raise the
110 invalid exception if so. */
111 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X) \
113 if (X##_e == _FP_EXPMAX_##fs \
114 && !_FP_FRAC_ZEROP_##wc(X) \
115 && !(_FP_FRAC_HIGH_##fs(X) & _FP_QNANBIT_SH_##fs)) \
116 FP_SET_EXCEPTION(FP_EX_INVALID); \
119 /* Choose a NaN result from an operation on two semi-raw NaN
121 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP) \
123 /* _FP_CHOOSENAN expects raw values, so shift as required. */ \
124 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
125 _FP_FRAC_SRL_##wc(Y, _FP_WORKBITS); \
126 _FP_CHOOSENAN(fs, wc, R, X, Y, OP); \
127 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
130 /* Test whether a biased exponent is normal (not zero or maximum). */
131 #define _FP_EXP_NORMAL(fs, wc, X) (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
133 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
134 rounded and shifted right, with the rounding possibly increasing
135 the exponent (including changing a finite value to infinity). */
136 #define _FP_PACK_SEMIRAW(fs, wc, X) \
139 if (X##_e == 0 && !_FP_FRAC_ZEROP_##wc(X)) \
141 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
142 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
143 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
145 if (_FP_FRAC_HIGH_##fs(X) \
146 & (_FP_OVERFLOW_##fs >> 1)) \
148 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1); \
150 if (X##_e == _FP_EXPMAX_##fs) \
151 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
153 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
154 if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
156 if (!_FP_KEEPNANFRACP) \
158 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
159 X##_s = _FP_NANSIGN_##fs; \
162 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
167 * Before packing the bits back into the native fp result, take care
168 * of such mundane things as rounding and overflow. Also, for some
169 * kinds of fp values, the original parts may not have been fully
170 * extracted -- but that is ok, we can regenerate them now.
173 #define _FP_PACK_CANONICAL(fs, wc, X) \
177 case FP_CLS_NORMAL: \
178 X##_e += _FP_EXPBIAS_##fs; \
182 if (_FP_FRAC_OVERP_##wc(fs, X)) \
184 _FP_FRAC_CLEAR_OVERP_##wc(fs, X); \
187 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
188 if (X##_e >= _FP_EXPMAX_##fs) \
191 switch (FP_ROUNDMODE) \
193 case FP_RND_NEAREST: \
194 X##_c = FP_CLS_INF; \
197 if (!X##_s) X##_c = FP_CLS_INF; \
200 if (X##_s) X##_c = FP_CLS_INF; \
203 if (X##_c == FP_CLS_INF) \
205 /* Overflow to infinity */ \
206 X##_e = _FP_EXPMAX_##fs; \
207 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
211 /* Overflow to maximum normal */ \
212 X##_e = _FP_EXPMAX_##fs - 1; \
213 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
215 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
216 FP_SET_EXCEPTION(FP_EX_INEXACT); \
221 /* we've got a denormalized number */ \
222 X##_e = -X##_e + 1; \
223 if (X##_e <= _FP_WFRACBITS_##fs) \
225 _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs); \
227 if (_FP_FRAC_HIGH_##fs(X) \
228 & (_FP_OVERFLOW_##fs >> 1)) \
231 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
232 FP_SET_EXCEPTION(FP_EX_INEXACT); \
237 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
239 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
240 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
241 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
245 /* underflow to zero */ \
247 if (!_FP_FRAC_ZEROP_##wc(X)) \
249 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
251 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
253 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
260 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
264 X##_e = _FP_EXPMAX_##fs; \
265 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
269 X##_e = _FP_EXPMAX_##fs; \
270 if (!_FP_KEEPNANFRACP) \
272 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
273 X##_s = _FP_NANSIGN_##fs; \
276 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
281 /* This one accepts raw argument and not cooked, returns
282 * 1 if X is a signaling NaN.
284 #define _FP_ISSIGNAN(fs, wc, X) \
287 if (X##_e == _FP_EXPMAX_##fs) \
289 if (!_FP_FRAC_ZEROP_##wc(X) \
290 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
300 /* Addition on semi-raw values. */
301 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
303 if (X##_s == Y##_s) \
307 int ediff = X##_e - Y##_e; \
313 /* Y is zero or denormalized. */ \
314 if (_FP_FRAC_ZEROP_##wc(Y)) \
316 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
317 _FP_FRAC_COPY_##wc(R, X); \
322 FP_SET_EXCEPTION(FP_EX_DENORM); \
326 _FP_FRAC_ADD_##wc(R, X, Y); \
329 if (X##_e == _FP_EXPMAX_##fs) \
331 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
332 _FP_FRAC_COPY_##wc(R, X); \
338 else if (X##_e == _FP_EXPMAX_##fs) \
340 /* X is NaN or Inf, Y is normal. */ \
341 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
342 _FP_FRAC_COPY_##wc(R, X); \
346 /* Insert implicit MSB of Y. */ \
347 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
350 /* Shift the mantissa of Y to the right EDIFF steps; \
351 remember to account later for the implicit MSB of X. */ \
352 if (ediff <= _FP_WFRACBITS_##fs) \
353 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
354 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
355 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
356 _FP_FRAC_ADD_##wc(R, X, Y); \
358 else if (ediff < 0) \
364 /* X is zero or denormalized. */ \
365 if (_FP_FRAC_ZEROP_##wc(X)) \
367 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
368 _FP_FRAC_COPY_##wc(R, Y); \
373 FP_SET_EXCEPTION(FP_EX_DENORM); \
377 _FP_FRAC_ADD_##wc(R, Y, X); \
380 if (Y##_e == _FP_EXPMAX_##fs) \
382 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
383 _FP_FRAC_COPY_##wc(R, Y); \
389 else if (Y##_e == _FP_EXPMAX_##fs) \
391 /* Y is NaN or Inf, X is normal. */ \
392 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
393 _FP_FRAC_COPY_##wc(R, Y); \
397 /* Insert implicit MSB of X. */ \
398 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
401 /* Shift the mantissa of X to the right EDIFF steps; \
402 remember to account later for the implicit MSB of Y. */ \
403 if (ediff <= _FP_WFRACBITS_##fs) \
404 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
405 else if (!_FP_FRAC_ZEROP_##wc(X)) \
406 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
407 _FP_FRAC_ADD_##wc(R, Y, X); \
412 if (!_FP_EXP_NORMAL(fs, wc, X)) \
416 /* X and Y are zero or denormalized. */ \
418 if (_FP_FRAC_ZEROP_##wc(X)) \
420 if (!_FP_FRAC_ZEROP_##wc(Y)) \
421 FP_SET_EXCEPTION(FP_EX_DENORM); \
422 _FP_FRAC_COPY_##wc(R, Y); \
425 else if (_FP_FRAC_ZEROP_##wc(Y)) \
427 FP_SET_EXCEPTION(FP_EX_DENORM); \
428 _FP_FRAC_COPY_##wc(R, X); \
433 FP_SET_EXCEPTION(FP_EX_DENORM); \
434 _FP_FRAC_ADD_##wc(R, X, Y); \
435 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
437 /* Normalized result. */ \
438 _FP_FRAC_HIGH_##fs(R) \
439 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
447 /* X and Y are NaN or Inf. */ \
448 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
449 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
450 R##_e = _FP_EXPMAX_##fs; \
451 if (_FP_FRAC_ZEROP_##wc(X)) \
452 _FP_FRAC_COPY_##wc(R, Y); \
453 else if (_FP_FRAC_ZEROP_##wc(Y)) \
454 _FP_FRAC_COPY_##wc(R, X); \
456 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
460 /* The exponents of X and Y, both normal, are equal. The \
461 implicit MSBs will always add to increase the \
463 _FP_FRAC_ADD_##wc(R, X, Y); \
465 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
466 if (R##_e == _FP_EXPMAX_##fs) \
467 /* Overflow to infinity (depending on rounding mode). */ \
468 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
472 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
475 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
477 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
478 if (R##_e == _FP_EXPMAX_##fs) \
479 /* Overflow to infinity (depending on rounding mode). */ \
480 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
487 int ediff = X##_e - Y##_e; \
494 /* Y is zero or denormalized. */ \
495 if (_FP_FRAC_ZEROP_##wc(Y)) \
497 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
498 _FP_FRAC_COPY_##wc(R, X); \
503 FP_SET_EXCEPTION(FP_EX_DENORM); \
507 _FP_FRAC_SUB_##wc(R, X, Y); \
510 if (X##_e == _FP_EXPMAX_##fs) \
512 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
513 _FP_FRAC_COPY_##wc(R, X); \
519 else if (X##_e == _FP_EXPMAX_##fs) \
521 /* X is NaN or Inf, Y is normal. */ \
522 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
523 _FP_FRAC_COPY_##wc(R, X); \
527 /* Insert implicit MSB of Y. */ \
528 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
531 /* Shift the mantissa of Y to the right EDIFF steps; \
532 remember to account later for the implicit MSB of X. */ \
533 if (ediff <= _FP_WFRACBITS_##fs) \
534 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
535 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
536 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
537 _FP_FRAC_SUB_##wc(R, X, Y); \
539 else if (ediff < 0) \
546 /* X is zero or denormalized. */ \
547 if (_FP_FRAC_ZEROP_##wc(X)) \
549 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
550 _FP_FRAC_COPY_##wc(R, Y); \
555 FP_SET_EXCEPTION(FP_EX_DENORM); \
559 _FP_FRAC_SUB_##wc(R, Y, X); \
562 if (Y##_e == _FP_EXPMAX_##fs) \
564 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
565 _FP_FRAC_COPY_##wc(R, Y); \
571 else if (Y##_e == _FP_EXPMAX_##fs) \
573 /* Y is NaN or Inf, X is normal. */ \
574 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
575 _FP_FRAC_COPY_##wc(R, Y); \
579 /* Insert implicit MSB of X. */ \
580 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
583 /* Shift the mantissa of X to the right EDIFF steps; \
584 remember to account later for the implicit MSB of Y. */ \
585 if (ediff <= _FP_WFRACBITS_##fs) \
586 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
587 else if (!_FP_FRAC_ZEROP_##wc(X)) \
588 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
589 _FP_FRAC_SUB_##wc(R, Y, X); \
594 if (!_FP_EXP_NORMAL(fs, wc, X)) \
598 /* X and Y are zero or denormalized. */ \
600 if (_FP_FRAC_ZEROP_##wc(X)) \
602 _FP_FRAC_COPY_##wc(R, Y); \
603 if (_FP_FRAC_ZEROP_##wc(Y)) \
604 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
607 FP_SET_EXCEPTION(FP_EX_DENORM); \
612 else if (_FP_FRAC_ZEROP_##wc(Y)) \
614 FP_SET_EXCEPTION(FP_EX_DENORM); \
615 _FP_FRAC_COPY_##wc(R, X); \
621 FP_SET_EXCEPTION(FP_EX_DENORM); \
622 _FP_FRAC_SUB_##wc(R, X, Y); \
624 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
626 /* |X| < |Y|, negate result. */ \
627 _FP_FRAC_SUB_##wc(R, Y, X); \
630 else if (_FP_FRAC_ZEROP_##wc(R)) \
631 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
637 /* X and Y are NaN or Inf, of opposite signs. */ \
638 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
639 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
640 R##_e = _FP_EXPMAX_##fs; \
641 if (_FP_FRAC_ZEROP_##wc(X)) \
643 if (_FP_FRAC_ZEROP_##wc(Y)) \
646 R##_s = _FP_NANSIGN_##fs; \
647 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
648 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
649 FP_SET_EXCEPTION(FP_EX_INVALID); \
655 _FP_FRAC_COPY_##wc(R, Y); \
660 if (_FP_FRAC_ZEROP_##wc(Y)) \
664 _FP_FRAC_COPY_##wc(R, X); \
669 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
675 /* The exponents of X and Y, both normal, are equal. The \
676 implicit MSBs cancel. */ \
678 _FP_FRAC_SUB_##wc(R, X, Y); \
680 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
682 /* |X| < |Y|, negate result. */ \
683 _FP_FRAC_SUB_##wc(R, Y, X); \
686 else if (_FP_FRAC_ZEROP_##wc(R)) \
689 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
695 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
698 /* Carry into most significant bit of larger one of X and Y, \
699 canceling it; renormalize. */ \
700 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \
702 _FP_FRAC_CLZ_##wc(diff, R); \
703 diff -= _FP_WFRACXBITS_##fs; \
704 _FP_FRAC_SLL_##wc(R, diff); \
707 /* R is denormalized. */ \
708 diff = diff - R##_e + 1; \
709 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \
715 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
722 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
723 #define _FP_SUB(fs, wc, R, X, Y) \
725 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
726 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \
731 * Main negation routine. FIXME -- when we care about setting exception
732 * bits reliably, this will not do. We should examine all of the fp classes.
735 #define _FP_NEG(fs, wc, R, X) \
737 _FP_FRAC_COPY_##wc(R, X); \
745 * Main multiplication routine. The input values should be cooked.
748 #define _FP_MUL(fs, wc, R, X, Y) \
750 R##_s = X##_s ^ Y##_s; \
751 R##_e = X##_e + Y##_e + 1; \
752 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
754 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
755 R##_c = FP_CLS_NORMAL; \
757 _FP_MUL_MEAT_##fs(R,X,Y); \
759 if (_FP_FRAC_OVERP_##wc(fs, R)) \
760 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
765 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
766 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \
769 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
770 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
771 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
774 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
775 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
776 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
777 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
778 _FP_FRAC_COPY_##wc(R, X); \
782 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
783 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
784 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
787 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
788 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
789 _FP_FRAC_COPY_##wc(R, Y); \
793 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
794 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
795 R##_s = _FP_NANSIGN_##fs; \
796 R##_c = FP_CLS_NAN; \
797 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
798 FP_SET_EXCEPTION(FP_EX_INVALID); \
808 * Main division routine. The input values should be cooked.
811 #define _FP_DIV(fs, wc, R, X, Y) \
813 R##_s = X##_s ^ Y##_s; \
814 R##_e = X##_e - Y##_e; \
815 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
817 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
818 R##_c = FP_CLS_NORMAL; \
820 _FP_DIV_MEAT_##fs(R,X,Y); \
823 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
824 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \
827 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
828 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
829 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
831 _FP_FRAC_COPY_##wc(R, X); \
835 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
836 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
837 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
839 _FP_FRAC_COPY_##wc(R, Y); \
843 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
844 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
845 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
846 R##_c = FP_CLS_ZERO; \
849 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
850 FP_SET_EXCEPTION(FP_EX_DIVZERO); \
851 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
852 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
853 R##_c = FP_CLS_INF; \
856 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
857 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
858 R##_s = _FP_NANSIGN_##fs; \
859 R##_c = FP_CLS_NAN; \
860 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
861 FP_SET_EXCEPTION(FP_EX_INVALID); \
871 * Main differential comparison routine. The inputs should be raw not
872 * cooked. The return is -1,0,1 for normal values, 2 otherwise.
875 #define _FP_CMP(fs, wc, ret, X, Y, un) \
877 /* NANs are unordered */ \
878 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
879 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
888 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \
889 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \
891 if (__is_zero_x && __is_zero_y) \
893 else if (__is_zero_x) \
894 ret = Y##_s ? 1 : -1; \
895 else if (__is_zero_y) \
896 ret = X##_s ? -1 : 1; \
897 else if (X##_s != Y##_s) \
898 ret = X##_s ? -1 : 1; \
899 else if (X##_e > Y##_e) \
900 ret = X##_s ? -1 : 1; \
901 else if (X##_e < Y##_e) \
902 ret = X##_s ? 1 : -1; \
903 else if (_FP_FRAC_GT_##wc(X, Y)) \
904 ret = X##_s ? -1 : 1; \
905 else if (_FP_FRAC_GT_##wc(Y, X)) \
906 ret = X##_s ? 1 : -1; \
913 /* Simplification for strict equality. */
915 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
917 /* NANs are unordered */ \
918 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
919 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
925 ret = !(X##_e == Y##_e \
926 && _FP_FRAC_EQ_##wc(X, Y) \
927 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
931 /* Version to test unordered. */
933 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
935 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
936 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \
940 * Main square root routine. The input value should be cooked.
943 #define _FP_SQRT(fs, wc, R, X) \
945 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \
950 _FP_FRAC_COPY_##wc(R, X); \
952 R##_c = FP_CLS_NAN; \
957 R##_s = _FP_NANSIGN_##fs; \
958 R##_c = FP_CLS_NAN; /* NAN */ \
959 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
960 FP_SET_EXCEPTION(FP_EX_INVALID); \
965 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
970 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
972 case FP_CLS_NORMAL: \
976 R##_c = FP_CLS_NAN; /* NAN */ \
977 R##_s = _FP_NANSIGN_##fs; \
978 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
979 FP_SET_EXCEPTION(FP_EX_INVALID); \
982 R##_c = FP_CLS_NORMAL; \
984 _FP_FRAC_SLL_##wc(X, 1); \
985 R##_e = X##_e >> 1; \
986 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \
987 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \
988 q = _FP_OVERFLOW_##fs >> 1; \
989 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \
994 * Convert from FP to integer. Input is raw.
997 /* RSIGNED can have following values:
998 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
999 * the result is either 0 or (2^rsize)-1 depending on the sign in such
1001 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1002 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1003 * depending on the sign in such case.
1004 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1005 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1006 * depending on the sign in such case.
1008 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1010 if (X##_e < _FP_EXPBIAS_##fs) \
1015 if (!_FP_FRAC_ZEROP_##wc(X)) \
1017 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1018 FP_SET_EXCEPTION(FP_EX_DENORM); \
1022 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1024 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1025 || (!rsigned && X##_s)) \
1027 /* Overflow or converting to the most negative integer. */ \
1039 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1041 /* Possibly converting to most negative integer; check the \
1044 (void)((_FP_FRACBITS_##fs > rsize) \
1045 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \
1046 _FP_FRACBITS_##fs - rsize, \
1047 _FP_FRACBITS_##fs); 0; }) \
1049 if (!_FP_FRAC_ZEROP_##wc(X)) \
1050 FP_SET_EXCEPTION(FP_EX_INVALID); \
1052 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1055 FP_SET_EXCEPTION(FP_EX_INVALID); \
1059 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
1060 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1062 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1063 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1068 _FP_FRAC_SRST_##wc(X, inexact, \
1069 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1071 _FP_FRACBITS_##fs); \
1073 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1074 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1076 if (rsigned && X##_s) \
1081 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1083 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1089 if ((X##_s = (r < 0))) \
1093 (void)((rsize <= _FP_W_TYPE_SIZE) \
1096 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \
1097 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1099 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1102 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
1104 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1109 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1110 && X##_e >= _FP_EXPMAX_##fs) \
1112 /* Exponent too big; overflow to infinity. (May also \
1113 happen after rounding below.) */ \
1114 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
1115 goto pack_semiraw; \
1118 if (rsize <= _FP_FRACBITS_##fs \
1119 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1121 /* Exactly representable; shift left. */ \
1122 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1123 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1124 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1128 /* More bits in integer than in floating type; need to \
1130 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1131 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \
1132 - _FP_WFRACBITS_##fs + 1)) \
1133 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1134 - _FP_WFRACBITS_##fs + 1))) \
1136 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1137 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1138 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1139 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1140 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
1142 _FP_PACK_SEMIRAW(fs, wc, X); \
1149 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
1154 /* Extend from a narrower floating-point format to a wider one. Input
1155 and output are raw. */
1156 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \
1158 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1159 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1160 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1161 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1162 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
1165 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1166 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1168 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1169 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1175 if (_FP_FRAC_ZEROP_##swc(S)) \
1177 else if (_FP_EXPBIAS_##dfs \
1178 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1180 FP_SET_EXCEPTION(FP_EX_DENORM); \
1181 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1182 - _FP_FRACBITS_##sfs)); \
1188 FP_SET_EXCEPTION(FP_EX_DENORM); \
1189 _FP_FRAC_CLZ_##swc(_lz, S); \
1190 _FP_FRAC_SLL_##dwc(D, \
1191 _lz + _FP_FRACBITS_##dfs \
1192 - _FP_FRACTBITS_##sfs); \
1193 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1194 + _FP_FRACXBITS_##sfs - _lz); \
1199 D##_e = _FP_EXPMAX_##dfs; \
1200 if (!_FP_FRAC_ZEROP_##swc(S)) \
1202 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1203 FP_SET_EXCEPTION(FP_EX_INVALID); \
1204 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1205 - _FP_FRACBITS_##sfs)); \
1211 /* Truncate from a wider floating-point format to a narrower one.
1212 Input and output are semi-raw. */
1213 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1215 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1216 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1217 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1220 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1222 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1223 if (D##_e >= _FP_EXPMAX_##dfs) \
1224 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1229 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1231 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1232 _FP_FRAC_LOW_##swc(S) |= 1; \
1236 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1237 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1238 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1239 _FP_WFRACBITS_##sfs); \
1244 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1245 - _FP_WFRACBITS_##dfs), \
1246 _FP_WFRACBITS_##sfs); \
1247 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1255 if (_FP_FRAC_ZEROP_##swc(S)) \
1256 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1259 FP_SET_EXCEPTION(FP_EX_DENORM); \
1260 if (_FP_EXPBIAS_##sfs \
1261 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1263 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1264 - _FP_WFRACBITS_##dfs), \
1265 _FP_WFRACBITS_##sfs); \
1266 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1270 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1271 _FP_FRAC_LOW_##dwc(D) |= 1; \
1277 D##_e = _FP_EXPMAX_##dfs; \
1278 if (_FP_FRAC_ZEROP_##swc(S)) \
1279 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1282 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1283 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1284 - _FP_WFRACBITS_##dfs)); \
1285 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1286 /* Semi-raw NaN must have all workbits cleared. */ \
1287 _FP_FRAC_LOW_##dwc(D) \
1288 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1289 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1296 * Helper primitives.
1299 /* Count leading zeros in a word. */
1302 /* GCC 3.4 and later provide the builtins for us. */
1303 #define __FP_CLZ(r, x) \
1305 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1306 r = __builtin_clz (x); \
1307 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1308 r = __builtin_clzl (x); \
1309 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1310 r = __builtin_clzll (x); \
1314 #endif /* ndef __FP_CLZ */
1316 #define _FP_DIV_HELP_imm(q, r, n, d) \
1318 q = n / d, r = n % d; \
1322 /* A restoring bit-by-bit division primitive. */
1324 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1326 int count = _FP_WFRACBITS_##fs; \
1327 _FP_FRAC_DECL_##wc (u); \
1328 _FP_FRAC_DECL_##wc (v); \
1329 _FP_FRAC_COPY_##wc (u, X); \
1330 _FP_FRAC_COPY_##wc (v, Y); \
1331 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1332 /* Normalize U and V. */ \
1333 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1334 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1335 /* First round. Since the operands are normalized, either the \
1336 first or second bit will be set in the fraction. Produce a \
1337 normalized result by checking which and adjusting the loop \
1338 count and exponent accordingly. */ \
1339 if (_FP_FRAC_GE_1 (u, v)) \
1341 _FP_FRAC_SUB_##wc (u, u, v); \
1342 _FP_FRAC_LOW_##wc (R) |= 1; \
1347 /* Subsequent rounds. */ \
1349 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1350 _FP_FRAC_SLL_##wc (u, 1); \
1351 _FP_FRAC_SLL_##wc (R, 1); \
1352 if (msb || _FP_FRAC_GE_1 (u, v)) \
1354 _FP_FRAC_SUB_##wc (u, u, v); \
1355 _FP_FRAC_LOW_##wc (R) |= 1; \
1357 } while (--count > 0); \
1358 /* If there's anything left in U, the result is inexact. */ \
1359 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1362 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1363 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1364 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)