xref: /netbsd-src/external/gpl3/gcc.old/dist/libgcc/soft-fp/op-common.h (revision bdc22b2e01993381dcefeff2bc9b56ca75a4235c)
1 /* Software floating-point emulation. Common operations.
2    Copyright (C) 1997-2014 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 #define _FP_DECL(wc, X)				\
33   _FP_I_TYPE X##_c __attribute__ ((unused));	\
34   _FP_I_TYPE X##_s __attribute__ ((unused));	\
35   _FP_I_TYPE X##_e __attribute__ ((unused));	\
36   _FP_FRAC_DECL_##wc (X)
37 
38 /* Test whether the qNaN bit denotes a signaling NaN.  */
39 #define _FP_FRAC_SNANP(fs, X)				\
40   ((_FP_QNANNEGATEDP)					\
41    ? (_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs)	\
42    : !(_FP_FRAC_HIGH_RAW_##fs (X) & _FP_QNANBIT_##fs))
43 #define _FP_FRAC_SNANP_SEMIRAW(fs, X)			\
44   ((_FP_QNANNEGATEDP)					\
45    ? (_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs)	\
46    : !(_FP_FRAC_HIGH_##fs (X) & _FP_QNANBIT_SH_##fs))
47 
48 /* Finish truly unpacking a native fp value by classifying the kind
49    of fp value and normalizing both the exponent and the fraction.  */
50 
51 #define _FP_UNPACK_CANONICAL(fs, wc, X)				\
52   do								\
53     {								\
54       switch (X##_e)						\
55 	{							\
56 	default:						\
57 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;	\
58 	  _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);			\
59 	  X##_e -= _FP_EXPBIAS_##fs;				\
60 	  X##_c = FP_CLS_NORMAL;				\
61 	  break;						\
62 								\
63 	case 0:							\
64 	  if (_FP_FRAC_ZEROP_##wc (X))				\
65 	    X##_c = FP_CLS_ZERO;				\
66 	  else if (FP_DENORM_ZERO)				\
67 	    {							\
68 	      X##_c = FP_CLS_ZERO;				\
69 	      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
70 	      FP_SET_EXCEPTION (FP_EX_DENORM);			\
71 	    }							\
72 	  else							\
73 	    {							\
74 	      /* A denormalized number.  */			\
75 	      _FP_I_TYPE _FP_UNPACK_CANONICAL_shift;		\
76 	      _FP_FRAC_CLZ_##wc (_FP_UNPACK_CANONICAL_shift,	\
77 				 X);				\
78 	      _FP_UNPACK_CANONICAL_shift -= _FP_FRACXBITS_##fs;	\
79 	      _FP_FRAC_SLL_##wc (X, (_FP_UNPACK_CANONICAL_shift \
80 				     + _FP_WORKBITS));		\
81 	      X##_e -= (_FP_EXPBIAS_##fs - 1			\
82 			+ _FP_UNPACK_CANONICAL_shift);		\
83 	      X##_c = FP_CLS_NORMAL;				\
84 	      FP_SET_EXCEPTION (FP_EX_DENORM);			\
85 	    }							\
86 	  break;						\
87 								\
88 	case _FP_EXPMAX_##fs:					\
89 	  if (_FP_FRAC_ZEROP_##wc (X))				\
90 	    X##_c = FP_CLS_INF;					\
91 	  else							\
92 	    {							\
93 	      X##_c = FP_CLS_NAN;				\
94 	      /* Check for signaling NaN.  */			\
95 	      if (_FP_FRAC_SNANP (fs, X))			\
96 		FP_SET_EXCEPTION (FP_EX_INVALID			\
97 				  | FP_EX_INVALID_SNAN);	\
98 	    }							\
99 	  break;						\
100 	}							\
101     }								\
102   while (0)
103 
104 /* Finish unpacking an fp value in semi-raw mode: the mantissa is
105    shifted by _FP_WORKBITS but the implicit MSB is not inserted and
106    other classification is not done.  */
107 #define _FP_UNPACK_SEMIRAW(fs, wc, X)	_FP_FRAC_SLL_##wc (X, _FP_WORKBITS)
108 
109 /* Check whether a raw or semi-raw input value should be flushed to
110    zero, and flush it to zero if so.  */
111 #define _FP_CHECK_FLUSH_ZERO(fs, wc, X)			\
112   do							\
113     {							\
114       if (FP_DENORM_ZERO				\
115 	  && X##_e == 0					\
116 	  && !_FP_FRAC_ZEROP_##wc (X))			\
117 	{						\
118 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);	\
119 	  FP_SET_EXCEPTION (FP_EX_DENORM);		\
120 	}						\
121     }							\
122   while (0)
123 
124 /* A semi-raw value has overflowed to infinity.  Adjust the mantissa
125    and exponent appropriately.  */
126 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X)			\
127   do							\
128     {							\
129       if (FP_ROUNDMODE == FP_RND_NEAREST		\
130 	  || (FP_ROUNDMODE == FP_RND_PINF && !X##_s)	\
131 	  || (FP_ROUNDMODE == FP_RND_MINF && X##_s))	\
132 	{						\
133 	  X##_e = _FP_EXPMAX_##fs;			\
134 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);	\
135 	}						\
136       else						\
137 	{						\
138 	  X##_e = _FP_EXPMAX_##fs - 1;			\
139 	  _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);	\
140 	}						\
141       FP_SET_EXCEPTION (FP_EX_INEXACT);			\
142       FP_SET_EXCEPTION (FP_EX_OVERFLOW);		\
143     }							\
144   while (0)
145 
146 /* Check for a semi-raw value being a signaling NaN and raise the
147    invalid exception if so.  */
148 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X)			\
149   do								\
150     {								\
151       if (X##_e == _FP_EXPMAX_##fs				\
152 	  && !_FP_FRAC_ZEROP_##wc (X)				\
153 	  && _FP_FRAC_SNANP_SEMIRAW (fs, X))			\
154 	FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SNAN);	\
155     }								\
156   while (0)
157 
158 /* Choose a NaN result from an operation on two semi-raw NaN
159    values.  */
160 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP)			\
161   do									\
162     {									\
163       /* _FP_CHOOSENAN expects raw values, so shift as required.  */	\
164       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);				\
165       _FP_FRAC_SRL_##wc (Y, _FP_WORKBITS);				\
166       _FP_CHOOSENAN (fs, wc, R, X, Y, OP);				\
167       _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);				\
168     }									\
169   while (0)
170 
171 /* Make the fractional part a quiet NaN, preserving the payload
172    if possible, otherwise make it the canonical quiet NaN and set
173    the sign bit accordingly.  */
174 #define _FP_SETQNAN(fs, wc, X)					\
175   do								\
176     {								\
177       if (_FP_QNANNEGATEDP)					\
178 	{							\
179 	  _FP_FRAC_HIGH_RAW_##fs (X) &= _FP_QNANBIT_##fs - 1;	\
180 	  if (_FP_FRAC_ZEROP_##wc (X))				\
181 	    {							\
182 	      X##_s = _FP_NANSIGN_##fs;				\
183 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
184 	    }							\
185 	}							\
186       else							\
187 	_FP_FRAC_HIGH_RAW_##fs (X) |= _FP_QNANBIT_##fs;		\
188     }								\
189   while (0)
190 #define _FP_SETQNAN_SEMIRAW(fs, wc, X)				\
191   do								\
192     {								\
193       if (_FP_QNANNEGATEDP)					\
194 	{							\
195 	  _FP_FRAC_HIGH_##fs (X) &= _FP_QNANBIT_SH_##fs - 1;	\
196 	  if (_FP_FRAC_ZEROP_##wc (X))				\
197 	    {							\
198 	      X##_s = _FP_NANSIGN_##fs;				\
199 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
200 	      _FP_FRAC_SLL_##wc (X, _FP_WORKBITS);		\
201 	    }							\
202 	}							\
203       else							\
204 	_FP_FRAC_HIGH_##fs (X) |= _FP_QNANBIT_SH_##fs;		\
205     }								\
206   while (0)
207 
208 /* Test whether a biased exponent is normal (not zero or maximum).  */
209 #define _FP_EXP_NORMAL(fs, wc, X)	(((X##_e + 1) & _FP_EXPMAX_##fs) > 1)
210 
211 /* Prepare to pack an fp value in semi-raw mode: the mantissa is
212    rounded and shifted right, with the rounding possibly increasing
213    the exponent (including changing a finite value to infinity).  */
214 #define _FP_PACK_SEMIRAW(fs, wc, X)				\
215   do								\
216     {								\
217       int _FP_PACK_SEMIRAW_is_tiny				\
218 	= X##_e == 0 && !_FP_FRAC_ZEROP_##wc (X);		\
219       if (_FP_TININESS_AFTER_ROUNDING				\
220 	  && _FP_PACK_SEMIRAW_is_tiny)				\
221 	{							\
222 	  FP_DECL_##fs (_FP_PACK_SEMIRAW_T);			\
223 	  _FP_FRAC_COPY_##wc (_FP_PACK_SEMIRAW_T, X);		\
224 	  _FP_PACK_SEMIRAW_T##_s = X##_s;			\
225 	  _FP_PACK_SEMIRAW_T##_e = X##_e;			\
226 	  _FP_FRAC_SLL_##wc (_FP_PACK_SEMIRAW_T, 1);		\
227 	  _FP_ROUND (wc, _FP_PACK_SEMIRAW_T);			\
228 	  if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_SEMIRAW_T))	\
229 	    _FP_PACK_SEMIRAW_is_tiny = 0;			\
230 	}							\
231       _FP_ROUND (wc, X);					\
232       if (_FP_PACK_SEMIRAW_is_tiny)				\
233 	{							\
234 	  if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)		\
235 	      || (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW))	\
236 	    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
237 	}							\
238       if (_FP_FRAC_HIGH_##fs (X)				\
239 	  & (_FP_OVERFLOW_##fs >> 1))				\
240 	{							\
241 	  _FP_FRAC_HIGH_##fs (X) &= ~(_FP_OVERFLOW_##fs >> 1);	\
242 	  X##_e++;						\
243 	  if (X##_e == _FP_EXPMAX_##fs)				\
244 	    _FP_OVERFLOW_SEMIRAW (fs, wc, X);			\
245 	}							\
246       _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
247       if (X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
248 	{							\
249 	  if (!_FP_KEEPNANFRACP)				\
250 	    {							\
251 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);		\
252 	      X##_s = _FP_NANSIGN_##fs;				\
253 	    }							\
254 	  else							\
255 	    _FP_SETQNAN (fs, wc, X);				\
256 	}							\
257     }								\
258   while (0)
259 
260 /* Before packing the bits back into the native fp result, take care
261    of such mundane things as rounding and overflow.  Also, for some
262    kinds of fp values, the original parts may not have been fully
263    extracted -- but that is ok, we can regenerate them now.  */
264 
265 #define _FP_PACK_CANONICAL(fs, wc, X)					\
266   do									\
267     {									\
268       switch (X##_c)							\
269 	{								\
270 	case FP_CLS_NORMAL:						\
271 	  X##_e += _FP_EXPBIAS_##fs;					\
272 	  if (X##_e > 0)						\
273 	    {								\
274 	      _FP_ROUND (wc, X);					\
275 	      if (_FP_FRAC_OVERP_##wc (fs, X))				\
276 		{							\
277 		  _FP_FRAC_CLEAR_OVERP_##wc (fs, X);			\
278 		  X##_e++;						\
279 		}							\
280 	      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);			\
281 	      if (X##_e >= _FP_EXPMAX_##fs)				\
282 		{							\
283 		  /* Overflow.  */					\
284 		  switch (FP_ROUNDMODE)					\
285 		    {							\
286 		    case FP_RND_NEAREST:				\
287 		      X##_c = FP_CLS_INF;				\
288 		      break;						\
289 		    case FP_RND_PINF:					\
290 		      if (!X##_s)					\
291 			X##_c = FP_CLS_INF;				\
292 		      break;						\
293 		    case FP_RND_MINF:					\
294 		      if (X##_s)					\
295 			X##_c = FP_CLS_INF;				\
296 		      break;						\
297 		    }							\
298 		  if (X##_c == FP_CLS_INF)				\
299 		    {							\
300 		      /* Overflow to infinity.  */			\
301 		      X##_e = _FP_EXPMAX_##fs;				\
302 		      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
303 		    }							\
304 		  else							\
305 		    {							\
306 		      /* Overflow to maximum normal.  */		\
307 		      X##_e = _FP_EXPMAX_##fs - 1;			\
308 		      _FP_FRAC_SET_##wc (X, _FP_MAXFRAC_##wc);		\
309 		    }							\
310 		  FP_SET_EXCEPTION (FP_EX_OVERFLOW);			\
311 		  FP_SET_EXCEPTION (FP_EX_INEXACT);			\
312 		}							\
313 	    }								\
314 	  else								\
315 	    {								\
316 	      /* We've got a denormalized number.  */			\
317 	      int _FP_PACK_CANONICAL_is_tiny = 1;			\
318 	      if (_FP_TININESS_AFTER_ROUNDING && X##_e == 0)		\
319 		{							\
320 		  FP_DECL_##fs (_FP_PACK_CANONICAL_T);			\
321 		  _FP_FRAC_COPY_##wc (_FP_PACK_CANONICAL_T, X);		\
322 		  _FP_PACK_CANONICAL_T##_s = X##_s;			\
323 		  _FP_PACK_CANONICAL_T##_e = X##_e;			\
324 		  _FP_ROUND (wc, _FP_PACK_CANONICAL_T);			\
325 		  if (_FP_FRAC_OVERP_##wc (fs, _FP_PACK_CANONICAL_T))	\
326 		    _FP_PACK_CANONICAL_is_tiny = 0;			\
327 		}							\
328 	      X##_e = -X##_e + 1;					\
329 	      if (X##_e <= _FP_WFRACBITS_##fs)				\
330 		{							\
331 		  _FP_FRAC_SRS_##wc (X, X##_e, _FP_WFRACBITS_##fs);	\
332 		  _FP_ROUND (wc, X);					\
333 		  if (_FP_FRAC_HIGH_##fs (X)				\
334 		      & (_FP_OVERFLOW_##fs >> 1))			\
335 		    {							\
336 		      X##_e = 1;					\
337 		      _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);		\
338 		      FP_SET_EXCEPTION (FP_EX_INEXACT);			\
339 		    }							\
340 		  else							\
341 		    {							\
342 		      X##_e = 0;					\
343 		      _FP_FRAC_SRL_##wc (X, _FP_WORKBITS);		\
344 		    }							\
345 		  if (_FP_PACK_CANONICAL_is_tiny			\
346 		      && ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT)		\
347 			  || (FP_TRAPPING_EXCEPTIONS			\
348 			      & FP_EX_UNDERFLOW)))			\
349 		    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
350 		}							\
351 	      else							\
352 		{							\
353 		  /* Underflow to zero.  */				\
354 		  X##_e = 0;						\
355 		  if (!_FP_FRAC_ZEROP_##wc (X))				\
356 		    {							\
357 		      _FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);		\
358 		      _FP_ROUND (wc, X);				\
359 		      _FP_FRAC_LOW_##wc (X) >>= (_FP_WORKBITS);		\
360 		    }							\
361 		  FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
362 		}							\
363 	    }								\
364 	  break;							\
365 									\
366 	case FP_CLS_ZERO:						\
367 	  X##_e = 0;							\
368 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
369 	  break;							\
370 									\
371 	case FP_CLS_INF:						\
372 	  X##_e = _FP_EXPMAX_##fs;					\
373 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
374 	  break;							\
375 									\
376 	case FP_CLS_NAN:						\
377 	  X##_e = _FP_EXPMAX_##fs;					\
378 	  if (!_FP_KEEPNANFRACP)					\
379 	    {								\
380 	      _FP_FRAC_SET_##wc (X, _FP_NANFRAC_##fs);			\
381 	      X##_s = _FP_NANSIGN_##fs;					\
382 	    }								\
383 	  else								\
384 	    _FP_SETQNAN (fs, wc, X);					\
385 	  break;							\
386 	}								\
387     }									\
388   while (0)
389 
390 /* This one accepts raw argument and not cooked,  returns
391    1 if X is a signaling NaN.  */
392 #define _FP_ISSIGNAN(fs, wc, X)			\
393   ({						\
394     int _FP_ISSIGNAN_ret = 0;			\
395     if (X##_e == _FP_EXPMAX_##fs)		\
396       {						\
397 	if (!_FP_FRAC_ZEROP_##wc (X)		\
398 	    && _FP_FRAC_SNANP (fs, X))		\
399 	  _FP_ISSIGNAN_ret = 1;			\
400       }						\
401     _FP_ISSIGNAN_ret;				\
402   })
403 
404 
405 
406 
407 
408 /* Addition on semi-raw values.  */
409 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP)				\
410   do									\
411     {									\
412       _FP_CHECK_FLUSH_ZERO (fs, wc, X);					\
413       _FP_CHECK_FLUSH_ZERO (fs, wc, Y);					\
414       if (X##_s == Y##_s)						\
415 	{								\
416 	  /* Addition.  */						\
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 	  int _FP_ADD_INTERNAL_ediff = X##_e - Y##_e;			\
603 	  if (_FP_ADD_INTERNAL_ediff > 0)				\
604 	    {								\
605 	      R##_e = X##_e;						\
606 	      R##_s = X##_s;						\
607 	      if (Y##_e == 0)						\
608 		{							\
609 		  /* Y is zero or denormalized.  */			\
610 		  if (_FP_FRAC_ZEROP_##wc (Y))				\
611 		    {							\
612 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
613 		      _FP_FRAC_COPY_##wc (R, X);			\
614 		      goto sub_done;					\
615 		    }							\
616 		  else							\
617 		    {							\
618 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
619 		      _FP_ADD_INTERNAL_ediff--;				\
620 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
621 			{						\
622 			  _FP_FRAC_SUB_##wc (R, X, Y);			\
623 			  goto sub3;					\
624 			}						\
625 		      if (X##_e == _FP_EXPMAX_##fs)			\
626 			{						\
627 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
628 			  _FP_FRAC_COPY_##wc (R, X);			\
629 			  goto sub_done;				\
630 			}						\
631 		      goto sub1;					\
632 		    }							\
633 		}							\
634 	      else if (X##_e == _FP_EXPMAX_##fs)			\
635 		{							\
636 		  /* X is NaN or Inf, Y is normal.  */			\
637 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);			\
638 		  _FP_FRAC_COPY_##wc (R, X);				\
639 		  goto sub_done;					\
640 		}							\
641 									\
642 	      /* Insert implicit MSB of Y.  */				\
643 	      _FP_FRAC_HIGH_##fs (Y) |= _FP_IMPLBIT_SH_##fs;		\
644 									\
645 	    sub1:							\
646 	      /* Shift the mantissa of Y to the right			\
647 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
648 		 later for the implicit MSB of X.  */			\
649 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
650 		_FP_FRAC_SRS_##wc (Y, _FP_ADD_INTERNAL_ediff,		\
651 				   _FP_WFRACBITS_##fs);			\
652 	      else if (!_FP_FRAC_ZEROP_##wc (Y))			\
653 		_FP_FRAC_SET_##wc (Y, _FP_MINFRAC_##wc);		\
654 	      _FP_FRAC_SUB_##wc (R, X, Y);				\
655 	    }								\
656 	  else if (_FP_ADD_INTERNAL_ediff < 0)				\
657 	    {								\
658 	      _FP_ADD_INTERNAL_ediff = -_FP_ADD_INTERNAL_ediff;		\
659 	      R##_e = Y##_e;						\
660 	      R##_s = Y##_s;						\
661 	      if (X##_e == 0)						\
662 		{							\
663 		  /* X is zero or denormalized.  */			\
664 		  if (_FP_FRAC_ZEROP_##wc (X))				\
665 		    {							\
666 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
667 		      _FP_FRAC_COPY_##wc (R, Y);			\
668 		      goto sub_done;					\
669 		    }							\
670 		  else							\
671 		    {							\
672 		      FP_SET_EXCEPTION (FP_EX_DENORM);			\
673 		      _FP_ADD_INTERNAL_ediff--;				\
674 		      if (_FP_ADD_INTERNAL_ediff == 0)			\
675 			{						\
676 			  _FP_FRAC_SUB_##wc (R, Y, X);			\
677 			  goto sub3;					\
678 			}						\
679 		      if (Y##_e == _FP_EXPMAX_##fs)			\
680 			{						\
681 			  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
682 			  _FP_FRAC_COPY_##wc (R, Y);			\
683 			  goto sub_done;				\
684 			}						\
685 		      goto sub2;					\
686 		    }							\
687 		}							\
688 	      else if (Y##_e == _FP_EXPMAX_##fs)			\
689 		{							\
690 		  /* Y is NaN or Inf, X is normal.  */			\
691 		  _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);			\
692 		  _FP_FRAC_COPY_##wc (R, Y);				\
693 		  goto sub_done;					\
694 		}							\
695 									\
696 	      /* Insert implicit MSB of X.  */				\
697 	      _FP_FRAC_HIGH_##fs (X) |= _FP_IMPLBIT_SH_##fs;		\
698 									\
699 	    sub2:							\
700 	      /* Shift the mantissa of X to the right			\
701 		 _FP_ADD_INTERNAL_EDIFF steps; remember to account	\
702 		 later for the implicit MSB of Y.  */			\
703 	      if (_FP_ADD_INTERNAL_ediff <= _FP_WFRACBITS_##fs)		\
704 		_FP_FRAC_SRS_##wc (X, _FP_ADD_INTERNAL_ediff,		\
705 				   _FP_WFRACBITS_##fs);			\
706 	      else if (!_FP_FRAC_ZEROP_##wc (X))			\
707 		_FP_FRAC_SET_##wc (X, _FP_MINFRAC_##wc);		\
708 	      _FP_FRAC_SUB_##wc (R, Y, X);				\
709 	    }								\
710 	  else								\
711 	    {								\
712 	      /* ediff == 0.  */					\
713 	      if (!_FP_EXP_NORMAL (fs, wc, X))				\
714 		{							\
715 		  if (X##_e == 0)					\
716 		    {							\
717 		      /* X and Y are zero or denormalized.  */		\
718 		      R##_e = 0;					\
719 		      if (_FP_FRAC_ZEROP_##wc (X))			\
720 			{						\
721 			  _FP_FRAC_COPY_##wc (R, Y);			\
722 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
723 			    R##_s = (FP_ROUNDMODE == FP_RND_MINF);	\
724 			  else						\
725 			    {						\
726 			      FP_SET_EXCEPTION (FP_EX_DENORM);		\
727 			      R##_s = Y##_s;				\
728 			    }						\
729 			  goto sub_done;				\
730 			}						\
731 		      else if (_FP_FRAC_ZEROP_##wc (Y))			\
732 			{						\
733 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
734 			  _FP_FRAC_COPY_##wc (R, X);			\
735 			  R##_s = X##_s;				\
736 			  goto sub_done;				\
737 			}						\
738 		      else						\
739 			{						\
740 			  FP_SET_EXCEPTION (FP_EX_DENORM);		\
741 			  _FP_FRAC_SUB_##wc (R, X, Y);			\
742 			  R##_s = X##_s;				\
743 			  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs) \
744 			    {						\
745 			      /* |X| < |Y|, negate result.  */		\
746 			      _FP_FRAC_SUB_##wc (R, Y, X);		\
747 			      R##_s = Y##_s;				\
748 			    }						\
749 			  else if (_FP_FRAC_ZEROP_##wc (R))		\
750 			    R##_s = (FP_ROUNDMODE == FP_RND_MINF);	\
751 			  goto sub_done;				\
752 			}						\
753 		    }							\
754 		  else							\
755 		    {							\
756 		      /* X and Y are NaN or Inf, of opposite signs.  */	\
757 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, X);		\
758 		      _FP_CHECK_SIGNAN_SEMIRAW (fs, wc, Y);		\
759 		      R##_e = _FP_EXPMAX_##fs;				\
760 		      if (_FP_FRAC_ZEROP_##wc (X))			\
761 			{						\
762 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
763 			    {						\
764 			      /* Inf - Inf.  */				\
765 			      R##_s = _FP_NANSIGN_##fs;			\
766 			      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);	\
767 			      _FP_FRAC_SLL_##wc (R, _FP_WORKBITS);	\
768 			      FP_SET_EXCEPTION (FP_EX_INVALID		\
769 						| FP_EX_INVALID_ISI);	\
770 			    }						\
771 			  else						\
772 			    {						\
773 			      /* Inf - NaN.  */				\
774 			      R##_s = Y##_s;				\
775 			      _FP_FRAC_COPY_##wc (R, Y);		\
776 			    }						\
777 			}						\
778 		      else						\
779 			{						\
780 			  if (_FP_FRAC_ZEROP_##wc (Y))			\
781 			    {						\
782 			      /* NaN - Inf.  */				\
783 			      R##_s = X##_s;				\
784 			      _FP_FRAC_COPY_##wc (R, X);		\
785 			    }						\
786 			  else						\
787 			    {						\
788 			      /* NaN - NaN.  */				\
789 			      _FP_CHOOSENAN_SEMIRAW (fs, wc, R, X, Y, OP); \
790 			    }						\
791 			}						\
792 		      goto sub_done;					\
793 		    }							\
794 		}							\
795 	      /* The exponents of X and Y, both normal, are equal.  The	\
796 		 implicit MSBs cancel.  */				\
797 	      R##_e = X##_e;						\
798 	      _FP_FRAC_SUB_##wc (R, X, Y);				\
799 	      R##_s = X##_s;						\
800 	      if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)		\
801 		{							\
802 		  /* |X| < |Y|, negate result.  */			\
803 		  _FP_FRAC_SUB_##wc (R, Y, X);				\
804 		  R##_s = Y##_s;					\
805 		}							\
806 	      else if (_FP_FRAC_ZEROP_##wc (R))				\
807 		{							\
808 		  R##_e = 0;						\
809 		  R##_s = (FP_ROUNDMODE == FP_RND_MINF);		\
810 		  goto sub_done;					\
811 		}							\
812 	      goto norm;						\
813 	    }								\
814 	sub3:								\
815 	  if (_FP_FRAC_HIGH_##fs (R) & _FP_IMPLBIT_SH_##fs)		\
816 	    {								\
817 	      int _FP_ADD_INTERNAL_diff;				\
818 	      /* Carry into most significant bit of larger one of X and Y, \
819 		 canceling it; renormalize.  */				\
820 	      _FP_FRAC_HIGH_##fs (R) &= _FP_IMPLBIT_SH_##fs - 1;	\
821 	    norm:							\
822 	      _FP_FRAC_CLZ_##wc (_FP_ADD_INTERNAL_diff, R);		\
823 	      _FP_ADD_INTERNAL_diff -= _FP_WFRACXBITS_##fs;		\
824 	      _FP_FRAC_SLL_##wc (R, _FP_ADD_INTERNAL_diff);		\
825 	      if (R##_e <= _FP_ADD_INTERNAL_diff)			\
826 		{							\
827 		  /* R is denormalized.  */				\
828 		  _FP_ADD_INTERNAL_diff					\
829 		    = _FP_ADD_INTERNAL_diff - R##_e + 1;		\
830 		  _FP_FRAC_SRS_##wc (R, _FP_ADD_INTERNAL_diff,		\
831 				     _FP_WFRACBITS_##fs);		\
832 		  R##_e = 0;						\
833 		}							\
834 	      else							\
835 		{							\
836 		  R##_e -= _FP_ADD_INTERNAL_diff;			\
837 		  _FP_FRAC_HIGH_##fs (R) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
838 		}							\
839 	    }								\
840 	sub_done: ;							\
841 	}								\
842     }									\
843   while (0)
844 
845 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL (fs, wc, R, X, Y, '+')
846 #define _FP_SUB(fs, wc, R, X, Y)					\
847   do									\
848     {									\
849       if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
850 	Y##_s ^= 1;							\
851       _FP_ADD_INTERNAL (fs, wc, R, X, Y, '-');				\
852     }									\
853   while (0)
854 
855 
856 /* Main negation routine.  The input value is raw.  */
857 
858 #define _FP_NEG(fs, wc, R, X)			\
859   do						\
860     {						\
861       _FP_FRAC_COPY_##wc (R, X);		\
862       R##_e = X##_e;				\
863       R##_s = 1 ^ X##_s;			\
864     }						\
865   while (0)
866 
867 
868 /* Main multiplication routine.  The input values should be cooked.  */
869 
870 #define _FP_MUL(fs, wc, R, X, Y)				\
871   do								\
872     {								\
873       R##_s = X##_s ^ Y##_s;					\
874       R##_e = X##_e + Y##_e + 1;				\
875       switch (_FP_CLS_COMBINE (X##_c, Y##_c))			\
876 	{							\
877 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):	\
878 	  R##_c = FP_CLS_NORMAL;				\
879 								\
880 	  _FP_MUL_MEAT_##fs (R, X, Y);				\
881 								\
882 	  if (_FP_FRAC_OVERP_##wc (fs, R))			\
883 	    _FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);	\
884 	  else							\
885 	    R##_e--;						\
886 	  break;						\
887 								\
888 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):		\
889 	  _FP_CHOOSENAN (fs, wc, R, X, Y, '*');			\
890 	  break;						\
891 								\
892 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):	\
893 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):		\
894 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):		\
895 	  R##_s = X##_s;					\
896 								\
897 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
898 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
899 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
900 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
901 	  _FP_FRAC_COPY_##wc (R, X);				\
902 	  R##_c = X##_c;					\
903 	  break;						\
904 								\
905 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
906 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
907 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
908 	  R##_s = Y##_s;					\
909 								\
910 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
911 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
912 	  _FP_FRAC_COPY_##wc (R, Y);				\
913 	  R##_c = Y##_c;					\
914 	  break;						\
915 								\
916 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
917 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
918 	  R##_s = _FP_NANSIGN_##fs;				\
919 	  R##_c = FP_CLS_NAN;					\
920 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
921 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ);	\
922 	  break;						\
923 								\
924 	default:						\
925 	  abort ();						\
926 	}							\
927     }								\
928   while (0)
929 
930 
931 /* Fused multiply-add.  The input values should be cooked.  */
932 
933 #define _FP_FMA(fs, wc, dwc, R, X, Y, Z)				\
934   do									\
935     {									\
936       FP_DECL_##fs (_FP_FMA_T);						\
937       _FP_FMA_T##_s = X##_s ^ Y##_s;					\
938       _FP_FMA_T##_e = X##_e + Y##_e + 1;				\
939       switch (_FP_CLS_COMBINE (X##_c, Y##_c))				\
940 	{								\
941 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):		\
942 	  switch (Z##_c)						\
943 	    {								\
944 	    case FP_CLS_INF:						\
945 	    case FP_CLS_NAN:						\
946 	      R##_s = Z##_s;						\
947 	      _FP_FRAC_COPY_##wc (R, Z);				\
948 	      R##_c = Z##_c;						\
949 	      break;							\
950 									\
951 	    case FP_CLS_ZERO:						\
952 	      R##_c = FP_CLS_NORMAL;					\
953 	      R##_s = _FP_FMA_T##_s;					\
954 	      R##_e = _FP_FMA_T##_e;					\
955 									\
956 	      _FP_MUL_MEAT_##fs (R, X, Y);				\
957 									\
958 	      if (_FP_FRAC_OVERP_##wc (fs, R))				\
959 		_FP_FRAC_SRS_##wc (R, 1, _FP_WFRACBITS_##fs);		\
960 	      else							\
961 		R##_e--;						\
962 	      break;							\
963 									\
964 	    case FP_CLS_NORMAL:;					\
965 	      _FP_FRAC_DECL_##dwc (_FP_FMA_TD);				\
966 	      _FP_FRAC_DECL_##dwc (_FP_FMA_ZD);				\
967 	      _FP_FRAC_DECL_##dwc (_FP_FMA_RD);				\
968 	      _FP_MUL_MEAT_DW_##fs (_FP_FMA_TD, X, Y);			\
969 	      R##_e = _FP_FMA_T##_e;					\
970 	      int _FP_FMA_tsh						\
971 		= _FP_FRAC_HIGHBIT_DW_##dwc (fs, _FP_FMA_TD) == 0;	\
972 	      _FP_FMA_T##_e -= _FP_FMA_tsh;				\
973 	      int _FP_FMA_ediff = _FP_FMA_T##_e - Z##_e;		\
974 	      if (_FP_FMA_ediff >= 0)					\
975 		{							\
976 		  int _FP_FMA_shift					\
977 		    = _FP_WFRACBITS_##fs - _FP_FMA_tsh - _FP_FMA_ediff;	\
978 		  if (_FP_FMA_shift <= -_FP_WFRACBITS_##fs)		\
979 		    _FP_FRAC_SET_##dwc (_FP_FMA_ZD, _FP_MINFRAC_##dwc);	\
980 		  else							\
981 		    {							\
982 		      _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);	\
983 		      if (_FP_FMA_shift < 0)				\
984 			_FP_FRAC_SRS_##dwc (_FP_FMA_ZD, -_FP_FMA_shift,	\
985 					    _FP_WFRACBITS_DW_##fs);	\
986 		      else if (_FP_FMA_shift > 0)			\
987 			_FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_FMA_shift);	\
988 		    }							\
989 		  R##_s = _FP_FMA_T##_s;				\
990 		  if (_FP_FMA_T##_s == Z##_s)				\
991 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_TD,		\
992 					_FP_FMA_ZD);			\
993 		  else							\
994 		    {							\
995 		      _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_TD,	\
996 					  _FP_FMA_ZD);			\
997 		      if (_FP_FRAC_NEGP_##dwc (_FP_FMA_RD))		\
998 			{						\
999 			  R##_s = Z##_s;				\
1000 			  _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,	\
1001 					      _FP_FMA_TD);		\
1002 			}						\
1003 		    }							\
1004 		}							\
1005 	      else							\
1006 		{							\
1007 		  R##_e = Z##_e;					\
1008 		  R##_s = Z##_s;					\
1009 		  _FP_FRAC_COPY_##dwc##_##wc (_FP_FMA_ZD, Z);		\
1010 		  _FP_FRAC_SLL_##dwc (_FP_FMA_ZD, _FP_WFRACBITS_##fs);	\
1011 		  int _FP_FMA_shift = -_FP_FMA_ediff - _FP_FMA_tsh;	\
1012 		  if (_FP_FMA_shift >= _FP_WFRACBITS_DW_##fs)		\
1013 		    _FP_FRAC_SET_##dwc (_FP_FMA_TD, _FP_MINFRAC_##dwc);	\
1014 		  else if (_FP_FMA_shift > 0)				\
1015 		    _FP_FRAC_SRS_##dwc (_FP_FMA_TD, _FP_FMA_shift,	\
1016 					_FP_WFRACBITS_DW_##fs);		\
1017 		  if (Z##_s == _FP_FMA_T##_s)				\
1018 		    _FP_FRAC_ADD_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1019 					_FP_FMA_TD);			\
1020 		  else							\
1021 		    _FP_FRAC_SUB_##dwc (_FP_FMA_RD, _FP_FMA_ZD,		\
1022 					_FP_FMA_TD);			\
1023 		}							\
1024 	      if (_FP_FRAC_ZEROP_##dwc (_FP_FMA_RD))			\
1025 		{							\
1026 		  if (_FP_FMA_T##_s == Z##_s)				\
1027 		    R##_s = Z##_s;					\
1028 		  else							\
1029 		    R##_s = (FP_ROUNDMODE == FP_RND_MINF);		\
1030 		  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);		\
1031 		  R##_c = FP_CLS_ZERO;					\
1032 		}							\
1033 	      else							\
1034 		{							\
1035 		  int _FP_FMA_rlz;					\
1036 		  _FP_FRAC_CLZ_##dwc (_FP_FMA_rlz, _FP_FMA_RD);		\
1037 		  _FP_FMA_rlz -= _FP_WFRACXBITS_DW_##fs;		\
1038 		  R##_e -= _FP_FMA_rlz;					\
1039 		  int _FP_FMA_shift = _FP_WFRACBITS_##fs - _FP_FMA_rlz;	\
1040 		  if (_FP_FMA_shift > 0)				\
1041 		    _FP_FRAC_SRS_##dwc (_FP_FMA_RD, _FP_FMA_shift,	\
1042 					_FP_WFRACBITS_DW_##fs);		\
1043 		  else if (_FP_FMA_shift < 0)				\
1044 		    _FP_FRAC_SLL_##dwc (_FP_FMA_RD, -_FP_FMA_shift);	\
1045 		  _FP_FRAC_COPY_##wc##_##dwc (R, _FP_FMA_RD);		\
1046 		  R##_c = FP_CLS_NORMAL;				\
1047 		}							\
1048 	      break;							\
1049 	    }								\
1050 	  goto done_fma;						\
1051 									\
1052 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1053 	  _FP_CHOOSENAN (fs, wc, _FP_FMA_T, X, Y, '*');			\
1054 	  break;							\
1055 									\
1056 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1057 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1058 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1059 	  _FP_FMA_T##_s = X##_s;					\
1060 									\
1061 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1062 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1063 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1064 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1065 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, X);				\
1066 	  _FP_FMA_T##_c = X##_c;					\
1067 	  break;							\
1068 									\
1069 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):		\
1070 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1071 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1072 	  _FP_FMA_T##_s = Y##_s;					\
1073 									\
1074 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):		\
1075 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):		\
1076 	  _FP_FRAC_COPY_##wc (_FP_FMA_T, Y);				\
1077 	  _FP_FMA_T##_c = Y##_c;					\
1078 	  break;							\
1079 									\
1080 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1081 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1082 	  _FP_FMA_T##_s = _FP_NANSIGN_##fs;				\
1083 	  _FP_FMA_T##_c = FP_CLS_NAN;					\
1084 	  _FP_FRAC_SET_##wc (_FP_FMA_T, _FP_NANFRAC_##fs);		\
1085 	  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_IMZ_FMA);	\
1086 	  break;							\
1087 									\
1088 	default:							\
1089 	  abort ();							\
1090 	}								\
1091 									\
1092       /* T = X * Y is zero, infinity or NaN.  */			\
1093       switch (_FP_CLS_COMBINE (_FP_FMA_T##_c, Z##_c))			\
1094 	{								\
1095 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):			\
1096 	  _FP_CHOOSENAN (fs, wc, R, _FP_FMA_T, Z, '+');			\
1097 	  break;							\
1098 									\
1099 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):		\
1100 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):			\
1101 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):			\
1102 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):		\
1103 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):			\
1104 	  R##_s = _FP_FMA_T##_s;					\
1105 	  _FP_FRAC_COPY_##wc (R, _FP_FMA_T);				\
1106 	  R##_c = _FP_FMA_T##_c;					\
1107 	  break;							\
1108 									\
1109 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):			\
1110 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):			\
1111 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):		\
1112 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):			\
1113 	  R##_s = Z##_s;						\
1114 	  _FP_FRAC_COPY_##wc (R, Z);					\
1115 	  R##_c = Z##_c;						\
1116 	  break;							\
1117 									\
1118 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):			\
1119 	  if (_FP_FMA_T##_s == Z##_s)					\
1120 	    {								\
1121 	      R##_s = Z##_s;						\
1122 	      _FP_FRAC_COPY_##wc (R, Z);				\
1123 	      R##_c = Z##_c;						\
1124 	    }								\
1125 	  else								\
1126 	    {								\
1127 	      R##_s = _FP_NANSIGN_##fs;					\
1128 	      R##_c = FP_CLS_NAN;					\
1129 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1130 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_ISI);	\
1131 	    }								\
1132 	  break;							\
1133 									\
1134 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):		\
1135 	  if (_FP_FMA_T##_s == Z##_s)					\
1136 	    R##_s = Z##_s;						\
1137 	  else								\
1138 	    R##_s = (FP_ROUNDMODE == FP_RND_MINF);			\
1139 	  _FP_FRAC_COPY_##wc (R, Z);					\
1140 	  R##_c = Z##_c;						\
1141 	  break;							\
1142 									\
1143 	default:							\
1144 	  abort ();							\
1145 	}								\
1146     done_fma: ;								\
1147     }									\
1148   while (0)
1149 
1150 
1151 /* Main division routine.  The input values should be cooked.  */
1152 
1153 #define _FP_DIV(fs, wc, R, X, Y)				\
1154   do								\
1155     {								\
1156       R##_s = X##_s ^ Y##_s;					\
1157       R##_e = X##_e - Y##_e;					\
1158       switch (_FP_CLS_COMBINE (X##_c, Y##_c))			\
1159 	{							\
1160 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NORMAL):	\
1161 	  R##_c = FP_CLS_NORMAL;				\
1162 								\
1163 	  _FP_DIV_MEAT_##fs (R, X, Y);				\
1164 	  break;						\
1165 								\
1166 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NAN):		\
1167 	  _FP_CHOOSENAN (fs, wc, R, X, Y, '/');			\
1168 	  break;						\
1169 								\
1170 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_NORMAL):	\
1171 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_INF):		\
1172 	case _FP_CLS_COMBINE (FP_CLS_NAN, FP_CLS_ZERO):		\
1173 	  R##_s = X##_s;					\
1174 	  _FP_FRAC_COPY_##wc (R, X);				\
1175 	  R##_c = X##_c;					\
1176 	  break;						\
1177 								\
1178 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_NAN):	\
1179 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NAN):		\
1180 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NAN):		\
1181 	  R##_s = Y##_s;					\
1182 	  _FP_FRAC_COPY_##wc (R, Y);				\
1183 	  R##_c = Y##_c;					\
1184 	  break;						\
1185 								\
1186 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_INF):	\
1187 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_INF):		\
1188 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_NORMAL):	\
1189 	  R##_c = FP_CLS_ZERO;					\
1190 	  break;						\
1191 								\
1192 	case _FP_CLS_COMBINE (FP_CLS_NORMAL, FP_CLS_ZERO):	\
1193 	  FP_SET_EXCEPTION (FP_EX_DIVZERO);			\
1194 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_ZERO):		\
1195 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_NORMAL):	\
1196 	  R##_c = FP_CLS_INF;					\
1197 	  break;						\
1198 								\
1199 	case _FP_CLS_COMBINE (FP_CLS_INF, FP_CLS_INF):		\
1200 	case _FP_CLS_COMBINE (FP_CLS_ZERO, FP_CLS_ZERO):	\
1201 	  R##_s = _FP_NANSIGN_##fs;				\
1202 	  R##_c = FP_CLS_NAN;					\
1203 	  _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);		\
1204 	  FP_SET_EXCEPTION (FP_EX_INVALID			\
1205 			    | (X##_c == FP_CLS_INF		\
1206 			       ? FP_EX_INVALID_IDI		\
1207 			       : FP_EX_INVALID_ZDZ));		\
1208 	  break;						\
1209 								\
1210 	default:						\
1211 	  abort ();						\
1212 	}							\
1213     }								\
1214   while (0)
1215 
1216 
1217 /* Helper for comparisons.  EX is 0 not to raise exceptions, 1 to
1218    raise exceptions for signaling NaN operands, 2 to raise exceptions
1219    for all NaN operands.  Conditionals are organized to allow the
1220    compiler to optimize away code based on the value of EX.  */
1221 
1222 #define _FP_CMP_CHECK_NAN(fs, wc, X, Y, ex)				\
1223   do									\
1224     {									\
1225       /* The arguments are unordered, which may or may not result in	\
1226 	 an exception.  */						\
1227       if (ex)								\
1228 	{								\
1229 	  /* At least some cases of unordered arguments result in	\
1230 	     exceptions; check whether this is one.  */			\
1231 	  if (FP_EX_INVALID_SNAN || FP_EX_INVALID_VC)			\
1232 	    {								\
1233 	      /* Check separately for each case of "invalid"		\
1234 		 exceptions.  */					\
1235 	      if ((ex) == 2)						\
1236 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_VC);	\
1237 	      if (_FP_ISSIGNAN (fs, wc, X)				\
1238 		  || _FP_ISSIGNAN (fs, wc, Y))				\
1239 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SNAN);	\
1240 	    }								\
1241 	  /* Otherwise, we only need to check whether to raise an	\
1242 	     exception, not which case or cases it is.  */		\
1243 	  else if ((ex) == 2						\
1244 		   || _FP_ISSIGNAN (fs, wc, X)				\
1245 		   || _FP_ISSIGNAN (fs, wc, Y))				\
1246 	    FP_SET_EXCEPTION (FP_EX_INVALID);				\
1247 	}								\
1248     }									\
1249   while (0)
1250 
1251 /* Main differential comparison routine.  The inputs should be raw not
1252    cooked.  The return is -1, 0, 1 for normal values, UN
1253    otherwise.  */
1254 
1255 #define _FP_CMP(fs, wc, ret, X, Y, un, ex)				\
1256   do									\
1257     {									\
1258       /* NANs are unordered.  */					\
1259       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1260 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1261 	{								\
1262 	  (ret) = (un);							\
1263 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1264 	}								\
1265       else								\
1266 	{								\
1267 	  int _FP_CMP_is_zero_x;					\
1268 	  int _FP_CMP_is_zero_y;					\
1269 									\
1270 	  _FP_CHECK_FLUSH_ZERO (fs, wc, X);				\
1271 	  _FP_CHECK_FLUSH_ZERO (fs, wc, Y);				\
1272 									\
1273 	  _FP_CMP_is_zero_x						\
1274 	    = (!X##_e && _FP_FRAC_ZEROP_##wc (X)) ? 1 : 0;		\
1275 	  _FP_CMP_is_zero_y						\
1276 	    = (!Y##_e && _FP_FRAC_ZEROP_##wc (Y)) ? 1 : 0;		\
1277 									\
1278 	  if (_FP_CMP_is_zero_x && _FP_CMP_is_zero_y)			\
1279 	    (ret) = 0;							\
1280 	  else if (_FP_CMP_is_zero_x)					\
1281 	    (ret) = Y##_s ? 1 : -1;					\
1282 	  else if (_FP_CMP_is_zero_y)					\
1283 	    (ret) = X##_s ? -1 : 1;					\
1284 	  else if (X##_s != Y##_s)					\
1285 	    (ret) = X##_s ? -1 : 1;					\
1286 	  else if (X##_e > Y##_e)					\
1287 	    (ret) = X##_s ? -1 : 1;					\
1288 	  else if (X##_e < Y##_e)					\
1289 	    (ret) = X##_s ? 1 : -1;					\
1290 	  else if (_FP_FRAC_GT_##wc (X, Y))				\
1291 	    (ret) = X##_s ? -1 : 1;					\
1292 	  else if (_FP_FRAC_GT_##wc (Y, X))				\
1293 	    (ret) = X##_s ? 1 : -1;					\
1294 	  else								\
1295 	    (ret) = 0;							\
1296 	}								\
1297     }									\
1298   while (0)
1299 
1300 
1301 /* Simplification for strict equality.  */
1302 
1303 #define _FP_CMP_EQ(fs, wc, ret, X, Y, ex)				\
1304   do									\
1305     {									\
1306       /* NANs are unordered.  */					\
1307       if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1308 	  || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y)))	\
1309 	{								\
1310 	  (ret) = 1;							\
1311 	  _FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));			\
1312 	}								\
1313       else								\
1314 	{								\
1315 	  _FP_CHECK_FLUSH_ZERO (fs, wc, X);				\
1316 	  _FP_CHECK_FLUSH_ZERO (fs, wc, Y);				\
1317 									\
1318 	  (ret) = !(X##_e == Y##_e					\
1319 		    && _FP_FRAC_EQ_##wc (X, Y)				\
1320 		    && (X##_s == Y##_s					\
1321 			|| (!X##_e && _FP_FRAC_ZEROP_##wc (X))));	\
1322 	}								\
1323     }									\
1324   while (0)
1325 
1326 /* Version to test unordered.  */
1327 
1328 #define _FP_CMP_UNORD(fs, wc, ret, X, Y, ex)				\
1329   do									\
1330     {									\
1331       (ret) = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (X))	\
1332 	       || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc (Y))); \
1333       if (ret)								\
1334 	_FP_CMP_CHECK_NAN (fs, wc, X, Y, (ex));				\
1335     }									\
1336   while (0)
1337 
1338 /* Main square root routine.  The input value should be cooked.  */
1339 
1340 #define _FP_SQRT(fs, wc, R, X)						\
1341   do									\
1342     {									\
1343       _FP_FRAC_DECL_##wc (_FP_SQRT_T);					\
1344       _FP_FRAC_DECL_##wc (_FP_SQRT_S);					\
1345       _FP_W_TYPE _FP_SQRT_q;						\
1346       switch (X##_c)							\
1347 	{								\
1348 	case FP_CLS_NAN:						\
1349 	  _FP_FRAC_COPY_##wc (R, X);					\
1350 	  R##_s = X##_s;						\
1351 	  R##_c = FP_CLS_NAN;						\
1352 	  break;							\
1353 	case FP_CLS_INF:						\
1354 	  if (X##_s)							\
1355 	    {								\
1356 	      R##_s = _FP_NANSIGN_##fs;					\
1357 	      R##_c = FP_CLS_NAN; /* NAN */				\
1358 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1359 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1360 	    }								\
1361 	  else								\
1362 	    {								\
1363 	      R##_s = 0;						\
1364 	      R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */		\
1365 	    }								\
1366 	  break;							\
1367 	case FP_CLS_ZERO:						\
1368 	  R##_s = X##_s;						\
1369 	  R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */			\
1370 	  break;							\
1371 	case FP_CLS_NORMAL:						\
1372 	  R##_s = 0;							\
1373 	  if (X##_s)							\
1374 	    {								\
1375 	      R##_c = FP_CLS_NAN; /* NAN */				\
1376 	      R##_s = _FP_NANSIGN_##fs;					\
1377 	      _FP_FRAC_SET_##wc (R, _FP_NANFRAC_##fs);			\
1378 	      FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_SQRT);	\
1379 	      break;							\
1380 	    }								\
1381 	  R##_c = FP_CLS_NORMAL;					\
1382 	  if (X##_e & 1)						\
1383 	    _FP_FRAC_SLL_##wc (X, 1);					\
1384 	  R##_e = X##_e >> 1;						\
1385 	  _FP_FRAC_SET_##wc (_FP_SQRT_S, _FP_ZEROFRAC_##wc);		\
1386 	  _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);			\
1387 	  _FP_SQRT_q = _FP_OVERFLOW_##fs >> 1;				\
1388 	  _FP_SQRT_MEAT_##wc (R, _FP_SQRT_S, _FP_SQRT_T, X,		\
1389 			      _FP_SQRT_q);				\
1390 	}								\
1391     }									\
1392   while (0)
1393 
1394 /* Convert from FP to integer.  Input is raw.  */
1395 
1396 /* RSIGNED can have following values:
1397    0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
1398        the result is either 0 or (2^rsize)-1 depending on the sign in such
1399        case.
1400    1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1401        NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1402        depending on the sign in such case.
1403    2:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not,
1404        NV is set plus the result is reduced modulo 2^rsize.
1405    -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
1406        set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1
1407        depending on the sign in such case.  */
1408 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)			\
1409   do									\
1410     {									\
1411       if (X##_e < _FP_EXPBIAS_##fs)					\
1412 	{								\
1413 	  (r) = 0;							\
1414 	  if (X##_e == 0)						\
1415 	    {								\
1416 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1417 		{							\
1418 		  if (!FP_DENORM_ZERO)					\
1419 		    FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1420 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1421 		}							\
1422 	    }								\
1423 	  else								\
1424 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1425 	}								\
1426       else if ((rsigned) == 2						\
1427 	       && (X##_e						\
1428 		   >= ((_FP_EXPMAX_##fs					\
1429 			< _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1) \
1430 		       ? _FP_EXPMAX_##fs				\
1431 		       : _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs + (rsize) - 1))) \
1432 	{								\
1433 	  /* Overflow resulting in 0.  */				\
1434 	  (r) = 0;							\
1435 	  FP_SET_EXCEPTION (FP_EX_INVALID				\
1436 			    | FP_EX_INVALID_CVI				\
1437 			    | ((FP_EX_INVALID_SNAN			\
1438 				&& _FP_ISSIGNAN (fs, wc, X))		\
1439 			       ? FP_EX_INVALID_SNAN			\
1440 			       : 0));					\
1441 	}								\
1442       else if ((rsigned) != 2						\
1443 	       && (X##_e >= (_FP_EXPMAX_##fs < _FP_EXPBIAS_##fs + (rsize) \
1444 			     ? _FP_EXPMAX_##fs				\
1445 			     : (_FP_EXPBIAS_##fs + (rsize)		\
1446 				- ((rsigned) > 0 || X##_s)))		\
1447 		   || (!(rsigned) && X##_s)))				\
1448 	{								\
1449 	  /* Overflow or converting to the most negative integer.  */	\
1450 	  if (rsigned)							\
1451 	    {								\
1452 	      (r) = 1;							\
1453 	      (r) <<= (rsize) - 1;					\
1454 	      (r) -= 1 - X##_s;						\
1455 	    }								\
1456 	  else								\
1457 	    {								\
1458 	      (r) = 0;							\
1459 	      if (!X##_s)						\
1460 		(r) = ~(r);						\
1461 	    }								\
1462 									\
1463 	  if (_FP_EXPBIAS_##fs + (rsize) - 1 < _FP_EXPMAX_##fs		\
1464 	      && (rsigned)						\
1465 	      && X##_s							\
1466 	      && X##_e == _FP_EXPBIAS_##fs + (rsize) - 1)		\
1467 	    {								\
1468 	      /* Possibly converting to most negative integer; check the \
1469 		 mantissa.  */						\
1470 	      int _FP_TO_INT_inexact = 0;				\
1471 	      (void) ((_FP_FRACBITS_##fs > (rsize))			\
1472 		      ? ({						\
1473 			  _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,	\
1474 					      _FP_FRACBITS_##fs - (rsize), \
1475 					      _FP_FRACBITS_##fs);	\
1476 			  0;						\
1477 			})						\
1478 		      : 0);						\
1479 	      if (!_FP_FRAC_ZEROP_##wc (X))				\
1480 		FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1481 	      else if (_FP_TO_INT_inexact)				\
1482 		FP_SET_EXCEPTION (FP_EX_INEXACT);			\
1483 	    }								\
1484 	  else								\
1485 	    FP_SET_EXCEPTION (FP_EX_INVALID				\
1486 			      | FP_EX_INVALID_CVI			\
1487 			      | ((FP_EX_INVALID_SNAN			\
1488 				  && _FP_ISSIGNAN (fs, wc, X))		\
1489 				 ? FP_EX_INVALID_SNAN			\
1490 				 : 0));					\
1491 	}								\
1492       else								\
1493 	{								\
1494 	  int _FP_TO_INT_inexact = 0;					\
1495 	  _FP_FRAC_HIGH_RAW_##fs (X) |= _FP_IMPLBIT_##fs;		\
1496 	  if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1)	\
1497 	    {								\
1498 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1499 	      (r) <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \
1500 	    }								\
1501 	  else								\
1502 	    {								\
1503 	      _FP_FRAC_SRST_##wc (X, _FP_TO_INT_inexact,		\
1504 				  (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \
1505 				   - X##_e),				\
1506 				  _FP_FRACBITS_##fs);			\
1507 	      _FP_FRAC_ASSEMBLE_##wc ((r), X, (rsize));			\
1508 	    }								\
1509 	  if ((rsigned) && X##_s)					\
1510 	    (r) = -(r);							\
1511 	  if ((rsigned) == 2 && X##_e >= _FP_EXPBIAS_##fs + (rsize) - 1) \
1512 	    {								\
1513 	      /* Overflow or converting to the most negative integer.  */ \
1514 	      if (X##_e > _FP_EXPBIAS_##fs + (rsize) - 1		\
1515 		  || !X##_s						\
1516 		  || (r) != (((typeof (r)) 1) << ((rsize) - 1)))	\
1517 		{							\
1518 		  _FP_TO_INT_inexact = 0;				\
1519 		  FP_SET_EXCEPTION (FP_EX_INVALID | FP_EX_INVALID_CVI);	\
1520 		}							\
1521 	    }								\
1522 	  if (_FP_TO_INT_inexact)					\
1523 	    FP_SET_EXCEPTION (FP_EX_INEXACT);				\
1524 	}								\
1525     }									\
1526   while (0)
1527 
1528 /* Convert integer to fp.  Output is raw.  RTYPE is unsigned even if
1529    input is signed.  */
1530 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)			\
1531   do									\
1532     {									\
1533       if (r)								\
1534 	{								\
1535 	  rtype _FP_FROM_INT_ur;					\
1536 									\
1537 	  if ((X##_s = ((r) < 0)))					\
1538 	    (r) = -(rtype) (r);						\
1539 									\
1540 	  _FP_FROM_INT_ur = (rtype) (r);				\
1541 	  (void) (((rsize) <= _FP_W_TYPE_SIZE)				\
1542 		  ? ({							\
1543 		      int _FP_FROM_INT_lz;				\
1544 		      __FP_CLZ (_FP_FROM_INT_lz,			\
1545 				(_FP_W_TYPE) _FP_FROM_INT_ur);		\
1546 		      X##_e = (_FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1	\
1547 			       - _FP_FROM_INT_lz);			\
1548 		    })							\
1549 		  : (((rsize) <= 2 * _FP_W_TYPE_SIZE)			\
1550 		     ? ({						\
1551 			 int _FP_FROM_INT_lz;				\
1552 			 __FP_CLZ_2 (_FP_FROM_INT_lz,			\
1553 				     (_FP_W_TYPE) (_FP_FROM_INT_ur	\
1554 						   >> _FP_W_TYPE_SIZE), \
1555 				     (_FP_W_TYPE) _FP_FROM_INT_ur);	\
1556 			 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \
1557 				  - _FP_FROM_INT_lz);			\
1558 		       })						\
1559 		     : (abort (), 0)));					\
1560 									\
1561 	  if ((rsize) - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs		\
1562 	      && X##_e >= _FP_EXPMAX_##fs)				\
1563 	    {								\
1564 	      /* Exponent too big; overflow to infinity.  (May also	\
1565 		 happen after rounding below.)  */			\
1566 	      _FP_OVERFLOW_SEMIRAW (fs, wc, X);				\
1567 	      goto pack_semiraw;					\
1568 	    }								\
1569 									\
1570 	  if ((rsize) <= _FP_FRACBITS_##fs				\
1571 	      || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs)		\
1572 	    {								\
1573 	      /* Exactly representable; shift left.  */			\
1574 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1575 	      if (_FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1 - X##_e > 0)	\
1576 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1577 				       + _FP_FRACBITS_##fs - 1 - X##_e)); \
1578 	    }								\
1579 	  else								\
1580 	    {								\
1581 	      /* More bits in integer than in floating type; need to	\
1582 		 round.  */						\
1583 	      if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e)	\
1584 		_FP_FROM_INT_ur						\
1585 		  = ((_FP_FROM_INT_ur >> (X##_e - _FP_EXPBIAS_##fs	\
1586 					  - _FP_WFRACBITS_##fs + 1))	\
1587 		     | ((_FP_FROM_INT_ur				\
1588 			 << ((rsize) - (X##_e - _FP_EXPBIAS_##fs	\
1589 					- _FP_WFRACBITS_##fs + 1)))	\
1590 			!= 0));						\
1591 	      _FP_FRAC_DISASSEMBLE_##wc (X, _FP_FROM_INT_ur, (rsize));	\
1592 	      if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \
1593 		_FP_FRAC_SLL_##wc (X, (_FP_EXPBIAS_##fs			\
1594 				       + _FP_WFRACBITS_##fs - 1 - X##_e)); \
1595 	      _FP_FRAC_HIGH_##fs (X) &= ~(_FP_W_TYPE) _FP_IMPLBIT_SH_##fs; \
1596 	    pack_semiraw:						\
1597 	      _FP_PACK_SEMIRAW (fs, wc, X);				\
1598 	    }								\
1599 	}								\
1600       else								\
1601 	{								\
1602 	  X##_s = 0;							\
1603 	  X##_e = 0;							\
1604 	  _FP_FRAC_SET_##wc (X, _FP_ZEROFRAC_##wc);			\
1605 	}								\
1606     }									\
1607   while (0)
1608 
1609 
1610 /* Extend from a narrower floating-point format to a wider one.  Input
1611    and output are raw.  */
1612 #define FP_EXTEND(dfs, sfs, dwc, swc, D, S)				\
1613   do									\
1614     {									\
1615       if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs			\
1616 	  || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs			\
1617 	      < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs)			\
1618 	  || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \
1619 	      && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs))		\
1620 	abort ();							\
1621       D##_s = S##_s;							\
1622       _FP_FRAC_COPY_##dwc##_##swc (D, S);				\
1623       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1624 	{								\
1625 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1626 	  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \
1627 	}								\
1628       else								\
1629 	{								\
1630 	  if (S##_e == 0)						\
1631 	    {								\
1632 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
1633 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1634 		D##_e = 0;						\
1635 	      else if (_FP_EXPBIAS_##dfs				\
1636 		       < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1)	\
1637 		{							\
1638 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1639 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1640 					  - _FP_FRACBITS_##sfs));	\
1641 		  D##_e = 0;						\
1642 		  if (FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW)		\
1643 		    FP_SET_EXCEPTION (FP_EX_UNDERFLOW);			\
1644 		}							\
1645 	      else							\
1646 		{							\
1647 		  int FP_EXTEND_lz;					\
1648 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1649 		  _FP_FRAC_CLZ_##swc (FP_EXTEND_lz, S);			\
1650 		  _FP_FRAC_SLL_##dwc (D,				\
1651 				      FP_EXTEND_lz + _FP_FRACBITS_##dfs	\
1652 				      - _FP_FRACTBITS_##sfs);		\
1653 		  D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1	\
1654 			   + _FP_FRACXBITS_##sfs - FP_EXTEND_lz);	\
1655 		}							\
1656 	    }								\
1657 	  else								\
1658 	    {								\
1659 	      D##_e = _FP_EXPMAX_##dfs;					\
1660 	      if (!_FP_FRAC_ZEROP_##swc (S))				\
1661 		{							\
1662 		  if (_FP_FRAC_SNANP (sfs, S))				\
1663 		    FP_SET_EXCEPTION (FP_EX_INVALID			\
1664 				      | FP_EX_INVALID_SNAN);		\
1665 		  _FP_FRAC_SLL_##dwc (D, (_FP_FRACBITS_##dfs		\
1666 					  - _FP_FRACBITS_##sfs));	\
1667 		  _FP_SETQNAN (dfs, dwc, D);				\
1668 		}							\
1669 	    }								\
1670 	}								\
1671     }									\
1672   while (0)
1673 
1674 /* Truncate from a wider floating-point format to a narrower one.
1675    Input and output are semi-raw.  */
1676 #define FP_TRUNC(dfs, sfs, dwc, swc, D, S)				\
1677   do									\
1678     {									\
1679       if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs			\
1680 	  || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \
1681 	      && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs))		\
1682 	abort ();							\
1683       D##_s = S##_s;							\
1684       if (_FP_EXP_NORMAL (sfs, swc, S))					\
1685 	{								\
1686 	  D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs;	\
1687 	  if (D##_e >= _FP_EXPMAX_##dfs)				\
1688 	    _FP_OVERFLOW_SEMIRAW (dfs, dwc, D);				\
1689 	  else								\
1690 	    {								\
1691 	      if (D##_e <= 0)						\
1692 		{							\
1693 		  if (D##_e < 1 - _FP_FRACBITS_##dfs)			\
1694 		    {							\
1695 		      _FP_FRAC_SET_##swc (S, _FP_ZEROFRAC_##swc);	\
1696 		      _FP_FRAC_LOW_##swc (S) |= 1;			\
1697 		    }							\
1698 		  else							\
1699 		    {							\
1700 		      _FP_FRAC_HIGH_##sfs (S) |= _FP_IMPLBIT_SH_##sfs;	\
1701 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
1702 					      - _FP_WFRACBITS_##dfs	\
1703 					      + 1 - D##_e),		\
1704 					  _FP_WFRACBITS_##sfs);		\
1705 		    }							\
1706 		  D##_e = 0;						\
1707 		}							\
1708 	      else							\
1709 		_FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs		\
1710 					- _FP_WFRACBITS_##dfs),		\
1711 				    _FP_WFRACBITS_##sfs);		\
1712 	      _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
1713 	    }								\
1714 	}								\
1715       else								\
1716 	{								\
1717 	  if (S##_e == 0)						\
1718 	    {								\
1719 	      _FP_CHECK_FLUSH_ZERO (sfs, swc, S);			\
1720 	      D##_e = 0;						\
1721 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1722 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
1723 	      else							\
1724 		{							\
1725 		  FP_SET_EXCEPTION (FP_EX_DENORM);			\
1726 		  if (_FP_EXPBIAS_##sfs					\
1727 		      < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1)	\
1728 		    {							\
1729 		      _FP_FRAC_SRS_##swc (S, (_FP_WFRACBITS_##sfs	\
1730 					      - _FP_WFRACBITS_##dfs),	\
1731 					  _FP_WFRACBITS_##sfs);		\
1732 		      _FP_FRAC_COPY_##dwc##_##swc (D, S);		\
1733 		    }							\
1734 		  else							\
1735 		    {							\
1736 		      _FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);	\
1737 		      _FP_FRAC_LOW_##dwc (D) |= 1;			\
1738 		    }							\
1739 		}							\
1740 	    }								\
1741 	  else								\
1742 	    {								\
1743 	      D##_e = _FP_EXPMAX_##dfs;					\
1744 	      if (_FP_FRAC_ZEROP_##swc (S))				\
1745 		_FP_FRAC_SET_##dwc (D, _FP_ZEROFRAC_##dwc);		\
1746 	      else							\
1747 		{							\
1748 		  _FP_CHECK_SIGNAN_SEMIRAW (sfs, swc, S);		\
1749 		  _FP_FRAC_SRL_##swc (S, (_FP_WFRACBITS_##sfs		\
1750 					  - _FP_WFRACBITS_##dfs));	\
1751 		  _FP_FRAC_COPY_##dwc##_##swc (D, S);			\
1752 		  /* Semi-raw NaN must have all workbits cleared.  */	\
1753 		  _FP_FRAC_LOW_##dwc (D)				\
1754 		    &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1);		\
1755 		  _FP_SETQNAN_SEMIRAW (dfs, dwc, D);			\
1756 		}							\
1757 	    }								\
1758 	}								\
1759     }									\
1760   while (0)
1761 
1762 /* Helper primitives.  */
1763 
1764 /* Count leading zeros in a word.  */
1765 
1766 #ifndef __FP_CLZ
1767 /* GCC 3.4 and later provide the builtins for us.  */
1768 # define __FP_CLZ(r, x)							\
1769   do									\
1770     {									\
1771       if (sizeof (_FP_W_TYPE) == sizeof (unsigned int))			\
1772 	(r) = __builtin_clz (x);					\
1773       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long))		\
1774 	(r) = __builtin_clzl (x);					\
1775       else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long))	\
1776 	(r) = __builtin_clzll (x);					\
1777       else								\
1778 	abort ();							\
1779     }									\
1780   while (0)
1781 #endif /* ndef __FP_CLZ */
1782 
1783 #define _FP_DIV_HELP_imm(q, r, n, d)		\
1784   do						\
1785     {						\
1786       (q) = (n) / (d), (r) = (n) % (d);		\
1787     }						\
1788   while (0)
1789 
1790 
1791 /* A restoring bit-by-bit division primitive.  */
1792 
1793 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y)				\
1794   do									\
1795     {									\
1796       int _FP_DIV_MEAT_N_loop_count = _FP_WFRACBITS_##fs;		\
1797       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_u);			\
1798       _FP_FRAC_DECL_##wc (_FP_DIV_MEAT_N_loop_v);			\
1799       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_u, X);			\
1800       _FP_FRAC_COPY_##wc (_FP_DIV_MEAT_N_loop_v, Y);			\
1801       _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc);				\
1802       /* Normalize _FP_DIV_MEAT_N_LOOP_U and _FP_DIV_MEAT_N_LOOP_V.  */	\
1803       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, _FP_WFRACXBITS_##fs);	\
1804       _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_v, _FP_WFRACXBITS_##fs);	\
1805       /* First round.  Since the operands are normalized, either the	\
1806 	 first or second bit will be set in the fraction.  Produce a	\
1807 	 normalized result by checking which and adjusting the loop	\
1808 	 count and exponent accordingly.  */				\
1809       if (_FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u, _FP_DIV_MEAT_N_loop_v))	\
1810 	{								\
1811 	  _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
1812 			     _FP_DIV_MEAT_N_loop_u,			\
1813 			     _FP_DIV_MEAT_N_loop_v);			\
1814 	  _FP_FRAC_LOW_##wc (R) |= 1;					\
1815 	  _FP_DIV_MEAT_N_loop_count--;					\
1816 	}								\
1817       else								\
1818 	R##_e--;							\
1819       /* Subsequent rounds.  */						\
1820       do								\
1821 	{								\
1822 	  int _FP_DIV_MEAT_N_loop_msb					\
1823 	    = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (_FP_DIV_MEAT_N_loop_u) < 0; \
1824 	  _FP_FRAC_SLL_##wc (_FP_DIV_MEAT_N_loop_u, 1);			\
1825 	  _FP_FRAC_SLL_##wc (R, 1);					\
1826 	  if (_FP_DIV_MEAT_N_loop_msb					\
1827 	      || _FP_FRAC_GE_1 (_FP_DIV_MEAT_N_loop_u,			\
1828 				_FP_DIV_MEAT_N_loop_v))			\
1829 	    {								\
1830 	      _FP_FRAC_SUB_##wc (_FP_DIV_MEAT_N_loop_u,			\
1831 				 _FP_DIV_MEAT_N_loop_u,			\
1832 				 _FP_DIV_MEAT_N_loop_v);		\
1833 	      _FP_FRAC_LOW_##wc (R) |= 1;				\
1834 	    }								\
1835 	}								\
1836       while (--_FP_DIV_MEAT_N_loop_count > 0);				\
1837       /* If there's anything left in _FP_DIV_MEAT_N_LOOP_U, the result	\
1838 	 is inexact.  */						\
1839       _FP_FRAC_LOW_##wc (R)						\
1840 	|= !_FP_FRAC_ZEROP_##wc (_FP_DIV_MEAT_N_loop_u);		\
1841     }									\
1842   while (0)
1843 
1844 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y)
1845 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y)
1846 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y)  _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y)
1847