xref: /netbsd-src/external/gpl3/gcc/dist/libgcc/soft-fp/op-common.h (revision bdc22b2e01993381dcefeff2bc9b56ca75a4235c)
1 /* Software floating-point emulation. Common operations.
2    Copyright (C) 1997-2015 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).
8 
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.
13 
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.)
22 
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.
27 
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/>.  */
31 
32 #ifndef SOFT_FP_OP_COMMON_H
33 #define SOFT_FP_OP_COMMON_H	1
34 
35 #define _FP_DECL(wc, X)						\
36   _FP_I_TYPE X##_c __attribute__ ((unused)) _FP_ZERO_INIT;	\
37   _FP_I_TYPE X##_s __attribute__ ((unused)) _FP_ZERO_INIT;	\
38   _FP_I_TYPE X##_e __attribute__ ((unused)) _FP_ZERO_INIT;	\
39   _FP_FRAC_DECL_##wc (X)
40 
41 /* Test whether the qNaN bit denotes a signaling NaN.  */
42 #define _FP_FRAC_SNANP(fs, X)				\
43   ((_FP_QNANNEGATEDP)					\
44    ? (_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs)	\
45    : !(_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs))
46 #define _FP_FRAC_SNANP_SEMIRAW(fs, X)			\
47   ((_FP_QNANNEGATEDP)					\
48    ? (_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs)	\
49    : !(_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs))
50 
51 /* Finish truly unpacking a native fp value by classifying the kind
52    of fp value and normalizing both the exponent and the fraction.  */
53 
54 #define _FP_UNPACK_CANONICAL(fs, wc, X)				\
55   do								\
56     {								\
57       switch (X##_e)						\
58 	{							\
59 	default:						\
60 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;	\
61 	  _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);			\
62 	  X##_e -= _FP_EXPBIAS_##fs;				\
63 	  X##_c = FP_CLS_NORMAL;				\
64 	  break;						\
65 								\
66 	case 0:							\
67 	  if (_FP_FRAC_ZEROP_##wc (X))				\
68 	    X##_c = FP_CLS_ZERO;				\
69 	  else if (FP_DENORM_ZERO)				\
70 	    {							\
71 	      X##_c = FP_CLS_ZERO;				\
72 	      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
73 	      FP_SET_EXCEPTION (FP_EX_DENORM);			\
74 	    }							\
75 	  else							\
76 	    {							\
77 	      /* A denormalized number.  */			\
78 	      _FP_I_TYPE _FP_UNPACK_CANONICAL_shift;		\
79 	      _FP_FRAC_CLZ_##wc (_FP_UNPACK_CANONICAL_shift,	\
80 				 X);				\
81 	      _FP_UNPACK_CANONICAL_shift -= _FP_FRACXBITS_##fs;	\
82 	      _FP_FRAC_SLL_##wc (X, (_FP_UNPACK_CANONICAL_shift \
83 				     + _FP_WORKBITS));		\
84 	      X##_e -= (_FP_EXPBIAS_##fs - 1			\
85 			+ _FP_UNPACK_CANONICAL_shift);		\
86 	      X##_c = FP_CLS_NORMAL;				\
87 	      FP_SET_EXCEPTION (FP_EX_DENORM);			\
88 	    }							\
89 	  break;						\
90 								\
91 	case _FP_EXPMAX_##fs:					\
92 	  if (_FP_FRAC_ZEROP_##wc (X))				\
93 	    X##_c = FP_CLS_INF;					\
94 	  else							\
95 	    {							\
96 	      X##_c = FP_CLS_NAN;				\
97 	      /* Check for signaling NaN.  */			\
98 	      if (_FP_FRAC_SNANP (fs, X))			\
99 		FP_SET_EXCEPTION (FP_EX_INVALID			\
100 				  | FP_EX_INVALID_SNAN);	\
101 	    }							\
102 	  break;						\
103 	}							\
104     }								\
105   while (0)
106 
107 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
108    shifted by _FP_WORKBITS but the implicit MSB is not inserted and
109    other classification is not done.  */
110 #define _FP_UNPACK_SEMIRAW(fs, wc, X)	_FP_FRAC_SLL_##wc (X, _FP_WORKBITS)
111 
112 /* Check whether a raw or semi-raw input value should be flushed to
113    zero, and flush it to zero if so.  */
114 #define _FP_CHECK_FLUSH_ZERO(fs, wc, X)			\
115   do							\
116     {							\
117       if (FP_DENORM_ZERO				\
118 	  && X##_e == 0					\
119 	  && !_FP_FRAC_ZEROP_##wc (X))			\
120 	{						\
121 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);	\
122 	  FP_SET_EXCEPTION (FP_EX_DENORM);		\
123 	}						\
124     }							\
125   while (0)
126 
127 /* A semi-raw value has overflowed to infinity.  Adjust the mantissa
128    and exponent appropriately.  */
129 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X)			\
130   do							\
131     {							\
132       if (FP_ROUNDMODE == FP_RND_NEAREST		\
133 	  || (FP_ROUNDMODE == FP_RND_PINF && !X##_s)	\
134 	  || (FP_ROUNDMODE == FP_RND_MINF && X##_s))	\
135 	{						\
136 	  X##_e = _FP_EXPMAX_##fs;			\
137 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);	\
138 	}						\
139       else						\
140 	{						\
141 	  X##_e = _FP_EXPMAX_##fs - 1;			\
142 	  _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);	\
143 	}						\
144       FP_SET_EXCEPTION (FP_EX_INEXACT);			\
145       FP_SET_EXCEPTION (FP_EX_OVERFLOW);		\
146     }							\
147   while (0)
148 
149 /* Check for a semi-raw value being a signaling NaN and raise the
150    invalid exception if so.  */
151 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X)			\
152   do								\
153     {								\
154       if (X##_e == _FP_EXPMAX_##fs				\
155 	  && !_FP_FRAC_ZEROP_##wc (X)				\
156 	  && _FP_FRAC_SNANP_SEMIRAW (fs, X))			\
157 	FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SNAN);	\
158     }								\
159   while (0)
160 
161 /* Choose a NaN result from an operation on two semi-raw NaN
162    values.  */
163 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP)			\
164   do									\
165     {									\
166       /* _FP_CHOOSENAN expects raw values, so shift as required.  */	\
167       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);				\
168       _FP_FRAC_SRL_##wc (Y, _FP_WORKBITS);				\
169       _FP_CHOOSENAN (fs, wc, R, X, Y, OP);				\
170       _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);				\
171     }									\
172   while (0)
173 
174 /* Make the fractional part a quiet NaN, preserving the payload
175    if possible, otherwise make it the canonical quiet NaN and set
176    the sign bit accordingly.  */
177 #define _FP_SETQNAN(fs, wc, X)					\
178   do								\
179     {								\
180       if (_FP_QNANNEGATEDP)					\
181 	{							\
182 	  _FP_FRAC_HIGH_RAW_##fs (X) &= _FP_QNANBIT_##fs - 1;	\
183 	  if (_FP_FRAC_ZEROP_##wc (X))				\
184 	    {							\
185 	      X##_s = _FP_NANSIGN_##fs;				\
186 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
187 	    }							\
188 	}							\
189       else							\
190 	_FP_FRAC_HIGH_RAW_##fs (X) |= _FP_QNANBIT_##fs;		\
191     }								\
192   while (0)
193 #define _FP_SETQNAN_SEMIRAW(fs, wc, X)				\
194   do								\
195     {								\
196       if (_FP_QNANNEGATEDP)					\
197 	{							\
198 	  _FP_FRAC_HIGH_##fs (X) &= _FP_QNANBIT_SH_##fs - 1;	\
199 	  if (_FP_FRAC_ZEROP_##wc (X))				\
200 	    {							\
201 	      X##_s = _FP_NANSIGN_##fs;				\
202 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
203 	      _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);		\
204 	    }							\
205 	}							\
206       else							\
207 	_FP_FRAC_HIGH_##fs (X) |= _FP_QNANBIT_SH_##fs;		\
208     }								\
209   while (0)
210 
211 /* Test whether a biased exponent is normal (not zero or maximum).  */
212 #define _FP_EXP_NORMAL(fs, wc, X)	(((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
213 
214 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
215    rounded and shifted right, with the rounding possibly increasing
216    the exponent (including changing a finite value to infinity).  */
217 #define _FP_PACK_SEMIRAW(fs, wc, X)				\
218   do								\
219     {								\
220       int _FP_PACK_SEMIRAW_is_tiny				\
221 	= X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X);		\
222       if (_FP_TININESS_AFTER_ROUNDING				\
223 	  && _FP_PACK_SEMIRAW_is_tiny)				\
224 	{							\
225 	  FP_DECL_##fs (_FP_PACK_SEMIRAW_T);			\
226 	  _FP_FRAC_COPY_##wc (_FP_PACK_SEMIRAW_T, X);		\
227 	  _FP_PACK_SEMIRAW_T##_s = X##_s;			\
228 	  _FP_PACK_SEMIRAW_T##_e = X##_e;			\
229 	  _FP_FRAC_SLL_##wc (_FP_PACK_SEMIRAW_T, 1);		\
230 	  _FP_ROUND (wc, _FP_PACK_SEMIRAW_T);			\
231 	  if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_SEMIRAW_T))	\
232 	    _FP_PACK_SEMIRAW_is_tiny = 0;			\
233 	}							\
234       _FP_ROUND (wc, X);					\
235       if (_FP_PACK_SEMIRAW_is_tiny)				\
236 	{							\
237 	  if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)		\
238 	      || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW))	\
239 	    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
240 	}							\
241       if (_FP_FRAC_HIGH_##fs (X)				\
242 	  & (_FP_OVERFLOW_##fs >> 1))				\
243 	{							\
244 	  _FP_FRAC_HIGH_##fs (X) &= ~(_FP_OVERFLOW_##fs >> 1);	\
245 	  X##_e++;						\
246 	  if (X##_e == _FP_EXPMAX_##fs)				\
247 	    _FP_OVERFLOW_SEMIRAW (fs, wc, X);			\
248 	}							\
249       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
250       if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
251 	{							\
252 	  if (!_FP_KEEPNANFRACP)				\
253 	    {							\
254 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
255 	      X##_s = _FP_NANSIGN_##fs;				\
256 	    }							\
257 	  else							\
258 	    _FP_SETQNAN (fs, wc, X);				\
259 	}							\
260     }								\
261   while (0)
262 
263 /* Before packing the bits back into the native fp result, take care
264    of such mundane things as rounding and overflow.  Also, for some
265    kinds of fp values, the original parts may not have been fully
266    extracted -- but that is ok, we can regenerate them now.  */
267 
268 #define _FP_PACK_CANONICAL(fs, wc, X)					\
269   do									\
270     {									\
271       switch (X##_c)							\
272 	{								\
273 	case FP_CLS_NORMAL:						\
274 	  X##_e += _FP_EXPBIAS_##fs;					\
275 	  if (X##_e > 0)						\
276 	    {								\
277 	      _FP_ROUND (wc, X);					\
278 	      if (_FP_FRAC_OVERP_##wc (fs, X))				\
279 		{							\
280 		  _FP_FRAC_CLEAR_OVERP_##wc (fs, X);			\
281 		  X##_e++;						\
282 		}							\
283 	      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
284 	      if (X##_e >= _FP_EXPMAX_##fs)				\
285 		{							\
286 		  /* Overflow.  */					\
287 		  switch (FP_ROUNDMODE)					\
288 		    {							\
289 		    case FP_RND_NEAREST:				\
290 		      X##_c = FP_CLS_INF;				\
291 		      break;						\
292 		    case FP_RND_PINF:					\
293 		      if (!X##_s)					\
294 			X##_c = FP_CLS_INF;				\
295 		      break;						\
296 		    case FP_RND_MINF:					\
297 		      if (X##_s)					\
298 			X##_c = FP_CLS_INF;				\
299 		      break;						\
300 		    }							\
301 		  if (X##_c == FP_CLS_INF)				\
302 		    {							\
303 		      /* Overflow to infinity.  */			\
304 		      X##_e = _FP_EXPMAX_##fs;				\
305 		      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
306 		    }							\
307 		  else							\
308 		    {							\
309 		      /* Overflow to maximum normal.  */		\
310 		      X##_e = _FP_EXPMAX_##fs - 1;			\
311 		      _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);		\
312 		    }							\
313 		  FP_SET_EXCEPTION (FP_EX_OVERFLOW);			\
314 		  FP_SET_EXCEPTION (FP_EX_INEXACT);			\
315 		}							\
316 	    }								\
317 	  else								\
318 	    {								\
319 	      /* We've got a denormalized number.  */			\
320 	      int _FP_PACK_CANONICAL_is_tiny = 1;			\
321 	      if (_FP_TININESS_AFTER_ROUNDING && X##_e == 0)		\
322 		{							\
323 		  FP_DECL_##fs (_FP_PACK_CANONICAL_T);			\
324 		  _FP_FRAC_COPY_##wc (_FP_PACK_CANONICAL_T, X);		\
325 		  _FP_PACK_CANONICAL_T##_s = X##_s;			\
326 		  _FP_PACK_CANONICAL_T##_e = X##_e;			\
327 		  _FP_ROUND (wc, _FP_PACK_CANONICAL_T);			\
328 		  if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_CANONICAL_T))	\
329 		    _FP_PACK_CANONICAL_is_tiny = 0;			\
330 		}							\
331 	      X##_e = -X##_e + 1;					\
332 	      if (X##_e <= _FP_WFRACBITS_##fs)				\
333 		{							\
334 		  _FP_FRAC_SRS_##wc (X, X##_e, _FP_WFRACBITS_##fs);	\
335 		  _FP_ROUND (wc, X);					\
336 		  if (_FP_FRAC_HIGH_##fs (X)				\
337 		      & (_FP_OVERFLOW_##fs >> 1))			\
338 		    {							\
339 		      X##_e = 1;					\
340 		      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
341 		      FP_SET_EXCEPTION (FP_EX_INEXACT);			\
342 		    }							\
343 		  else							\
344 		    {							\
345 		      X##_e = 0;					\
346 		      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);		\
347 		    }							\
348 		  if (_FP_PACK_CANONICAL_is_tiny			\
349 		      && ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)		\
350 			  || (FP_TRAPPING_EXCEPTIONS			\
351 			      & FP_EX_UNDERFLOW)))			\
352 		    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
353 		}							\
354 	      else							\
355 		{							\
356 		  /* Underflow to zero.  */				\
357 		  X##_e = 0;						\
358 		  if (!_FP_FRAC_ZEROP_##wc (X))				\
359 		    {							\
360 		      _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);		\
361 		      _FP_ROUND (wc, X);				\
362 		      _FP_FRAC_LOW_##wc (X) >>= (_FP_WORKBITS);		\
363 		    }							\
364 		  FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
365 		}							\
366 	    }								\
367 	  break;							\
368 									\
369 	case FP_CLS_ZERO:						\
370 	  X##_e = 0;							\
371 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
372 	  break;							\
373 									\
374 	case FP_CLS_INF:						\
375 	  X##_e = _FP_EXPMAX_##fs;					\
376 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
377 	  break;							\
378 									\
379 	case FP_CLS_NAN:						\
380 	  X##_e = _FP_EXPMAX_##fs;					\
381 	  if (!_FP_KEEPNANFRACP)					\
382 	    {								\
383 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);			\
384 	      X##_s = _FP_NANSIGN_##fs;					\
385 	    }								\
386 	  else								\
387 	    _FP_SETQNAN (fs, wc, X);					\
388 	  break;							\
389 	}								\
390     }									\
391   while (0)
392 
393 /* This one accepts raw argument and not cooked,  returns
394    1 if X is a signaling NaN.  */
395 #define _FP_ISSIGNAN(fs, wc, X)			\
396   ({						\
397     int _FP_ISSIGNAN_ret = 0;			\
398     if (X##_e == _FP_EXPMAX_##fs)		\
399       {						\
400 	if (!_FP_FRAC_ZEROP_##wc (X)		\
401 	    && _FP_FRAC_SNANP (fs, X))		\
402 	  _FP_ISSIGNAN_ret = 1;			\
403       }						\
404     _FP_ISSIGNAN_ret;				\
405   })
406 
407 
408 
409 
410 
411 /* Addition on semi-raw values.  */
412 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP)				\
413   do									\
414     {									\
415       _FP_CHECK_FLUSH_ZERO (fs, wc, X);					\
416       _FP_CHECK_FLUSH_ZERO (fs, wc, Y);					\
417       if (X##_s == Y##_s)						\
418 	{								\
419 	  /* Addition.  */						\
420 	  __label__ add1, add2, add3, add_done;				\
421 	  R##_s = X##_s;						\
422 	  int _FP_ADD_INTERNAL_ediff = X##_e - Y##_e;			\
423 	  if (_FP_ADD_INTERNAL_ediff > 0)				\
424 	    {								\
425 	      R##_e = X##_e;						\
426 	      if (Y##_e == 0)						\
427 		{							\
428 		  /* Y is zero or denormalized.  */			\
429 		  if (_FP_FRAC_ZEROP_##wc (Y))				\
430 		    {							\
431 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
432 		      _FP_FRAC_COPY_##wc (R, X);			\
433 		      goto add_done;					\
434 		    }							\
435 		  else							\
436 		    {							\
437 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
438 		      _FP_ADD_INTERNAL_ediff--;				\
439 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
440 			{						\
441 			  _FP_FRAC_ADD_##wc (R, X, Y);			\
442 			  goto add3;					\
443 			}						\
444 		      if (X##_e == _FP_EXPMAX_##fs)			\
445 			{						\
446 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
447 			  _FP_FRAC_COPY_##wc (R, X);			\
448 			  goto add_done;				\
449 			}						\
450 		      goto add1;					\
451 		    }							\
452 		}							\
453 	      else if (X##_e == _FP_EXPMAX_##fs)			\
454 		{							\
455 		  /* X is NaN or Inf, Y is normal.  */			\
456 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);			\
457 		  _FP_FRAC_COPY_##wc (R, X);				\
458 		  goto add_done;					\
459 		}							\
460 									\
461 	      /* Insert implicit MSB of Y.  */				\
462 	      _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs;		\
463 									\
464 	    add1:							\
465 	      /* Shift the mantissa of Y to the right			\
466 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
467 		 later for the implicit MSB of X.  */			\
468 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
469 		_FP_FRAC_SRS_##wc (Y, _FP_ADD_INTERNAL_ediff,		\
470 				   _FP_WFRACBITS_##fs);			\
471 	      else if (!_FP_FRAC_ZEROP_##wc (Y))			\
472 		_FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc);		\
473 	      _FP_FRAC_ADD_##wc (R, X, Y);				\
474 	    }								\
475 	  else if (_FP_ADD_INTERNAL_ediff < 0)				\
476 	    {								\
477 	      _FP_ADD_INTERNAL_ediff = -_FP_ADD_INTERNAL_ediff;		\
478 	      R##_e = Y##_e;						\
479 	      if (X##_e == 0)						\
480 		{							\
481 		  /* X is zero or denormalized.  */			\
482 		  if (_FP_FRAC_ZEROP_##wc (X))				\
483 		    {							\
484 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
485 		      _FP_FRAC_COPY_##wc (R, Y);			\
486 		      goto add_done;					\
487 		    }							\
488 		  else							\
489 		    {							\
490 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
491 		      _FP_ADD_INTERNAL_ediff--;				\
492 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
493 			{						\
494 			  _FP_FRAC_ADD_##wc (R, Y, X);			\
495 			  goto add3;					\
496 			}						\
497 		      if (Y##_e == _FP_EXPMAX_##fs)			\
498 			{						\
499 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
500 			  _FP_FRAC_COPY_##wc (R, Y);			\
501 			  goto add_done;				\
502 			}						\
503 		      goto add2;					\
504 		    }							\
505 		}							\
506 	      else if (Y##_e == _FP_EXPMAX_##fs)			\
507 		{							\
508 		  /* Y is NaN or Inf, X is normal.  */			\
509 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);			\
510 		  _FP_FRAC_COPY_##wc (R, Y);				\
511 		  goto add_done;					\
512 		}							\
513 									\
514 	      /* Insert implicit MSB of X.  */				\
515 	      _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs;		\
516 									\
517 	    add2:							\
518 	      /* Shift the mantissa of X to the right			\
519 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
520 		 later for the implicit MSB of Y.  */			\
521 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
522 		_FP_FRAC_SRS_##wc (X, _FP_ADD_INTERNAL_ediff,		\
523 				   _FP_WFRACBITS_##fs);			\
524 	      else if (!_FP_FRAC_ZEROP_##wc (X))			\
525 		_FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);		\
526 	      _FP_FRAC_ADD_##wc (R, Y, X);				\
527 	    }								\
528 	  else								\
529 	    {								\
530 	      /* _FP_ADD_INTERNAL_ediff == 0.  */			\
531 	      if (!_FP_EXP_NORMAL (fs, wc, X))				\
532 		{							\
533 		  if (X##_e == 0)					\
534 		    {							\
535 		      /* X and Y are zero or denormalized.  */		\
536 		      R##_e = 0;					\
537 		      if (_FP_FRAC_ZEROP_##wc (X))			\
538 			{						\
539 			  if (!_FP_FRAC_ZEROP_##wc (Y))			\
540 			    FP_SET_EXCEPTION (FP_EX_DENORM);		\
541 			  _FP_FRAC_COPY_##wc (R, Y);			\
542 			  goto add_done;				\
543 			}						\
544 		      else if (_FP_FRAC_ZEROP_##wc (Y))			\
545 			{						\
546 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
547 			  _FP_FRAC_COPY_##wc (R, X);			\
548 			  goto add_done;				\
549 			}						\
550 		      else						\
551 			{						\
552 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
553 			  _FP_FRAC_ADD_##wc (R, X, Y);			\
554 			  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
555 			    {						\
556 			      /* Normalized result.  */			\
557 			      _FP_FRAC_HIGH_##fs (R)			\
558 				&= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs;	\
559 			      R##_e = 1;				\
560 			    }						\
561 			  goto add_done;				\
562 			}						\
563 		    }							\
564 		  else							\
565 		    {							\
566 		      /* X and Y are NaN or Inf.  */			\
567 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
568 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
569 		      R##_e = _FP_EXPMAX_##fs;				\
570 		      if (_FP_FRAC_ZEROP_##wc (X))			\
571 			_FP_FRAC_COPY_##wc (R, Y);			\
572 		      else if (_FP_FRAC_ZEROP_##wc (Y))			\
573 			_FP_FRAC_COPY_##wc (R, X);			\
574 		      else						\
575 			_FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP);	\
576 		      goto add_done;					\
577 		    }							\
578 		}							\
579 	      /* The exponents of X and Y, both normal, are equal.  The	\
580 		 implicit MSBs will always add to increase the		\
581 		 exponent.  */						\
582 	      _FP_FRAC_ADD_##wc (R, X, Y);				\
583 	      R##_e = X##_e + 1;					\
584 	      _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);		\
585 	      if (R##_e == _FP_EXPMAX_##fs)				\
586 		/* Overflow to infinity (depending on rounding mode).  */ \
587 		_FP_OVERFLOW_SEMIRAW (fs, wc, R);			\
588 	      goto add_done;						\
589 	    }								\
590 	add3:								\
591 	  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)		\
592 	    {								\
593 	      /* Overflow.  */						\
594 	      _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
595 	      R##_e++;							\
596 	      _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);		\
597 	      if (R##_e == _FP_EXPMAX_##fs)				\
598 		/* Overflow to infinity (depending on rounding mode).  */ \
599 		_FP_OVERFLOW_SEMIRAW (fs, wc, R);			\
600 	    }								\
601 	add_done: ;							\
602 	}								\
603       else								\
604 	{								\
605 	  /* Subtraction.  */						\
606 	  __label__ sub1, sub2, sub3, norm, sub_done;			\
607 	  int _FP_ADD_INTERNAL_ediff = X##_e - Y##_e;			\
608 	  if (_FP_ADD_INTERNAL_ediff > 0)				\
609 	    {								\
610 	      R##_e = X##_e;						\
611 	      R##_s = X##_s;						\
612 	      if (Y##_e == 0)						\
613 		{							\
614 		  /* Y is zero or denormalized.  */			\
615 		  if (_FP_FRAC_ZEROP_##wc (Y))				\
616 		    {							\
617 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
618 		      _FP_FRAC_COPY_##wc (R, X);			\
619 		      goto sub_done;					\
620 		    }							\
621 		  else							\
622 		    {							\
623 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
624 		      _FP_ADD_INTERNAL_ediff--;				\
625 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
626 			{						\
627 			  _FP_FRAC_SUB_##wc (R, X, Y);			\
628 			  goto sub3;					\
629 			}						\
630 		      if (X##_e == _FP_EXPMAX_##fs)			\
631 			{						\
632 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
633 			  _FP_FRAC_COPY_##wc (R, X);			\
634 			  goto sub_done;				\
635 			}						\
636 		      goto sub1;					\
637 		    }							\
638 		}							\
639 	      else if (X##_e == _FP_EXPMAX_##fs)			\
640 		{							\
641 		  /* X is NaN or Inf, Y is normal.  */			\
642 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);			\
643 		  _FP_FRAC_COPY_##wc (R, X);				\
644 		  goto sub_done;					\
645 		}							\
646 									\
647 	      /* Insert implicit MSB of Y.  */				\
648 	      _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs;		\
649 									\
650 	    sub1:							\
651 	      /* Shift the mantissa of Y to the right			\
652 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
653 		 later for the implicit MSB of X.  */			\
654 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
655 		_FP_FRAC_SRS_##wc (Y, _FP_ADD_INTERNAL_ediff,		\
656 				   _FP_WFRACBITS_##fs);			\
657 	      else if (!_FP_FRAC_ZEROP_##wc (Y))			\
658 		_FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc);		\
659 	      _FP_FRAC_SUB_##wc (R, X, Y);				\
660 	    }								\
661 	  else if (_FP_ADD_INTERNAL_ediff < 0)				\
662 	    {								\
663 	      _FP_ADD_INTERNAL_ediff = -_FP_ADD_INTERNAL_ediff;		\
664 	      R##_e = Y##_e;						\
665 	      R##_s = Y##_s;						\
666 	      if (X##_e == 0)						\
667 		{							\
668 		  /* X is zero or denormalized.  */			\
669 		  if (_FP_FRAC_ZEROP_##wc (X))				\
670 		    {							\
671 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
672 		      _FP_FRAC_COPY_##wc (R, Y);			\
673 		      goto sub_done;					\
674 		    }							\
675 		  else							\
676 		    {							\
677 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
678 		      _FP_ADD_INTERNAL_ediff--;				\
679 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
680 			{						\
681 			  _FP_FRAC_SUB_##wc (R, Y, X);			\
682 			  goto sub3;					\
683 			}						\
684 		      if (Y##_e == _FP_EXPMAX_##fs)			\
685 			{						\
686 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
687 			  _FP_FRAC_COPY_##wc (R, Y);			\
688 			  goto sub_done;				\
689 			}						\
690 		      goto sub2;					\
691 		    }							\
692 		}							\
693 	      else if (Y##_e == _FP_EXPMAX_##fs)			\
694 		{							\
695 		  /* Y is NaN or Inf, X is normal.  */			\
696 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);			\
697 		  _FP_FRAC_COPY_##wc (R, Y);				\
698 		  goto sub_done;					\
699 		}							\
700 									\
701 	      /* Insert implicit MSB of X.  */				\
702 	      _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs;		\
703 									\
704 	    sub2:							\
705 	      /* Shift the mantissa of X to the right			\
706 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
707 		 later for the implicit MSB of Y.  */			\
708 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
709 		_FP_FRAC_SRS_##wc (X, _FP_ADD_INTERNAL_ediff,		\
710 				   _FP_WFRACBITS_##fs);			\
711 	      else if (!_FP_FRAC_ZEROP_##wc (X))			\
712 		_FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);		\
713 	      _FP_FRAC_SUB_##wc (R, Y, X);				\
714 	    }								\
715 	  else								\
716 	    {								\
717 	      /* ediff == 0.  */					\
718 	      if (!_FP_EXP_NORMAL (fs, wc, X))				\
719 		{							\
720 		  if (X##_e == 0)					\
721 		    {							\
722 		      /* X and Y are zero or denormalized.  */		\
723 		      R##_e = 0;					\
724 		      if (_FP_FRAC_ZEROP_##wc (X))			\
725 			{						\
726 			  _FP_FRAC_COPY_##wc (R, Y);			\
727 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
728 			    R##_s = (FP_ROUNDMODE == FP_RND_MINF);	\
729 			  else						\
730 			    {						\
731 			      FP_SET_EXCEPTION (FP_EX_DENORM);		\
732 			      R##_s = Y##_s;				\
733 			    }						\
734 			  goto sub_done;				\
735 			}						\
736 		      else if (_FP_FRAC_ZEROP_##wc (Y))			\
737 			{						\
738 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
739 			  _FP_FRAC_COPY_##wc (R, X);			\
740 			  R##_s = X##_s;				\
741 			  goto sub_done;				\
742 			}						\
743 		      else						\
744 			{						\
745 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
746 			  _FP_FRAC_SUB_##wc (R, X, Y);			\
747 			  R##_s = X##_s;				\
748 			  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
749 			    {						\
750 			      /* |X| < |Y|, negate result.  */		\
751 			      _FP_FRAC_SUB_##wc (R, Y, X);		\
752 			      R##_s = Y##_s;				\
753 			    }						\
754 			  else if (_FP_FRAC_ZEROP_##wc (R))		\
755 			    R##_s = (FP_ROUNDMODE == FP_RND_MINF);	\
756 			  goto sub_done;				\
757 			}						\
758 		    }							\
759 		  else							\
760 		    {							\
761 		      /* X and Y are NaN or Inf, of opposite signs.  */	\
762 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
763 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
764 		      R##_e = _FP_EXPMAX_##fs;				\
765 		      if (_FP_FRAC_ZEROP_##wc (X))			\
766 			{						\
767 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
768 			    {						\
769 			      /* Inf - Inf.  */				\
770 			      R##_s = _FP_NANSIGN_##fs;			\
771 			      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);	\
772 			      _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);	\
773 			      FP_SET_EXCEPTION (FP_EX_INVALID		\
774 						| FP_EX_INVALID_ISI);	\
775 			    }						\
776 			  else						\
777 			    {						\
778 			      /* Inf - NaN.  */				\
779 			      R##_s = Y##_s;				\
780 			      _FP_FRAC_COPY_##wc (R, Y);		\
781 			    }						\
782 			}						\
783 		      else						\
784 			{						\
785 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
786 			    {						\
787 			      /* NaN - Inf.  */				\
788 			      R##_s = X##_s;				\
789 			      _FP_FRAC_COPY_##wc (R, X);		\
790 			    }						\
791 			  else						\
792 			    {						\
793 			      /* NaN - NaN.  */				\
794 			      _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP); \
795 			    }						\
796 			}						\
797 		      goto sub_done;					\
798 		    }							\
799 		}							\
800 	      /* The exponents of X and Y, both normal, are equal.  The	\
801 		 implicit MSBs cancel.  */				\
802 	      R##_e = X##_e;						\
803 	      _FP_FRAC_SUB_##wc (R, X, Y);				\
804 	      R##_s = X##_s;						\
805 	      if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)		\
806 		{							\
807 		  /* |X| < |Y|, negate result.  */			\
808 		  _FP_FRAC_SUB_##wc (R, Y, X);				\
809 		  R##_s = Y##_s;					\
810 		}							\
811 	      else if (_FP_FRAC_ZEROP_##wc (R))				\
812 		{							\
813 		  R##_e = 0;						\
814 		  R##_s = (FP_ROUNDMODE == FP_RND_MINF);		\
815 		  goto sub_done;					\
816 		}							\
817 	      goto norm;						\
818 	    }								\
819 	sub3:								\
820 	  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)		\
821 	    {								\
822 	      int _FP_ADD_INTERNAL_diff;				\
823 	      /* Carry into most significant bit of larger one of X and Y, \
824 		 canceling it; renormalize.  */				\
825 	      _FP_FRAC_HIGH_##fs (R) &= _FP_IMPLBIT_SH_##fs - 1;	\
826 	    norm:							\
827 	      _FP_FRAC_CLZ_##wc (_FP_ADD_INTERNAL_diff, R);		\
828 	      _FP_ADD_INTERNAL_diff -= _FP_WFRACXBITS_##fs;		\
829 	      _FP_FRAC_SLL_##wc (R, _FP_ADD_INTERNAL_diff);		\
830 	      if (R##_e <= _FP_ADD_INTERNAL_diff)			\
831 		{							\
832 		  /* R is denormalized.  */				\
833 		  _FP_ADD_INTERNAL_diff					\
834 		    = _FP_ADD_INTERNAL_diff - R##_e + 1;		\
835 		  _FP_FRAC_SRS_##wc (R, _FP_ADD_INTERNAL_diff,		\
836 				     _FP_WFRACBITS_##fs);		\
837 		  R##_e = 0;						\
838 		}							\
839 	      else							\
840 		{							\
841 		  R##_e -= _FP_ADD_INTERNAL_diff;			\
842 		  _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
843 		}							\
844 	    }								\
845 	sub_done: ;							\
846 	}								\
847     }									\
848   while (0)
849 
850 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL (fs, wc, R, X, Y, '+')
851 #define _FP_SUB(fs, wc, R, X, Y)					\
852   do									\
853     {									\
854       if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
855 	Y##_s ^= 1;							\
856       _FP_ADD_INTERNAL (fs, wc, R, X, Y, '-');				\
857     }									\
858   while (0)
859 
860 
861 /* Main negation routine.  The input value is raw.  */
862 
863 #define _FP_NEG(fs, wc, R, X)			\
864   do						\
865     {						\
866       _FP_FRAC_COPY_##wc (R, X);		\
867       R##_e = X##_e;				\
868       R##_s = 1 ^ X##_s;			\
869     }						\
870   while (0)
871 
872 
873 /* Main multiplication routine.  The input values should be cooked.  */
874 
875 #define _FP_MUL(fs, wc, R, X, Y)				\
876   do								\
877     {								\
878       R##_s = X##_s ^ Y##_s;					\
879       R##_e = X##_e + Y##_e + 1;				\
880       switch (_FP_CLS_COMBINE (X##_c, Y##_c))			\
881 	{							\
882 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):	\
883 	  R##_c = FP_CLS_NORMAL;				\
884 								\
885 	  _FP_MUL_MEAT_##fs (R, X, Y);				\
886 								\
887 	  if (_FP_FRAC_OVERP_##wc (fs, R))			\
888 	    _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);	\
889 	  else							\
890 	    R##_e--;						\
891 	  break;						\
892 								\
893 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):		\
894 	  _FP_CHOOSENAN (fs, wc, R, X, Y, '*');			\
895 	  break;						\
896 								\
897 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):	\
898 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):		\
899 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):		\
900 	  R##_s = X##_s;					\
901 								\
902 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
903 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
904 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
905 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
906 	  _FP_FRAC_COPY_##wc (R, X);				\
907 	  R##_c = X##_c;					\
908 	  break;						\
909 								\
910 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
911 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
912 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
913 	  R##_s = Y##_s;					\
914 								\
915 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
916 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
917 	  _FP_FRAC_COPY_##wc (R, Y);				\
918 	  R##_c = Y##_c;					\
919 	  break;						\
920 								\
921 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
922 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
923 	  R##_s = _FP_NANSIGN_##fs;				\
924 	  R##_c = FP_CLS_NAN;					\
925 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
926 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ);	\
927 	  break;						\
928 								\
929 	default:						\
930 	  _FP_UNREACHABLE;					\
931 	}							\
932     }								\
933   while (0)
934 
935 
936 /* Fused multiply-add.  The input values should be cooked.  */
937 
938 #define _FP_FMA(fs, wc, dwc, R, X, Y, Z)				\
939   do									\
940     {									\
941       __label__ done_fma;						\
942       FP_DECL_##fs (_FP_FMA_T);						\
943       _FP_FMA_T##_s = X##_s ^ Y##_s;					\
944       _FP_FMA_T##_e = X##_e + Y##_e + 1;				\
945       switch (_FP_CLS_COMBINE (X##_c, Y##_c))				\
946 	{								\
947 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):		\
948 	  switch (Z##_c)						\
949 	    {								\
950 	    case FP_CLS_INF:						\
951 	    case FP_CLS_NAN:						\
952 	      R##_s = Z##_s;						\
953 	      _FP_FRAC_COPY_##wc (R, Z);				\
954 	      R##_c = Z##_c;						\
955 	      break;							\
956 									\
957 	    case FP_CLS_ZERO:						\
958 	      R##_c = FP_CLS_NORMAL;					\
959 	      R##_s = _FP_FMA_T##_s;					\
960 	      R##_e = _FP_FMA_T##_e;					\
961 									\
962 	      _FP_MUL_MEAT_##fs (R, X, Y);				\
963 									\
964 	      if (_FP_FRAC_OVERP_##wc (fs, R))				\
965 		_FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);		\
966 	      else							\
967 		R##_e--;						\
968 	      break;							\
969 									\
970 	    case FP_CLS_NORMAL:;					\
971 	      _FP_FRAC_DECL_##dwc (_FP_FMA_TD);				\
972 	      _FP_FRAC_DECL_##dwc (_FP_FMA_ZD);				\
973 	      _FP_FRAC_DECL_##dwc (_FP_FMA_RD);				\
974 	      _FP_MUL_MEAT_DW_##fs (_FP_FMA_TD, X, Y);			\
975 	      R##_e = _FP_FMA_T##_e;					\
976 	      int _FP_FMA_tsh						\
977 		= _FP_FRAC_HIGHBIT_DW_##dwc (fs, _FP_FMA_TD) == 0;	\
978 	      _FP_FMA_T##_e -= _FP_FMA_tsh;				\
979 	      int _FP_FMA_ediff = _FP_FMA_T##_e - Z##_e;		\
980 	      if (_FP_FMA_ediff >= 0)					\
981 		{							\
982 		  int _FP_FMA_shift					\
983 		    = _FP_WFRACBITS_##fs - _FP_FMA_tsh - _FP_FMA_ediff;	\
984 		  if (_FP_FMA_shift <= -_FP_WFRACBITS_##fs)		\
985 		    _FP_FRAC_SET_##dwc (_FP_FMA_ZD, _FP_MINFRAC_##dwc);	\
986 		  else							\
987 		    {							\
988 		      _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);	\
989 		      if (_FP_FMA_shift < 0)				\
990 			_FP_FRAC_SRS_##dwc (_FP_FMA_ZD, -_FP_FMA_shift,	\
991 					    _FP_WFRACBITS_DW_##fs);	\
992 		      else if (_FP_FMA_shift > 0)			\
993 			_FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_FMA_shift);	\
994 		    }							\
995 		  R##_s = _FP_FMA_T##_s;				\
996 		  if (_FP_FMA_T##_s == Z##_s)				\
997 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_TD,		\
998 					_FP_FMA_ZD);			\
999 		  else							\
1000 		    {							\
1001 		      _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_TD,	\
1002 					  _FP_FMA_ZD);			\
1003 		      if (_FP_FRAC_NEGP_##dwc (_FP_FMA_RD))		\
1004 			{						\
1005 			  R##_s = Z##_s;				\
1006 			  _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,	\
1007 					      _FP_FMA_TD);		\
1008 			}						\
1009 		    }							\
1010 		}							\
1011 	      else							\
1012 		{							\
1013 		  R##_e = Z##_e;					\
1014 		  R##_s = Z##_s;					\
1015 		  _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);		\
1016 		  _FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_WFRACBITS_##fs);	\
1017 		  int _FP_FMA_shift = -_FP_FMA_ediff - _FP_FMA_tsh;	\
1018 		  if (_FP_FMA_shift >= _FP_WFRACBITS_DW_##fs)		\
1019 		    _FP_FRAC_SET_##dwc (_FP_FMA_TD, _FP_MINFRAC_##dwc);	\
1020 		  else if (_FP_FMA_shift > 0)				\
1021 		    _FP_FRAC_SRS_##dwc (_FP_FMA_TD, _FP_FMA_shift,	\
1022 					_FP_WFRACBITS_DW_##fs);		\
1023 		  if (Z##_s == _FP_FMA_T##_s)				\
1024 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1025 					_FP_FMA_TD);			\
1026 		  else							\
1027 		    _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1028 					_FP_FMA_TD);			\
1029 		}							\
1030 	      if (_FP_FRAC_ZEROP_##dwc (_FP_FMA_RD))			\
1031 		{							\
1032 		  if (_FP_FMA_T##_s == Z##_s)				\
1033 		    R##_s = Z##_s;					\
1034 		  else							\
1035 		    R##_s = (FP_ROUNDMODE == FP_RND_MINF);		\
1036 		  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);		\
1037 		  R##_c = FP_CLS_ZERO;					\
1038 		}							\
1039 	      else							\
1040 		{							\
1041 		  int _FP_FMA_rlz;					\
1042 		  _FP_FRAC_CLZ_##dwc (_FP_FMA_rlz, _FP_FMA_RD);		\
1043 		  _FP_FMA_rlz -= _FP_WFRACXBITS_DW_##fs;		\
1044 		  R##_e -= _FP_FMA_rlz;					\
1045 		  int _FP_FMA_shift = _FP_WFRACBITS_##fs - _FP_FMA_rlz;	\
1046 		  if (_FP_FMA_shift > 0)				\
1047 		    _FP_FRAC_SRS_##dwc (_FP_FMA_RD, _FP_FMA_shift,	\
1048 					_FP_WFRACBITS_DW_##fs);		\
1049 		  else if (_FP_FMA_shift < 0)				\
1050 		    _FP_FRAC_SLL_##dwc (_FP_FMA_RD, -_FP_FMA_shift);	\
1051 		  _FP_FRAC_COPY_##wc##_##dwc (R, _FP_FMA_RD);		\
1052 		  R##_c = FP_CLS_NORMAL;				\
1053 		}							\
1054 	      break;							\
1055 	    }								\
1056 	  goto done_fma;						\
1057 									\
1058 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1059 	  _FP_CHOOSENAN (fs, wc, _FP_FMA_T, X, Y, '*');			\
1060 	  break;							\
1061 									\
1062 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1063 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1064 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1065 	  _FP_FMA_T##_s = X##_s;					\
1066 									\
1067 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1068 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1069 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1070 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1071 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, X);				\
1072 	  _FP_FMA_T##_c = X##_c;					\
1073 	  break;							\
1074 									\
1075 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):		\
1076 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1077 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1078 	  _FP_FMA_T##_s = Y##_s;					\
1079 									\
1080 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):		\
1081 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):		\
1082 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, Y);				\
1083 	  _FP_FMA_T##_c = Y##_c;					\
1084 	  break;							\
1085 									\
1086 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1087 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1088 	  _FP_FMA_T##_s = _FP_NANSIGN_##fs;				\
1089 	  _FP_FMA_T##_c = FP_CLS_NAN;					\
1090 	  _FP_FRAC_SET_##wc (_FP_FMA_T, _FP_NANFRAC_##fs);		\
1091 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ_FMA);	\
1092 	  break;							\
1093 									\
1094 	default:							\
1095 	  _FP_UNREACHABLE;						\
1096 	}								\
1097 									\
1098       /* T = X * Y is zero, infinity or NaN.  */			\
1099       switch (_FP_CLS_COMBINE (_FP_FMA_T##_c, Z##_c))			\
1100 	{								\
1101 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1102 	  _FP_CHOOSENAN (fs, wc, R, _FP_FMA_T, Z, '+');			\
1103 	  break;							\
1104 									\
1105 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1106 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1107 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1108 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1109 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1110 	  R##_s = _FP_FMA_T##_s;					\
1111 	  _FP_FRAC_COPY_##wc (R, _FP_FMA_T);				\
1112 	  R##_c = _FP_FMA_T##_c;					\
1113 	  break;							\
1114 									\
1115 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1116 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1117 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1118 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1119 	  R##_s = Z##_s;						\
1120 	  _FP_FRAC_COPY_##wc (R, Z);					\
1121 	  R##_c = Z##_c;						\
1122 	  R##_e = Z##_e;						\
1123 	  break;							\
1124 									\
1125 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1126 	  if (_FP_FMA_T##_s == Z##_s)					\
1127 	    {								\
1128 	      R##_s = Z##_s;						\
1129 	      _FP_FRAC_COPY_##wc (R, Z);				\
1130 	      R##_c = Z##_c;						\
1131 	    }								\
1132 	  else								\
1133 	    {								\
1134 	      R##_s = _FP_NANSIGN_##fs;					\
1135 	      R##_c = FP_CLS_NAN;					\
1136 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1137 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_ISI);	\
1138 	    }								\
1139 	  break;							\
1140 									\
1141 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1142 	  if (_FP_FMA_T##_s == Z##_s)					\
1143 	    R##_s = Z##_s;						\
1144 	  else								\
1145 	    R##_s = (FP_ROUNDMODE == FP_RND_MINF);			\
1146 	  _FP_FRAC_COPY_##wc (R, Z);					\
1147 	  R##_c = Z##_c;						\
1148 	  break;							\
1149 									\
1150 	default:							\
1151 	  _FP_UNREACHABLE;						\
1152 	}								\
1153     done_fma: ;								\
1154     }									\
1155   while (0)
1156 
1157 
1158 /* Main division routine.  The input values should be cooked.  */
1159 
1160 #define _FP_DIV(fs, wc, R, X, Y)				\
1161   do								\
1162     {								\
1163       R##_s = X##_s ^ Y##_s;					\
1164       R##_e = X##_e - Y##_e;					\
1165       switch (_FP_CLS_COMBINE (X##_c, Y##_c))			\
1166 	{							\
1167 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):	\
1168 	  R##_c = FP_CLS_NORMAL;				\
1169 								\
1170 	  _FP_DIV_MEAT_##fs (R, X, Y);				\
1171 	  break;						\
1172 								\
1173 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):		\
1174 	  _FP_CHOOSENAN (fs, wc, R, X, Y, '/');			\
1175 	  break;						\
1176 								\
1177 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):	\
1178 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):		\
1179 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):		\
1180 	  R##_s = X##_s;					\
1181 	  _FP_FRAC_COPY_##wc (R, X);				\
1182 	  R##_c = X##_c;					\
1183 	  break;						\
1184 								\
1185 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
1186 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
1187 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
1188 	  R##_s = Y##_s;					\
1189 	  _FP_FRAC_COPY_##wc (R, Y);				\
1190 	  R##_c = Y##_c;					\
1191 	  break;						\
1192 								\
1193 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
1194 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
1195 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
1196 	  R##_c = FP_CLS_ZERO;					\
1197 	  break;						\
1198 								\
1199 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
1200 	  FP_SET_EXCEPTION (FP_EX_DIVZERO);			\
1201 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
1202 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
1203 	  R##_c = FP_CLS_INF;					\
1204 	  break;						\
1205 								\
1206 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
1207 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
1208 	  R##_s = _FP_NANSIGN_##fs;				\
1209 	  R##_c = FP_CLS_NAN;					\
1210 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
1211 	  FP_SET_EXCEPTION (FP_EX_INVALID			\
1212 			    | (X##_c == FP_CLS_INF		\
1213 			       ? FP_EX_INVALID_IDI		\
1214 			       : FP_EX_INVALID_ZDZ));		\
1215 	  break;						\
1216 								\
1217 	default:						\
1218 	  _FP_UNREACHABLE;					\
1219 	}							\
1220     }								\
1221   while (0)
1222 
1223 
1224 /* Helper for comparisons.  EX is 0 not to raise exceptions, 1 to
1225    raise exceptions for signaling NaN operands, 2 to raise exceptions
1226    for all NaN operands.  Conditionals are organized to allow the
1227    compiler to optimize away code based on the value of EX.  */
1228 
1229 #define _FP_CMP_CHECK_NAN(fs, wc, X, Y, ex)				\
1230   do									\
1231     {									\
1232       /* The arguments are unordered, which may or may not result in	\
1233 	 an exception.  */						\
1234       if (ex)								\
1235 	{								\
1236 	  /* At least some cases of unordered arguments result in	\
1237 	     exceptions; check whether this is one.  */			\
1238 	  if (FP_EX_INVALID_SNAN || FP_EX_INVALID_VC)			\
1239 	    {								\
1240 	      /* Check separately for each case of "invalid"		\
1241 		 exceptions.  */					\
1242 	      if ((ex) == 2)						\
1243 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_VC);	\
1244 	      if (_FP_ISSIGNAN (fs, wc, X)				\
1245 		  || _FP_ISSIGNAN (fs, wc, Y))				\
1246 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SNAN);	\
1247 	    }								\
1248 	  /* Otherwise, we only need to check whether to raise an	\
1249 	     exception, not which case or cases it is.  */		\
1250 	  else if ((ex) == 2						\
1251 		   || _FP_ISSIGNAN (fs, wc, X)				\
1252 		   || _FP_ISSIGNAN (fs, wc, Y))				\
1253 	    FP_SET_EXCEPTION (FP_EX_INVALID);				\
1254 	}								\
1255     }									\
1256   while (0)
1257 
1258 /* Helper for comparisons.  If denormal operands would raise an
1259    exception, check for them, and flush to zero as appropriate
1260    (otherwise, we need only check and flush to zero if it might affect
1261    the result, which is done later with _FP_CMP_CHECK_FLUSH_ZERO).  */
1262 #define _FP_CMP_CHECK_DENORM(fs, wc, X, Y)				\
1263   do									\
1264     {									\
1265       if (FP_EX_DENORM != 0)						\
1266 	{								\
1267 	  /* We must ensure the correct exceptions are raised for	\
1268 	     denormal operands, even though this may not affect the	\
1269 	     result of the comparison.  */				\
1270 	  if (FP_DENORM_ZERO)						\
1271 	    {								\
1272 	      _FP_CHECK_FLUSH_ZERO (fs, wc, X);				\
1273 	      _FP_CHECK_FLUSH_ZERO (fs, wc, Y);				\
1274 	    }								\
1275 	  else								\
1276 	    {								\
1277 	      if ((X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X))		\
1278 		  || (Y##_e == 0 && !_FP_FRAC_ZEROP_##wc (Y)))		\
1279 		FP_SET_EXCEPTION (FP_EX_DENORM);			\
1280 	    }								\
1281 	}								\
1282     }									\
1283   while (0)
1284 
1285 /* Helper for comparisons.  Check for flushing denormals for zero if
1286    we didn't need to check earlier for any denormal operands.  */
1287 #define _FP_CMP_CHECK_FLUSH_ZERO(fs, wc, X, Y)	\
1288   do						\
1289     {						\
1290       if (FP_EX_DENORM == 0)			\
1291 	{					\
1292 	  _FP_CHECK_FLUSH_ZERO (fs, wc, X);	\
1293 	  _FP_CHECK_FLUSH_ZERO (fs, wc, Y);	\
1294 	}					\
1295     }						\
1296   while (0)
1297 
1298 /* Main differential comparison routine.  The inputs should be raw not
1299    cooked.  The return is -1, 0, 1 for normal values, UN
1300    otherwise.  */
1301 
1302 #define _FP_CMP(fs, wc, ret, X, Y, un, ex)				\
1303   do									\
1304     {									\
1305       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1306       /* NANs are unordered.  */					\
1307       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1308 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1309 	{								\
1310 	  (ret) = (un);							\
1311 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1312 	}								\
1313       else								\
1314 	{								\
1315 	  int _FP_CMP_is_zero_x;					\
1316 	  int _FP_CMP_is_zero_y;					\
1317 									\
1318 	  _FP_CMP_CHECK_FLUSH_ZERO (fs, wc, X, Y);			\
1319 									\
1320 	  _FP_CMP_is_zero_x						\
1321 	    = (!X##_e && _FP_FRAC_ZEROP_##wc (X)) ? 1 : 0;		\
1322 	  _FP_CMP_is_zero_y						\
1323 	    = (!Y##_e && _FP_FRAC_ZEROP_##wc (Y)) ? 1 : 0;		\
1324 									\
1325 	  if (_FP_CMP_is_zero_x && _FP_CMP_is_zero_y)			\
1326 	    (ret) = 0;							\
1327 	  else if (_FP_CMP_is_zero_x)					\
1328 	    (ret) = Y##_s ? 1 : -1;					\
1329 	  else if (_FP_CMP_is_zero_y)					\
1330 	    (ret) = X##_s ? -1 : 1;					\
1331 	  else if (X##_s != Y##_s)					\
1332 	    (ret) = X##_s ? -1 : 1;					\
1333 	  else if (X##_e > Y##_e)					\
1334 	    (ret) = X##_s ? -1 : 1;					\
1335 	  else if (X##_e < Y##_e)					\
1336 	    (ret) = X##_s ? 1 : -1;					\
1337 	  else if (_FP_FRAC_GT_##wc (X, Y))				\
1338 	    (ret) = X##_s ? -1 : 1;					\
1339 	  else if (_FP_FRAC_GT_##wc (Y, X))				\
1340 	    (ret) = X##_s ? 1 : -1;					\
1341 	  else								\
1342 	    (ret) = 0;							\
1343 	}								\
1344     }									\
1345   while (0)
1346 
1347 
1348 /* Simplification for strict equality.  */
1349 
1350 #define _FP_CMP_EQ(fs, wc, ret, X, Y, ex)				\
1351   do									\
1352     {									\
1353       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1354       /* NANs are unordered.  */					\
1355       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1356 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1357 	{								\
1358 	  (ret) = 1;							\
1359 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1360 	}								\
1361       else								\
1362 	{								\
1363 	  _FP_CMP_CHECK_FLUSH_ZERO (fs, wc, X, Y);			\
1364 									\
1365 	  (ret) = !(X##_e == Y##_e					\
1366 		    && _FP_FRAC_EQ_##wc (X, Y)				\
1367 		    && (X##_s == Y##_s					\
1368 			|| (!X##_e && _FP_FRAC_ZEROP_##wc (X))));	\
1369 	}								\
1370     }									\
1371   while (0)
1372 
1373 /* Version to test unordered.  */
1374 
1375 #define _FP_CMP_UNORD(fs, wc, ret, X, Y, ex)				\
1376   do									\
1377     {									\
1378       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1379       (ret) = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1380 	       || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))); \
1381       if (ret)								\
1382 	_FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));				\
1383     }									\
1384   while (0)
1385 
1386 /* Main square root routine.  The input value should be cooked.  */
1387 
1388 #define _FP_SQRT(fs, wc, R, X)						\
1389   do									\
1390     {									\
1391       _FP_FRAC_DECL_##wc (_FP_SQRT_T);					\
1392       _FP_FRAC_DECL_##wc (_FP_SQRT_S);					\
1393       _FP_W_TYPE _FP_SQRT_q;						\
1394       switch (X##_c)							\
1395 	{								\
1396 	case FP_CLS_NAN:						\
1397 	  _FP_FRAC_COPY_##wc (R, X);					\
1398 	  R##_s = X##_s;						\
1399 	  R##_c = FP_CLS_NAN;						\
1400 	  break;							\
1401 	case FP_CLS_INF:						\
1402 	  if (X##_s)							\
1403 	    {								\
1404 	      R##_s = _FP_NANSIGN_##fs;					\
1405 	      R##_c = FP_CLS_NAN; /* NAN */				\
1406 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1407 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1408 	    }								\
1409 	  else								\
1410 	    {								\
1411 	      R##_s = 0;						\
1412 	      R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */		\
1413 	    }								\
1414 	  break;							\
1415 	case FP_CLS_ZERO:						\
1416 	  R##_s = X##_s;						\
1417 	  R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */			\
1418 	  break;							\
1419 	case FP_CLS_NORMAL:						\
1420 	  R##_s = 0;							\
1421 	  if (X##_s)							\
1422 	    {								\
1423 	      R##_c = FP_CLS_NAN; /* NAN */				\
1424 	      R##_s = _FP_NANSIGN_##fs;					\
1425 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1426 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1427 	      break;							\
1428 	    }								\
1429 	  R##_c = FP_CLS_NORMAL;					\
1430 	  if (X##_e & 1)						\
1431 	    _FP_FRAC_SLL_##wc (X, 1);					\
1432 	  R##_e = X##_e >> 1;						\
1433 	  _FP_FRAC_SET_##wc (_FP_SQRT_S, _FP_ZEROFRAC_##wc);		\
1434 	  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);			\
1435 	  _FP_SQRT_q = _FP_OVERFLOW_##fs >> 1;				\
1436 	  _FP_SQRT_MEAT_##wc (R, _FP_SQRT_S, _FP_SQRT_T, X,		\
1437 			      _FP_SQRT_q);				\
1438 	}								\
1439     }									\
1440   while (0)
1441 
1442 /* Convert from FP to integer.  Input is raw.  */
1443 
1444 /* RSIGNED can have following values:
1445    0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1446        the result is either 0 or (2^rsize)-1 depending on the sign in such
1447        case.
1448    1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1449        NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1450        depending on the sign in such case.
1451    2:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1452        NV is set plus the result is reduced modulo 2^rsize.
1453    -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1454        set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1455        depending on the sign in such case.  */
1456 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)			\
1457   do									\
1458     {									\
1459       if (X##_e < _FP_EXPBIAS_##fs)					\
1460 	{								\
1461 	  (r) = 0;							\
1462 	  if (X##_e == 0)						\
1463 	    {								\
1464 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1465 		{							\
1466 		  if (!FP_DENORM_ZERO)					\
1467 		    FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1468 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1469 		}							\
1470 	    }								\
1471 	  else								\
1472 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1473 	}								\
1474       else if ((rsigned) == 2						\
1475 	       && (X##_e						\
1476 		   >= ((_FP_EXPMAX_##fs					\
1477 			< _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1) \
1478 		       ? _FP_EXPMAX_##fs				\
1479 		       : _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1))) \
1480 	{								\
1481 	  /* Overflow resulting in 0.  */				\
1482 	  (r) = 0;							\
1483 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1484 			    | FP_EX_INVALID_CVI				\
1485 			    | ((FP_EX_INVALID_SNAN			\
1486 				&& _FP_ISSIGNAN (fs, wc, X))		\
1487 			       ? FP_EX_INVALID_SNAN			\
1488 			       : 0));					\
1489 	}								\
1490       else if ((rsigned) != 2						\
1491 	       && (X##_e >= (_FP_EXPMAX_##fs < _FP_EXPBIAS_##fs + (rsize) \
1492 			     ? _FP_EXPMAX_##fs				\
1493 			     : (_FP_EXPBIAS_##fs + (rsize)		\
1494 				- ((rsigned) > 0 || X##_s)))		\
1495 		   || (!(rsigned) && X##_s)))				\
1496 	{								\
1497 	  /* Overflow or converting to the most negative integer.  */	\
1498 	  if (rsigned)							\
1499 	    {								\
1500 	      (r) = 1;							\
1501 	      (r) <<= (rsize) - 1;					\
1502 	      (r) -= 1 - X##_s;						\
1503 	    }								\
1504 	  else								\
1505 	    {								\
1506 	      (r) = 0;							\
1507 	      if (!X##_s)						\
1508 		(r) = ~(r);						\
1509 	    }								\
1510 									\
1511 	  if (_FP_EXPBIAS_##fs + (rsize) - 1 < _FP_EXPMAX_##fs		\
1512 	      && (rsigned)						\
1513 	      && X##_s							\
1514 	      && X##_e == _FP_EXPBIAS_##fs + (rsize) - 1)		\
1515 	    {								\
1516 	      /* Possibly converting to most negative integer; check the \
1517 		 mantissa.  */						\
1518 	      int _FP_TO_INT_inexact = 0;				\
1519 	      (void) ((_FP_FRACBITS_##fs > (rsize))			\
1520 		      ? ({						\
1521 			  _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,	\
1522 					      _FP_FRACBITS_##fs - (rsize), \
1523 					      _FP_FRACBITS_##fs);	\
1524 			  0;						\
1525 			})						\
1526 		      : 0);						\
1527 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1528 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1529 	      else if (_FP_TO_INT_inexact)				\
1530 		FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1531 	    }								\
1532 	  else								\
1533 	    FP_SET_EXCEPTION (FP_EX_INVALID				\
1534 			      | FP_EX_INVALID_CVI			\
1535 			      | ((FP_EX_INVALID_SNAN			\
1536 				  && _FP_ISSIGNAN (fs, wc, X))		\
1537 				 ? FP_EX_INVALID_SNAN			\
1538 				 : 0));					\
1539 	}								\
1540       else								\
1541 	{								\
1542 	  int _FP_TO_INT_inexact = 0;					\
1543 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;		\
1544 	  if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)	\
1545 	    {								\
1546 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1547 	      (r) <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1548 	    }								\
1549 	  else								\
1550 	    {								\
1551 	      _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,		\
1552 				  (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1553 				   - X##_e),				\
1554 				  _FP_FRACBITS_##fs);			\
1555 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1556 	    }								\
1557 	  if ((rsigned) && X##_s)					\
1558 	    (r) = -(r);							\
1559 	  if ((rsigned) == 2 && X##_e >= _FP_EXPBIAS_##fs + (rsize) - 1) \
1560 	    {								\
1561 	      /* Overflow or converting to the most negative integer.  */ \
1562 	      if (X##_e > _FP_EXPBIAS_##fs + (rsize) - 1		\
1563 		  || !X##_s						\
1564 		  || (r) != (((typeof (r)) 1) << ((rsize) - 1)))	\
1565 		{							\
1566 		  _FP_TO_INT_inexact = 0;				\
1567 		  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1568 		}							\
1569 	    }								\
1570 	  if (_FP_TO_INT_inexact)					\
1571 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1572 	}								\
1573     }									\
1574   while (0)
1575 
1576 /* Convert from floating point to integer, rounding according to the
1577    current rounding direction.  Input is raw.  RSIGNED is as for
1578    _FP_TO_INT.  */
1579 #define _FP_TO_INT_ROUND(fs, wc, r, X, rsize, rsigned)			\
1580   do									\
1581     {									\
1582       __label__ _FP_TO_INT_ROUND_done;					\
1583       if (X##_e < _FP_EXPBIAS_##fs)					\
1584 	{								\
1585 	  int _FP_TO_INT_ROUND_rounds_away = 0;				\
1586 	  if (X##_e == 0)						\
1587 	    {								\
1588 	      if (_FP_FRAC_ZEROP_##wc (X))				\
1589 		{							\
1590 		  (r) = 0;						\
1591 		  goto _FP_TO_INT_ROUND_done;				\
1592 		}							\
1593 	      else							\
1594 		{							\
1595 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1596 		  if (FP_DENORM_ZERO)					\
1597 		    {							\
1598 		      (r) = 0;						\
1599 		      goto _FP_TO_INT_ROUND_done;			\
1600 		    }							\
1601 		}							\
1602 	    }								\
1603 	  /* The result is 0, 1 or -1 depending on the rounding mode;	\
1604 	     -1 may cause overflow in the unsigned case.  */		\
1605 	  switch (FP_ROUNDMODE)						\
1606 	    {								\
1607 	    case FP_RND_NEAREST:					\
1608 	      _FP_TO_INT_ROUND_rounds_away				\
1609 		= (X##_e == _FP_EXPBIAS_##fs - 1			\
1610 		   && !_FP_FRAC_ZEROP_##wc (X));			\
1611 	      break;							\
1612 	    case FP_RND_ZERO:						\
1613 	      /* _FP_TO_INT_ROUND_rounds_away is already 0.  */		\
1614 	      break;							\
1615 	    case FP_RND_PINF:						\
1616 	      _FP_TO_INT_ROUND_rounds_away = !X##_s;			\
1617 	      break;							\
1618 	    case FP_RND_MINF:						\
1619 	      _FP_TO_INT_ROUND_rounds_away = X##_s;			\
1620 	      break;							\
1621 	    }								\
1622 	  if ((rsigned) == 0 && _FP_TO_INT_ROUND_rounds_away && X##_s)	\
1623 	    {								\
1624 	      /* Result of -1 for an unsigned conversion.  */		\
1625 	      (r) = 0;							\
1626 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1627 	    }								\
1628 	  else if ((rsize) == 1 && (rsigned) > 0			\
1629 		   && _FP_TO_INT_ROUND_rounds_away && !X##_s)		\
1630 	    {								\
1631 	      /* Converting to a 1-bit signed bit-field, which cannot	\
1632 		 represent +1.  */					\
1633 	      (r) = ((rsigned) == 2 ? -1 : 0);				\
1634 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1635 	    }								\
1636 	  else								\
1637 	    {								\
1638 	      (r) = (_FP_TO_INT_ROUND_rounds_away			\
1639 		     ? (X##_s ? -1 : 1)					\
1640 		     : 0);						\
1641 	      FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1642 	    }								\
1643 	}								\
1644       else if ((rsigned) == 2						\
1645 	       && (X##_e						\
1646 		   >= ((_FP_EXPMAX_##fs					\
1647 			< _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1) \
1648 		       ? _FP_EXPMAX_##fs				\
1649 		       : _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1))) \
1650 	{								\
1651 	  /* Overflow resulting in 0.  */				\
1652 	  (r) = 0;							\
1653 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1654 			    | FP_EX_INVALID_CVI				\
1655 			    | ((FP_EX_INVALID_SNAN			\
1656 				&& _FP_ISSIGNAN (fs, wc, X))		\
1657 			       ? FP_EX_INVALID_SNAN			\
1658 			       : 0));					\
1659 	}								\
1660       else if ((rsigned) != 2						\
1661 	       && (X##_e >= (_FP_EXPMAX_##fs < _FP_EXPBIAS_##fs + (rsize) \
1662 			     ? _FP_EXPMAX_##fs				\
1663 			     : (_FP_EXPBIAS_##fs + (rsize)		\
1664 				- ((rsigned) > 0 && !X##_s)))		\
1665 		   || ((rsigned) == 0 && X##_s)))			\
1666 	{								\
1667 	  /* Definite overflow (does not require rounding to tell).  */	\
1668 	  if ((rsigned) != 0)						\
1669 	    {								\
1670 	      (r) = 1;							\
1671 	      (r) <<= (rsize) - 1;					\
1672 	      (r) -= 1 - X##_s;						\
1673 	    }								\
1674 	  else								\
1675 	    {								\
1676 	      (r) = 0;							\
1677 	      if (!X##_s)						\
1678 		(r) = ~(r);						\
1679 	    }								\
1680 									\
1681 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1682 			    | FP_EX_INVALID_CVI				\
1683 			    | ((FP_EX_INVALID_SNAN			\
1684 				&& _FP_ISSIGNAN (fs, wc, X))		\
1685 			       ? FP_EX_INVALID_SNAN			\
1686 			       : 0));					\
1687 	}								\
1688       else								\
1689 	{								\
1690 	  /* The value is finite, with magnitude at least 1.  If	\
1691 	     the conversion is unsigned, the value is positive.		\
1692 	     If RSIGNED is not 2, the value does not definitely		\
1693 	     overflow by virtue of its exponent, but may still turn	\
1694 	     out to overflow after rounding; if RSIGNED is 2, the	\
1695 	     exponent may be such that the value definitely overflows,	\
1696 	     but at least one mantissa bit will not be shifted out.  */ \
1697 	  int _FP_TO_INT_ROUND_inexact = 0;				\
1698 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;		\
1699 	  if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)	\
1700 	    {								\
1701 	      /* The value is an integer, no rounding needed.  */	\
1702 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1703 	      (r) <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1704 	    }								\
1705 	  else								\
1706 	    {								\
1707 	      /* May need to shift in order to round (unless there	\
1708 		 are exactly _FP_WORKBITS fractional bits already).  */	\
1709 	      int _FP_TO_INT_ROUND_rshift				\
1710 		= (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs			\
1711 		   - 1 - _FP_WORKBITS - X##_e);				\
1712 	      if (_FP_TO_INT_ROUND_rshift > 0)				\
1713 		_FP_FRAC_SRS_##wc (X, _FP_TO_INT_ROUND_rshift,		\
1714 				   _FP_WFRACBITS_##fs);			\
1715 	      else if (_FP_TO_INT_ROUND_rshift < 0)			\
1716 		_FP_FRAC_SLL_##wc (X, -_FP_TO_INT_ROUND_rshift);	\
1717 	      /* Round like _FP_ROUND, but setting			\
1718 		 _FP_TO_INT_ROUND_inexact instead of directly setting	\
1719 		 the "inexact" exception, since it may turn out we	\
1720 		 should set "invalid" instead.  */			\
1721 	      if (_FP_FRAC_LOW_##wc (X) & 7)				\
1722 		{							\
1723 		  _FP_TO_INT_ROUND_inexact = 1;				\
1724 		  switch (FP_ROUNDMODE)					\
1725 		    {							\
1726 		    case FP_RND_NEAREST:				\
1727 		      _FP_ROUND_NEAREST (wc, X);			\
1728 		      break;						\
1729 		    case FP_RND_ZERO:					\
1730 		      _FP_ROUND_ZERO (wc, X);				\
1731 		      break;						\
1732 		    case FP_RND_PINF:					\
1733 		      _FP_ROUND_PINF (wc, X);				\
1734 		      break;						\
1735 		    case FP_RND_MINF:					\
1736 		      _FP_ROUND_MINF (wc, X);				\
1737 		      break;						\
1738 		    }							\
1739 		}							\
1740 	      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
1741 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1742 	    }								\
1743 	  if ((rsigned) != 0 && X##_s)					\
1744 	    (r) = -(r);							\
1745 	  /* An exponent of RSIZE - 1 always needs testing for		\
1746 	     overflow (either directly overflowing, or overflowing	\
1747 	     when rounding up results in 2^RSIZE).  An exponent of	\
1748 	     RSIZE - 2 can overflow for positive values when rounding	\
1749 	     up to 2^(RSIZE-1), but cannot overflow for negative	\
1750 	     values.  Smaller exponents cannot overflow.  */		\
1751 	  if (X##_e >= (_FP_EXPBIAS_##fs + (rsize) - 1			\
1752 			- ((rsigned) > 0 && !X##_s)))			\
1753 	    {								\
1754 	      if (X##_e > _FP_EXPBIAS_##fs + (rsize) - 1		\
1755 		  || (X##_e == _FP_EXPBIAS_##fs + (rsize) - 1		\
1756 		      && (X##_s						\
1757 			  ? (r) != (((typeof (r)) 1) << ((rsize) - 1))	\
1758 			  : ((rsigned) > 0 || (r) == 0)))		\
1759 		  || ((rsigned) > 0					\
1760 		      && !X##_s						\
1761 		      && X##_e == _FP_EXPBIAS_##fs + (rsize) - 2	\
1762 		      && (r) == (((typeof (r)) 1) << ((rsize) - 1))))	\
1763 		{							\
1764 		  if ((rsigned) != 2)					\
1765 		    {							\
1766 		      if ((rsigned) != 0)				\
1767 			{						\
1768 			  (r) = 1;					\
1769 			  (r) <<= (rsize) - 1;				\
1770 			  (r) -= 1 - X##_s;				\
1771 			}						\
1772 		      else						\
1773 			{						\
1774 			  (r) = 0;					\
1775 			  (r) = ~(r);					\
1776 			}						\
1777 		    }							\
1778 		  _FP_TO_INT_ROUND_inexact = 0;				\
1779 		  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1780 		}							\
1781 	    }								\
1782 	  if (_FP_TO_INT_ROUND_inexact)					\
1783 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1784 	}								\
1785     _FP_TO_INT_ROUND_done: ;						\
1786     }									\
1787   while (0)
1788 
1789 /* Convert integer to fp.  Output is raw.  RTYPE is unsigned even if
1790    input is signed.  */
1791 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)			\
1792   do									\
1793     {									\
1794       __label__ pack_semiraw;						\
1795       if (r)								\
1796 	{								\
1797 	  rtype _FP_FROM_INT_ur = (r);					\
1798 									\
1799 	  if ((X##_s = ((r) < 0)))					\
1800 	    _FP_FROM_INT_ur = -_FP_FROM_INT_ur;				\
1801 									\
1802 	  _FP_STATIC_ASSERT ((rsize) <= 2 * _FP_W_TYPE_SIZE,		\
1803 			     "rsize too large");			\
1804 	  (void) (((rsize) <= _FP_W_TYPE_SIZE)				\
1805 		  ? ({							\
1806 		      int _FP_FROM_INT_lz;				\
1807 		      __FP_CLZ (_FP_FROM_INT_lz,			\
1808 				(_FP_W_TYPE) _FP_FROM_INT_ur);		\
1809 		      X##_e = (_FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1	\
1810 			       - _FP_FROM_INT_lz);			\
1811 		    })							\
1812 		  : ({						\
1813 		      int _FP_FROM_INT_lz;				\
1814 		      __FP_CLZ_2 (_FP_FROM_INT_lz,			\
1815 				  (_FP_W_TYPE) (_FP_FROM_INT_ur		\
1816 						>> _FP_W_TYPE_SIZE),	\
1817 				  (_FP_W_TYPE) _FP_FROM_INT_ur);	\
1818 		      X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1819 			       - _FP_FROM_INT_lz);			\
1820 		    }));						\
1821 									\
1822 	  if ((rsize) - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs		\
1823 	      && X##_e >= _FP_EXPMAX_##fs)				\
1824 	    {								\
1825 	      /* Exponent too big; overflow to infinity.  (May also	\
1826 		 happen after rounding below.)  */			\
1827 	      _FP_OVERFLOW_SEMIRAW (fs, wc, X);				\
1828 	      goto pack_semiraw;					\
1829 	    }								\
1830 									\
1831 	  if ((rsize) <= _FP_FRACBITS_##fs				\
1832 	      || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs)		\
1833 	    {								\
1834 	      /* Exactly representable; shift left.  */			\
1835 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1836 	      if (_FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1 - X##_e > 0)	\
1837 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1838 				       + _FP_FRACBITS_##fs - 1 - X##_e)); \
1839 	    }								\
1840 	  else								\
1841 	    {								\
1842 	      /* More bits in integer than in floating type; need to	\
1843 		 round.  */						\
1844 	      if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e)	\
1845 		_FP_FROM_INT_ur						\
1846 		  = ((_FP_FROM_INT_ur >> (X##_e - _FP_EXPBIAS_##fs	\
1847 					  - _FP_WFRACBITS_##fs + 1))	\
1848 		     | ((_FP_FROM_INT_ur				\
1849 			 << ((rsize) - (X##_e - _FP_EXPBIAS_##fs	\
1850 					- _FP_WFRACBITS_##fs + 1)))	\
1851 			!= 0));						\
1852 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1853 	      if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1854 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1855 				       + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1856 	      _FP_FRAC_HIGH_##fs (X) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
1857 	    pack_semiraw:						\
1858 	      _FP_PACK_SEMIRAW (fs, wc, X);				\
1859 	    }								\
1860 	}								\
1861       else								\
1862 	{								\
1863 	  X##_s = 0;							\
1864 	  X##_e = 0;							\
1865 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
1866 	}								\
1867     }									\
1868   while (0)
1869 
1870 
1871 /* Extend from a narrower floating-point format to a wider one.  Input
1872    and output are raw.  If CHECK_NAN, then signaling NaNs are
1873    converted to quiet with the "invalid" exception raised; otherwise
1874    signaling NaNs remain signaling with no exception.  */
1875 #define _FP_EXTEND_CNAN(dfs, sfs, dwc, swc, D, S, check_nan)		\
1876   do									\
1877     {									\
1878       _FP_STATIC_ASSERT (_FP_FRACBITS_##dfs >= _FP_FRACBITS_##sfs,	\
1879 			 "destination mantissa narrower than source");	\
1880       _FP_STATIC_ASSERT ((_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs		\
1881 			  >= _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs),	\
1882 			 "destination max exponent smaller"		\
1883 			 " than source");				\
1884       _FP_STATIC_ASSERT (((_FP_EXPBIAS_##dfs				\
1885 			   >= (_FP_EXPBIAS_##sfs			\
1886 			       + _FP_FRACBITS_##sfs - 1))		\
1887 			  || (_FP_EXPBIAS_##dfs == _FP_EXPBIAS_##sfs)), \
1888 			 "source subnormals do not all become normal,"	\
1889 			 " but bias not the same");			\
1890       D##_s = S##_s;							\
1891       _FP_FRAC_COPY_##dwc##_##swc (D, S);				\
1892       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1893 	{								\
1894 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1895 	  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1896 	}								\
1897       else								\
1898 	{								\
1899 	  if (S##_e == 0)						\
1900 	    {								\
1901 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
1902 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1903 		D##_e = 0;						\
1904 	      else if (_FP_EXPBIAS_##dfs				\
1905 		       < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1)	\
1906 		{							\
1907 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1908 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1909 					  - _FP_FRACBITS_##sfs));	\
1910 		  D##_e = 0;						\
1911 		  if (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)		\
1912 		    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
1913 		}							\
1914 	      else							\
1915 		{							\
1916 		  int FP_EXTEND_lz;					\
1917 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1918 		  _FP_FRAC_CLZ_##swc (FP_EXTEND_lz, S);			\
1919 		  _FP_FRAC_SLL_##dwc (D,				\
1920 				      FP_EXTEND_lz + _FP_FRACBITS_##dfs	\
1921 				      - _FP_FRACTBITS_##sfs);		\
1922 		  D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1	\
1923 			   + _FP_FRACXBITS_##sfs - FP_EXTEND_lz);	\
1924 		}							\
1925 	    }								\
1926 	  else								\
1927 	    {								\
1928 	      D##_e = _FP_EXPMAX_##dfs;					\
1929 	      if (!_FP_FRAC_ZEROP_##swc (S))				\
1930 		{							\
1931 		  if (check_nan && _FP_FRAC_SNANP (sfs, S))		\
1932 		    FP_SET_EXCEPTION (FP_EX_INVALID			\
1933 				      | FP_EX_INVALID_SNAN);		\
1934 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1935 					  - _FP_FRACBITS_##sfs));	\
1936 		  if (check_nan)					\
1937 		    _FP_SETQNAN (dfs, dwc, D);				\
1938 		}							\
1939 	    }								\
1940 	}								\
1941     }									\
1942   while (0)
1943 
1944 #define FP_EXTEND(dfs, sfs, dwc, swc, D, S)		\
1945     _FP_EXTEND_CNAN (dfs, sfs, dwc, swc, D, S, 1)
1946 
1947 /* Truncate from a wider floating-point format to a narrower one.
1948    Input and output are semi-raw.  */
1949 #define FP_TRUNC(dfs, sfs, dwc, swc, D, S)				\
1950   do									\
1951     {									\
1952       _FP_STATIC_ASSERT (_FP_FRACBITS_##sfs >= _FP_FRACBITS_##dfs,	\
1953 			 "destination mantissa wider than source");	\
1954       _FP_STATIC_ASSERT (((_FP_EXPBIAS_##sfs				\
1955 			   >= (_FP_EXPBIAS_##dfs			\
1956 			       + _FP_FRACBITS_##dfs - 1))		\
1957 			  || _FP_EXPBIAS_##sfs == _FP_EXPBIAS_##dfs),	\
1958 			 "source subnormals do not all become same,"	\
1959 			 " but bias not the same");			\
1960       D##_s = S##_s;							\
1961       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1962 	{								\
1963 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1964 	  if (D##_e >= _FP_EXPMAX_##dfs)				\
1965 	    _FP_OVERFLOW_SEMIRAW (dfs, dwc, D);				\
1966 	  else								\
1967 	    {								\
1968 	      if (D##_e <= 0)						\
1969 		{							\
1970 		  if (D##_e < 1 - _FP_FRACBITS_##dfs)			\
1971 		    {							\
1972 		      _FP_FRAC_SET_##swc (S, _FP_ZEROFRAC_##swc);	\
1973 		      _FP_FRAC_LOW_##swc (S) |= 1;			\
1974 		    }							\
1975 		  else							\
1976 		    {							\
1977 		      _FP_FRAC_HIGH_##sfs (S) |= _FP_IMPLBIT_SH_##sfs;	\
1978 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
1979 					      - _FP_WFRACBITS_##dfs	\
1980 					      + 1 - D##_e),		\
1981 					  _FP_WFRACBITS_##sfs);		\
1982 		    }							\
1983 		  D##_e = 0;						\
1984 		}							\
1985 	      else							\
1986 		_FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs		\
1987 					- _FP_WFRACBITS_##dfs),		\
1988 				    _FP_WFRACBITS_##sfs);		\
1989 	      _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
1990 	    }								\
1991 	}								\
1992       else								\
1993 	{								\
1994 	  if (S##_e == 0)						\
1995 	    {								\
1996 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
1997 	      D##_e = 0;						\
1998 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1999 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
2000 	      else							\
2001 		{							\
2002 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
2003 		  if (_FP_EXPBIAS_##sfs					\
2004 		      < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1)	\
2005 		    {							\
2006 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
2007 					      - _FP_WFRACBITS_##dfs),	\
2008 					  _FP_WFRACBITS_##sfs);		\
2009 		      _FP_FRAC_COPY_##dwc##_##swc (D, S);		\
2010 		    }							\
2011 		  else							\
2012 		    {							\
2013 		      _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);	\
2014 		      _FP_FRAC_LOW_##dwc (D) |= 1;			\
2015 		    }							\
2016 		}							\
2017 	    }								\
2018 	  else								\
2019 	    {								\
2020 	      D##_e = _FP_EXPMAX_##dfs;					\
2021 	      if (_FP_FRAC_ZEROP_##swc (S))				\
2022 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
2023 	      else							\
2024 		{							\
2025 		  _FP_CHECK_SIGNAN_SEMIRAW (sfs, swc, S);		\
2026 		  _FP_FRAC_SRL_##swc (S, (_FP_WFRACBITS_##sfs		\
2027 					  - _FP_WFRACBITS_##dfs));	\
2028 		  _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
2029 		  /* Semi-raw NaN must have all workbits cleared.  */	\
2030 		  _FP_FRAC_LOW_##dwc (D)				\
2031 		    &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1);		\
2032 		  _FP_SETQNAN_SEMIRAW (dfs, dwc, D);			\
2033 		}							\
2034 	    }								\
2035 	}								\
2036     }									\
2037   while (0)
2038 
2039 /* Helper primitives.  */
2040 
2041 /* Count leading zeros in a word.  */
2042 
2043 #ifndef __FP_CLZ
2044 /* GCC 3.4 and later provide the builtins for us.  */
2045 # define __FP_CLZ(r, x)							\
2046   do									\
2047     {									\
2048       _FP_STATIC_ASSERT ((sizeof (_FP_W_TYPE) == sizeof (unsigned int)	\
2049 			  || (sizeof (_FP_W_TYPE)			\
2050 			      == sizeof (unsigned long))		\
2051 			  || (sizeof (_FP_W_TYPE)			\
2052 			      == sizeof (unsigned long long))),		\
2053 			 "_FP_W_TYPE size unsupported for clz");	\
2054       if (sizeof (_FP_W_TYPE) == sizeof (unsigned int))			\
2055 	(r) = __builtin_clz (x);					\
2056       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long))		\
2057 	(r) = __builtin_clzl (x);					\
2058       else /* sizeof (_FP_W_TYPE) == sizeof (unsigned long long).  */	\
2059 	(r) = __builtin_clzll (x);					\
2060     }									\
2061   while (0)
2062 #endif /* ndef __FP_CLZ */
2063 
2064 #define _FP_DIV_HELP_imm(q, r, n, d)		\
2065   do						\
2066     {						\
2067       (q) = (n) / (d), (r) = (n) % (d);		\
2068     }						\
2069   while (0)
2070 
2071 
2072 /* A restoring bit-by-bit division primitive.  */
2073 
2074 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y)				\
2075   do									\
2076     {									\
2077       int _FP_DIV_MEAT_N_loop_count = _FP_WFRACBITS_##fs;		\
2078       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_u);			\
2079       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_v);			\
2080       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_u, X);			\
2081       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_v, Y);			\
2082       _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);				\
2083       /* Normalize _FP_DIV_MEAT_N_LOOP_U and _FP_DIV_MEAT_N_LOOP_V.  */	\
2084       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, _FP_WFRACXBITS_##fs);	\
2085       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_v, _FP_WFRACXBITS_##fs);	\
2086       /* First round.  Since the operands are normalized, either the	\
2087 	 first or second bit will be set in the fraction.  Produce a	\
2088 	 normalized result by checking which and adjusting the loop	\
2089 	 count and exponent accordingly.  */				\
2090       if (_FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u, _FP_DIV_MEAT_N_loop_v))	\
2091 	{								\
2092 	  _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
2093 			     _FP_DIV_MEAT_N_loop_u,			\
2094 			     _FP_DIV_MEAT_N_loop_v);			\
2095 	  _FP_FRAC_LOW_##wc (R) |= 1;					\
2096 	  _FP_DIV_MEAT_N_loop_count--;					\
2097 	}								\
2098       else								\
2099 	R##_e--;							\
2100       /* Subsequent rounds.  */						\
2101       do								\
2102 	{								\
2103 	  int _FP_DIV_MEAT_N_loop_msb					\
2104 	    = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (_FP_DIV_MEAT_N_loop_u) < 0; \
2105 	  _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, 1);			\
2106 	  _FP_FRAC_SLL_##wc (R, 1);					\
2107 	  if (_FP_DIV_MEAT_N_loop_msb					\
2108 	      || _FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u,			\
2109 				_FP_DIV_MEAT_N_loop_v))			\
2110 	    {								\
2111 	      _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
2112 				 _FP_DIV_MEAT_N_loop_u,			\
2113 				 _FP_DIV_MEAT_N_loop_v);		\
2114 	      _FP_FRAC_LOW_##wc (R) |= 1;				\
2115 	    }								\
2116 	}								\
2117       while (--_FP_DIV_MEAT_N_loop_count > 0);				\
2118       /* If there's anything left in _FP_DIV_MEAT_N_LOOP_U, the result	\
2119 	 is inexact.  */						\
2120       _FP_FRAC_LOW_##wc (R)						\
2121 	|= !_FP_FRAC_ZEROP_##wc (_FP_DIV_MEAT_N_loop_u);		\
2122     }									\
2123   while (0)
2124 
2125 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
2126 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
2127 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)
2128 
2129 #endif /* !SOFT_FP_OP_COMMON_H */
2130