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 (_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)) \
150 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
151 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
152 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
156 if (!_FP_KEEPNANFRACP) \
158 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
159 X##_s = _FP_NANSIGN_##fs; \
162 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
168 * Before packing the bits back into the native fp result, take care
169 * of such mundane things as rounding and overflow. Also, for some
170 * kinds of fp values, the original parts may not have been fully
171 * extracted -- but that is ok, we can regenerate them now.
174 #define _FP_PACK_CANONICAL(fs, wc, X) \
178 case FP_CLS_NORMAL: \
179 X##_e += _FP_EXPBIAS_##fs; \
183 if (_FP_FRAC_OVERP_##wc(fs, X)) \
185 _FP_FRAC_CLEAR_OVERP_##wc(fs, X); \
188 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
189 if (X##_e >= _FP_EXPMAX_##fs) \
192 switch (FP_ROUNDMODE) \
194 case FP_RND_NEAREST: \
195 X##_c = FP_CLS_INF; \
198 if (!X##_s) X##_c = FP_CLS_INF; \
201 if (X##_s) X##_c = FP_CLS_INF; \
204 if (X##_c == FP_CLS_INF) \
206 /* Overflow to infinity */ \
207 X##_e = _FP_EXPMAX_##fs; \
208 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
212 /* Overflow to maximum normal */ \
213 X##_e = _FP_EXPMAX_##fs - 1; \
214 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \
216 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \
217 FP_SET_EXCEPTION(FP_EX_INEXACT); \
222 /* we've got a denormalized number */ \
223 X##_e = -X##_e + 1; \
224 if (X##_e <= _FP_WFRACBITS_##fs) \
226 _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs); \
228 if (_FP_FRAC_HIGH_##fs(X) \
229 & (_FP_OVERFLOW_##fs >> 1)) \
232 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
233 FP_SET_EXCEPTION(FP_EX_INEXACT); \
238 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \
240 if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) \
241 || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)) \
242 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
246 /* underflow to zero */ \
248 if (!_FP_FRAC_ZEROP_##wc(X)) \
250 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
252 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \
254 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \
261 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
265 X##_e = _FP_EXPMAX_##fs; \
266 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
270 X##_e = _FP_EXPMAX_##fs; \
271 if (!_FP_KEEPNANFRACP) \
273 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \
274 X##_s = _FP_NANSIGN_##fs; \
277 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \
282 /* This one accepts raw argument and not cooked, returns
283 * 1 if X is a signaling NaN.
285 #define _FP_ISSIGNAN(fs, wc, X) \
288 if (X##_e == _FP_EXPMAX_##fs) \
290 if (!_FP_FRAC_ZEROP_##wc(X) \
291 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \
301 /* Addition on semi-raw values. */
302 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \
304 if (X##_s == Y##_s) \
308 int ediff = X##_e - Y##_e; \
314 /* Y is zero or denormalized. */ \
315 if (_FP_FRAC_ZEROP_##wc(Y)) \
317 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
318 _FP_FRAC_COPY_##wc(R, X); \
323 FP_SET_EXCEPTION(FP_EX_DENORM); \
327 _FP_FRAC_ADD_##wc(R, X, Y); \
330 if (X##_e == _FP_EXPMAX_##fs) \
332 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
333 _FP_FRAC_COPY_##wc(R, X); \
339 else if (X##_e == _FP_EXPMAX_##fs) \
341 /* X is NaN or Inf, Y is normal. */ \
342 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
343 _FP_FRAC_COPY_##wc(R, X); \
347 /* Insert implicit MSB of Y. */ \
348 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
351 /* Shift the mantissa of Y to the right EDIFF steps; \
352 remember to account later for the implicit MSB of X. */ \
353 if (ediff <= _FP_WFRACBITS_##fs) \
354 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
355 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
356 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
357 _FP_FRAC_ADD_##wc(R, X, Y); \
359 else if (ediff < 0) \
365 /* X is zero or denormalized. */ \
366 if (_FP_FRAC_ZEROP_##wc(X)) \
368 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
369 _FP_FRAC_COPY_##wc(R, Y); \
374 FP_SET_EXCEPTION(FP_EX_DENORM); \
378 _FP_FRAC_ADD_##wc(R, Y, X); \
381 if (Y##_e == _FP_EXPMAX_##fs) \
383 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
384 _FP_FRAC_COPY_##wc(R, Y); \
390 else if (Y##_e == _FP_EXPMAX_##fs) \
392 /* Y is NaN or Inf, X is normal. */ \
393 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
394 _FP_FRAC_COPY_##wc(R, Y); \
398 /* Insert implicit MSB of X. */ \
399 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
402 /* Shift the mantissa of X to the right EDIFF steps; \
403 remember to account later for the implicit MSB of Y. */ \
404 if (ediff <= _FP_WFRACBITS_##fs) \
405 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
406 else if (!_FP_FRAC_ZEROP_##wc(X)) \
407 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
408 _FP_FRAC_ADD_##wc(R, Y, X); \
413 if (!_FP_EXP_NORMAL(fs, wc, X)) \
417 /* X and Y are zero or denormalized. */ \
419 if (_FP_FRAC_ZEROP_##wc(X)) \
421 if (!_FP_FRAC_ZEROP_##wc(Y)) \
422 FP_SET_EXCEPTION(FP_EX_DENORM); \
423 _FP_FRAC_COPY_##wc(R, Y); \
426 else if (_FP_FRAC_ZEROP_##wc(Y)) \
428 FP_SET_EXCEPTION(FP_EX_DENORM); \
429 _FP_FRAC_COPY_##wc(R, X); \
434 FP_SET_EXCEPTION(FP_EX_DENORM); \
435 _FP_FRAC_ADD_##wc(R, X, Y); \
436 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
438 /* Normalized result. */ \
439 _FP_FRAC_HIGH_##fs(R) \
440 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
448 /* X and Y are NaN or Inf. */ \
449 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
450 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
451 R##_e = _FP_EXPMAX_##fs; \
452 if (_FP_FRAC_ZEROP_##wc(X)) \
453 _FP_FRAC_COPY_##wc(R, Y); \
454 else if (_FP_FRAC_ZEROP_##wc(Y)) \
455 _FP_FRAC_COPY_##wc(R, X); \
457 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
461 /* The exponents of X and Y, both normal, are equal. The \
462 implicit MSBs will always add to increase the \
464 _FP_FRAC_ADD_##wc(R, X, Y); \
466 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
467 if (R##_e == _FP_EXPMAX_##fs) \
468 /* Overflow to infinity (depending on rounding mode). */ \
469 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
473 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
476 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
478 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
479 if (R##_e == _FP_EXPMAX_##fs) \
480 /* Overflow to infinity (depending on rounding mode). */ \
481 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \
488 int ediff = X##_e - Y##_e; \
495 /* Y is zero or denormalized. */ \
496 if (_FP_FRAC_ZEROP_##wc(Y)) \
498 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
499 _FP_FRAC_COPY_##wc(R, X); \
504 FP_SET_EXCEPTION(FP_EX_DENORM); \
508 _FP_FRAC_SUB_##wc(R, X, Y); \
511 if (X##_e == _FP_EXPMAX_##fs) \
513 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
514 _FP_FRAC_COPY_##wc(R, X); \
520 else if (X##_e == _FP_EXPMAX_##fs) \
522 /* X is NaN or Inf, Y is normal. */ \
523 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
524 _FP_FRAC_COPY_##wc(R, X); \
528 /* Insert implicit MSB of Y. */ \
529 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \
532 /* Shift the mantissa of Y to the right EDIFF steps; \
533 remember to account later for the implicit MSB of X. */ \
534 if (ediff <= _FP_WFRACBITS_##fs) \
535 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \
536 else if (!_FP_FRAC_ZEROP_##wc(Y)) \
537 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \
538 _FP_FRAC_SUB_##wc(R, X, Y); \
540 else if (ediff < 0) \
547 /* X is zero or denormalized. */ \
548 if (_FP_FRAC_ZEROP_##wc(X)) \
550 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
551 _FP_FRAC_COPY_##wc(R, Y); \
556 FP_SET_EXCEPTION(FP_EX_DENORM); \
560 _FP_FRAC_SUB_##wc(R, Y, X); \
563 if (Y##_e == _FP_EXPMAX_##fs) \
565 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
566 _FP_FRAC_COPY_##wc(R, Y); \
572 else if (Y##_e == _FP_EXPMAX_##fs) \
574 /* Y is NaN or Inf, X is normal. */ \
575 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
576 _FP_FRAC_COPY_##wc(R, Y); \
580 /* Insert implicit MSB of X. */ \
581 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \
584 /* Shift the mantissa of X to the right EDIFF steps; \
585 remember to account later for the implicit MSB of Y. */ \
586 if (ediff <= _FP_WFRACBITS_##fs) \
587 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \
588 else if (!_FP_FRAC_ZEROP_##wc(X)) \
589 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \
590 _FP_FRAC_SUB_##wc(R, Y, X); \
595 if (!_FP_EXP_NORMAL(fs, wc, X)) \
599 /* X and Y are zero or denormalized. */ \
601 if (_FP_FRAC_ZEROP_##wc(X)) \
603 _FP_FRAC_COPY_##wc(R, Y); \
604 if (_FP_FRAC_ZEROP_##wc(Y)) \
605 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
608 FP_SET_EXCEPTION(FP_EX_DENORM); \
613 else if (_FP_FRAC_ZEROP_##wc(Y)) \
615 FP_SET_EXCEPTION(FP_EX_DENORM); \
616 _FP_FRAC_COPY_##wc(R, X); \
622 FP_SET_EXCEPTION(FP_EX_DENORM); \
623 _FP_FRAC_SUB_##wc(R, X, Y); \
625 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
627 /* |X| < |Y|, negate result. */ \
628 _FP_FRAC_SUB_##wc(R, Y, X); \
631 else if (_FP_FRAC_ZEROP_##wc(R)) \
632 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
638 /* X and Y are NaN or Inf, of opposite signs. */ \
639 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \
640 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \
641 R##_e = _FP_EXPMAX_##fs; \
642 if (_FP_FRAC_ZEROP_##wc(X)) \
644 if (_FP_FRAC_ZEROP_##wc(Y)) \
647 R##_s = _FP_NANSIGN_##fs; \
648 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
649 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \
650 FP_SET_EXCEPTION(FP_EX_INVALID); \
656 _FP_FRAC_COPY_##wc(R, Y); \
661 if (_FP_FRAC_ZEROP_##wc(Y)) \
665 _FP_FRAC_COPY_##wc(R, X); \
670 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \
676 /* The exponents of X and Y, both normal, are equal. The \
677 implicit MSBs cancel. */ \
679 _FP_FRAC_SUB_##wc(R, X, Y); \
681 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
683 /* |X| < |Y|, negate result. */ \
684 _FP_FRAC_SUB_##wc(R, Y, X); \
687 else if (_FP_FRAC_ZEROP_##wc(R)) \
690 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \
696 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \
699 /* Carry into most significant bit of larger one of X and Y, \
700 canceling it; renormalize. */ \
701 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \
703 _FP_FRAC_CLZ_##wc(diff, R); \
704 diff -= _FP_WFRACXBITS_##fs; \
705 _FP_FRAC_SLL_##wc(R, diff); \
708 /* R is denormalized. */ \
709 diff = diff - R##_e + 1; \
710 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \
716 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
723 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
724 #define _FP_SUB(fs, wc, R, X, Y) \
726 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \
727 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \
732 * Main negation routine. FIXME -- when we care about setting exception
733 * bits reliably, this will not do. We should examine all of the fp classes.
736 #define _FP_NEG(fs, wc, R, X) \
738 _FP_FRAC_COPY_##wc(R, X); \
746 * Main multiplication routine. The input values should be cooked.
749 #define _FP_MUL(fs, wc, R, X, Y) \
751 R##_s = X##_s ^ Y##_s; \
752 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
754 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
755 R##_c = FP_CLS_NORMAL; \
756 R##_e = X##_e + Y##_e + 1; \
758 _FP_MUL_MEAT_##fs(R,X,Y); \
760 if (_FP_FRAC_OVERP_##wc(fs, R)) \
761 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \
766 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
767 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \
770 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
771 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
772 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
775 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
776 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
777 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
778 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
779 _FP_FRAC_COPY_##wc(R, X); \
783 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
784 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
785 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
788 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
789 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
790 _FP_FRAC_COPY_##wc(R, Y); \
794 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
795 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
796 R##_s = _FP_NANSIGN_##fs; \
797 R##_c = FP_CLS_NAN; \
798 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
799 FP_SET_EXCEPTION(FP_EX_INVALID); \
809 * Main division routine. The input values should be cooked.
812 #define _FP_DIV(fs, wc, R, X, Y) \
814 R##_s = X##_s ^ Y##_s; \
815 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \
817 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \
818 R##_c = FP_CLS_NORMAL; \
819 R##_e = X##_e - Y##_e; \
821 _FP_DIV_MEAT_##fs(R,X,Y); \
824 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \
825 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \
828 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \
829 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \
830 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \
832 _FP_FRAC_COPY_##wc(R, X); \
836 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \
837 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \
838 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \
840 _FP_FRAC_COPY_##wc(R, Y); \
844 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \
845 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \
846 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \
847 R##_c = FP_CLS_ZERO; \
850 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \
851 FP_SET_EXCEPTION(FP_EX_DIVZERO); \
852 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \
853 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \
854 R##_c = FP_CLS_INF; \
857 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \
858 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \
859 R##_s = _FP_NANSIGN_##fs; \
860 R##_c = FP_CLS_NAN; \
861 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
862 FP_SET_EXCEPTION(FP_EX_INVALID); \
872 * Main differential comparison routine. The inputs should be raw not
873 * cooked. The return is -1,0,1 for normal values, 2 otherwise.
876 #define _FP_CMP(fs, wc, ret, X, Y, un) \
878 /* NANs are unordered */ \
879 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
880 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
889 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \
890 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \
892 if (__is_zero_x && __is_zero_y) \
894 else if (__is_zero_x) \
895 ret = Y##_s ? 1 : -1; \
896 else if (__is_zero_y) \
897 ret = X##_s ? -1 : 1; \
898 else if (X##_s != Y##_s) \
899 ret = X##_s ? -1 : 1; \
900 else if (X##_e > Y##_e) \
901 ret = X##_s ? -1 : 1; \
902 else if (X##_e < Y##_e) \
903 ret = X##_s ? 1 : -1; \
904 else if (_FP_FRAC_GT_##wc(X, Y)) \
905 ret = X##_s ? -1 : 1; \
906 else if (_FP_FRAC_GT_##wc(Y, X)) \
907 ret = X##_s ? 1 : -1; \
914 /* Simplification for strict equality. */
916 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \
918 /* NANs are unordered */ \
919 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
920 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \
926 ret = !(X##_e == Y##_e \
927 && _FP_FRAC_EQ_##wc(X, Y) \
928 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \
932 /* Version to test unordered. */
934 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \
936 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \
937 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \
941 * Main square root routine. The input value should be cooked.
944 #define _FP_SQRT(fs, wc, R, X) \
946 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \
951 _FP_FRAC_COPY_##wc(R, X); \
953 R##_c = FP_CLS_NAN; \
958 R##_s = _FP_NANSIGN_##fs; \
959 R##_c = FP_CLS_NAN; /* NAN */ \
960 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
961 FP_SET_EXCEPTION(FP_EX_INVALID); \
966 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \
971 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \
973 case FP_CLS_NORMAL: \
977 R##_c = FP_CLS_NAN; /* sNAN */ \
978 R##_s = _FP_NANSIGN_##fs; \
979 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \
980 FP_SET_EXCEPTION(FP_EX_INVALID); \
983 R##_c = FP_CLS_NORMAL; \
985 _FP_FRAC_SLL_##wc(X, 1); \
986 R##_e = X##_e >> 1; \
987 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \
988 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \
989 q = _FP_OVERFLOW_##fs >> 1; \
990 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \
995 * Convert from FP to integer. Input is raw.
998 /* RSIGNED can have following values:
999 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1000 * the result is either 0 or (2^rsize)-1 depending on the sign in such
1002 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1003 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1004 * depending on the sign in such case.
1005 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1006 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1007 * depending on the sign in such case.
1009 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \
1011 if (X##_e < _FP_EXPBIAS_##fs) \
1016 if (!_FP_FRAC_ZEROP_##wc(X)) \
1018 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1019 FP_SET_EXCEPTION(FP_EX_DENORM); \
1023 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1025 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \
1026 || (!rsigned && X##_s)) \
1028 /* Overflow or converting to the most negative integer. */ \
1040 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \
1042 /* Possibly converting to most negative integer; check the \
1045 (void)((_FP_FRACBITS_##fs > rsize) \
1046 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \
1047 _FP_FRACBITS_##fs - rsize, \
1048 _FP_FRACBITS_##fs); 0; }) \
1050 if (!_FP_FRAC_ZEROP_##wc(X)) \
1051 FP_SET_EXCEPTION(FP_EX_INVALID); \
1053 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1056 FP_SET_EXCEPTION(FP_EX_INVALID); \
1060 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \
1061 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \
1063 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1064 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1069 _FP_FRAC_SRST_##wc(X, inexact, \
1070 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1072 _FP_FRACBITS_##fs); \
1074 FP_SET_EXCEPTION(FP_EX_INEXACT); \
1075 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \
1077 if (rsigned && X##_s) \
1082 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if
1084 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \
1090 if ((X##_s = (r < 0))) \
1094 (void)((rsize <= _FP_W_TYPE_SIZE) \
1097 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \
1098 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \
1100 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \
1103 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \
1105 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1110 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \
1111 && X##_e >= _FP_EXPMAX_##fs) \
1113 /* Exponent too big; overflow to infinity. (May also \
1114 happen after rounding below.) */ \
1115 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \
1116 goto pack_semiraw; \
1119 if (rsize <= _FP_FRACBITS_##fs \
1120 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \
1122 /* Exactly representable; shift left. */ \
1123 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1124 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1125 + _FP_FRACBITS_##fs - 1 - X##_e)); \
1129 /* More bits in integer than in floating type; need to \
1131 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \
1132 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \
1133 - _FP_WFRACBITS_##fs + 1)) \
1134 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \
1135 - _FP_WFRACBITS_##fs + 1))) \
1137 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \
1138 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1139 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \
1140 + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1141 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \
1143 _FP_PACK_SEMIRAW(fs, wc, X); \
1150 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \
1155 /* Extend from a narrower floating-point format to a wider one. Input
1156 and output are raw. */
1157 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \
1159 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \
1160 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \
1161 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \
1162 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1163 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \
1166 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1167 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1169 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1170 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1176 if (_FP_FRAC_ZEROP_##swc(S)) \
1178 else if (_FP_EXPBIAS_##dfs \
1179 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \
1181 FP_SET_EXCEPTION(FP_EX_DENORM); \
1182 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1183 - _FP_FRACBITS_##sfs)); \
1189 FP_SET_EXCEPTION(FP_EX_DENORM); \
1190 _FP_FRAC_CLZ_##swc(_lz, S); \
1191 _FP_FRAC_SLL_##dwc(D, \
1192 _lz + _FP_FRACBITS_##dfs \
1193 - _FP_FRACTBITS_##sfs); \
1194 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \
1195 + _FP_FRACXBITS_##sfs - _lz); \
1200 D##_e = _FP_EXPMAX_##dfs; \
1201 if (!_FP_FRAC_ZEROP_##swc(S)) \
1203 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \
1204 FP_SET_EXCEPTION(FP_EX_INVALID); \
1205 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \
1206 - _FP_FRACBITS_##sfs)); \
1212 /* Truncate from a wider floating-point format to a narrower one.
1213 Input and output are semi-raw. */
1214 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \
1216 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \
1217 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1218 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \
1221 if (_FP_EXP_NORMAL(sfs, swc, S)) \
1223 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \
1224 if (D##_e >= _FP_EXPMAX_##dfs) \
1225 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \
1230 if (D##_e < 1 - _FP_FRACBITS_##dfs) \
1232 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \
1233 _FP_FRAC_LOW_##swc(S) |= 1; \
1237 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \
1238 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1239 - _FP_WFRACBITS_##dfs + 1 - D##_e), \
1240 _FP_WFRACBITS_##sfs); \
1245 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1246 - _FP_WFRACBITS_##dfs), \
1247 _FP_WFRACBITS_##sfs); \
1248 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1256 if (_FP_FRAC_ZEROP_##swc(S)) \
1257 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1260 FP_SET_EXCEPTION(FP_EX_DENORM); \
1261 if (_FP_EXPBIAS_##sfs \
1262 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \
1264 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \
1265 - _FP_WFRACBITS_##dfs), \
1266 _FP_WFRACBITS_##sfs); \
1267 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1271 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1272 _FP_FRAC_LOW_##dwc(D) |= 1; \
1278 D##_e = _FP_EXPMAX_##dfs; \
1279 if (_FP_FRAC_ZEROP_##swc(S)) \
1280 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \
1283 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \
1284 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \
1285 - _FP_WFRACBITS_##dfs)); \
1286 _FP_FRAC_COPY_##dwc##_##swc(D, S); \
1287 /* Semi-raw NaN must have all workbits cleared. */ \
1288 _FP_FRAC_LOW_##dwc(D) \
1289 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \
1290 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \
1297 * Helper primitives.
1300 /* Count leading zeros in a word. */
1303 /* GCC 3.4 and later provide the builtins for us. */
1304 #define __FP_CLZ(r, x) \
1306 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \
1307 r = __builtin_clz (x); \
1308 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \
1309 r = __builtin_clzl (x); \
1310 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \
1311 r = __builtin_clzll (x); \
1315 #endif /* ndef __FP_CLZ */
1317 #define _FP_DIV_HELP_imm(q, r, n, d) \
1319 q = n / d, r = n % d; \
1323 /* A restoring bit-by-bit division primitive. */
1325 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \
1327 int count = _FP_WFRACBITS_##fs; \
1328 _FP_FRAC_DECL_##wc (u); \
1329 _FP_FRAC_DECL_##wc (v); \
1330 _FP_FRAC_COPY_##wc (u, X); \
1331 _FP_FRAC_COPY_##wc (v, Y); \
1332 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \
1333 /* Normalize U and V. */ \
1334 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \
1335 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \
1336 /* First round. Since the operands are normalized, either the \
1337 first or second bit will be set in the fraction. Produce a \
1338 normalized result by checking which and adjusting the loop \
1339 count and exponent accordingly. */ \
1340 if (_FP_FRAC_GE_1 (u, v)) \
1342 _FP_FRAC_SUB_##wc (u, u, v); \
1343 _FP_FRAC_LOW_##wc (R) |= 1; \
1348 /* Subsequent rounds. */ \
1350 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \
1351 _FP_FRAC_SLL_##wc (u, 1); \
1352 _FP_FRAC_SLL_##wc (R, 1); \
1353 if (msb || _FP_FRAC_GE_1 (u, v)) \
1355 _FP_FRAC_SUB_##wc (u, u, v); \
1356 _FP_FRAC_LOW_##wc (R) |= 1; \
1358 } while (--count > 0); \
1359 /* If there's anything left in U, the result is inexact. */ \
1360 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \
1363 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1364 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1365 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)