1 /* Software floating-point emulation. Common operations.
2 Copyright (C) 1997,1998,1999,2006,2007,2012 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)), X##_s, X##_e; \
37 * Finish truely unpacking a native fp value by classifying the kind
38 * of fp value and normalizing both the exponent and the fraction.
41 #define _FP_UNPACK_CANONICAL(fs, wc, X) \
46 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
47 _FP_FRAC_SLL_##wc(X, _FP_WORKBITS); \
48 X##_e -= _FP_EXPBIAS_##fs; \
49 X##_c = FP_CLS_NORMAL; \
53 if (_FP_FRAC_ZEROP_##wc(X)) \
54 X##_c = FP_CLS_ZERO; \
57 /* a denormalized number */ \
59 _FP_FRAC_CLZ_##wc(_shift, X); \
60 _shift -= _FP_FRACXBITS_##fs; \
61 _FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS)); \
62 X##_e -= _FP_EXPBIAS_##fs - 1 + _shift; \
63 X##_c = FP_CLS_NORMAL; \
64 FP_SET_EXCEPTION(FP_EX_DENORM); \
68 case _FP_EXPMAX_##fs: \
69 if (_FP_FRAC_ZEROP_##wc(X)) \
74 /* Check for signaling NaN */ \
75 if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
76 FP_SET_EXCEPTION(FP_EX_INVALID); \
82 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
83 shifted by _FP_WORKBITS but the implicit MSB is not inserted and
84 other classification is not done. */
85 #define _FP_UNPACK_SEMIRAW(fs, wc, X) _FP_FRAC_SLL_##wc(X, _FP_WORKBITS)
87 /* A semi-raw value has overflowed to infinity. Adjust the mantissa
88 and exponent appropriately. */
89 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X) \
91 if (FP_ROUNDMODE == FP_RND_NEAREST \
92 || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \
93 || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \
95 X##_e = _FP_EXPMAX_##fs; \
96 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
100 X##_e = _FP_EXPMAX_##fs - 1; \
101 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
103 FP_SET_EXCEPTION(FP_EX_INEXACT); \
104 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
107 /* Check for a semi-raw value being a signaling NaN and raise the
108 invalid exception if so. */
109 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X) \
111 if (X##_e == _FP_EXPMAX_##fs \
112 && !_FP_FRAC_ZEROP_##wc(X) \
113 && !(_FP_FRAC_HIGH_##fs(X) & _FP_QNANBIT_SH_##fs)) \
114 FP_SET_EXCEPTION(FP_EX_INVALID); \
117 /* Choose a NaN result from an operation on two semi-raw NaN
119 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP) \
121 /* _FP_CHOOSENAN expects raw values, so shift as required. */ \
122 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
123 _FP_FRAC_SRL_##wc(Y, _FP_WORKBITS); \
124 _FP_CHOOSENAN(fs, wc, R, X, Y, OP); \
125 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
128 /* Test whether a biased exponent is normal (not zero or maximum). */
129 #define _FP_EXP_NORMAL(fs, wc, X) (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
131 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
132 rounded and shifted right, with the rounding possibly increasing
133 the exponent (including changing a finite value to infinity). */
134 #define _FP_PACK_SEMIRAW(fs, wc, X) \
137 if (X##_e == 0 && !_FP_FRAC_ZEROP_##wc(X)) \
139 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
140 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
141 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
143 if (_FP_FRAC_HIGH_##fs(X) \
144 & (_FP_OVERFLOW_##fs >> 1)) \
146 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1); \
148 if (X##_e == _FP_EXPMAX_##fs) \
149 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
151 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
152 if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
154 if (!_FP_KEEPNANFRACP) \
156 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
157 X##_s = _FP_NANSIGN_##fs; \
160 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
165 * Before packing the bits back into the native fp result, take care
166 * of such mundane things as rounding and overflow. Also, for some
167 * kinds of fp values, the original parts may not have been fully
168 * extracted -- but that is ok, we can regenerate them now.
171 #define _FP_PACK_CANONICAL(fs, wc, X) \
175 case FP_CLS_NORMAL: \
176 X##_e += _FP_EXPBIAS_##fs; \
180 if (_FP_FRAC_OVERP_##wc(fs, X)) \
182 _FP_FRAC_CLEAR_OVERP_##wc(fs, X); \
185 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
186 if (X##_e >= _FP_EXPMAX_##fs) \
189 switch (FP_ROUNDMODE) \
191 case FP_RND_NEAREST: \
192 X##_c = FP_CLS_INF; \
195 if (!X##_s) X##_c = FP_CLS_INF; \
198 if (X##_s) X##_c = FP_CLS_INF; \
201 if (X##_c == FP_CLS_INF) \
203 /* Overflow to infinity */ \
204 X##_e = _FP_EXPMAX_##fs; \
205 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
209 /* Overflow to maximum normal */ \
210 X##_e = _FP_EXPMAX_##fs - 1; \
211 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
213 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
214 FP_SET_EXCEPTION(FP_EX_INEXACT); \
219 /* we've got a denormalized number */ \
220 X##_e = -X##_e + 1; \
221 if (X##_e <= _FP_WFRACBITS_##fs) \
223 _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs); \
225 if (_FP_FRAC_HIGH_##fs(X) \
226 & (_FP_OVERFLOW_##fs >> 1)) \
229 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
230 FP_SET_EXCEPTION(FP_EX_INEXACT); \
235 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
237 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
238 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
239 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
243 /* underflow to zero */ \
245 if (!_FP_FRAC_ZEROP_##wc(X)) \
247 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
249 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
251 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
258 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
262 X##_e = _FP_EXPMAX_##fs; \
263 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
267 X##_e = _FP_EXPMAX_##fs; \
268 if (!_FP_KEEPNANFRACP) \
270 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
271 X##_s = _FP_NANSIGN_##fs; \
274 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
279 /* This one accepts raw argument and not cooked, returns
280 * 1 if X is a signaling NaN.
282 #define _FP_ISSIGNAN(fs, wc, X) \
285 if (X##_e == _FP_EXPMAX_##fs) \
287 if (!_FP_FRAC_ZEROP_##wc(X) \
288 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
298 /* Addition on semi-raw values. */
299 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
301 if (X##_s == Y##_s) \
305 int ediff = X##_e - Y##_e; \
311 /* Y is zero or denormalized. */ \
312 if (_FP_FRAC_ZEROP_##wc(Y)) \
314 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
315 _FP_FRAC_COPY_##wc(R, X); \
320 FP_SET_EXCEPTION(FP_EX_DENORM); \
324 _FP_FRAC_ADD_##wc(R, X, Y); \
327 if (X##_e == _FP_EXPMAX_##fs) \
329 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
330 _FP_FRAC_COPY_##wc(R, X); \
336 else if (X##_e == _FP_EXPMAX_##fs) \
338 /* X is NaN or Inf, Y is normal. */ \
339 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
340 _FP_FRAC_COPY_##wc(R, X); \
344 /* Insert implicit MSB of Y. */ \
345 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
348 /* Shift the mantissa of Y to the right EDIFF steps; \
349 remember to account later for the implicit MSB of X. */ \
350 if (ediff <= _FP_WFRACBITS_##fs) \
351 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
352 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
353 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
354 _FP_FRAC_ADD_##wc(R, X, Y); \
356 else if (ediff < 0) \
362 /* X is zero or denormalized. */ \
363 if (_FP_FRAC_ZEROP_##wc(X)) \
365 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
366 _FP_FRAC_COPY_##wc(R, Y); \
371 FP_SET_EXCEPTION(FP_EX_DENORM); \
375 _FP_FRAC_ADD_##wc(R, Y, X); \
378 if (Y##_e == _FP_EXPMAX_##fs) \
380 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
381 _FP_FRAC_COPY_##wc(R, Y); \
387 else if (Y##_e == _FP_EXPMAX_##fs) \
389 /* Y is NaN or Inf, X is normal. */ \
390 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
391 _FP_FRAC_COPY_##wc(R, Y); \
395 /* Insert implicit MSB of X. */ \
396 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
399 /* Shift the mantissa of X to the right EDIFF steps; \
400 remember to account later for the implicit MSB of Y. */ \
401 if (ediff <= _FP_WFRACBITS_##fs) \
402 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
403 else if (!_FP_FRAC_ZEROP_##wc(X)) \
404 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
405 _FP_FRAC_ADD_##wc(R, Y, X); \
410 if (!_FP_EXP_NORMAL(fs, wc, X)) \
414 /* X and Y are zero or denormalized. */ \
416 if (_FP_FRAC_ZEROP_##wc(X)) \
418 if (!_FP_FRAC_ZEROP_##wc(Y)) \
419 FP_SET_EXCEPTION(FP_EX_DENORM); \
420 _FP_FRAC_COPY_##wc(R, Y); \
423 else if (_FP_FRAC_ZEROP_##wc(Y)) \
425 FP_SET_EXCEPTION(FP_EX_DENORM); \
426 _FP_FRAC_COPY_##wc(R, X); \
431 FP_SET_EXCEPTION(FP_EX_DENORM); \
432 _FP_FRAC_ADD_##wc(R, X, Y); \
433 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
435 /* Normalized result. */ \
436 _FP_FRAC_HIGH_##fs(R) \
437 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
445 /* X and Y are NaN or Inf. */ \
446 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
447 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
448 R##_e = _FP_EXPMAX_##fs; \
449 if (_FP_FRAC_ZEROP_##wc(X)) \
450 _FP_FRAC_COPY_##wc(R, Y); \
451 else if (_FP_FRAC_ZEROP_##wc(Y)) \
452 _FP_FRAC_COPY_##wc(R, X); \
454 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
458 /* The exponents of X and Y, both normal, are equal. The \
459 implicit MSBs will always add to increase the \
461 _FP_FRAC_ADD_##wc(R, X, Y); \
463 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
464 if (R##_e == _FP_EXPMAX_##fs) \
465 /* Overflow to infinity (depending on rounding mode). */ \
466 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
470 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
473 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
475 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
476 if (R##_e == _FP_EXPMAX_##fs) \
477 /* Overflow to infinity (depending on rounding mode). */ \
478 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
485 int ediff = X##_e - Y##_e; \
492 /* Y is zero or denormalized. */ \
493 if (_FP_FRAC_ZEROP_##wc(Y)) \
495 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
496 _FP_FRAC_COPY_##wc(R, X); \
501 FP_SET_EXCEPTION(FP_EX_DENORM); \
505 _FP_FRAC_SUB_##wc(R, X, Y); \
508 if (X##_e == _FP_EXPMAX_##fs) \
510 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
511 _FP_FRAC_COPY_##wc(R, X); \
517 else if (X##_e == _FP_EXPMAX_##fs) \
519 /* X is NaN or Inf, Y is normal. */ \
520 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
521 _FP_FRAC_COPY_##wc(R, X); \
525 /* Insert implicit MSB of Y. */ \
526 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
529 /* Shift the mantissa of Y to the right EDIFF steps; \
530 remember to account later for the implicit MSB of X. */ \
531 if (ediff <= _FP_WFRACBITS_##fs) \
532 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
533 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
534 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
535 _FP_FRAC_SUB_##wc(R, X, Y); \
537 else if (ediff < 0) \
544 /* X is zero or denormalized. */ \
545 if (_FP_FRAC_ZEROP_##wc(X)) \
547 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
548 _FP_FRAC_COPY_##wc(R, Y); \
553 FP_SET_EXCEPTION(FP_EX_DENORM); \
557 _FP_FRAC_SUB_##wc(R, Y, X); \
560 if (Y##_e == _FP_EXPMAX_##fs) \
562 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
563 _FP_FRAC_COPY_##wc(R, Y); \
569 else if (Y##_e == _FP_EXPMAX_##fs) \
571 /* Y is NaN or Inf, X is normal. */ \
572 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
573 _FP_FRAC_COPY_##wc(R, Y); \
577 /* Insert implicit MSB of X. */ \
578 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
581 /* Shift the mantissa of X to the right EDIFF steps; \
582 remember to account later for the implicit MSB of Y. */ \
583 if (ediff <= _FP_WFRACBITS_##fs) \
584 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
585 else if (!_FP_FRAC_ZEROP_##wc(X)) \
586 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
587 _FP_FRAC_SUB_##wc(R, Y, X); \
592 if (!_FP_EXP_NORMAL(fs, wc, X)) \
596 /* X and Y are zero or denormalized. */ \
598 if (_FP_FRAC_ZEROP_##wc(X)) \
600 _FP_FRAC_COPY_##wc(R, Y); \
601 if (_FP_FRAC_ZEROP_##wc(Y)) \
602 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
605 FP_SET_EXCEPTION(FP_EX_DENORM); \
610 else if (_FP_FRAC_ZEROP_##wc(Y)) \
612 FP_SET_EXCEPTION(FP_EX_DENORM); \
613 _FP_FRAC_COPY_##wc(R, X); \
619 FP_SET_EXCEPTION(FP_EX_DENORM); \
620 _FP_FRAC_SUB_##wc(R, X, Y); \
622 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
624 /* |X| < |Y|, negate result. */ \
625 _FP_FRAC_SUB_##wc(R, Y, X); \
628 else if (_FP_FRAC_ZEROP_##wc(R)) \
629 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
635 /* X and Y are NaN or Inf, of opposite signs. */ \
636 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
637 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
638 R##_e = _FP_EXPMAX_##fs; \
639 if (_FP_FRAC_ZEROP_##wc(X)) \
641 if (_FP_FRAC_ZEROP_##wc(Y)) \
644 R##_s = _FP_NANSIGN_##fs; \
645 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
646 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
647 FP_SET_EXCEPTION(FP_EX_INVALID); \
653 _FP_FRAC_COPY_##wc(R, Y); \
658 if (_FP_FRAC_ZEROP_##wc(Y)) \
662 _FP_FRAC_COPY_##wc(R, X); \
667 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
673 /* The exponents of X and Y, both normal, are equal. The \
674 implicit MSBs cancel. */ \
676 _FP_FRAC_SUB_##wc(R, X, Y); \
678 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
680 /* |X| < |Y|, negate result. */ \
681 _FP_FRAC_SUB_##wc(R, Y, X); \
684 else if (_FP_FRAC_ZEROP_##wc(R)) \
687 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
693 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
696 /* Carry into most significant bit of larger one of X and Y, \
697 canceling it; renormalize. */ \
698 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \
700 _FP_FRAC_CLZ_##wc(diff, R); \
701 diff -= _FP_WFRACXBITS_##fs; \
702 _FP_FRAC_SLL_##wc(R, diff); \
705 /* R is denormalized. */ \
706 diff = diff - R##_e + 1; \
707 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \
713 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
720 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
721 #define _FP_SUB(fs, wc, R, X, Y) \
723 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
724 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \
729 * Main negation routine. FIXME -- when we care about setting exception
730 * bits reliably, this will not do. We should examine all of the fp classes.
733 #define _FP_NEG(fs, wc, R, X) \
735 _FP_FRAC_COPY_##wc(R, X); \
743 * Main multiplication routine. The input values should be cooked.
746 #define _FP_MUL(fs, wc, R, X, Y) \
748 R##_s = X##_s ^ Y##_s; \
749 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
751 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
752 R##_c = FP_CLS_NORMAL; \
753 R##_e = X##_e + Y##_e + 1; \
755 _FP_MUL_MEAT_##fs(R,X,Y); \
757 if (_FP_FRAC_OVERP_##wc(fs, R)) \
758 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
763 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
764 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \
767 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
768 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
769 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
772 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
773 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
774 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
775 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
776 _FP_FRAC_COPY_##wc(R, X); \
780 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
781 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
782 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
785 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
786 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
787 _FP_FRAC_COPY_##wc(R, Y); \
791 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
792 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
793 R##_s = _FP_NANSIGN_##fs; \
794 R##_c = FP_CLS_NAN; \
795 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
796 FP_SET_EXCEPTION(FP_EX_INVALID); \
806 * Main division routine. The input values should be cooked.
809 #define _FP_DIV(fs, wc, R, X, Y) \
811 R##_s = X##_s ^ Y##_s; \
812 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
814 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
815 R##_c = FP_CLS_NORMAL; \
816 R##_e = X##_e - Y##_e; \
818 _FP_DIV_MEAT_##fs(R,X,Y); \
821 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
822 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \
825 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
826 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
827 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
829 _FP_FRAC_COPY_##wc(R, X); \
833 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
834 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
835 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
837 _FP_FRAC_COPY_##wc(R, Y); \
841 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
842 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
843 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
844 R##_c = FP_CLS_ZERO; \
847 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
848 FP_SET_EXCEPTION(FP_EX_DIVZERO); \
849 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
850 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
851 R##_c = FP_CLS_INF; \
854 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
855 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
856 R##_s = _FP_NANSIGN_##fs; \
857 R##_c = FP_CLS_NAN; \
858 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
859 FP_SET_EXCEPTION(FP_EX_INVALID); \
869 * Main differential comparison routine. The inputs should be raw not
870 * cooked. The return is -1,0,1 for normal values, 2 otherwise.
873 #define _FP_CMP(fs, wc, ret, X, Y, un) \
875 /* NANs are unordered */ \
876 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
877 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
886 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \
887 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \
889 if (__is_zero_x && __is_zero_y) \
891 else if (__is_zero_x) \
892 ret = Y##_s ? 1 : -1; \
893 else if (__is_zero_y) \
894 ret = X##_s ? -1 : 1; \
895 else if (X##_s != Y##_s) \
896 ret = X##_s ? -1 : 1; \
897 else if (X##_e > Y##_e) \
898 ret = X##_s ? -1 : 1; \
899 else if (X##_e < Y##_e) \
900 ret = X##_s ? 1 : -1; \
901 else if (_FP_FRAC_GT_##wc(X, Y)) \
902 ret = X##_s ? -1 : 1; \
903 else if (_FP_FRAC_GT_##wc(Y, X)) \
904 ret = X##_s ? 1 : -1; \
911 /* Simplification for strict equality. */
913 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
915 /* NANs are unordered */ \
916 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
917 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
923 ret = !(X##_e == Y##_e \
924 && _FP_FRAC_EQ_##wc(X, Y) \
925 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
929 /* Version to test unordered. */
931 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
933 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
934 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \
938 * Main square root routine. The input value should be cooked.
941 #define _FP_SQRT(fs, wc, R, X) \
943 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \
948 _FP_FRAC_COPY_##wc(R, X); \
950 R##_c = FP_CLS_NAN; \
955 R##_s = _FP_NANSIGN_##fs; \
956 R##_c = FP_CLS_NAN; /* NAN */ \
957 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
958 FP_SET_EXCEPTION(FP_EX_INVALID); \
963 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
968 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
970 case FP_CLS_NORMAL: \
974 R##_c = FP_CLS_NAN; /* sNAN */ \
975 R##_s = _FP_NANSIGN_##fs; \
976 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
977 FP_SET_EXCEPTION(FP_EX_INVALID); \
980 R##_c = FP_CLS_NORMAL; \
982 _FP_FRAC_SLL_##wc(X, 1); \
983 R##_e = X##_e >> 1; \
984 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \
985 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \
986 q = _FP_OVERFLOW_##fs >> 1; \
987 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \
992 * Convert from FP to integer. Input is raw.
995 /* RSIGNED can have following values:
996 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
997 * the result is either 0 or (2^rsize)-1 depending on the sign in such
999 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1000 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1001 * depending on the sign in such case.
1002 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1003 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1004 * depending on the sign in such case.
1006 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1008 if (X##_e < _FP_EXPBIAS_##fs) \
1013 if (!_FP_FRAC_ZEROP_##wc(X)) \
1015 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1016 FP_SET_EXCEPTION(FP_EX_DENORM); \
1020 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1022 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1023 || (!rsigned && X##_s)) \
1025 /* Overflow or converting to the most negative integer. */ \
1037 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1039 /* Possibly converting to most negative integer; check the \
1042 (void)((_FP_FRACBITS_##fs > rsize) \
1043 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \
1044 _FP_FRACBITS_##fs - rsize, \
1045 _FP_FRACBITS_##fs); 0; }) \
1047 if (!_FP_FRAC_ZEROP_##wc(X)) \
1048 FP_SET_EXCEPTION(FP_EX_INVALID); \
1050 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1053 FP_SET_EXCEPTION(FP_EX_INVALID); \
1057 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
1058 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1060 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1061 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1066 _FP_FRAC_SRST_##wc(X, inexact, \
1067 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1069 _FP_FRACBITS_##fs); \
1071 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1072 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1074 if (rsigned && X##_s) \
1079 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1081 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1087 if ((X##_s = (r < 0))) \
1091 (void)((rsize <= _FP_W_TYPE_SIZE) \
1094 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \
1095 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1097 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1100 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
1102 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1107 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1108 && X##_e >= _FP_EXPMAX_##fs) \
1110 /* Exponent too big; overflow to infinity. (May also \
1111 happen after rounding below.) */ \
1112 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
1113 goto pack_semiraw; \
1116 if (rsize <= _FP_FRACBITS_##fs \
1117 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1119 /* Exactly representable; shift left. */ \
1120 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1121 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1122 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1126 /* More bits in integer than in floating type; need to \
1128 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1129 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \
1130 - _FP_WFRACBITS_##fs + 1)) \
1131 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1132 - _FP_WFRACBITS_##fs + 1))) \
1134 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1135 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1136 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1137 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1138 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
1140 _FP_PACK_SEMIRAW(fs, wc, X); \
1147 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
1152 /* Extend from a narrower floating-point format to a wider one. Input
1153 and output are raw. */
1154 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \
1156 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1157 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1158 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1159 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1160 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
1163 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1164 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1166 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1167 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1173 if (_FP_FRAC_ZEROP_##swc(S)) \
1175 else if (_FP_EXPBIAS_##dfs \
1176 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1178 FP_SET_EXCEPTION(FP_EX_DENORM); \
1179 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1180 - _FP_FRACBITS_##sfs)); \
1186 FP_SET_EXCEPTION(FP_EX_DENORM); \
1187 _FP_FRAC_CLZ_##swc(_lz, S); \
1188 _FP_FRAC_SLL_##dwc(D, \
1189 _lz + _FP_FRACBITS_##dfs \
1190 - _FP_FRACTBITS_##sfs); \
1191 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1192 + _FP_FRACXBITS_##sfs - _lz); \
1197 D##_e = _FP_EXPMAX_##dfs; \
1198 if (!_FP_FRAC_ZEROP_##swc(S)) \
1200 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1201 FP_SET_EXCEPTION(FP_EX_INVALID); \
1202 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1203 - _FP_FRACBITS_##sfs)); \
1209 /* Truncate from a wider floating-point format to a narrower one.
1210 Input and output are semi-raw. */
1211 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1213 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1214 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1215 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1218 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1220 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1221 if (D##_e >= _FP_EXPMAX_##dfs) \
1222 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1227 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1229 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1230 _FP_FRAC_LOW_##swc(S) |= 1; \
1234 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1235 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1236 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1237 _FP_WFRACBITS_##sfs); \
1242 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1243 - _FP_WFRACBITS_##dfs), \
1244 _FP_WFRACBITS_##sfs); \
1245 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1253 if (_FP_FRAC_ZEROP_##swc(S)) \
1254 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1257 FP_SET_EXCEPTION(FP_EX_DENORM); \
1258 if (_FP_EXPBIAS_##sfs \
1259 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1261 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1262 - _FP_WFRACBITS_##dfs), \
1263 _FP_WFRACBITS_##sfs); \
1264 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1268 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1269 _FP_FRAC_LOW_##dwc(D) |= 1; \
1275 D##_e = _FP_EXPMAX_##dfs; \
1276 if (_FP_FRAC_ZEROP_##swc(S)) \
1277 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1280 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1281 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1282 - _FP_WFRACBITS_##dfs)); \
1283 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1284 /* Semi-raw NaN must have all workbits cleared. */ \
1285 _FP_FRAC_LOW_##dwc(D) \
1286 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1287 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1294 * Helper primitives.
1297 /* Count leading zeros in a word. */
1300 /* GCC 3.4 and later provide the builtins for us. */
1301 #define __FP_CLZ(r, x) \
1303 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1304 r = __builtin_clz (x); \
1305 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1306 r = __builtin_clzl (x); \
1307 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1308 r = __builtin_clzll (x); \
1312 #endif /* ndef __FP_CLZ */
1314 #define _FP_DIV_HELP_imm(q, r, n, d) \
1316 q = n / d, r = n % d; \
1320 /* A restoring bit-by-bit division primitive. */
1322 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1324 int count = _FP_WFRACBITS_##fs; \
1325 _FP_FRAC_DECL_##wc (u); \
1326 _FP_FRAC_DECL_##wc (v); \
1327 _FP_FRAC_COPY_##wc (u, X); \
1328 _FP_FRAC_COPY_##wc (v, Y); \
1329 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1330 /* Normalize U and V. */ \
1331 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1332 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1333 /* First round. Since the operands are normalized, either the \
1334 first or second bit will be set in the fraction. Produce a \
1335 normalized result by checking which and adjusting the loop \
1336 count and exponent accordingly. */ \
1337 if (_FP_FRAC_GE_1 (u, v)) \
1339 _FP_FRAC_SUB_##wc (u, u, v); \
1340 _FP_FRAC_LOW_##wc (R) |= 1; \
1345 /* Subsequent rounds. */ \
1347 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1348 _FP_FRAC_SLL_##wc (u, 1); \
1349 _FP_FRAC_SLL_##wc (R, 1); \
1350 if (msb || _FP_FRAC_GE_1 (u, v)) \
1352 _FP_FRAC_SUB_##wc (u, u, v); \
1353 _FP_FRAC_LOW_##wc (R) |= 1; \
1355 } while (--count > 0); \
1356 /* If there's anything left in U, the result is inexact. */ \
1357 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1360 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1361 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1362 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)