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