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, 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 (_FP_FRAC_HIGH_##fs(X) \
138 & (_FP_OVERFLOW_##fs >> 1)) \
140 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1); \
142 if (X##_e == _FP_EXPMAX_##fs) \
143 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
145 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
146 if (!_FP_EXP_NORMAL(fs, wc, X) && !_FP_FRAC_ZEROP_##wc(X)) \
149 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
152 if (!_FP_KEEPNANFRACP) \
154 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
155 X##_s = _FP_NANSIGN_##fs; \
158 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
164 * Before packing the bits back into the native fp result, take care
165 * of such mundane things as rounding and overflow. Also, for some
166 * kinds of fp values, the original parts may not have been fully
167 * extracted -- but that is ok, we can regenerate them now.
170 #define _FP_PACK_CANONICAL(fs, wc, X) \
174 case FP_CLS_NORMAL: \
175 X##_e += _FP_EXPBIAS_##fs; \
179 if (_FP_FRAC_OVERP_##wc(fs, X)) \
181 _FP_FRAC_CLEAR_OVERP_##wc(fs, X); \
184 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
185 if (X##_e >= _FP_EXPMAX_##fs) \
188 switch (FP_ROUNDMODE) \
190 case FP_RND_NEAREST: \
191 X##_c = FP_CLS_INF; \
194 if (!X##_s) X##_c = FP_CLS_INF; \
197 if (X##_s) X##_c = FP_CLS_INF; \
200 if (X##_c == FP_CLS_INF) \
202 /* Overflow to infinity */ \
203 X##_e = _FP_EXPMAX_##fs; \
204 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
208 /* Overflow to maximum normal */ \
209 X##_e = _FP_EXPMAX_##fs - 1; \
210 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
212 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
213 FP_SET_EXCEPTION(FP_EX_INEXACT); \
218 /* we've got a denormalized number */ \
219 X##_e = -X##_e + 1; \
220 if (X##_e <= _FP_WFRACBITS_##fs) \
222 _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs); \
224 if (_FP_FRAC_HIGH_##fs(X) \
225 & (_FP_OVERFLOW_##fs >> 1)) \
228 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
233 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
234 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
239 /* underflow to zero */ \
241 if (!_FP_FRAC_ZEROP_##wc(X)) \
243 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
245 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
247 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
254 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
258 X##_e = _FP_EXPMAX_##fs; \
259 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
263 X##_e = _FP_EXPMAX_##fs; \
264 if (!_FP_KEEPNANFRACP) \
266 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
267 X##_s = _FP_NANSIGN_##fs; \
270 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
275 /* This one accepts raw argument and not cooked, returns
276 * 1 if X is a signaling NaN.
278 #define _FP_ISSIGNAN(fs, wc, X) \
281 if (X##_e == _FP_EXPMAX_##fs) \
283 if (!_FP_FRAC_ZEROP_##wc(X) \
284 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
294 /* Addition on semi-raw values. */
295 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
297 if (X##_s == Y##_s) \
301 int ediff = X##_e - Y##_e; \
307 /* Y is zero or denormalized. */ \
308 if (_FP_FRAC_ZEROP_##wc(Y)) \
310 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
311 _FP_FRAC_COPY_##wc(R, X); \
316 FP_SET_EXCEPTION(FP_EX_DENORM); \
320 _FP_FRAC_ADD_##wc(R, X, Y); \
323 if (X##_e == _FP_EXPMAX_##fs) \
325 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
326 _FP_FRAC_COPY_##wc(R, X); \
332 else if (X##_e == _FP_EXPMAX_##fs) \
334 /* X is NaN or Inf, Y is normal. */ \
335 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
336 _FP_FRAC_COPY_##wc(R, X); \
340 /* Insert implicit MSB of Y. */ \
341 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
344 /* Shift the mantissa of Y to the right EDIFF steps; \
345 remember to account later for the implicit MSB of X. */ \
346 if (ediff <= _FP_WFRACBITS_##fs) \
347 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
348 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
349 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
350 _FP_FRAC_ADD_##wc(R, X, Y); \
352 else if (ediff < 0) \
358 /* X is zero or denormalized. */ \
359 if (_FP_FRAC_ZEROP_##wc(X)) \
361 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
362 _FP_FRAC_COPY_##wc(R, Y); \
367 FP_SET_EXCEPTION(FP_EX_DENORM); \
371 _FP_FRAC_ADD_##wc(R, Y, X); \
374 if (Y##_e == _FP_EXPMAX_##fs) \
376 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
377 _FP_FRAC_COPY_##wc(R, Y); \
383 else if (Y##_e == _FP_EXPMAX_##fs) \
385 /* Y is NaN or Inf, X is normal. */ \
386 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
387 _FP_FRAC_COPY_##wc(R, Y); \
391 /* Insert implicit MSB of X. */ \
392 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
395 /* Shift the mantissa of X to the right EDIFF steps; \
396 remember to account later for the implicit MSB of Y. */ \
397 if (ediff <= _FP_WFRACBITS_##fs) \
398 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
399 else if (!_FP_FRAC_ZEROP_##wc(X)) \
400 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
401 _FP_FRAC_ADD_##wc(R, Y, X); \
406 if (!_FP_EXP_NORMAL(fs, wc, X)) \
410 /* X and Y are zero or denormalized. */ \
412 if (_FP_FRAC_ZEROP_##wc(X)) \
414 if (!_FP_FRAC_ZEROP_##wc(Y)) \
415 FP_SET_EXCEPTION(FP_EX_DENORM); \
416 _FP_FRAC_COPY_##wc(R, Y); \
419 else if (_FP_FRAC_ZEROP_##wc(Y)) \
421 FP_SET_EXCEPTION(FP_EX_DENORM); \
422 _FP_FRAC_COPY_##wc(R, X); \
427 FP_SET_EXCEPTION(FP_EX_DENORM); \
428 _FP_FRAC_ADD_##wc(R, X, Y); \
429 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
431 /* Normalized result. */ \
432 _FP_FRAC_HIGH_##fs(R) \
433 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
441 /* X and Y are NaN or Inf. */ \
442 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
443 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
444 R##_e = _FP_EXPMAX_##fs; \
445 if (_FP_FRAC_ZEROP_##wc(X)) \
446 _FP_FRAC_COPY_##wc(R, Y); \
447 else if (_FP_FRAC_ZEROP_##wc(Y)) \
448 _FP_FRAC_COPY_##wc(R, X); \
450 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
454 /* The exponents of X and Y, both normal, are equal. The \
455 implicit MSBs will always add to increase the \
457 _FP_FRAC_ADD_##wc(R, X, Y); \
459 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
460 if (R##_e == _FP_EXPMAX_##fs) \
461 /* Overflow to infinity (depending on rounding mode). */ \
462 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
466 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
469 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
471 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
472 if (R##_e == _FP_EXPMAX_##fs) \
473 /* Overflow to infinity (depending on rounding mode). */ \
474 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
481 int ediff = X##_e - Y##_e; \
488 /* Y is zero or denormalized. */ \
489 if (_FP_FRAC_ZEROP_##wc(Y)) \
491 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
492 _FP_FRAC_COPY_##wc(R, X); \
497 FP_SET_EXCEPTION(FP_EX_DENORM); \
501 _FP_FRAC_SUB_##wc(R, X, Y); \
504 if (X##_e == _FP_EXPMAX_##fs) \
506 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
507 _FP_FRAC_COPY_##wc(R, X); \
513 else if (X##_e == _FP_EXPMAX_##fs) \
515 /* X is NaN or Inf, Y is normal. */ \
516 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
517 _FP_FRAC_COPY_##wc(R, X); \
521 /* Insert implicit MSB of Y. */ \
522 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
525 /* Shift the mantissa of Y to the right EDIFF steps; \
526 remember to account later for the implicit MSB of X. */ \
527 if (ediff <= _FP_WFRACBITS_##fs) \
528 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
529 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
530 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
531 _FP_FRAC_SUB_##wc(R, X, Y); \
533 else if (ediff < 0) \
540 /* X is zero or denormalized. */ \
541 if (_FP_FRAC_ZEROP_##wc(X)) \
543 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
544 _FP_FRAC_COPY_##wc(R, Y); \
549 FP_SET_EXCEPTION(FP_EX_DENORM); \
553 _FP_FRAC_SUB_##wc(R, Y, X); \
556 if (Y##_e == _FP_EXPMAX_##fs) \
558 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
559 _FP_FRAC_COPY_##wc(R, Y); \
565 else if (Y##_e == _FP_EXPMAX_##fs) \
567 /* Y is NaN or Inf, X is normal. */ \
568 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
569 _FP_FRAC_COPY_##wc(R, Y); \
573 /* Insert implicit MSB of X. */ \
574 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
577 /* Shift the mantissa of X to the right EDIFF steps; \
578 remember to account later for the implicit MSB of Y. */ \
579 if (ediff <= _FP_WFRACBITS_##fs) \
580 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
581 else if (!_FP_FRAC_ZEROP_##wc(X)) \
582 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
583 _FP_FRAC_SUB_##wc(R, Y, X); \
588 if (!_FP_EXP_NORMAL(fs, wc, X)) \
592 /* X and Y are zero or denormalized. */ \
594 if (_FP_FRAC_ZEROP_##wc(X)) \
596 _FP_FRAC_COPY_##wc(R, Y); \
597 if (_FP_FRAC_ZEROP_##wc(Y)) \
598 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
601 FP_SET_EXCEPTION(FP_EX_DENORM); \
606 else if (_FP_FRAC_ZEROP_##wc(Y)) \
608 FP_SET_EXCEPTION(FP_EX_DENORM); \
609 _FP_FRAC_COPY_##wc(R, X); \
615 FP_SET_EXCEPTION(FP_EX_DENORM); \
616 _FP_FRAC_SUB_##wc(R, X, Y); \
618 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
620 /* |X| < |Y|, negate result. */ \
621 _FP_FRAC_SUB_##wc(R, Y, X); \
624 else if (_FP_FRAC_ZEROP_##wc(R)) \
625 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
631 /* X and Y are NaN or Inf, of opposite signs. */ \
632 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
633 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
634 R##_e = _FP_EXPMAX_##fs; \
635 if (_FP_FRAC_ZEROP_##wc(X)) \
637 if (_FP_FRAC_ZEROP_##wc(Y)) \
640 R##_s = _FP_NANSIGN_##fs; \
641 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
642 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
643 FP_SET_EXCEPTION(FP_EX_INVALID); \
649 _FP_FRAC_COPY_##wc(R, Y); \
654 if (_FP_FRAC_ZEROP_##wc(Y)) \
658 _FP_FRAC_COPY_##wc(R, X); \
663 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
669 /* The exponents of X and Y, both normal, are equal. The \
670 implicit MSBs cancel. */ \
672 _FP_FRAC_SUB_##wc(R, X, Y); \
674 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
676 /* |X| < |Y|, negate result. */ \
677 _FP_FRAC_SUB_##wc(R, Y, X); \
680 else if (_FP_FRAC_ZEROP_##wc(R)) \
683 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
689 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
692 /* Carry into most significant bit of larger one of X and Y, \
693 canceling it; renormalize. */ \
694 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \
696 _FP_FRAC_CLZ_##wc(diff, R); \
697 diff -= _FP_WFRACXBITS_##fs; \
698 _FP_FRAC_SLL_##wc(R, diff); \
701 /* R is denormalized. */ \
702 diff = diff - R##_e + 1; \
703 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \
709 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
716 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
717 #define _FP_SUB(fs, wc, R, X, Y) \
719 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
720 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \
725 * Main negation routine. FIXME -- when we care about setting exception
726 * bits reliably, this will not do. We should examine all of the fp classes.
729 #define _FP_NEG(fs, wc, R, X) \
731 _FP_FRAC_COPY_##wc(R, X); \
739 * Main multiplication routine. The input values should be cooked.
742 #define _FP_MUL(fs, wc, R, X, Y) \
744 R##_s = X##_s ^ Y##_s; \
745 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
747 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
748 R##_c = FP_CLS_NORMAL; \
749 R##_e = X##_e + Y##_e + 1; \
751 _FP_MUL_MEAT_##fs(R,X,Y); \
753 if (_FP_FRAC_OVERP_##wc(fs, R)) \
754 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
759 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
760 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \
763 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
764 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
765 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
768 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
769 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
770 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
771 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
772 _FP_FRAC_COPY_##wc(R, X); \
776 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
777 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
778 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
781 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
782 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
783 _FP_FRAC_COPY_##wc(R, Y); \
787 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
788 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
789 R##_s = _FP_NANSIGN_##fs; \
790 R##_c = FP_CLS_NAN; \
791 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
792 FP_SET_EXCEPTION(FP_EX_INVALID); \
802 * Main division routine. The input values should be cooked.
805 #define _FP_DIV(fs, wc, R, X, Y) \
807 R##_s = X##_s ^ Y##_s; \
808 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
810 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
811 R##_c = FP_CLS_NORMAL; \
812 R##_e = X##_e - Y##_e; \
814 _FP_DIV_MEAT_##fs(R,X,Y); \
817 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
818 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \
821 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
822 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
823 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
825 _FP_FRAC_COPY_##wc(R, X); \
829 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
830 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
831 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
833 _FP_FRAC_COPY_##wc(R, Y); \
837 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
838 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
839 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
840 R##_c = FP_CLS_ZERO; \
843 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
844 FP_SET_EXCEPTION(FP_EX_DIVZERO); \
845 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
846 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
847 R##_c = FP_CLS_INF; \
850 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
851 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
852 R##_s = _FP_NANSIGN_##fs; \
853 R##_c = FP_CLS_NAN; \
854 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
855 FP_SET_EXCEPTION(FP_EX_INVALID); \
865 * Main differential comparison routine. The inputs should be raw not
866 * cooked. The return is -1,0,1 for normal values, 2 otherwise.
869 #define _FP_CMP(fs, wc, ret, X, Y, un) \
871 /* NANs are unordered */ \
872 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
873 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
882 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \
883 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \
885 if (__is_zero_x && __is_zero_y) \
887 else if (__is_zero_x) \
888 ret = Y##_s ? 1 : -1; \
889 else if (__is_zero_y) \
890 ret = X##_s ? -1 : 1; \
891 else if (X##_s != Y##_s) \
892 ret = X##_s ? -1 : 1; \
893 else if (X##_e > Y##_e) \
894 ret = X##_s ? -1 : 1; \
895 else if (X##_e < Y##_e) \
896 ret = X##_s ? 1 : -1; \
897 else if (_FP_FRAC_GT_##wc(X, Y)) \
898 ret = X##_s ? -1 : 1; \
899 else if (_FP_FRAC_GT_##wc(Y, X)) \
900 ret = X##_s ? 1 : -1; \
907 /* Simplification for strict equality. */
909 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
911 /* NANs are unordered */ \
912 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
913 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
919 ret = !(X##_e == Y##_e \
920 && _FP_FRAC_EQ_##wc(X, Y) \
921 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
925 /* Version to test unordered. */
927 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
929 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
930 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \
934 * Main square root routine. The input value should be cooked.
937 #define _FP_SQRT(fs, wc, R, X) \
939 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \
944 _FP_FRAC_COPY_##wc(R, X); \
946 R##_c = FP_CLS_NAN; \
951 R##_s = _FP_NANSIGN_##fs; \
952 R##_c = FP_CLS_NAN; /* NAN */ \
953 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
954 FP_SET_EXCEPTION(FP_EX_INVALID); \
959 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
964 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
966 case FP_CLS_NORMAL: \
970 R##_c = FP_CLS_NAN; /* sNAN */ \
971 R##_s = _FP_NANSIGN_##fs; \
972 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
973 FP_SET_EXCEPTION(FP_EX_INVALID); \
976 R##_c = FP_CLS_NORMAL; \
978 _FP_FRAC_SLL_##wc(X, 1); \
979 R##_e = X##_e >> 1; \
980 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \
981 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \
982 q = _FP_OVERFLOW_##fs >> 1; \
983 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \
988 * Convert from FP to integer. Input is raw.
991 /* RSIGNED can have following values:
992 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
993 * the result is either 0 or (2^rsize)-1 depending on the sign in such
995 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
996 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
997 * depending on the sign in such case.
998 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
999 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1000 * depending on the sign in such case.
1002 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1004 if (X##_e < _FP_EXPBIAS_##fs) \
1009 if (!_FP_FRAC_ZEROP_##wc(X)) \
1011 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1012 FP_SET_EXCEPTION(FP_EX_DENORM); \
1016 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1018 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1019 || (!rsigned && X##_s)) \
1021 /* Overflow or converting to the most negative integer. */ \
1033 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1035 /* Possibly converting to most negative integer; check the \
1038 (void)((_FP_FRACBITS_##fs > rsize) \
1039 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \
1040 _FP_FRACBITS_##fs - rsize, \
1041 _FP_FRACBITS_##fs); 0; }) \
1043 if (!_FP_FRAC_ZEROP_##wc(X)) \
1044 FP_SET_EXCEPTION(FP_EX_INVALID); \
1046 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1049 FP_SET_EXCEPTION(FP_EX_INVALID); \
1053 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
1054 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1056 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1057 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1062 _FP_FRAC_SRST_##wc(X, inexact, \
1063 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1065 _FP_FRACBITS_##fs); \
1067 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1068 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1070 if (rsigned && X##_s) \
1075 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1077 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1083 if ((X##_s = (r < 0))) \
1087 (void)((rsize <= _FP_W_TYPE_SIZE) \
1090 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \
1091 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1093 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1096 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
1098 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1103 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1104 && X##_e >= _FP_EXPMAX_##fs) \
1106 /* Exponent too big; overflow to infinity. (May also \
1107 happen after rounding below.) */ \
1108 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
1109 goto pack_semiraw; \
1112 if (rsize <= _FP_FRACBITS_##fs \
1113 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1115 /* Exactly representable; shift left. */ \
1116 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1117 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1118 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1122 /* More bits in integer than in floating type; need to \
1124 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1125 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \
1126 - _FP_WFRACBITS_##fs + 1)) \
1127 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1128 - _FP_WFRACBITS_##fs + 1))) \
1130 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1131 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1132 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1133 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1134 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
1136 _FP_PACK_SEMIRAW(fs, wc, X); \
1143 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
1148 /* Extend from a narrower floating-point format to a wider one. Input
1149 and output are raw. */
1150 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \
1152 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1153 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1154 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1155 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1156 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
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)) \
1171 else if (_FP_EXPBIAS_##dfs \
1172 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1174 FP_SET_EXCEPTION(FP_EX_DENORM); \
1175 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1176 - _FP_FRACBITS_##sfs)); \
1182 FP_SET_EXCEPTION(FP_EX_DENORM); \
1183 _FP_FRAC_CLZ_##swc(_lz, S); \
1184 _FP_FRAC_SLL_##dwc(D, \
1185 _lz + _FP_FRACBITS_##dfs \
1186 - _FP_FRACTBITS_##sfs); \
1187 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1188 + _FP_FRACXBITS_##sfs - _lz); \
1193 D##_e = _FP_EXPMAX_##dfs; \
1194 if (!_FP_FRAC_ZEROP_##swc(S)) \
1196 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1197 FP_SET_EXCEPTION(FP_EX_INVALID); \
1198 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1199 - _FP_FRACBITS_##sfs)); \
1205 /* Truncate from a wider floating-point format to a narrower one.
1206 Input and output are semi-raw. */
1207 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1209 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1210 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1211 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1214 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1216 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1217 if (D##_e >= _FP_EXPMAX_##dfs) \
1218 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1223 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1225 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1226 _FP_FRAC_LOW_##swc(S) |= 1; \
1230 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1231 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1232 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1233 _FP_WFRACBITS_##sfs); \
1238 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1239 - _FP_WFRACBITS_##dfs), \
1240 _FP_WFRACBITS_##sfs); \
1241 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1249 if (_FP_FRAC_ZEROP_##swc(S)) \
1250 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1253 FP_SET_EXCEPTION(FP_EX_DENORM); \
1254 if (_FP_EXPBIAS_##sfs \
1255 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1257 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1258 - _FP_WFRACBITS_##dfs), \
1259 _FP_WFRACBITS_##sfs); \
1260 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1264 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1265 _FP_FRAC_LOW_##dwc(D) |= 1; \
1271 D##_e = _FP_EXPMAX_##dfs; \
1272 if (_FP_FRAC_ZEROP_##swc(S)) \
1273 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1276 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1277 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1278 - _FP_WFRACBITS_##dfs)); \
1279 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1280 /* Semi-raw NaN must have all workbits cleared. */ \
1281 _FP_FRAC_LOW_##dwc(D) \
1282 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1283 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1290 * Helper primitives.
1293 /* Count leading zeros in a word. */
1296 /* GCC 3.4 and later provide the builtins for us. */
1297 #define __FP_CLZ(r, x) \
1299 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1300 r = __builtin_clz (x); \
1301 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1302 r = __builtin_clzl (x); \
1303 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1304 r = __builtin_clzll (x); \
1308 #endif /* ndef __FP_CLZ */
1310 #define _FP_DIV_HELP_imm(q, r, n, d) \
1312 q = n / d, r = n % d; \
1316 /* A restoring bit-by-bit division primitive. */
1318 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1320 int count = _FP_WFRACBITS_##fs; \
1321 _FP_FRAC_DECL_##wc (u); \
1322 _FP_FRAC_DECL_##wc (v); \
1323 _FP_FRAC_COPY_##wc (u, X); \
1324 _FP_FRAC_COPY_##wc (v, Y); \
1325 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1326 /* Normalize U and V. */ \
1327 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1328 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1329 /* First round. Since the operands are normalized, either the \
1330 first or second bit will be set in the fraction. Produce a \
1331 normalized result by checking which and adjusting the loop \
1332 count and exponent accordingly. */ \
1333 if (_FP_FRAC_GE_1 (u, v)) \
1335 _FP_FRAC_SUB_##wc (u, u, v); \
1336 _FP_FRAC_LOW_##wc (R) |= 1; \
1341 /* Subsequent rounds. */ \
1343 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1344 _FP_FRAC_SLL_##wc (u, 1); \
1345 _FP_FRAC_SLL_##wc (R, 1); \
1346 if (msb || _FP_FRAC_GE_1 (u, v)) \
1348 _FP_FRAC_SUB_##wc (u, u, v); \
1349 _FP_FRAC_LOW_##wc (R) |= 1; \
1351 } while (--count > 0); \
1352 /* If there's anything left in U, the result is inexact. */ \
1353 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1356 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1357 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1358 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)