1 /* Software floating-point emulation. Common operations.
2 Copyright (C) 1997,1998,1999,2006,2007 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_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
104 FP_SET_EXCEPTION(FP_EX_INEXACT); \
105 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
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 \
1157 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
1160 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1161 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1163 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1164 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1170 if (_FP_FRAC_ZEROP_##swc(S)) \
1172 else if (_FP_EXPBIAS_##dfs \
1173 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1175 FP_SET_EXCEPTION(FP_EX_DENORM); \
1176 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1177 - _FP_FRACBITS_##sfs)); \
1183 FP_SET_EXCEPTION(FP_EX_DENORM); \
1184 _FP_FRAC_CLZ_##swc(_lz, S); \
1185 _FP_FRAC_SLL_##dwc(D, \
1186 _lz + _FP_FRACBITS_##dfs \
1187 - _FP_FRACTBITS_##sfs); \
1188 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1189 + _FP_FRACXBITS_##sfs - _lz); \
1194 D##_e = _FP_EXPMAX_##dfs; \
1195 if (!_FP_FRAC_ZEROP_##swc(S)) \
1197 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1198 FP_SET_EXCEPTION(FP_EX_INVALID); \
1199 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1200 - _FP_FRACBITS_##sfs)); \
1206 /* Truncate from a wider floating-point format to a narrower one.
1207 Input and output are semi-raw. */
1208 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1210 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1211 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1212 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1215 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1217 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1218 if (D##_e >= _FP_EXPMAX_##dfs) \
1219 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1224 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1226 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1227 _FP_FRAC_LOW_##swc(S) |= 1; \
1231 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1232 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1233 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1234 _FP_WFRACBITS_##sfs); \
1239 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1240 - _FP_WFRACBITS_##dfs), \
1241 _FP_WFRACBITS_##sfs); \
1242 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1250 if (_FP_FRAC_ZEROP_##swc(S)) \
1251 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1254 FP_SET_EXCEPTION(FP_EX_DENORM); \
1255 if (_FP_EXPBIAS_##sfs \
1256 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1258 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1259 - _FP_WFRACBITS_##dfs), \
1260 _FP_WFRACBITS_##sfs); \
1261 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1265 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1266 _FP_FRAC_LOW_##dwc(D) |= 1; \
1272 D##_e = _FP_EXPMAX_##dfs; \
1273 if (_FP_FRAC_ZEROP_##swc(S)) \
1274 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1277 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1278 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1279 - _FP_WFRACBITS_##dfs)); \
1280 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1281 /* Semi-raw NaN must have all workbits cleared. */ \
1282 _FP_FRAC_LOW_##dwc(D) \
1283 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1284 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1291 * Helper primitives.
1294 /* Count leading zeros in a word. */
1297 /* GCC 3.4 and later provide the builtins for us. */
1298 #define __FP_CLZ(r, x) \
1300 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1301 r = __builtin_clz (x); \
1302 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1303 r = __builtin_clzl (x); \
1304 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1305 r = __builtin_clzll (x); \
1309 #endif /* ndef __FP_CLZ */
1311 #define _FP_DIV_HELP_imm(q, r, n, d) \
1313 q = n / d, r = n % d; \
1317 /* A restoring bit-by-bit division primitive. */
1319 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1321 int count = _FP_WFRACBITS_##fs; \
1322 _FP_FRAC_DECL_##wc (u); \
1323 _FP_FRAC_DECL_##wc (v); \
1324 _FP_FRAC_COPY_##wc (u, X); \
1325 _FP_FRAC_COPY_##wc (v, Y); \
1326 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1327 /* Normalize U and V. */ \
1328 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1329 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1330 /* First round. Since the operands are normalized, either the \
1331 first or second bit will be set in the fraction. Produce a \
1332 normalized result by checking which and adjusting the loop \
1333 count and exponent accordingly. */ \
1334 if (_FP_FRAC_GE_1 (u, v)) \
1336 _FP_FRAC_SUB_##wc (u, u, v); \
1337 _FP_FRAC_LOW_##wc (R) |= 1; \
1342 /* Subsequent rounds. */ \
1344 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1345 _FP_FRAC_SLL_##wc (u, 1); \
1346 _FP_FRAC_SLL_##wc (R, 1); \
1347 if (msb || _FP_FRAC_GE_1 (u, v)) \
1349 _FP_FRAC_SUB_##wc (u, u, v); \
1350 _FP_FRAC_LOW_##wc (R) |= 1; \
1352 } while (--count > 0); \
1353 /* If there's anything left in U, the result is inexact. */ \
1354 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1357 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1358 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1359 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)