1 /* Software floating-point emulation. Common operations.
2    Copyright (C) 1997-2017 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 	  /* FALLTHRU */					\
902 								\
903 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
904 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
905 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
906 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
907 	  _FP_FRAC_COPY_##wc (R, X);				\
908 	  R##_c = X##_c;					\
909 	  break;						\
910 								\
911 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
912 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
913 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
914 	  R##_s = Y##_s;					\
915 	  /* FALLTHRU */					\
916 								\
917 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
918 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
919 	  _FP_FRAC_COPY_##wc (R, Y);				\
920 	  R##_c = Y##_c;					\
921 	  break;						\
922 								\
923 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
924 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
925 	  R##_s = _FP_NANSIGN_##fs;				\
926 	  R##_c = FP_CLS_NAN;					\
927 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
928 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ);	\
929 	  break;						\
930 								\
931 	default:						\
932 	  _FP_UNREACHABLE;					\
933 	}							\
934     }								\
935   while (0)
936 
937 
938 /* Fused multiply-add.  The input values should be cooked.  */
939 
940 #define _FP_FMA(fs, wc, dwc, R, X, Y, Z)				\
941   do									\
942     {									\
943       __label__ done_fma;						\
944       FP_DECL_##fs (_FP_FMA_T);						\
945       _FP_FMA_T##_s = X##_s ^ Y##_s;					\
946       _FP_FMA_T##_e = X##_e + Y##_e + 1;				\
947       switch (_FP_CLS_COMBINE (X##_c, Y##_c))				\
948 	{								\
949 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):		\
950 	  switch (Z##_c)						\
951 	    {								\
952 	    case FP_CLS_INF:						\
953 	    case FP_CLS_NAN:						\
954 	      R##_s = Z##_s;						\
955 	      _FP_FRAC_COPY_##wc (R, Z);				\
956 	      R##_c = Z##_c;						\
957 	      break;							\
958 									\
959 	    case FP_CLS_ZERO:						\
960 	      R##_c = FP_CLS_NORMAL;					\
961 	      R##_s = _FP_FMA_T##_s;					\
962 	      R##_e = _FP_FMA_T##_e;					\
963 									\
964 	      _FP_MUL_MEAT_##fs (R, X, Y);				\
965 									\
966 	      if (_FP_FRAC_OVERP_##wc (fs, R))				\
967 		_FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);		\
968 	      else							\
969 		R##_e--;						\
970 	      break;							\
971 									\
972 	    case FP_CLS_NORMAL:;					\
973 	      _FP_FRAC_DECL_##dwc (_FP_FMA_TD);				\
974 	      _FP_FRAC_DECL_##dwc (_FP_FMA_ZD);				\
975 	      _FP_FRAC_DECL_##dwc (_FP_FMA_RD);				\
976 	      _FP_MUL_MEAT_DW_##fs (_FP_FMA_TD, X, Y);			\
977 	      R##_e = _FP_FMA_T##_e;					\
978 	      int _FP_FMA_tsh						\
979 		= _FP_FRAC_HIGHBIT_DW_##dwc (fs, _FP_FMA_TD) == 0;	\
980 	      _FP_FMA_T##_e -= _FP_FMA_tsh;				\
981 	      int _FP_FMA_ediff = _FP_FMA_T##_e - Z##_e;		\
982 	      if (_FP_FMA_ediff >= 0)					\
983 		{							\
984 		  int _FP_FMA_shift					\
985 		    = _FP_WFRACBITS_##fs - _FP_FMA_tsh - _FP_FMA_ediff;	\
986 		  if (_FP_FMA_shift <= -_FP_WFRACBITS_##fs)		\
987 		    _FP_FRAC_SET_##dwc (_FP_FMA_ZD, _FP_MINFRAC_##dwc);	\
988 		  else							\
989 		    {							\
990 		      _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);	\
991 		      if (_FP_FMA_shift < 0)				\
992 			_FP_FRAC_SRS_##dwc (_FP_FMA_ZD, -_FP_FMA_shift,	\
993 					    _FP_WFRACBITS_DW_##fs);	\
994 		      else if (_FP_FMA_shift > 0)			\
995 			_FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_FMA_shift);	\
996 		    }							\
997 		  R##_s = _FP_FMA_T##_s;				\
998 		  if (_FP_FMA_T##_s == Z##_s)				\
999 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_TD,		\
1000 					_FP_FMA_ZD);			\
1001 		  else							\
1002 		    {							\
1003 		      _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_TD,	\
1004 					  _FP_FMA_ZD);			\
1005 		      if (_FP_FRAC_NEGP_##dwc (_FP_FMA_RD))		\
1006 			{						\
1007 			  R##_s = Z##_s;				\
1008 			  _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,	\
1009 					      _FP_FMA_TD);		\
1010 			}						\
1011 		    }							\
1012 		}							\
1013 	      else							\
1014 		{							\
1015 		  R##_e = Z##_e;					\
1016 		  R##_s = Z##_s;					\
1017 		  _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);		\
1018 		  _FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_WFRACBITS_##fs);	\
1019 		  int _FP_FMA_shift = -_FP_FMA_ediff - _FP_FMA_tsh;	\
1020 		  if (_FP_FMA_shift >= _FP_WFRACBITS_DW_##fs)		\
1021 		    _FP_FRAC_SET_##dwc (_FP_FMA_TD, _FP_MINFRAC_##dwc);	\
1022 		  else if (_FP_FMA_shift > 0)				\
1023 		    _FP_FRAC_SRS_##dwc (_FP_FMA_TD, _FP_FMA_shift,	\
1024 					_FP_WFRACBITS_DW_##fs);		\
1025 		  if (Z##_s == _FP_FMA_T##_s)				\
1026 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1027 					_FP_FMA_TD);			\
1028 		  else							\
1029 		    _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1030 					_FP_FMA_TD);			\
1031 		}							\
1032 	      if (_FP_FRAC_ZEROP_##dwc (_FP_FMA_RD))			\
1033 		{							\
1034 		  if (_FP_FMA_T##_s == Z##_s)				\
1035 		    R##_s = Z##_s;					\
1036 		  else							\
1037 		    R##_s = (FP_ROUNDMODE == FP_RND_MINF);		\
1038 		  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);		\
1039 		  R##_c = FP_CLS_ZERO;					\
1040 		}							\
1041 	      else							\
1042 		{							\
1043 		  int _FP_FMA_rlz;					\
1044 		  _FP_FRAC_CLZ_##dwc (_FP_FMA_rlz, _FP_FMA_RD);		\
1045 		  _FP_FMA_rlz -= _FP_WFRACXBITS_DW_##fs;		\
1046 		  R##_e -= _FP_FMA_rlz;					\
1047 		  int _FP_FMA_shift = _FP_WFRACBITS_##fs - _FP_FMA_rlz;	\
1048 		  if (_FP_FMA_shift > 0)				\
1049 		    _FP_FRAC_SRS_##dwc (_FP_FMA_RD, _FP_FMA_shift,	\
1050 					_FP_WFRACBITS_DW_##fs);		\
1051 		  else if (_FP_FMA_shift < 0)				\
1052 		    _FP_FRAC_SLL_##dwc (_FP_FMA_RD, -_FP_FMA_shift);	\
1053 		  _FP_FRAC_COPY_##wc##_##dwc (R, _FP_FMA_RD);		\
1054 		  R##_c = FP_CLS_NORMAL;				\
1055 		}							\
1056 	      break;							\
1057 	    }								\
1058 	  goto done_fma;						\
1059 									\
1060 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1061 	  _FP_CHOOSENAN (fs, wc, _FP_FMA_T, X, Y, '*');			\
1062 	  break;							\
1063 									\
1064 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1065 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1066 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1067 	  _FP_FMA_T##_s = X##_s;					\
1068 	  /* FALLTHRU */						\
1069 									\
1070 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1071 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1072 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1073 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1074 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, X);				\
1075 	  _FP_FMA_T##_c = X##_c;					\
1076 	  break;							\
1077 									\
1078 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):		\
1079 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1080 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1081 	  _FP_FMA_T##_s = Y##_s;					\
1082 	  /* FALLTHRU */						\
1083 									\
1084 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):		\
1085 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):		\
1086 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, Y);				\
1087 	  _FP_FMA_T##_c = Y##_c;					\
1088 	  break;							\
1089 									\
1090 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1091 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1092 	  _FP_FMA_T##_s = _FP_NANSIGN_##fs;				\
1093 	  _FP_FMA_T##_c = FP_CLS_NAN;					\
1094 	  _FP_FRAC_SET_##wc (_FP_FMA_T, _FP_NANFRAC_##fs);		\
1095 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ_FMA);	\
1096 	  break;							\
1097 									\
1098 	default:							\
1099 	  _FP_UNREACHABLE;						\
1100 	}								\
1101 									\
1102       /* T = X * Y is zero, infinity or NaN.  */			\
1103       switch (_FP_CLS_COMBINE (_FP_FMA_T##_c, Z##_c))			\
1104 	{								\
1105 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1106 	  _FP_CHOOSENAN (fs, wc, R, _FP_FMA_T, Z, '+');			\
1107 	  break;							\
1108 									\
1109 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1110 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1111 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1112 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1113 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1114 	  R##_s = _FP_FMA_T##_s;					\
1115 	  _FP_FRAC_COPY_##wc (R, _FP_FMA_T);				\
1116 	  R##_c = _FP_FMA_T##_c;					\
1117 	  break;							\
1118 									\
1119 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1120 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1121 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1122 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1123 	  R##_s = Z##_s;						\
1124 	  _FP_FRAC_COPY_##wc (R, Z);					\
1125 	  R##_c = Z##_c;						\
1126 	  R##_e = Z##_e;						\
1127 	  break;							\
1128 									\
1129 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1130 	  if (_FP_FMA_T##_s == Z##_s)					\
1131 	    {								\
1132 	      R##_s = Z##_s;						\
1133 	      _FP_FRAC_COPY_##wc (R, Z);				\
1134 	      R##_c = Z##_c;						\
1135 	    }								\
1136 	  else								\
1137 	    {								\
1138 	      R##_s = _FP_NANSIGN_##fs;					\
1139 	      R##_c = FP_CLS_NAN;					\
1140 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1141 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_ISI);	\
1142 	    }								\
1143 	  break;							\
1144 									\
1145 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1146 	  if (_FP_FMA_T##_s == Z##_s)					\
1147 	    R##_s = Z##_s;						\
1148 	  else								\
1149 	    R##_s = (FP_ROUNDMODE == FP_RND_MINF);			\
1150 	  _FP_FRAC_COPY_##wc (R, Z);					\
1151 	  R##_c = Z##_c;						\
1152 	  break;							\
1153 									\
1154 	default:							\
1155 	  _FP_UNREACHABLE;						\
1156 	}								\
1157     done_fma: ;								\
1158     }									\
1159   while (0)
1160 
1161 
1162 /* Main division routine.  The input values should be cooked.  */
1163 
1164 #define _FP_DIV(fs, wc, R, X, Y)				\
1165   do								\
1166     {								\
1167       R##_s = X##_s ^ Y##_s;					\
1168       R##_e = X##_e - Y##_e;					\
1169       switch (_FP_CLS_COMBINE (X##_c, Y##_c))			\
1170 	{							\
1171 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):	\
1172 	  R##_c = FP_CLS_NORMAL;				\
1173 								\
1174 	  _FP_DIV_MEAT_##fs (R, X, Y);				\
1175 	  break;						\
1176 								\
1177 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):		\
1178 	  _FP_CHOOSENAN (fs, wc, R, X, Y, '/');			\
1179 	  break;						\
1180 								\
1181 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):	\
1182 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):		\
1183 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):		\
1184 	  R##_s = X##_s;					\
1185 	  _FP_FRAC_COPY_##wc (R, X);				\
1186 	  R##_c = X##_c;					\
1187 	  break;						\
1188 								\
1189 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
1190 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
1191 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
1192 	  R##_s = Y##_s;					\
1193 	  _FP_FRAC_COPY_##wc (R, Y);				\
1194 	  R##_c = Y##_c;					\
1195 	  break;						\
1196 								\
1197 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
1198 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
1199 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
1200 	  R##_c = FP_CLS_ZERO;					\
1201 	  break;						\
1202 								\
1203 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
1204 	  FP_SET_EXCEPTION (FP_EX_DIVZERO);			\
1205 	  /* FALLTHRU */					\
1206 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
1207 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
1208 	  R##_c = FP_CLS_INF;					\
1209 	  break;						\
1210 								\
1211 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
1212 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
1213 	  R##_s = _FP_NANSIGN_##fs;				\
1214 	  R##_c = FP_CLS_NAN;					\
1215 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
1216 	  FP_SET_EXCEPTION (FP_EX_INVALID			\
1217 			    | (X##_c == FP_CLS_INF		\
1218 			       ? FP_EX_INVALID_IDI		\
1219 			       : FP_EX_INVALID_ZDZ));		\
1220 	  break;						\
1221 								\
1222 	default:						\
1223 	  _FP_UNREACHABLE;					\
1224 	}							\
1225     }								\
1226   while (0)
1227 
1228 
1229 /* Helper for comparisons.  EX is 0 not to raise exceptions, 1 to
1230    raise exceptions for signaling NaN operands, 2 to raise exceptions
1231    for all NaN operands.  Conditionals are organized to allow the
1232    compiler to optimize away code based on the value of EX.  */
1233 
1234 #define _FP_CMP_CHECK_NAN(fs, wc, X, Y, ex)				\
1235   do									\
1236     {									\
1237       /* The arguments are unordered, which may or may not result in	\
1238 	 an exception.  */						\
1239       if (ex)								\
1240 	{								\
1241 	  /* At least some cases of unordered arguments result in	\
1242 	     exceptions; check whether this is one.  */			\
1243 	  if (FP_EX_INVALID_SNAN || FP_EX_INVALID_VC)			\
1244 	    {								\
1245 	      /* Check separately for each case of "invalid"		\
1246 		 exceptions.  */					\
1247 	      if ((ex) == 2)						\
1248 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_VC);	\
1249 	      if (_FP_ISSIGNAN (fs, wc, X)				\
1250 		  || _FP_ISSIGNAN (fs, wc, Y))				\
1251 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SNAN);	\
1252 	    }								\
1253 	  /* Otherwise, we only need to check whether to raise an	\
1254 	     exception, not which case or cases it is.  */		\
1255 	  else if ((ex) == 2						\
1256 		   || _FP_ISSIGNAN (fs, wc, X)				\
1257 		   || _FP_ISSIGNAN (fs, wc, Y))				\
1258 	    FP_SET_EXCEPTION (FP_EX_INVALID);				\
1259 	}								\
1260     }									\
1261   while (0)
1262 
1263 /* Helper for comparisons.  If denormal operands would raise an
1264    exception, check for them, and flush to zero as appropriate
1265    (otherwise, we need only check and flush to zero if it might affect
1266    the result, which is done later with _FP_CMP_CHECK_FLUSH_ZERO).  */
1267 #define _FP_CMP_CHECK_DENORM(fs, wc, X, Y)				\
1268   do									\
1269     {									\
1270       if (FP_EX_DENORM != 0)						\
1271 	{								\
1272 	  /* We must ensure the correct exceptions are raised for	\
1273 	     denormal operands, even though this may not affect the	\
1274 	     result of the comparison.  */				\
1275 	  if (FP_DENORM_ZERO)						\
1276 	    {								\
1277 	      _FP_CHECK_FLUSH_ZERO (fs, wc, X);				\
1278 	      _FP_CHECK_FLUSH_ZERO (fs, wc, Y);				\
1279 	    }								\
1280 	  else								\
1281 	    {								\
1282 	      if ((X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X))		\
1283 		  || (Y##_e == 0 && !_FP_FRAC_ZEROP_##wc (Y)))		\
1284 		FP_SET_EXCEPTION (FP_EX_DENORM);			\
1285 	    }								\
1286 	}								\
1287     }									\
1288   while (0)
1289 
1290 /* Helper for comparisons.  Check for flushing denormals for zero if
1291    we didn't need to check earlier for any denormal operands.  */
1292 #define _FP_CMP_CHECK_FLUSH_ZERO(fs, wc, X, Y)	\
1293   do						\
1294     {						\
1295       if (FP_EX_DENORM == 0)			\
1296 	{					\
1297 	  _FP_CHECK_FLUSH_ZERO (fs, wc, X);	\
1298 	  _FP_CHECK_FLUSH_ZERO (fs, wc, Y);	\
1299 	}					\
1300     }						\
1301   while (0)
1302 
1303 /* Main differential comparison routine.  The inputs should be raw not
1304    cooked.  The return is -1, 0, 1 for normal values, UN
1305    otherwise.  */
1306 
1307 #define _FP_CMP(fs, wc, ret, X, Y, un, ex)				\
1308   do									\
1309     {									\
1310       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1311       /* NANs are unordered.  */					\
1312       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1313 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1314 	{								\
1315 	  (ret) = (un);							\
1316 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1317 	}								\
1318       else								\
1319 	{								\
1320 	  int _FP_CMP_is_zero_x;					\
1321 	  int _FP_CMP_is_zero_y;					\
1322 									\
1323 	  _FP_CMP_CHECK_FLUSH_ZERO (fs, wc, X, Y);			\
1324 									\
1325 	  _FP_CMP_is_zero_x						\
1326 	    = (!X##_e && _FP_FRAC_ZEROP_##wc (X)) ? 1 : 0;		\
1327 	  _FP_CMP_is_zero_y						\
1328 	    = (!Y##_e && _FP_FRAC_ZEROP_##wc (Y)) ? 1 : 0;		\
1329 									\
1330 	  if (_FP_CMP_is_zero_x && _FP_CMP_is_zero_y)			\
1331 	    (ret) = 0;							\
1332 	  else if (_FP_CMP_is_zero_x)					\
1333 	    (ret) = Y##_s ? 1 : -1;					\
1334 	  else if (_FP_CMP_is_zero_y)					\
1335 	    (ret) = X##_s ? -1 : 1;					\
1336 	  else if (X##_s != Y##_s)					\
1337 	    (ret) = X##_s ? -1 : 1;					\
1338 	  else if (X##_e > Y##_e)					\
1339 	    (ret) = X##_s ? -1 : 1;					\
1340 	  else if (X##_e < Y##_e)					\
1341 	    (ret) = X##_s ? 1 : -1;					\
1342 	  else if (_FP_FRAC_GT_##wc (X, Y))				\
1343 	    (ret) = X##_s ? -1 : 1;					\
1344 	  else if (_FP_FRAC_GT_##wc (Y, X))				\
1345 	    (ret) = X##_s ? 1 : -1;					\
1346 	  else								\
1347 	    (ret) = 0;							\
1348 	}								\
1349     }									\
1350   while (0)
1351 
1352 
1353 /* Simplification for strict equality.  */
1354 
1355 #define _FP_CMP_EQ(fs, wc, ret, X, Y, ex)				\
1356   do									\
1357     {									\
1358       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1359       /* NANs are unordered.  */					\
1360       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1361 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1362 	{								\
1363 	  (ret) = 1;							\
1364 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1365 	}								\
1366       else								\
1367 	{								\
1368 	  _FP_CMP_CHECK_FLUSH_ZERO (fs, wc, X, Y);			\
1369 									\
1370 	  (ret) = !(X##_e == Y##_e					\
1371 		    && _FP_FRAC_EQ_##wc (X, Y)				\
1372 		    && (X##_s == Y##_s					\
1373 			|| (!X##_e && _FP_FRAC_ZEROP_##wc (X))));	\
1374 	}								\
1375     }									\
1376   while (0)
1377 
1378 /* Version to test unordered.  */
1379 
1380 #define _FP_CMP_UNORD(fs, wc, ret, X, Y, ex)				\
1381   do									\
1382     {									\
1383       _FP_CMP_CHECK_DENORM (fs, wc, X, Y);				\
1384       (ret) = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1385 	       || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))); \
1386       if (ret)								\
1387 	_FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));				\
1388     }									\
1389   while (0)
1390 
1391 /* Main square root routine.  The input value should be cooked.  */
1392 
1393 #define _FP_SQRT(fs, wc, R, X)						\
1394   do									\
1395     {									\
1396       _FP_FRAC_DECL_##wc (_FP_SQRT_T);					\
1397       _FP_FRAC_DECL_##wc (_FP_SQRT_S);					\
1398       _FP_W_TYPE _FP_SQRT_q;						\
1399       switch (X##_c)							\
1400 	{								\
1401 	case FP_CLS_NAN:						\
1402 	  _FP_FRAC_COPY_##wc (R, X);					\
1403 	  R##_s = X##_s;						\
1404 	  R##_c = FP_CLS_NAN;						\
1405 	  break;							\
1406 	case FP_CLS_INF:						\
1407 	  if (X##_s)							\
1408 	    {								\
1409 	      R##_s = _FP_NANSIGN_##fs;					\
1410 	      R##_c = FP_CLS_NAN; /* NAN */				\
1411 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1412 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1413 	    }								\
1414 	  else								\
1415 	    {								\
1416 	      R##_s = 0;						\
1417 	      R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */		\
1418 	    }								\
1419 	  break;							\
1420 	case FP_CLS_ZERO:						\
1421 	  R##_s = X##_s;						\
1422 	  R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */			\
1423 	  break;							\
1424 	case FP_CLS_NORMAL:						\
1425 	  R##_s = 0;							\
1426 	  if (X##_s)							\
1427 	    {								\
1428 	      R##_c = FP_CLS_NAN; /* NAN */				\
1429 	      R##_s = _FP_NANSIGN_##fs;					\
1430 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1431 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1432 	      break;							\
1433 	    }								\
1434 	  R##_c = FP_CLS_NORMAL;					\
1435 	  if (X##_e & 1)						\
1436 	    _FP_FRAC_SLL_##wc (X, 1);					\
1437 	  R##_e = X##_e >> 1;						\
1438 	  _FP_FRAC_SET_##wc (_FP_SQRT_S, _FP_ZEROFRAC_##wc);		\
1439 	  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);			\
1440 	  _FP_SQRT_q = _FP_OVERFLOW_##fs >> 1;				\
1441 	  _FP_SQRT_MEAT_##wc (R, _FP_SQRT_S, _FP_SQRT_T, X,		\
1442 			      _FP_SQRT_q);				\
1443 	}								\
1444     }									\
1445   while (0)
1446 
1447 /* Convert from FP to integer.  Input is raw.  */
1448 
1449 /* RSIGNED can have following values:
1450    0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1451        the result is either 0 or (2^rsize)-1 depending on the sign in such
1452        case.
1453    1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1454        NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1455        depending on the sign in such case.
1456    2:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1457        NV is set plus the result is reduced modulo 2^rsize.
1458    -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1459        set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1460        depending on the sign in such case.  */
1461 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)			\
1462   do									\
1463     {									\
1464       if (X##_e < _FP_EXPBIAS_##fs)					\
1465 	{								\
1466 	  (r) = 0;							\
1467 	  if (X##_e == 0)						\
1468 	    {								\
1469 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1470 		{							\
1471 		  if (!FP_DENORM_ZERO)					\
1472 		    FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1473 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1474 		}							\
1475 	    }								\
1476 	  else								\
1477 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1478 	}								\
1479       else if ((rsigned) == 2						\
1480 	       && (X##_e						\
1481 		   >= ((_FP_EXPMAX_##fs					\
1482 			< _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1) \
1483 		       ? _FP_EXPMAX_##fs				\
1484 		       : _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1))) \
1485 	{								\
1486 	  /* Overflow resulting in 0.  */				\
1487 	  (r) = 0;							\
1488 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1489 			    | FP_EX_INVALID_CVI				\
1490 			    | ((FP_EX_INVALID_SNAN			\
1491 				&& _FP_ISSIGNAN (fs, wc, X))		\
1492 			       ? FP_EX_INVALID_SNAN			\
1493 			       : 0));					\
1494 	}								\
1495       else if ((rsigned) != 2						\
1496 	       && (X##_e >= (_FP_EXPMAX_##fs < _FP_EXPBIAS_##fs + (rsize) \
1497 			     ? _FP_EXPMAX_##fs				\
1498 			     : (_FP_EXPBIAS_##fs + (rsize)		\
1499 				- ((rsigned) > 0 || X##_s)))		\
1500 		   || (!(rsigned) && X##_s)))				\
1501 	{								\
1502 	  /* Overflow or converting to the most negative integer.  */	\
1503 	  if (rsigned)							\
1504 	    {								\
1505 	      (r) = 1;							\
1506 	      (r) <<= (rsize) - 1;					\
1507 	      (r) -= 1 - X##_s;						\
1508 	    }								\
1509 	  else								\
1510 	    {								\
1511 	      (r) = 0;							\
1512 	      if (!X##_s)						\
1513 		(r) = ~(r);						\
1514 	    }								\
1515 									\
1516 	  if (_FP_EXPBIAS_##fs + (rsize) - 1 < _FP_EXPMAX_##fs		\
1517 	      && (rsigned)						\
1518 	      && X##_s							\
1519 	      && X##_e == _FP_EXPBIAS_##fs + (rsize) - 1)		\
1520 	    {								\
1521 	      /* Possibly converting to most negative integer; check the \
1522 		 mantissa.  */						\
1523 	      int _FP_TO_INT_inexact = 0;				\
1524 	      (void) ((_FP_FRACBITS_##fs > (rsize))			\
1525 		      ? ({						\
1526 			  _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,	\
1527 					      _FP_FRACBITS_##fs - (rsize), \
1528 					      _FP_FRACBITS_##fs);	\
1529 			  0;						\
1530 			})						\
1531 		      : 0);						\
1532 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1533 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1534 	      else if (_FP_TO_INT_inexact)				\
1535 		FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1536 	    }								\
1537 	  else								\
1538 	    FP_SET_EXCEPTION (FP_EX_INVALID				\
1539 			      | FP_EX_INVALID_CVI			\
1540 			      | ((FP_EX_INVALID_SNAN			\
1541 				  && _FP_ISSIGNAN (fs, wc, X))		\
1542 				 ? FP_EX_INVALID_SNAN			\
1543 				 : 0));					\
1544 	}								\
1545       else								\
1546 	{								\
1547 	  int _FP_TO_INT_inexact = 0;					\
1548 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;		\
1549 	  if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)	\
1550 	    {								\
1551 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1552 	      (r) <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1553 	    }								\
1554 	  else								\
1555 	    {								\
1556 	      _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,		\
1557 				  (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1558 				   - X##_e),				\
1559 				  _FP_FRACBITS_##fs);			\
1560 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1561 	    }								\
1562 	  if ((rsigned) && X##_s)					\
1563 	    (r) = -(r);							\
1564 	  if ((rsigned) == 2 && X##_e >= _FP_EXPBIAS_##fs + (rsize) - 1) \
1565 	    {								\
1566 	      /* Overflow or converting to the most negative integer.  */ \
1567 	      if (X##_e > _FP_EXPBIAS_##fs + (rsize) - 1		\
1568 		  || !X##_s						\
1569 		  || (r) != (((typeof (r)) 1) << ((rsize) - 1)))	\
1570 		{							\
1571 		  _FP_TO_INT_inexact = 0;				\
1572 		  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1573 		}							\
1574 	    }								\
1575 	  if (_FP_TO_INT_inexact)					\
1576 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1577 	}								\
1578     }									\
1579   while (0)
1580 
1581 /* Convert from floating point to integer, rounding according to the
1582    current rounding direction.  Input is raw.  RSIGNED is as for
1583    _FP_TO_INT.  */
1584 #define _FP_TO_INT_ROUND(fs, wc, r, X, rsize, rsigned)			\
1585   do									\
1586     {									\
1587       __label__ _FP_TO_INT_ROUND_done;					\
1588       if (X##_e < _FP_EXPBIAS_##fs)					\
1589 	{								\
1590 	  int _FP_TO_INT_ROUND_rounds_away = 0;				\
1591 	  if (X##_e == 0)						\
1592 	    {								\
1593 	      if (_FP_FRAC_ZEROP_##wc (X))				\
1594 		{							\
1595 		  (r) = 0;						\
1596 		  goto _FP_TO_INT_ROUND_done;				\
1597 		}							\
1598 	      else							\
1599 		{							\
1600 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1601 		  if (FP_DENORM_ZERO)					\
1602 		    {							\
1603 		      (r) = 0;						\
1604 		      goto _FP_TO_INT_ROUND_done;			\
1605 		    }							\
1606 		}							\
1607 	    }								\
1608 	  /* The result is 0, 1 or -1 depending on the rounding mode;	\
1609 	     -1 may cause overflow in the unsigned case.  */		\
1610 	  switch (FP_ROUNDMODE)						\
1611 	    {								\
1612 	    case FP_RND_NEAREST:					\
1613 	      _FP_TO_INT_ROUND_rounds_away				\
1614 		= (X##_e == _FP_EXPBIAS_##fs - 1			\
1615 		   && !_FP_FRAC_ZEROP_##wc (X));			\
1616 	      break;							\
1617 	    case FP_RND_ZERO:						\
1618 	      /* _FP_TO_INT_ROUND_rounds_away is already 0.  */		\
1619 	      break;							\
1620 	    case FP_RND_PINF:						\
1621 	      _FP_TO_INT_ROUND_rounds_away = !X##_s;			\
1622 	      break;							\
1623 	    case FP_RND_MINF:						\
1624 	      _FP_TO_INT_ROUND_rounds_away = X##_s;			\
1625 	      break;							\
1626 	    }								\
1627 	  if ((rsigned) == 0 && _FP_TO_INT_ROUND_rounds_away && X##_s)	\
1628 	    {								\
1629 	      /* Result of -1 for an unsigned conversion.  */		\
1630 	      (r) = 0;							\
1631 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1632 	    }								\
1633 	  else if ((rsize) == 1 && (rsigned) > 0			\
1634 		   && _FP_TO_INT_ROUND_rounds_away && !X##_s)		\
1635 	    {								\
1636 	      /* Converting to a 1-bit signed bit-field, which cannot	\
1637 		 represent +1.  */					\
1638 	      (r) = ((rsigned) == 2 ? -1 : 0);				\
1639 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1640 	    }								\
1641 	  else								\
1642 	    {								\
1643 	      (r) = (_FP_TO_INT_ROUND_rounds_away			\
1644 		     ? (X##_s ? -1 : 1)					\
1645 		     : 0);						\
1646 	      FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1647 	    }								\
1648 	}								\
1649       else if ((rsigned) == 2						\
1650 	       && (X##_e						\
1651 		   >= ((_FP_EXPMAX_##fs					\
1652 			< _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1) \
1653 		       ? _FP_EXPMAX_##fs				\
1654 		       : _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1))) \
1655 	{								\
1656 	  /* Overflow resulting in 0.  */				\
1657 	  (r) = 0;							\
1658 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1659 			    | FP_EX_INVALID_CVI				\
1660 			    | ((FP_EX_INVALID_SNAN			\
1661 				&& _FP_ISSIGNAN (fs, wc, X))		\
1662 			       ? FP_EX_INVALID_SNAN			\
1663 			       : 0));					\
1664 	}								\
1665       else if ((rsigned) != 2						\
1666 	       && (X##_e >= (_FP_EXPMAX_##fs < _FP_EXPBIAS_##fs + (rsize) \
1667 			     ? _FP_EXPMAX_##fs				\
1668 			     : (_FP_EXPBIAS_##fs + (rsize)		\
1669 				- ((rsigned) > 0 && !X##_s)))		\
1670 		   || ((rsigned) == 0 && X##_s)))			\
1671 	{								\
1672 	  /* Definite overflow (does not require rounding to tell).  */	\
1673 	  if ((rsigned) != 0)						\
1674 	    {								\
1675 	      (r) = 1;							\
1676 	      (r) <<= (rsize) - 1;					\
1677 	      (r) -= 1 - X##_s;						\
1678 	    }								\
1679 	  else								\
1680 	    {								\
1681 	      (r) = 0;							\
1682 	      if (!X##_s)						\
1683 		(r) = ~(r);						\
1684 	    }								\
1685 									\
1686 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1687 			    | FP_EX_INVALID_CVI				\
1688 			    | ((FP_EX_INVALID_SNAN			\
1689 				&& _FP_ISSIGNAN (fs, wc, X))		\
1690 			       ? FP_EX_INVALID_SNAN			\
1691 			       : 0));					\
1692 	}								\
1693       else								\
1694 	{								\
1695 	  /* The value is finite, with magnitude at least 1.  If	\
1696 	     the conversion is unsigned, the value is positive.		\
1697 	     If RSIGNED is not 2, the value does not definitely		\
1698 	     overflow by virtue of its exponent, but may still turn	\
1699 	     out to overflow after rounding; if RSIGNED is 2, the	\
1700 	     exponent may be such that the value definitely overflows,	\
1701 	     but at least one mantissa bit will not be shifted out.  */ \
1702 	  int _FP_TO_INT_ROUND_inexact = 0;				\
1703 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;		\
1704 	  if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)	\
1705 	    {								\
1706 	      /* The value is an integer, no rounding needed.  */	\
1707 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1708 	      (r) <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1709 	    }								\
1710 	  else								\
1711 	    {								\
1712 	      /* May need to shift in order to round (unless there	\
1713 		 are exactly _FP_WORKBITS fractional bits already).  */	\
1714 	      int _FP_TO_INT_ROUND_rshift				\
1715 		= (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs			\
1716 		   - 1 - _FP_WORKBITS - X##_e);				\
1717 	      if (_FP_TO_INT_ROUND_rshift > 0)				\
1718 		_FP_FRAC_SRS_##wc (X, _FP_TO_INT_ROUND_rshift,		\
1719 				   _FP_WFRACBITS_##fs);			\
1720 	      else if (_FP_TO_INT_ROUND_rshift < 0)			\
1721 		_FP_FRAC_SLL_##wc (X, -_FP_TO_INT_ROUND_rshift);	\
1722 	      /* Round like _FP_ROUND, but setting			\
1723 		 _FP_TO_INT_ROUND_inexact instead of directly setting	\
1724 		 the "inexact" exception, since it may turn out we	\
1725 		 should set "invalid" instead.  */			\
1726 	      if (_FP_FRAC_LOW_##wc (X) & 7)				\
1727 		{							\
1728 		  _FP_TO_INT_ROUND_inexact = 1;				\
1729 		  switch (FP_ROUNDMODE)					\
1730 		    {							\
1731 		    case FP_RND_NEAREST:				\
1732 		      _FP_ROUND_NEAREST (wc, X);			\
1733 		      break;						\
1734 		    case FP_RND_ZERO:					\
1735 		      _FP_ROUND_ZERO (wc, X);				\
1736 		      break;						\
1737 		    case FP_RND_PINF:					\
1738 		      _FP_ROUND_PINF (wc, X);				\
1739 		      break;						\
1740 		    case FP_RND_MINF:					\
1741 		      _FP_ROUND_MINF (wc, X);				\
1742 		      break;						\
1743 		    }							\
1744 		}							\
1745 	      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
1746 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1747 	    }								\
1748 	  if ((rsigned) != 0 && X##_s)					\
1749 	    (r) = -(r);							\
1750 	  /* An exponent of RSIZE - 1 always needs testing for		\
1751 	     overflow (either directly overflowing, or overflowing	\
1752 	     when rounding up results in 2^RSIZE).  An exponent of	\
1753 	     RSIZE - 2 can overflow for positive values when rounding	\
1754 	     up to 2^(RSIZE-1), but cannot overflow for negative	\
1755 	     values.  Smaller exponents cannot overflow.  */		\
1756 	  if (X##_e >= (_FP_EXPBIAS_##fs + (rsize) - 1			\
1757 			- ((rsigned) > 0 && !X##_s)))			\
1758 	    {								\
1759 	      if (X##_e > _FP_EXPBIAS_##fs + (rsize) - 1		\
1760 		  || (X##_e == _FP_EXPBIAS_##fs + (rsize) - 1		\
1761 		      && (X##_s						\
1762 			  ? (r) != (((typeof (r)) 1) << ((rsize) - 1))	\
1763 			  : ((rsigned) > 0 || (r) == 0)))		\
1764 		  || ((rsigned) > 0					\
1765 		      && !X##_s						\
1766 		      && X##_e == _FP_EXPBIAS_##fs + (rsize) - 2	\
1767 		      && (r) == (((typeof (r)) 1) << ((rsize) - 1))))	\
1768 		{							\
1769 		  if ((rsigned) != 2)					\
1770 		    {							\
1771 		      if ((rsigned) != 0)				\
1772 			{						\
1773 			  (r) = 1;					\
1774 			  (r) <<= (rsize) - 1;				\
1775 			  (r) -= 1 - X##_s;				\
1776 			}						\
1777 		      else						\
1778 			{						\
1779 			  (r) = 0;					\
1780 			  (r) = ~(r);					\
1781 			}						\
1782 		    }							\
1783 		  _FP_TO_INT_ROUND_inexact = 0;				\
1784 		  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1785 		}							\
1786 	    }								\
1787 	  if (_FP_TO_INT_ROUND_inexact)					\
1788 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1789 	}								\
1790     _FP_TO_INT_ROUND_done: ;						\
1791     }									\
1792   while (0)
1793 
1794 /* Convert integer to fp.  Output is raw.  RTYPE is unsigned even if
1795    input is signed.  */
1796 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)			\
1797   do									\
1798     {									\
1799       __label__ pack_semiraw;						\
1800       if (r)								\
1801 	{								\
1802 	  rtype _FP_FROM_INT_ur = (r);					\
1803 									\
1804 	  if ((X##_s = ((r) < 0)))					\
1805 	    _FP_FROM_INT_ur = -_FP_FROM_INT_ur;				\
1806 									\
1807 	  _FP_STATIC_ASSERT ((rsize) <= 2 * _FP_W_TYPE_SIZE,		\
1808 			     "rsize too large");			\
1809 	  (void) (((rsize) <= _FP_W_TYPE_SIZE)				\
1810 		  ? ({							\
1811 		      int _FP_FROM_INT_lz;				\
1812 		      __FP_CLZ (_FP_FROM_INT_lz,			\
1813 				(_FP_W_TYPE) _FP_FROM_INT_ur);		\
1814 		      X##_e = (_FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1	\
1815 			       - _FP_FROM_INT_lz);			\
1816 		    })							\
1817 		  : ({						\
1818 		      int _FP_FROM_INT_lz;				\
1819 		      __FP_CLZ_2 (_FP_FROM_INT_lz,			\
1820 				  (_FP_W_TYPE) (_FP_FROM_INT_ur		\
1821 						>> _FP_W_TYPE_SIZE),	\
1822 				  (_FP_W_TYPE) _FP_FROM_INT_ur);	\
1823 		      X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1824 			       - _FP_FROM_INT_lz);			\
1825 		    }));						\
1826 									\
1827 	  if ((rsize) - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs		\
1828 	      && X##_e >= _FP_EXPMAX_##fs)				\
1829 	    {								\
1830 	      /* Exponent too big; overflow to infinity.  (May also	\
1831 		 happen after rounding below.)  */			\
1832 	      _FP_OVERFLOW_SEMIRAW (fs, wc, X);				\
1833 	      goto pack_semiraw;					\
1834 	    }								\
1835 									\
1836 	  if ((rsize) <= _FP_FRACBITS_##fs				\
1837 	      || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs)		\
1838 	    {								\
1839 	      /* Exactly representable; shift left.  */			\
1840 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1841 	      if (_FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1 - X##_e > 0)	\
1842 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1843 				       + _FP_FRACBITS_##fs - 1 - X##_e)); \
1844 	    }								\
1845 	  else								\
1846 	    {								\
1847 	      /* More bits in integer than in floating type; need to	\
1848 		 round.  */						\
1849 	      if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e)	\
1850 		_FP_FROM_INT_ur						\
1851 		  = ((_FP_FROM_INT_ur >> (X##_e - _FP_EXPBIAS_##fs	\
1852 					  - _FP_WFRACBITS_##fs + 1))	\
1853 		     | ((_FP_FROM_INT_ur				\
1854 			 << ((rsize) - (X##_e - _FP_EXPBIAS_##fs	\
1855 					- _FP_WFRACBITS_##fs + 1)))	\
1856 			!= 0));						\
1857 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1858 	      if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1859 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1860 				       + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1861 	      _FP_FRAC_HIGH_##fs (X) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
1862 	    pack_semiraw:						\
1863 	      _FP_PACK_SEMIRAW (fs, wc, X);				\
1864 	    }								\
1865 	}								\
1866       else								\
1867 	{								\
1868 	  X##_s = 0;							\
1869 	  X##_e = 0;							\
1870 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
1871 	}								\
1872     }									\
1873   while (0)
1874 
1875 
1876 /* Extend from a narrower floating-point format to a wider one.  Input
1877    and output are raw.  If CHECK_NAN, then signaling NaNs are
1878    converted to quiet with the "invalid" exception raised; otherwise
1879    signaling NaNs remain signaling with no exception.  */
1880 #define _FP_EXTEND_CNAN(dfs, sfs, dwc, swc, D, S, check_nan)		\
1881   do									\
1882     {									\
1883       _FP_STATIC_ASSERT (_FP_FRACBITS_##dfs >= _FP_FRACBITS_##sfs,	\
1884 			 "destination mantissa narrower than source");	\
1885       _FP_STATIC_ASSERT ((_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs		\
1886 			  >= _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs),	\
1887 			 "destination max exponent smaller"		\
1888 			 " than source");				\
1889       _FP_STATIC_ASSERT (((_FP_EXPBIAS_##dfs				\
1890 			   >= (_FP_EXPBIAS_##sfs			\
1891 			       + _FP_FRACBITS_##sfs - 1))		\
1892 			  || (_FP_EXPBIAS_##dfs == _FP_EXPBIAS_##sfs)), \
1893 			 "source subnormals do not all become normal,"	\
1894 			 " but bias not the same");			\
1895       D##_s = S##_s;							\
1896       _FP_FRAC_COPY_##dwc##_##swc (D, S);				\
1897       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1898 	{								\
1899 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1900 	  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1901 	}								\
1902       else								\
1903 	{								\
1904 	  if (S##_e == 0)						\
1905 	    {								\
1906 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
1907 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1908 		D##_e = 0;						\
1909 	      else if (_FP_EXPBIAS_##dfs				\
1910 		       < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1)	\
1911 		{							\
1912 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1913 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1914 					  - _FP_FRACBITS_##sfs));	\
1915 		  D##_e = 0;						\
1916 		  if (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)		\
1917 		    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
1918 		}							\
1919 	      else							\
1920 		{							\
1921 		  int FP_EXTEND_lz;					\
1922 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1923 		  _FP_FRAC_CLZ_##swc (FP_EXTEND_lz, S);			\
1924 		  _FP_FRAC_SLL_##dwc (D,				\
1925 				      FP_EXTEND_lz + _FP_FRACBITS_##dfs	\
1926 				      - _FP_FRACTBITS_##sfs);		\
1927 		  D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1	\
1928 			   + _FP_FRACXBITS_##sfs - FP_EXTEND_lz);	\
1929 		}							\
1930 	    }								\
1931 	  else								\
1932 	    {								\
1933 	      D##_e = _FP_EXPMAX_##dfs;					\
1934 	      if (!_FP_FRAC_ZEROP_##swc (S))				\
1935 		{							\
1936 		  if (check_nan && _FP_FRAC_SNANP (sfs, S))		\
1937 		    FP_SET_EXCEPTION (FP_EX_INVALID			\
1938 				      | FP_EX_INVALID_SNAN);		\
1939 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1940 					  - _FP_FRACBITS_##sfs));	\
1941 		  if (check_nan)					\
1942 		    _FP_SETQNAN (dfs, dwc, D);				\
1943 		}							\
1944 	    }								\
1945 	}								\
1946     }									\
1947   while (0)
1948 
1949 #define FP_EXTEND(dfs, sfs, dwc, swc, D, S)		\
1950     _FP_EXTEND_CNAN (dfs, sfs, dwc, swc, D, S, 1)
1951 
1952 /* Truncate from a wider floating-point format to a narrower one.
1953    Input and output are semi-raw.  */
1954 #define FP_TRUNC(dfs, sfs, dwc, swc, D, S)				\
1955   do									\
1956     {									\
1957       _FP_STATIC_ASSERT (_FP_FRACBITS_##sfs >= _FP_FRACBITS_##dfs,	\
1958 			 "destination mantissa wider than source");	\
1959       _FP_STATIC_ASSERT (((_FP_EXPBIAS_##sfs				\
1960 			   >= (_FP_EXPBIAS_##dfs			\
1961 			       + _FP_FRACBITS_##dfs - 1))		\
1962 			  || _FP_EXPBIAS_##sfs == _FP_EXPBIAS_##dfs),	\
1963 			 "source subnormals do not all become same,"	\
1964 			 " but bias not the same");			\
1965       D##_s = S##_s;							\
1966       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1967 	{								\
1968 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1969 	  if (D##_e >= _FP_EXPMAX_##dfs)				\
1970 	    _FP_OVERFLOW_SEMIRAW (dfs, dwc, D);				\
1971 	  else								\
1972 	    {								\
1973 	      if (D##_e <= 0)						\
1974 		{							\
1975 		  if (D##_e < 1 - _FP_FRACBITS_##dfs)			\
1976 		    {							\
1977 		      _FP_FRAC_SET_##swc (S, _FP_ZEROFRAC_##swc);	\
1978 		      _FP_FRAC_LOW_##swc (S) |= 1;			\
1979 		    }							\
1980 		  else							\
1981 		    {							\
1982 		      _FP_FRAC_HIGH_##sfs (S) |= _FP_IMPLBIT_SH_##sfs;	\
1983 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
1984 					      - _FP_WFRACBITS_##dfs	\
1985 					      + 1 - D##_e),		\
1986 					  _FP_WFRACBITS_##sfs);		\
1987 		    }							\
1988 		  D##_e = 0;						\
1989 		}							\
1990 	      else							\
1991 		_FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs		\
1992 					- _FP_WFRACBITS_##dfs),		\
1993 				    _FP_WFRACBITS_##sfs);		\
1994 	      _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
1995 	    }								\
1996 	}								\
1997       else								\
1998 	{								\
1999 	  if (S##_e == 0)						\
2000 	    {								\
2001 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
2002 	      D##_e = 0;						\
2003 	      if (_FP_FRAC_ZEROP_##swc (S))				\
2004 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
2005 	      else							\
2006 		{							\
2007 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
2008 		  if (_FP_EXPBIAS_##sfs					\
2009 		      < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1)	\
2010 		    {							\
2011 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
2012 					      - _FP_WFRACBITS_##dfs),	\
2013 					  _FP_WFRACBITS_##sfs);		\
2014 		      _FP_FRAC_COPY_##dwc##_##swc (D, S);		\
2015 		    }							\
2016 		  else							\
2017 		    {							\
2018 		      _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);	\
2019 		      _FP_FRAC_LOW_##dwc (D) |= 1;			\
2020 		    }							\
2021 		}							\
2022 	    }								\
2023 	  else								\
2024 	    {								\
2025 	      D##_e = _FP_EXPMAX_##dfs;					\
2026 	      if (_FP_FRAC_ZEROP_##swc (S))				\
2027 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
2028 	      else							\
2029 		{							\
2030 		  _FP_CHECK_SIGNAN_SEMIRAW (sfs, swc, S);		\
2031 		  _FP_FRAC_SRL_##swc (S, (_FP_WFRACBITS_##sfs		\
2032 					  - _FP_WFRACBITS_##dfs));	\
2033 		  _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
2034 		  /* Semi-raw NaN must have all workbits cleared.  */	\
2035 		  _FP_FRAC_LOW_##dwc (D)				\
2036 		    &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1);		\
2037 		  _FP_SETQNAN_SEMIRAW (dfs, dwc, D);			\
2038 		}							\
2039 	    }								\
2040 	}								\
2041     }									\
2042   while (0)
2043 
2044 /* Helper primitives.  */
2045 
2046 /* Count leading zeros in a word.  */
2047 
2048 #ifndef __FP_CLZ
2049 /* GCC 3.4 and later provide the builtins for us.  */
2050 # define __FP_CLZ(r, x)							\
2051   do									\
2052     {									\
2053       _FP_STATIC_ASSERT ((sizeof (_FP_W_TYPE) == sizeof (unsigned int)	\
2054 			  || (sizeof (_FP_W_TYPE)			\
2055 			      == sizeof (unsigned long))		\
2056 			  || (sizeof (_FP_W_TYPE)			\
2057 			      == sizeof (unsigned long long))),		\
2058 			 "_FP_W_TYPE size unsupported for clz");	\
2059       if (sizeof (_FP_W_TYPE) == sizeof (unsigned int))			\
2060 	(r) = __builtin_clz (x);					\
2061       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long))		\
2062 	(r) = __builtin_clzl (x);					\
2063       else /* sizeof (_FP_W_TYPE) == sizeof (unsigned long long).  */	\
2064 	(r) = __builtin_clzll (x);					\
2065     }									\
2066   while (0)
2067 #endif /* ndef __FP_CLZ */
2068 
2069 #define _FP_DIV_HELP_imm(q, r, n, d)		\
2070   do						\
2071     {						\
2072       (q) = (n) / (d), (r) = (n) % (d);		\
2073     }						\
2074   while (0)
2075 
2076 
2077 /* A restoring bit-by-bit division primitive.  */
2078 
2079 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y)				\
2080   do									\
2081     {									\
2082       int _FP_DIV_MEAT_N_loop_count = _FP_WFRACBITS_##fs;		\
2083       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_u);			\
2084       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_v);			\
2085       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_u, X);			\
2086       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_v, Y);			\
2087       _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);				\
2088       /* Normalize _FP_DIV_MEAT_N_LOOP_U and _FP_DIV_MEAT_N_LOOP_V.  */	\
2089       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, _FP_WFRACXBITS_##fs);	\
2090       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_v, _FP_WFRACXBITS_##fs);	\
2091       /* First round.  Since the operands are normalized, either the	\
2092 	 first or second bit will be set in the fraction.  Produce a	\
2093 	 normalized result by checking which and adjusting the loop	\
2094 	 count and exponent accordingly.  */				\
2095       if (_FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u, _FP_DIV_MEAT_N_loop_v))	\
2096 	{								\
2097 	  _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
2098 			     _FP_DIV_MEAT_N_loop_u,			\
2099 			     _FP_DIV_MEAT_N_loop_v);			\
2100 	  _FP_FRAC_LOW_##wc (R) |= 1;					\
2101 	  _FP_DIV_MEAT_N_loop_count--;					\
2102 	}								\
2103       else								\
2104 	R##_e--;							\
2105       /* Subsequent rounds.  */						\
2106       do								\
2107 	{								\
2108 	  int _FP_DIV_MEAT_N_loop_msb					\
2109 	    = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (_FP_DIV_MEAT_N_loop_u) < 0; \
2110 	  _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, 1);			\
2111 	  _FP_FRAC_SLL_##wc (R, 1);					\
2112 	  if (_FP_DIV_MEAT_N_loop_msb					\
2113 	      || _FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u,			\
2114 				_FP_DIV_MEAT_N_loop_v))			\
2115 	    {								\
2116 	      _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
2117 				 _FP_DIV_MEAT_N_loop_u,			\
2118 				 _FP_DIV_MEAT_N_loop_v);		\
2119 	      _FP_FRAC_LOW_##wc (R) |= 1;				\
2120 	    }								\
2121 	}								\
2122       while (--_FP_DIV_MEAT_N_loop_count > 0);				\
2123       /* If there's anything left in _FP_DIV_MEAT_N_LOOP_U, the result	\
2124 	 is inexact.  */						\
2125       _FP_FRAC_LOW_##wc (R)						\
2126 	|= !_FP_FRAC_ZEROP_##wc (_FP_DIV_MEAT_N_loop_u);		\
2127     }									\
2128   while (0)
2129 
2130 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
2131 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
2132 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)
2133 
2134 #endif /* !SOFT_FP_OP_COMMON_H */
2135