1 /* Software floating-point emulation. Common operations.
2 Copyright (C) 1997,1998,1999,2006 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, write to the Free
30 Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
31 MA 02110-1301, USA. */
33 #define _FP_DECL(wc, X) \
34 _FP_I_TYPE X##_c __attribute__((unused)), X##_s, X##_e; \
38 * Finish truely unpacking a native fp value by classifying the kind
39 * of fp value and normalizing both the exponent and the fraction.
42 #define _FP_UNPACK_CANONICAL(fs, wc, X) \
47 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
48 _FP_FRAC_SLL_##wc(X, _FP_WORKBITS); \
49 X##_e -= _FP_EXPBIAS_##fs; \
50 X##_c = FP_CLS_NORMAL; \
54 if (_FP_FRAC_ZEROP_##wc(X)) \
55 X##_c = FP_CLS_ZERO; \
58 /* a denormalized number */ \
60 _FP_FRAC_CLZ_##wc(_shift, X); \
61 _shift -= _FP_FRACXBITS_##fs; \
62 _FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS)); \
63 X##_e -= _FP_EXPBIAS_##fs - 1 + _shift; \
64 X##_c = FP_CLS_NORMAL; \
65 FP_SET_EXCEPTION(FP_EX_DENORM); \
69 case _FP_EXPMAX_##fs: \
70 if (_FP_FRAC_ZEROP_##wc(X)) \
75 /* Check for signaling NaN */ \
76 if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
77 FP_SET_EXCEPTION(FP_EX_INVALID); \
83 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
84 shifted by _FP_WORKBITS but the implicit MSB is not inserted and
85 other classification is not done. */
86 #define _FP_UNPACK_SEMIRAW(fs, wc, X) _FP_FRAC_SLL_##wc(X, _FP_WORKBITS)
88 /* A semi-raw value has overflowed to infinity. Adjust the mantissa
89 and exponent appropriately. */
90 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X) \
92 if (FP_ROUNDMODE == FP_RND_NEAREST \
93 || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \
94 || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \
96 X##_e = _FP_EXPMAX_##fs; \
97 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
101 X##_e = _FP_EXPMAX_##fs - 1; \
102 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
103 FP_SET_EXCEPTION(FP_EX_INEXACT); \
104 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
108 /* Check for a semi-raw value being a signaling NaN and raise the
109 invalid exception if so. */
110 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X) \
112 if (X##_e == _FP_EXPMAX_##fs \
113 && !_FP_FRAC_ZEROP_##wc(X) \
114 && !(_FP_FRAC_HIGH_##fs(X) & _FP_QNANBIT_SH_##fs)) \
115 FP_SET_EXCEPTION(FP_EX_INVALID); \
118 /* Choose a NaN result from an operation on two semi-raw NaN
120 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP) \
122 /* _FP_CHOOSENAN expects raw values, so shift as required. */ \
123 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
124 _FP_FRAC_SRL_##wc(Y, _FP_WORKBITS); \
125 _FP_CHOOSENAN(fs, wc, R, X, Y, OP); \
126 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
129 /* Test whether a biased exponent is normal (not zero or maximum). */
130 #define _FP_EXP_NORMAL(fs, wc, X) (((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
132 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
133 rounded and shifted right, with the rounding possibly increasing
134 the exponent (including changing a finite value to infinity). */
135 #define _FP_PACK_SEMIRAW(fs, wc, X) \
138 if (_FP_FRAC_HIGH_##fs(X) \
139 & (_FP_OVERFLOW_##fs >> 1)) \
141 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1); \
143 if (X##_e == _FP_EXPMAX_##fs) \
144 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
146 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
147 if (!_FP_EXP_NORMAL(fs, wc, X) && !_FP_FRAC_ZEROP_##wc(X)) \
150 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
153 if (!_FP_KEEPNANFRACP) \
155 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
156 X##_s = _FP_NANSIGN_##fs; \
159 _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); \
234 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
235 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
240 /* underflow to zero */ \
242 if (!_FP_FRAC_ZEROP_##wc(X)) \
244 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
246 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
248 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
255 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
259 X##_e = _FP_EXPMAX_##fs; \
260 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
264 X##_e = _FP_EXPMAX_##fs; \
265 if (!_FP_KEEPNANFRACP) \
267 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
268 X##_s = _FP_NANSIGN_##fs; \
271 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
276 /* This one accepts raw argument and not cooked, returns
277 * 1 if X is a signaling NaN.
279 #define _FP_ISSIGNAN(fs, wc, X) \
282 if (X##_e == _FP_EXPMAX_##fs) \
284 if (!_FP_FRAC_ZEROP_##wc(X) \
285 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
295 /* Addition on semi-raw values. */
296 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
298 if (X##_s == Y##_s) \
302 int ediff = X##_e - Y##_e; \
308 /* Y is zero or denormalized. */ \
309 if (_FP_FRAC_ZEROP_##wc(Y)) \
311 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
312 _FP_FRAC_COPY_##wc(R, X); \
317 FP_SET_EXCEPTION(FP_EX_DENORM); \
321 _FP_FRAC_ADD_##wc(R, X, Y); \
324 if (X##_e == _FP_EXPMAX_##fs) \
326 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
327 _FP_FRAC_COPY_##wc(R, X); \
333 else if (X##_e == _FP_EXPMAX_##fs) \
335 /* X is NaN or Inf, Y is normal. */ \
336 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
337 _FP_FRAC_COPY_##wc(R, X); \
341 /* Insert implicit MSB of Y. */ \
342 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
345 /* Shift the mantissa of Y to the right EDIFF steps; \
346 remember to account later for the implicit MSB of X. */ \
347 if (ediff <= _FP_WFRACBITS_##fs) \
348 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
349 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
350 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
351 _FP_FRAC_ADD_##wc(R, X, Y); \
353 else if (ediff < 0) \
359 /* X is zero or denormalized. */ \
360 if (_FP_FRAC_ZEROP_##wc(X)) \
362 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
363 _FP_FRAC_COPY_##wc(R, Y); \
368 FP_SET_EXCEPTION(FP_EX_DENORM); \
372 _FP_FRAC_ADD_##wc(R, Y, X); \
375 if (Y##_e == _FP_EXPMAX_##fs) \
377 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
378 _FP_FRAC_COPY_##wc(R, Y); \
384 else if (Y##_e == _FP_EXPMAX_##fs) \
386 /* Y is NaN or Inf, X is normal. */ \
387 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
388 _FP_FRAC_COPY_##wc(R, Y); \
392 /* Insert implicit MSB of X. */ \
393 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
396 /* Shift the mantissa of X to the right EDIFF steps; \
397 remember to account later for the implicit MSB of Y. */ \
398 if (ediff <= _FP_WFRACBITS_##fs) \
399 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
400 else if (!_FP_FRAC_ZEROP_##wc(X)) \
401 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
402 _FP_FRAC_ADD_##wc(R, Y, X); \
407 if (!_FP_EXP_NORMAL(fs, wc, X)) \
411 /* X and Y are zero or denormalized. */ \
413 if (_FP_FRAC_ZEROP_##wc(X)) \
415 if (!_FP_FRAC_ZEROP_##wc(Y)) \
416 FP_SET_EXCEPTION(FP_EX_DENORM); \
417 _FP_FRAC_COPY_##wc(R, Y); \
420 else if (_FP_FRAC_ZEROP_##wc(Y)) \
422 FP_SET_EXCEPTION(FP_EX_DENORM); \
423 _FP_FRAC_COPY_##wc(R, X); \
428 FP_SET_EXCEPTION(FP_EX_DENORM); \
429 _FP_FRAC_ADD_##wc(R, X, Y); \
430 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
432 /* Normalized result. */ \
433 _FP_FRAC_HIGH_##fs(R) \
434 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
442 /* X and Y are NaN or Inf. */ \
443 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
444 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
445 R##_e = _FP_EXPMAX_##fs; \
446 if (_FP_FRAC_ZEROP_##wc(X)) \
447 _FP_FRAC_COPY_##wc(R, Y); \
448 else if (_FP_FRAC_ZEROP_##wc(Y)) \
449 _FP_FRAC_COPY_##wc(R, X); \
451 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
455 /* The exponents of X and Y, both normal, are equal. The \
456 implicit MSBs will always add to increase the \
458 _FP_FRAC_ADD_##wc(R, X, Y); \
460 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
461 if (R##_e == _FP_EXPMAX_##fs) \
462 /* Overflow to infinity (depending on rounding mode). */ \
463 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
467 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
470 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
472 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
473 if (R##_e == _FP_EXPMAX_##fs) \
474 /* Overflow to infinity (depending on rounding mode). */ \
475 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
482 int ediff = X##_e - Y##_e; \
489 /* Y is zero or denormalized. */ \
490 if (_FP_FRAC_ZEROP_##wc(Y)) \
492 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
493 _FP_FRAC_COPY_##wc(R, X); \
498 FP_SET_EXCEPTION(FP_EX_DENORM); \
502 _FP_FRAC_SUB_##wc(R, X, Y); \
505 if (X##_e == _FP_EXPMAX_##fs) \
507 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
508 _FP_FRAC_COPY_##wc(R, X); \
514 else if (X##_e == _FP_EXPMAX_##fs) \
516 /* X is NaN or Inf, Y is normal. */ \
517 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
518 _FP_FRAC_COPY_##wc(R, X); \
522 /* Insert implicit MSB of Y. */ \
523 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
526 /* Shift the mantissa of Y to the right EDIFF steps; \
527 remember to account later for the implicit MSB of X. */ \
528 if (ediff <= _FP_WFRACBITS_##fs) \
529 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
530 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
531 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
532 _FP_FRAC_SUB_##wc(R, X, Y); \
534 else if (ediff < 0) \
541 /* X is zero or denormalized. */ \
542 if (_FP_FRAC_ZEROP_##wc(X)) \
544 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
545 _FP_FRAC_COPY_##wc(R, Y); \
550 FP_SET_EXCEPTION(FP_EX_DENORM); \
554 _FP_FRAC_SUB_##wc(R, Y, X); \
557 if (Y##_e == _FP_EXPMAX_##fs) \
559 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
560 _FP_FRAC_COPY_##wc(R, Y); \
566 else if (Y##_e == _FP_EXPMAX_##fs) \
568 /* Y is NaN or Inf, X is normal. */ \
569 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
570 _FP_FRAC_COPY_##wc(R, Y); \
574 /* Insert implicit MSB of X. */ \
575 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
578 /* Shift the mantissa of X to the right EDIFF steps; \
579 remember to account later for the implicit MSB of Y. */ \
580 if (ediff <= _FP_WFRACBITS_##fs) \
581 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
582 else if (!_FP_FRAC_ZEROP_##wc(X)) \
583 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
584 _FP_FRAC_SUB_##wc(R, Y, X); \
589 if (!_FP_EXP_NORMAL(fs, wc, X)) \
593 /* X and Y are zero or denormalized. */ \
595 if (_FP_FRAC_ZEROP_##wc(X)) \
597 _FP_FRAC_COPY_##wc(R, Y); \
598 if (_FP_FRAC_ZEROP_##wc(Y)) \
599 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
602 FP_SET_EXCEPTION(FP_EX_DENORM); \
607 else if (_FP_FRAC_ZEROP_##wc(Y)) \
609 FP_SET_EXCEPTION(FP_EX_DENORM); \
610 _FP_FRAC_COPY_##wc(R, X); \
616 FP_SET_EXCEPTION(FP_EX_DENORM); \
617 _FP_FRAC_SUB_##wc(R, X, Y); \
619 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
621 /* |X| < |Y|, negate result. */ \
622 _FP_FRAC_SUB_##wc(R, Y, X); \
625 else if (_FP_FRAC_ZEROP_##wc(R)) \
626 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
632 /* X and Y are NaN or Inf, of opposite signs. */ \
633 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
634 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
635 R##_e = _FP_EXPMAX_##fs; \
636 if (_FP_FRAC_ZEROP_##wc(X)) \
638 if (_FP_FRAC_ZEROP_##wc(Y)) \
641 R##_s = _FP_NANSIGN_##fs; \
642 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
643 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
644 FP_SET_EXCEPTION(FP_EX_INVALID); \
650 _FP_FRAC_COPY_##wc(R, Y); \
655 if (_FP_FRAC_ZEROP_##wc(Y)) \
659 _FP_FRAC_COPY_##wc(R, X); \
664 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
670 /* The exponents of X and Y, both normal, are equal. The \
671 implicit MSBs cancel. */ \
673 _FP_FRAC_SUB_##wc(R, X, Y); \
675 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
677 /* |X| < |Y|, negate result. */ \
678 _FP_FRAC_SUB_##wc(R, Y, X); \
681 else if (_FP_FRAC_ZEROP_##wc(R)) \
684 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
690 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
693 /* Carry into most significant bit of larger one of X and Y, \
694 canceling it; renormalize. */ \
695 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \
697 _FP_FRAC_CLZ_##wc(diff, R); \
698 diff -= _FP_WFRACXBITS_##fs; \
699 _FP_FRAC_SLL_##wc(R, diff); \
702 /* R is denormalized. */ \
703 diff = diff - R##_e + 1; \
704 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \
710 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
717 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
718 #define _FP_SUB(fs, wc, R, X, Y) \
720 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
721 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \
726 * Main negation routine. FIXME -- when we care about setting exception
727 * bits reliably, this will not do. We should examine all of the fp classes.
730 #define _FP_NEG(fs, wc, R, X) \
732 _FP_FRAC_COPY_##wc(R, X); \
740 * Main multiplication routine. The input values should be cooked.
743 #define _FP_MUL(fs, wc, R, X, Y) \
745 R##_s = X##_s ^ Y##_s; \
746 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
748 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
749 R##_c = FP_CLS_NORMAL; \
750 R##_e = X##_e + Y##_e + 1; \
752 _FP_MUL_MEAT_##fs(R,X,Y); \
754 if (_FP_FRAC_OVERP_##wc(fs, R)) \
755 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
760 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
761 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \
764 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
765 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
766 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
769 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
770 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
771 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
772 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
773 _FP_FRAC_COPY_##wc(R, X); \
777 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
778 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
779 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
782 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
783 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
784 _FP_FRAC_COPY_##wc(R, Y); \
788 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
789 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
790 R##_s = _FP_NANSIGN_##fs; \
791 R##_c = FP_CLS_NAN; \
792 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
793 FP_SET_EXCEPTION(FP_EX_INVALID); \
803 * Main division routine. The input values should be cooked.
806 #define _FP_DIV(fs, wc, R, X, Y) \
808 R##_s = X##_s ^ Y##_s; \
809 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
811 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
812 R##_c = FP_CLS_NORMAL; \
813 R##_e = X##_e - Y##_e; \
815 _FP_DIV_MEAT_##fs(R,X,Y); \
818 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
819 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \
822 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
823 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
824 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
826 _FP_FRAC_COPY_##wc(R, X); \
830 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
831 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
832 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
834 _FP_FRAC_COPY_##wc(R, Y); \
838 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
839 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
840 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
841 R##_c = FP_CLS_ZERO; \
844 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
845 FP_SET_EXCEPTION(FP_EX_DIVZERO); \
846 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
847 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
848 R##_c = FP_CLS_INF; \
851 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
852 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
853 R##_s = _FP_NANSIGN_##fs; \
854 R##_c = FP_CLS_NAN; \
855 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
856 FP_SET_EXCEPTION(FP_EX_INVALID); \
866 * Main differential comparison routine. The inputs should be raw not
867 * cooked. The return is -1,0,1 for normal values, 2 otherwise.
870 #define _FP_CMP(fs, wc, ret, X, Y, un) \
872 /* NANs are unordered */ \
873 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
874 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
883 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \
884 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \
886 if (__is_zero_x && __is_zero_y) \
888 else if (__is_zero_x) \
889 ret = Y##_s ? 1 : -1; \
890 else if (__is_zero_y) \
891 ret = X##_s ? -1 : 1; \
892 else if (X##_s != Y##_s) \
893 ret = X##_s ? -1 : 1; \
894 else if (X##_e > Y##_e) \
895 ret = X##_s ? -1 : 1; \
896 else if (X##_e < Y##_e) \
897 ret = X##_s ? 1 : -1; \
898 else if (_FP_FRAC_GT_##wc(X, Y)) \
899 ret = X##_s ? -1 : 1; \
900 else if (_FP_FRAC_GT_##wc(Y, X)) \
901 ret = X##_s ? 1 : -1; \
908 /* Simplification for strict equality. */
910 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
912 /* NANs are unordered */ \
913 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
914 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
920 ret = !(X##_e == Y##_e \
921 && _FP_FRAC_EQ_##wc(X, Y) \
922 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
926 /* Version to test unordered. */
928 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
930 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
931 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \
935 * Main square root routine. The input value should be cooked.
938 #define _FP_SQRT(fs, wc, R, X) \
940 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \
945 _FP_FRAC_COPY_##wc(R, X); \
947 R##_c = FP_CLS_NAN; \
952 R##_s = _FP_NANSIGN_##fs; \
953 R##_c = FP_CLS_NAN; /* NAN */ \
954 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
955 FP_SET_EXCEPTION(FP_EX_INVALID); \
960 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
965 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
967 case FP_CLS_NORMAL: \
971 R##_c = FP_CLS_NAN; /* sNAN */ \
972 R##_s = _FP_NANSIGN_##fs; \
973 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
974 FP_SET_EXCEPTION(FP_EX_INVALID); \
977 R##_c = FP_CLS_NORMAL; \
979 _FP_FRAC_SLL_##wc(X, 1); \
980 R##_e = X##_e >> 1; \
981 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \
982 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \
983 q = _FP_OVERFLOW_##fs >> 1; \
984 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \
989 * Convert from FP to integer. Input is raw.
992 /* RSIGNED can have following values:
993 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
994 * the result is either 0 or (2^rsize)-1 depending on the sign in such
996 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
997 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
998 * depending on the sign in such case.
999 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1000 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1001 * depending on the sign in such case.
1003 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1005 if (X##_e < _FP_EXPBIAS_##fs) \
1010 if (!_FP_FRAC_ZEROP_##wc(X)) \
1012 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1013 FP_SET_EXCEPTION(FP_EX_DENORM); \
1017 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1019 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1020 || (!rsigned && X##_s)) \
1022 /* Overflow or converting to the most negative integer. */ \
1034 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1036 /* Possibly converting to most negative integer; check the \
1039 (void)((_FP_FRACBITS_##fs > rsize) \
1040 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \
1041 _FP_FRACBITS_##fs - rsize, \
1042 _FP_FRACBITS_##fs); 0; }) \
1044 if (!_FP_FRAC_ZEROP_##wc(X)) \
1045 FP_SET_EXCEPTION(FP_EX_INVALID); \
1047 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1050 FP_SET_EXCEPTION(FP_EX_INVALID); \
1054 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
1055 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1057 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1058 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1063 _FP_FRAC_SRST_##wc(X, inexact, \
1064 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1066 _FP_FRACBITS_##fs); \
1068 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1069 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1071 if (rsigned && X##_s) \
1076 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1078 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1084 if ((X##_s = (r < 0))) \
1088 (void)((rsize <= _FP_W_TYPE_SIZE) \
1091 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \
1092 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1094 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1097 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
1099 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1104 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1105 && X##_e >= _FP_EXPMAX_##fs) \
1107 /* Exponent too big; overflow to infinity. (May also \
1108 happen after rounding below.) */ \
1109 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
1110 goto pack_semiraw; \
1113 if (rsize <= _FP_FRACBITS_##fs \
1114 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1116 /* Exactly representable; shift left. */ \
1117 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1118 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1119 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1123 /* More bits in integer than in floating type; need to \
1125 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1126 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \
1127 - _FP_WFRACBITS_##fs + 1)) \
1128 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1129 - _FP_WFRACBITS_##fs + 1))) \
1131 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1132 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1133 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1134 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1135 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
1137 _FP_PACK_SEMIRAW(fs, wc, X); \
1144 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
1149 /* Extend from a narrower floating-point format to a wider one. Input
1150 and output are raw. */
1151 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \
1153 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1154 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1155 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1156 || _FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1159 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1160 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1162 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1163 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1169 if (_FP_FRAC_ZEROP_##swc(S)) \
1174 FP_SET_EXCEPTION(FP_EX_DENORM); \
1175 _FP_FRAC_CLZ_##swc(_lz, S); \
1176 _FP_FRAC_SLL_##dwc(D, \
1177 _lz + _FP_FRACBITS_##dfs \
1178 - _FP_FRACTBITS_##sfs); \
1179 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1180 + _FP_FRACXBITS_##sfs - _lz); \
1185 D##_e = _FP_EXPMAX_##dfs; \
1186 if (!_FP_FRAC_ZEROP_##swc(S)) \
1188 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1189 FP_SET_EXCEPTION(FP_EX_INVALID); \
1190 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1191 - _FP_FRACBITS_##sfs)); \
1197 /* Truncate from a wider floating-point format to a narrower one.
1198 Input and output are semi-raw. */
1199 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1201 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1202 || _FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1205 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1207 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1208 if (D##_e >= _FP_EXPMAX_##dfs) \
1209 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1214 if (D##_e <= 1 - _FP_FRACBITS_##dfs) \
1215 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1218 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1219 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1220 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1221 _FP_WFRACBITS_##sfs); \
1226 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1227 - _FP_WFRACBITS_##dfs), \
1228 _FP_WFRACBITS_##sfs); \
1229 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1237 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1238 if (!_FP_FRAC_ZEROP_##swc(S)) \
1240 FP_SET_EXCEPTION(FP_EX_DENORM); \
1241 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1246 D##_e = _FP_EXPMAX_##dfs; \
1247 if (_FP_FRAC_ZEROP_##swc(S)) \
1248 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1251 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1252 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1253 - _FP_WFRACBITS_##dfs)); \
1254 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1255 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1262 * Helper primitives.
1265 /* Count leading zeros in a word. */
1268 /* GCC 3.4 and later provide the builtins for us. */
1269 #define __FP_CLZ(r, x) \
1271 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1272 r = __builtin_clz (x); \
1273 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1274 r = __builtin_clzl (x); \
1275 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1276 r = __builtin_clzll (x); \
1280 #endif /* ndef __FP_CLZ */
1282 #define _FP_DIV_HELP_imm(q, r, n, d) \
1284 q = n / d, r = n % d; \
1288 /* A restoring bit-by-bit division primitive. */
1290 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1292 int count = _FP_WFRACBITS_##fs; \
1293 _FP_FRAC_DECL_##wc (u); \
1294 _FP_FRAC_DECL_##wc (v); \
1295 _FP_FRAC_COPY_##wc (u, X); \
1296 _FP_FRAC_COPY_##wc (v, Y); \
1297 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1298 /* Normalize U and V. */ \
1299 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1300 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1301 /* First round. Since the operands are normalized, either the \
1302 first or second bit will be set in the fraction. Produce a \
1303 normalized result by checking which and adjusting the loop \
1304 count and exponent accordingly. */ \
1305 if (_FP_FRAC_GE_1 (u, v)) \
1307 _FP_FRAC_SUB_##wc (u, u, v); \
1308 _FP_FRAC_LOW_##wc (R) |= 1; \
1313 /* Subsequent rounds. */ \
1315 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1316 _FP_FRAC_SLL_##wc (u, 1); \
1317 _FP_FRAC_SLL_##wc (R, 1); \
1318 if (msb || _FP_FRAC_GE_1 (u, v)) \
1320 _FP_FRAC_SUB_##wc (u, u, v); \
1321 _FP_FRAC_LOW_##wc (R) |= 1; \
1323 } while (--count > 0); \
1324 /* If there's anything left in U, the result is inexact. */ \
1325 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1328 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1329 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1330 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)