BXE: min->MIN, plus an spatch
[akaros.git] / kern / arch / riscv / softfloat.c
1 \r
2 /*============================================================================\r
3 \r
4 This C source file is part of the SoftFloat IEC/IEEE Floating-point Arithmetic\r
5 Package, Release 2b.\r
6 \r
7 Written by John R. Hauser.  This work was made possible in part by the\r
8 International Computer Science Institute, located at Suite 600, 1947 Center\r
9 Street, Berkeley, California 94704.  Funding was partially provided by the\r
10 National Science Foundation under grant MIP-9311980.  The original version\r
11 of this code was written as part of a project to build a fixed-point vector\r
12 processor in collaboration with the University of California at Berkeley,\r
13 overseen by Profs. Nelson Morgan and John Wawrzynek.  More information\r
14 is available through the Web page `http://www.cs.berkeley.edu/~jhauser/\r
15 arithmetic/SoftFloat.html'.\r
16 \r
17 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has\r
18 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES\r
19 RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS\r
20 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,\r
21 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE\r
22 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE\r
23 INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR\r
24 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.\r
25 \r
26 Derivative works are acceptable, even for commercial purposes, so long as\r
27 (1) the source code for the derivative work includes prominent notice that\r
28 the work is derivative, and (2) the source code includes prominent notice with\r
29 these four paragraphs for those parts of this code that are retained.\r
30 \r
31 =============================================================================*/\r
32 \r
33 //#include "milieu.h"\r
34 #include "softfloat.h"\r
35 \r
36 void softfloat_init(softfloat_t* sf)\r
37 {\r
38   sf->float_detect_tininess = float_tininess_before_rounding;\r
39   sf->float_rounding_mode = float_round_nearest_even;\r
40   sf->float_exception_flags = 0;\r
41   #ifdef FLOATX80\r
42     sf->floatx80_rounding_precision = 80;\r
43   #endif\r
44 };\r
45 \r
46 \r
47 /*----------------------------------------------------------------------------\r
48 | Floating-point rounding mode, extended double-precision rounding precision,\r
49 | and exception flags.\r
50 *----------------------------------------------------------------------------*/\r
51 //int8_t sf->float_rounding_mode = float_round_nearest_even;\r
52 //int8_t sf->float_exception_flags = 0;\r
53 #ifdef FLOATX80\r
54 //int sf->floatx80_rounding_precision = 80;\r
55 #endif\r
56 \r
57 /*----------------------------------------------------------------------------\r
58 | Primitive arithmetic functions, including multi-word arithmetic, and\r
59 | division and square root approximations.  (Can be specialized to target if\r
60 | desired.)\r
61 *----------------------------------------------------------------------------*/\r
62 #include "softfloat-macros.h"\r
63 \r
64 /*----------------------------------------------------------------------------\r
65 | Functions and definitions to determine:  (1) whether tininess for underflow\r
66 | is detected before or after rounding by default, (2) what (if anything)\r
67 | happens when exceptions are raised, (3) how signaling NaNs are distinguished\r
68 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs\r
69 | are propagated from function inputs to output.  These details are target-\r
70 | specific.\r
71 *----------------------------------------------------------------------------*/\r
72 #include "softfloat-specialize.h"\r
73 \r
74 /*----------------------------------------------------------------------------\r
75 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6\r
76 | and 7, and returns the properly rounded 32-bit integer corresponding to the\r
77 | input.  If `zSign' is 1, the input is negated before being converted to an\r
78 | integer.  Bit 63 of `absZ' must be zero.  Ordinarily, the fixed-point input\r
79 | is simply rounded to an integer, with the inexact exception raised if the\r
80 | input cannot be represented exactly as an integer.  However, if the fixed-\r
81 | point input is too large, the invalid exception is raised and the largest\r
82 | positive or negative integer is returned.\r
83 *----------------------------------------------------------------------------*/\r
84 \r
85 int32_t roundAndPackInt32( softfloat_t* sf, flag zSign, bits64 absZ )\r
86 {\r
87     int8_t roundingMode;\r
88     flag roundNearestEven;\r
89     int8_t roundIncrement, roundBits;\r
90     int32_t z;\r
91 \r
92     roundingMode = sf->float_rounding_mode;\r
93     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
94     roundIncrement = 0x40;\r
95     if ( ! roundNearestEven ) {\r
96         if ( roundingMode == float_round_to_zero ) {\r
97             roundIncrement = 0;\r
98         }\r
99         else {\r
100             roundIncrement = 0x7F;\r
101             if ( zSign ) {\r
102                 if ( roundingMode == float_round_up ) roundIncrement = 0;\r
103             }\r
104             else {\r
105                 if ( roundingMode == float_round_down ) roundIncrement = 0;\r
106             }\r
107         }\r
108     }\r
109     roundBits = absZ & 0x7F;\r
110     absZ = ( absZ + roundIncrement )>>7;\r
111     absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );\r
112     z = absZ;\r
113     if ( zSign ) z = - z;\r
114     if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {\r
115         float_raise( sf, float_flag_invalid );\r
116         return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;\r
117     }\r
118     if ( roundBits ) sf->float_exception_flags |= float_flag_inexact;\r
119     return z;\r
120 \r
121 }\r
122 \r
123 /*----------------------------------------------------------------------------\r
124 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and\r
125 | `absZ1', with binary point between bits 63 and 64 (between the input words),\r
126 | and returns the properly rounded 64-bit integer corresponding to the input.\r
127 | If `zSign' is 1, the input is negated before being converted to an integer.\r
128 | Ordinarily, the fixed-point input is simply rounded to an integer, with\r
129 | the inexact exception raised if the input cannot be represented exactly as\r
130 | an integer.  However, if the fixed-point input is too large, the invalid\r
131 | exception is raised and the largest positive or negative integer is\r
132 | returned.\r
133 *----------------------------------------------------------------------------*/\r
134 \r
135 int64_t roundAndPackInt64( softfloat_t* sf, flag zSign, bits64 absZ0, bits64 absZ1 )\r
136 {\r
137     int8_t roundingMode;\r
138     flag roundNearestEven, increment;\r
139     int64_t z;\r
140 \r
141     roundingMode = sf->float_rounding_mode;\r
142     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
143     increment = ( (sbits64) absZ1 < 0 );\r
144     if ( ! roundNearestEven ) {\r
145         if ( roundingMode == float_round_to_zero ) {\r
146             increment = 0;\r
147         }\r
148         else {\r
149             if ( zSign ) {\r
150                 increment = ( roundingMode == float_round_down ) && absZ1;\r
151             }\r
152             else {\r
153                 increment = ( roundingMode == float_round_up ) && absZ1;\r
154             }\r
155         }\r
156     }\r
157     if ( increment ) {\r
158         ++absZ0;\r
159         if ( absZ0 == 0 ) goto overflow;\r
160         absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );\r
161     }\r
162     z = absZ0;\r
163     if ( zSign ) z = - z;\r
164     if ( z && ( ( z < 0 ) ^ zSign ) ) {\r
165  overflow:\r
166         float_raise( sf, float_flag_invalid );\r
167         return\r
168               zSign ? (sbits64) LIT64( 0x8000000000000000 )\r
169             : LIT64( 0x7FFFFFFFFFFFFFFF );\r
170     }\r
171     if ( absZ1 ) sf->float_exception_flags |= float_flag_inexact;\r
172     return z;\r
173 \r
174 }\r
175 \r
176 /*----------------------------------------------------------------------------\r
177 | Returns the fraction bits of the single-precision floating-point value `a'.\r
178 *----------------------------------------------------------------------------*/\r
179 \r
180 INLINE bits32 extractFloat32Frac( float32 a )\r
181 {\r
182 \r
183     return a & 0x007FFFFF;\r
184 \r
185 }\r
186 \r
187 /*----------------------------------------------------------------------------\r
188 | Returns the exponent bits of the single-precision floating-point value `a'.\r
189 *----------------------------------------------------------------------------*/\r
190 \r
191 INLINE int16_t extractFloat32Exp( float32 a )\r
192 {\r
193 \r
194     return ( a>>23 ) & 0xFF;\r
195 \r
196 }\r
197 \r
198 /*----------------------------------------------------------------------------\r
199 | Returns the sign bit of the single-precision floating-point value `a'.\r
200 *----------------------------------------------------------------------------*/\r
201 \r
202 INLINE flag extractFloat32Sign( float32 a )\r
203 {\r
204 \r
205     return a>>31;\r
206 \r
207 }\r
208 \r
209 /*----------------------------------------------------------------------------\r
210 | Normalizes the subnormal single-precision floating-point value represented\r
211 | by the denormalized significand `aSig'.  The normalized exponent and\r
212 | significand are stored at the locations pointed to by `zExpPtr' and\r
213 | `zSigPtr', respectively.\r
214 *----------------------------------------------------------------------------*/\r
215 \r
216 void normalizeFloat32Subnormal( bits32 aSig, int16_t *zExpPtr, bits32 *zSigPtr )\r
217 {\r
218     int8_t shiftCount;\r
219 \r
220     shiftCount = countLeadingZeros32( aSig ) - 8;\r
221     *zSigPtr = aSig<<shiftCount;\r
222     *zExpPtr = 1 - shiftCount;\r
223 \r
224 }\r
225 \r
226 /*----------------------------------------------------------------------------\r
227 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a\r
228 | single-precision floating-point value, returning the result.  After being\r
229 | shifted into the proper positions, the three fields are simply added\r
230 | together to form the result.  This means that any integer portion of `zSig'\r
231 | will be added into the exponent.  Since a properly normalized significand\r
232 | will have an integer portion equal to 1, the `zExp' input should be 1 less\r
233 | than the desired result exponent whenever `zSig' is a complete, normalized\r
234 | significand.\r
235 *----------------------------------------------------------------------------*/\r
236 \r
237 INLINE float32 packFloat32( flag zSign, int16_t zExp, bits32 zSig )\r
238 {\r
239 \r
240     return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;\r
241 \r
242 }\r
243 \r
244 /*----------------------------------------------------------------------------\r
245 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
246 | and significand `zSig', and returns the proper single-precision floating-\r
247 | point value corresponding to the abstract input.  Ordinarily, the abstract\r
248 | value is simply rounded and packed into the single-precision format, with\r
249 | the inexact exception raised if the abstract input cannot be represented\r
250 | exactly.  However, if the abstract value is too large, the overflow and\r
251 | inexact exceptions are raised and an infinity or maximal finite value is\r
252 | returned.  If the abstract value is too small, the input value is rounded to\r
253 | a subnormal number, and the underflow and inexact exceptions are raised if\r
254 | the abstract input cannot be represented exactly as a subnormal single-\r
255 | precision floating-point number.\r
256 |     The input significand `zSig' has its binary point between bits 30\r
257 | and 29, which is 7 bits to the left of the usual location.  This shifted\r
258 | significand must be normalized or smaller.  If `zSig' is not normalized,\r
259 | `zExp' must be 0; in that case, the result returned is a subnormal number,\r
260 | and it must not require rounding.  In the usual case that `zSig' is\r
261 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.\r
262 | The handling of underflow and overflow follows the IEC/IEEE Standard for\r
263 | Binary Floating-Point Arithmetic.\r
264 *----------------------------------------------------------------------------*/\r
265 \r
266 float32 roundAndPackFloat32( softfloat_t* sf, flag zSign, int16_t zExp, bits32 zSig )\r
267 {\r
268     int8_t roundingMode;\r
269     flag roundNearestEven;\r
270     int8_t roundIncrement, roundBits;\r
271     flag isTiny;\r
272 \r
273     roundingMode = sf->float_rounding_mode;\r
274     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
275     roundIncrement = 0x40;\r
276     if ( ! roundNearestEven ) {\r
277         if ( roundingMode == float_round_to_zero ) {\r
278             roundIncrement = 0;\r
279         }\r
280         else {\r
281             roundIncrement = 0x7F;\r
282             if ( zSign ) {\r
283                 if ( roundingMode == float_round_up ) roundIncrement = 0;\r
284             }\r
285             else {\r
286                 if ( roundingMode == float_round_down ) roundIncrement = 0;\r
287             }\r
288         }\r
289     }\r
290     roundBits = zSig & 0x7F;\r
291     if ( 0xFD <= (bits16) zExp ) {\r
292         if (    ( 0xFD < zExp )\r
293              || (    ( zExp == 0xFD )\r
294                   && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )\r
295            ) {\r
296             float_raise( sf, float_flag_overflow | float_flag_inexact );\r
297             return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );\r
298         }\r
299         if ( zExp < 0 ) {\r
300             isTiny =\r
301                    ( sf->float_detect_tininess == float_tininess_before_rounding )\r
302                 || ( zExp < -1 )\r
303                 || ( zSig + roundIncrement < 0x80000000 );\r
304             shift32RightJamming( zSig, - zExp, &zSig );\r
305             zExp = 0;\r
306             roundBits = zSig & 0x7F;\r
307             if ( isTiny && roundBits ) float_raise( sf, float_flag_underflow );\r
308         }\r
309     }\r
310     if ( roundBits ) sf->float_exception_flags |= float_flag_inexact;\r
311     zSig = ( zSig + roundIncrement )>>7;\r
312     zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );\r
313     if ( zSig == 0 ) zExp = 0;\r
314     return packFloat32( zSign, zExp, zSig );\r
315 \r
316 }\r
317 \r
318 /*----------------------------------------------------------------------------\r
319 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
320 | and significand `zSig', and returns the proper single-precision floating-\r
321 | point value corresponding to the abstract input.  This routine is just like\r
322 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.\r
323 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''\r
324 | floating-point exponent.\r
325 *----------------------------------------------------------------------------*/\r
326 \r
327 float32 normalizeRoundAndPackFloat32( softfloat_t* sf, flag zSign, int16_t zExp, bits32 zSig )\r
328 {\r
329     int8_t shiftCount;\r
330 \r
331     shiftCount = countLeadingZeros32( zSig ) - 1;\r
332     return roundAndPackFloat32( sf, zSign, zExp - shiftCount, zSig<<shiftCount );\r
333 \r
334 }\r
335 \r
336 /*----------------------------------------------------------------------------\r
337 | Returns the fraction bits of the double-precision floating-point value `a'.\r
338 *----------------------------------------------------------------------------*/\r
339 \r
340 INLINE bits64 extractFloat64Frac( float64 a )\r
341 {\r
342 \r
343     return a & LIT64( 0x000FFFFFFFFFFFFF );\r
344 \r
345 }\r
346 \r
347 /*----------------------------------------------------------------------------\r
348 | Returns the exponent bits of the double-precision floating-point value `a'.\r
349 *----------------------------------------------------------------------------*/\r
350 \r
351 INLINE int16_t extractFloat64Exp( float64 a )\r
352 {\r
353 \r
354     return ( a>>52 ) & 0x7FF;\r
355 \r
356 }\r
357 \r
358 /*----------------------------------------------------------------------------\r
359 | Returns the sign bit of the double-precision floating-point value `a'.\r
360 *----------------------------------------------------------------------------*/\r
361 \r
362 INLINE flag extractFloat64Sign( float64 a )\r
363 {\r
364 \r
365     return a>>63;\r
366 \r
367 }\r
368 \r
369 /*----------------------------------------------------------------------------\r
370 | Normalizes the subnormal double-precision floating-point value represented\r
371 | by the denormalized significand `aSig'.  The normalized exponent and\r
372 | significand are stored at the locations pointed to by `zExpPtr' and\r
373 | `zSigPtr', respectively.\r
374 *----------------------------------------------------------------------------*/\r
375 \r
376 void normalizeFloat64Subnormal( bits64 aSig, int16_t *zExpPtr, bits64 *zSigPtr )\r
377 {\r
378     int8_t shiftCount;\r
379 \r
380     shiftCount = countLeadingZeros64( aSig ) - 11;\r
381     *zSigPtr = aSig<<shiftCount;\r
382     *zExpPtr = 1 - shiftCount;\r
383 \r
384 }\r
385 \r
386 /*----------------------------------------------------------------------------\r
387 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a\r
388 | double-precision floating-point value, returning the result.  After being\r
389 | shifted into the proper positions, the three fields are simply added\r
390 | together to form the result.  This means that any integer portion of `zSig'\r
391 | will be added into the exponent.  Since a properly normalized significand\r
392 | will have an integer portion equal to 1, the `zExp' input should be 1 less\r
393 | than the desired result exponent whenever `zSig' is a complete, normalized\r
394 | significand.\r
395 *----------------------------------------------------------------------------*/\r
396 \r
397 INLINE float64 packFloat64( flag zSign, int16_t zExp, bits64 zSig )\r
398 {\r
399 \r
400     return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;\r
401 \r
402 }\r
403 \r
404 /*----------------------------------------------------------------------------\r
405 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
406 | and significand `zSig', and returns the proper double-precision floating-\r
407 | point value corresponding to the abstract input.  Ordinarily, the abstract\r
408 | value is simply rounded and packed into the double-precision format, with\r
409 | the inexact exception raised if the abstract input cannot be represented\r
410 | exactly.  However, if the abstract value is too large, the overflow and\r
411 | inexact exceptions are raised and an infinity or maximal finite value is\r
412 | returned.  If the abstract value is too small, the input value is rounded\r
413 | to a subnormal number, and the underflow and inexact exceptions are raised\r
414 | if the abstract input cannot be represented exactly as a subnormal double-\r
415 | precision floating-point number.\r
416 |     The input significand `zSig' has its binary point between bits 62\r
417 | and 61, which is 10 bits to the left of the usual location.  This shifted\r
418 | significand must be normalized or smaller.  If `zSig' is not normalized,\r
419 | `zExp' must be 0; in that case, the result returned is a subnormal number,\r
420 | and it must not require rounding.  In the usual case that `zSig' is\r
421 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.\r
422 | The handling of underflow and overflow follows the IEC/IEEE Standard for\r
423 | Binary Floating-Point Arithmetic.\r
424 *----------------------------------------------------------------------------*/\r
425 \r
426 float64 roundAndPackFloat64( softfloat_t* sf, flag zSign, int16_t zExp, bits64 zSig )\r
427 {\r
428     int8_t roundingMode;\r
429     flag roundNearestEven;\r
430     int16_t roundIncrement, roundBits;\r
431     flag isTiny;\r
432 \r
433     roundingMode = sf->float_rounding_mode;\r
434     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
435     roundIncrement = 0x200;\r
436     if ( ! roundNearestEven ) {\r
437         if ( roundingMode == float_round_to_zero ) {\r
438             roundIncrement = 0;\r
439         }\r
440         else {\r
441             roundIncrement = 0x3FF;\r
442             if ( zSign ) {\r
443                 if ( roundingMode == float_round_up ) roundIncrement = 0;\r
444             }\r
445             else {\r
446                 if ( roundingMode == float_round_down ) roundIncrement = 0;\r
447             }\r
448         }\r
449     }\r
450     roundBits = zSig & 0x3FF;\r
451     if ( 0x7FD <= (bits16) zExp ) {\r
452         if (    ( 0x7FD < zExp )\r
453              || (    ( zExp == 0x7FD )\r
454                   && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )\r
455            ) {\r
456             float_raise( sf, float_flag_overflow | float_flag_inexact );\r
457             return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );\r
458         }\r
459         if ( zExp < 0 ) {\r
460             isTiny =\r
461                    ( sf->float_detect_tininess == float_tininess_before_rounding )\r
462                 || ( zExp < -1 )\r
463                 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );\r
464             shift64RightJamming( zSig, - zExp, &zSig );\r
465             zExp = 0;\r
466             roundBits = zSig & 0x3FF;\r
467             if ( isTiny && roundBits ) float_raise( sf, float_flag_underflow );\r
468         }\r
469     }\r
470     if ( roundBits ) sf->float_exception_flags |= float_flag_inexact;\r
471     zSig = ( zSig + roundIncrement )>>10;\r
472     zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );\r
473     if ( zSig == 0 ) zExp = 0;\r
474     return packFloat64( zSign, zExp, zSig );\r
475 \r
476 }\r
477 \r
478 /*----------------------------------------------------------------------------\r
479 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
480 | and significand `zSig', and returns the proper double-precision floating-\r
481 | point value corresponding to the abstract input.  This routine is just like\r
482 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.\r
483 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''\r
484 | floating-point exponent.\r
485 *----------------------------------------------------------------------------*/\r
486 \r
487 float64 normalizeRoundAndPackFloat64( softfloat_t* sf, flag zSign, int16_t zExp, bits64 zSig )\r
488 {\r
489     int8_t shiftCount;\r
490 \r
491     shiftCount = countLeadingZeros64( zSig ) - 1;\r
492     return roundAndPackFloat64( sf, zSign, zExp - shiftCount, zSig<<shiftCount );\r
493 \r
494 }\r
495 \r
496 #ifdef FLOATX80\r
497 \r
498 /*----------------------------------------------------------------------------\r
499 | Returns the fraction bits of the extended double-precision floating-point\r
500 | value `a'.\r
501 *----------------------------------------------------------------------------*/\r
502 \r
503 INLINE bits64 extractFloatx80Frac( floatx80 a )\r
504 {\r
505 \r
506     return a.low;\r
507 \r
508 }\r
509 \r
510 /*----------------------------------------------------------------------------\r
511 | Returns the exponent bits of the extended double-precision floating-point\r
512 | value `a'.\r
513 *----------------------------------------------------------------------------*/\r
514 \r
515 INLINE int32_t extractFloatx80Exp( floatx80 a )\r
516 {\r
517 \r
518     return a.high & 0x7FFF;\r
519 \r
520 }\r
521 \r
522 /*----------------------------------------------------------------------------\r
523 | Returns the sign bit of the extended double-precision floating-point value\r
524 | `a'.\r
525 *----------------------------------------------------------------------------*/\r
526 \r
527 INLINE flag extractFloatx80Sign( floatx80 a )\r
528 {\r
529 \r
530     return a.high>>15;\r
531 \r
532 }\r
533 \r
534 /*----------------------------------------------------------------------------\r
535 | Normalizes the subnormal extended double-precision floating-point value\r
536 | represented by the denormalized significand `aSig'.  The normalized exponent\r
537 | and significand are stored at the locations pointed to by `zExpPtr' and\r
538 | `zSigPtr', respectively.\r
539 *----------------------------------------------------------------------------*/\r
540 \r
541 void normalizeFloatx80Subnormal( bits64 aSig, int32_t *zExpPtr, bits64 *zSigPtr )\r
542 {\r
543     int8_t shiftCount;\r
544 \r
545     shiftCount = countLeadingZeros64( aSig );\r
546     *zSigPtr = aSig<<shiftCount;\r
547     *zExpPtr = 1 - shiftCount;\r
548 \r
549 }\r
550 \r
551 /*----------------------------------------------------------------------------\r
552 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an\r
553 | extended double-precision floating-point value, returning the result.\r
554 *----------------------------------------------------------------------------*/\r
555 \r
556 INLINE floatx80 packFloatx80( flag zSign, int32_t zExp, bits64 zSig )\r
557 {\r
558     floatx80 z;\r
559 \r
560     z.low = zSig;\r
561     z.high = ( ( (bits16) zSign )<<15 ) + zExp;\r
562     return z;\r
563 \r
564 }\r
565 \r
566 /*----------------------------------------------------------------------------\r
567 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
568 | and extended significand formed by the concatenation of `zSig0' and `zSig1',\r
569 | and returns the proper extended double-precision floating-point value\r
570 | corresponding to the abstract input.  Ordinarily, the abstract value is\r
571 | rounded and packed into the extended double-precision format, with the\r
572 | inexact exception raised if the abstract input cannot be represented\r
573 | exactly.  However, if the abstract value is too large, the overflow and\r
574 | inexact exceptions are raised and an infinity or maximal finite value is\r
575 | returned.  If the abstract value is too small, the input value is rounded to\r
576 | a subnormal number, and the underflow and inexact exceptions are raised if\r
577 | the abstract input cannot be represented exactly as a subnormal extended\r
578 | double-precision floating-point number.\r
579 |     If `roundingPrecision' is 32 or 64, the result is rounded to the same\r
580 | number of bits as single or double precision, respectively.  Otherwise, the\r
581 | result is rounded to the full precision of the extended double-precision\r
582 | format.\r
583 |     The input significand must be normalized or smaller.  If the input\r
584 | significand is not normalized, `zExp' must be 0; in that case, the result\r
585 | returned is a subnormal number, and it must not require rounding.  The\r
586 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary\r
587 | Floating-Point Arithmetic.\r
588 *----------------------------------------------------------------------------*/\r
589 \r
590 floatx80 roundAndPackFloatx80( softfloat_t* sf,\r
591      int8_t roundingPrecision, flag zSign, int32_t zExp, bits64 zSig0, bits64 zSig1\r
592  )\r
593 {\r
594     int8_t roundingMode;\r
595     flag roundNearestEven, increment, isTiny;\r
596     int64_t roundIncrement, roundMask, roundBits;\r
597 \r
598     roundingMode = sf->float_rounding_mode;\r
599     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
600     if ( roundingPrecision == 80 ) goto precision80;\r
601     if ( roundingPrecision == 64 ) {\r
602         roundIncrement = LIT64( 0x0000000000000400 );\r
603         roundMask = LIT64( 0x00000000000007FF );\r
604     }\r
605     else if ( roundingPrecision == 32 ) {\r
606         roundIncrement = LIT64( 0x0000008000000000 );\r
607         roundMask = LIT64( 0x000000FFFFFFFFFF );\r
608     }\r
609     else {\r
610         goto precision80;\r
611     }\r
612     zSig0 |= ( zSig1 != 0 );\r
613     if ( ! roundNearestEven ) {\r
614         if ( roundingMode == float_round_to_zero ) {\r
615             roundIncrement = 0;\r
616         }\r
617         else {\r
618             roundIncrement = roundMask;\r
619             if ( zSign ) {\r
620                 if ( roundingMode == float_round_up ) roundIncrement = 0;\r
621             }\r
622             else {\r
623                 if ( roundingMode == float_round_down ) roundIncrement = 0;\r
624             }\r
625         }\r
626     }\r
627     roundBits = zSig0 & roundMask;\r
628     if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {\r
629         if (    ( 0x7FFE < zExp )\r
630              || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )\r
631            ) {\r
632             goto overflow;\r
633         }\r
634         if ( zExp <= 0 ) {\r
635             isTiny =\r
636                    ( sf->float_detect_tininess == float_tininess_before_rounding )\r
637                 || ( zExp < 0 )\r
638                 || ( zSig0 <= zSig0 + roundIncrement );\r
639             shift64RightJamming( zSig0, 1 - zExp, &zSig0 );\r
640             zExp = 0;\r
641             roundBits = zSig0 & roundMask;\r
642             if ( isTiny && roundBits ) float_raise( sf, float_flag_underflow );\r
643             if ( roundBits ) sf->float_exception_flags |= float_flag_inexact;\r
644             zSig0 += roundIncrement;\r
645             if ( (sbits64) zSig0 < 0 ) zExp = 1;\r
646             roundIncrement = roundMask + 1;\r
647             if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {\r
648                 roundMask |= roundIncrement;\r
649             }\r
650             zSig0 &= ~ roundMask;\r
651             return packFloatx80( zSign, zExp, zSig0 );\r
652         }\r
653     }\r
654     if ( roundBits ) sf->float_exception_flags |= float_flag_inexact;\r
655     zSig0 += roundIncrement;\r
656     if ( zSig0 < roundIncrement ) {\r
657         ++zExp;\r
658         zSig0 = LIT64( 0x8000000000000000 );\r
659     }\r
660     roundIncrement = roundMask + 1;\r
661     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {\r
662         roundMask |= roundIncrement;\r
663     }\r
664     zSig0 &= ~ roundMask;\r
665     if ( zSig0 == 0 ) zExp = 0;\r
666     return packFloatx80( zSign, zExp, zSig0 );\r
667  precision80:\r
668     increment = ( (sbits64) zSig1 < 0 );\r
669     if ( ! roundNearestEven ) {\r
670         if ( roundingMode == float_round_to_zero ) {\r
671             increment = 0;\r
672         }\r
673         else {\r
674             if ( zSign ) {\r
675                 increment = ( roundingMode == float_round_down ) && zSig1;\r
676             }\r
677             else {\r
678                 increment = ( roundingMode == float_round_up ) && zSig1;\r
679             }\r
680         }\r
681     }\r
682     if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {\r
683         if (    ( 0x7FFE < zExp )\r
684              || (    ( zExp == 0x7FFE )\r
685                   && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )\r
686                   && increment\r
687                 )\r
688            ) {\r
689             roundMask = 0;\r
690  overflow:\r
691             float_raise( sf, float_flag_overflow | float_flag_inexact );\r
692             if (    ( roundingMode == float_round_to_zero )\r
693                  || ( zSign && ( roundingMode == float_round_up ) )\r
694                  || ( ! zSign && ( roundingMode == float_round_down ) )\r
695                ) {\r
696                 return packFloatx80( zSign, 0x7FFE, ~ roundMask );\r
697             }\r
698             return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );\r
699         }\r
700         if ( zExp <= 0 ) {\r
701             isTiny =\r
702                    ( sf->float_detect_tininess == float_tininess_before_rounding )\r
703                 || ( zExp < 0 )\r
704                 || ! increment\r
705                 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );\r
706             shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );\r
707             zExp = 0;\r
708             if ( isTiny && zSig1 ) float_raise( sf, float_flag_underflow );\r
709             if ( zSig1 ) sf->float_exception_flags |= float_flag_inexact;\r
710             if ( roundNearestEven ) {\r
711                 increment = ( (sbits64) zSig1 < 0 );\r
712             }\r
713             else {\r
714                 if ( zSign ) {\r
715                     increment = ( roundingMode == float_round_down ) && zSig1;\r
716                 }\r
717                 else {\r
718                     increment = ( roundingMode == float_round_up ) && zSig1;\r
719                 }\r
720             }\r
721             if ( increment ) {\r
722                 ++zSig0;\r
723                 zSig0 &=\r
724                     ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );\r
725                 if ( (sbits64) zSig0 < 0 ) zExp = 1;\r
726             }\r
727             return packFloatx80( zSign, zExp, zSig0 );\r
728         }\r
729     }\r
730     if ( zSig1 ) sf->float_exception_flags |= float_flag_inexact;\r
731     if ( increment ) {\r
732         ++zSig0;\r
733         if ( zSig0 == 0 ) {\r
734             ++zExp;\r
735             zSig0 = LIT64( 0x8000000000000000 );\r
736         }\r
737         else {\r
738             zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );\r
739         }\r
740     }\r
741     else {\r
742         if ( zSig0 == 0 ) zExp = 0;\r
743     }\r
744     return packFloatx80( zSign, zExp, zSig0 );\r
745 \r
746 }\r
747 \r
748 /*----------------------------------------------------------------------------\r
749 | Takes an abstract floating-point value having sign `zSign', exponent\r
750 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',\r
751 | and returns the proper extended double-precision floating-point value\r
752 | corresponding to the abstract input.  This routine is just like\r
753 | `roundAndPackFloatx80' except that the input significand does not have to be\r
754 | normalized.\r
755 *----------------------------------------------------------------------------*/\r
756 \r
757 floatx80 normalizeRoundAndPackFloatx80( softfloat_t* sf, \r
758      int8_t roundingPrecision, flag zSign, int32_t zExp, bits64 zSig0, bits64 zSig1\r
759  )\r
760 {\r
761     int8_t shiftCount;\r
762 \r
763     if ( zSig0 == 0 ) {\r
764         zSig0 = zSig1;\r
765         zSig1 = 0;\r
766         zExp -= 64;\r
767     }\r
768     shiftCount = countLeadingZeros64( zSig0 );\r
769     shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );\r
770     zExp -= shiftCount;\r
771     return\r
772         roundAndPackFloatx80( sf,  roundingPrecision, zSign, zExp, zSig0, zSig1 );\r
773 \r
774 }\r
775 \r
776 #endif\r
777 \r
778 #ifdef FLOAT128\r
779 \r
780 /*----------------------------------------------------------------------------\r
781 | Returns the least-significant 64 fraction bits of the quadruple-precision\r
782 | floating-point value `a'.\r
783 *----------------------------------------------------------------------------*/\r
784 \r
785 INLINE bits64 extractFloat128Frac1( float128 a )\r
786 {\r
787 \r
788     return a.low;\r
789 \r
790 }\r
791 \r
792 /*----------------------------------------------------------------------------\r
793 | Returns the most-significant 48 fraction bits of the quadruple-precision\r
794 | floating-point value `a'.\r
795 *----------------------------------------------------------------------------*/\r
796 \r
797 INLINE bits64 extractFloat128Frac0( float128 a )\r
798 {\r
799 \r
800     return a.high & LIT64( 0x0000FFFFFFFFFFFF );\r
801 \r
802 }\r
803 \r
804 /*----------------------------------------------------------------------------\r
805 | Returns the exponent bits of the quadruple-precision floating-point value\r
806 | `a'.\r
807 *----------------------------------------------------------------------------*/\r
808 \r
809 INLINE int32_t extractFloat128Exp( float128 a )\r
810 {\r
811 \r
812     return ( a.high>>48 ) & 0x7FFF;\r
813 \r
814 }\r
815 \r
816 /*----------------------------------------------------------------------------\r
817 | Returns the sign bit of the quadruple-precision floating-point value `a'.\r
818 *----------------------------------------------------------------------------*/\r
819 \r
820 INLINE flag extractFloat128Sign( float128 a )\r
821 {\r
822 \r
823     return a.high>>63;\r
824 \r
825 }\r
826 \r
827 /*----------------------------------------------------------------------------\r
828 | Normalizes the subnormal quadruple-precision floating-point value\r
829 | represented by the denormalized significand formed by the concatenation of\r
830 | `aSig0' and `aSig1'.  The normalized exponent is stored at the location\r
831 | pointed to by `zExpPtr'.  The most significant 49 bits of the normalized\r
832 | significand are stored at the location pointed to by `zSig0Ptr', and the\r
833 | least significant 64 bits of the normalized significand are stored at the\r
834 | location pointed to by `zSig1Ptr'.\r
835 *----------------------------------------------------------------------------*/\r
836 \r
837 void normalizeFloat128Subnormal(\r
838      bits64 aSig0,\r
839      bits64 aSig1,\r
840      int32_t *zExpPtr,\r
841      bits64 *zSig0Ptr,\r
842      bits64 *zSig1Ptr\r
843  )\r
844 {\r
845     int8_t shiftCount;\r
846 \r
847     if ( aSig0 == 0 ) {\r
848         shiftCount = countLeadingZeros64( aSig1 ) - 15;\r
849         if ( shiftCount < 0 ) {\r
850             *zSig0Ptr = aSig1>>( - shiftCount );\r
851             *zSig1Ptr = aSig1<<( shiftCount & 63 );\r
852         }\r
853         else {\r
854             *zSig0Ptr = aSig1<<shiftCount;\r
855             *zSig1Ptr = 0;\r
856         }\r
857         *zExpPtr = - shiftCount - 63;\r
858     }\r
859     else {\r
860         shiftCount = countLeadingZeros64( aSig0 ) - 15;\r
861         shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );\r
862         *zExpPtr = 1 - shiftCount;\r
863     }\r
864 \r
865 }\r
866 \r
867 /*----------------------------------------------------------------------------\r
868 | Packs the sign `zSign', the exponent `zExp', and the significand formed\r
869 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision\r
870 | floating-point value, returning the result.  After being shifted into the\r
871 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply\r
872 | added together to form the most significant 32 bits of the result.  This\r
873 | means that any integer portion of `zSig0' will be added into the exponent.\r
874 | Since a properly normalized significand will have an integer portion equal\r
875 | to 1, the `zExp' input should be 1 less than the desired result exponent\r
876 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized\r
877 | significand.\r
878 *----------------------------------------------------------------------------*/\r
879 \r
880 INLINE float128 packFloat128( flag zSign, int32_t zExp, bits64 zSig0, bits64 zSig1 )\r
881 {\r
882     float128 z;\r
883 \r
884     z.low = zSig1;\r
885     z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;\r
886     return z;\r
887 \r
888 }\r
889 \r
890 /*----------------------------------------------------------------------------\r
891 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
892 | and extended significand formed by the concatenation of `zSig0', `zSig1',\r
893 | and `zSig2', and returns the proper quadruple-precision floating-point value\r
894 | corresponding to the abstract input.  Ordinarily, the abstract value is\r
895 | simply rounded and packed into the quadruple-precision format, with the\r
896 | inexact exception raised if the abstract input cannot be represented\r
897 | exactly.  However, if the abstract value is too large, the overflow and\r
898 | inexact exceptions are raised and an infinity or maximal finite value is\r
899 | returned.  If the abstract value is too small, the input value is rounded to\r
900 | a subnormal number, and the underflow and inexact exceptions are raised if\r
901 | the abstract input cannot be represented exactly as a subnormal quadruple-\r
902 | precision floating-point number.\r
903 |     The input significand must be normalized or smaller.  If the input\r
904 | significand is not normalized, `zExp' must be 0; in that case, the result\r
905 | returned is a subnormal number, and it must not require rounding.  In the\r
906 | usual case that the input significand is normalized, `zExp' must be 1 less\r
907 | than the ``true'' floating-point exponent.  The handling of underflow and\r
908 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
909 *----------------------------------------------------------------------------*/\r
910 \r
911 float128 roundAndPackFloat128( softfloat_t* sf,\r
912      flag zSign, int32_t zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 )\r
913 {\r
914     int8_t roundingMode;\r
915     flag roundNearestEven, increment, isTiny;\r
916 \r
917     roundingMode = sf->float_rounding_mode;\r
918     roundNearestEven = ( roundingMode == float_round_nearest_even );\r
919     increment = ( (sbits64) zSig2 < 0 );\r
920     if ( ! roundNearestEven ) {\r
921         if ( roundingMode == float_round_to_zero ) {\r
922             increment = 0;\r
923         }\r
924         else {\r
925             if ( zSign ) {\r
926                 increment = ( roundingMode == float_round_down ) && zSig2;\r
927             }\r
928             else {\r
929                 increment = ( roundingMode == float_round_up ) && zSig2;\r
930             }\r
931         }\r
932     }\r
933     if ( 0x7FFD <= (bits32) zExp ) {\r
934         if (    ( 0x7FFD < zExp )\r
935              || (    ( zExp == 0x7FFD )\r
936                   && eq128(\r
937                          LIT64( 0x0001FFFFFFFFFFFF ),\r
938                          LIT64( 0xFFFFFFFFFFFFFFFF ),\r
939                          zSig0,\r
940                          zSig1\r
941                      )\r
942                   && increment\r
943                 )\r
944            ) {\r
945             float_raise( sf, float_flag_overflow | float_flag_inexact );\r
946             if (    ( roundingMode == float_round_to_zero )\r
947                  || ( zSign && ( roundingMode == float_round_up ) )\r
948                  || ( ! zSign && ( roundingMode == float_round_down ) )\r
949                ) {\r
950                 return\r
951                     packFloat128(\r
952                         zSign,\r
953                         0x7FFE,\r
954                         LIT64( 0x0000FFFFFFFFFFFF ),\r
955                         LIT64( 0xFFFFFFFFFFFFFFFF )\r
956                     );\r
957             }\r
958             return packFloat128( zSign, 0x7FFF, 0, 0 );\r
959         }\r
960         if ( zExp < 0 ) {\r
961             isTiny =\r
962                    ( sf->float_detect_tininess == float_tininess_before_rounding )\r
963                 || ( zExp < -1 )\r
964                 || ! increment\r
965                 || lt128(\r
966                        zSig0,\r
967                        zSig1,\r
968                        LIT64( 0x0001FFFFFFFFFFFF ),\r
969                        LIT64( 0xFFFFFFFFFFFFFFFF )\r
970                    );\r
971             shift128ExtraRightJamming(\r
972                 zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );\r
973             zExp = 0;\r
974             if ( isTiny && zSig2 ) float_raise( sf, float_flag_underflow );\r
975             if ( roundNearestEven ) {\r
976                 increment = ( (sbits64) zSig2 < 0 );\r
977             }\r
978             else {\r
979                 if ( zSign ) {\r
980                     increment = ( roundingMode == float_round_down ) && zSig2;\r
981                 }\r
982                 else {\r
983                     increment = ( roundingMode == float_round_up ) && zSig2;\r
984                 }\r
985             }\r
986         }\r
987     }\r
988     if ( zSig2 ) sf->float_exception_flags |= float_flag_inexact;\r
989     if ( increment ) {\r
990         add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );\r
991         zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );\r
992     }\r
993     else {\r
994         if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;\r
995     }\r
996     return packFloat128( zSign, zExp, zSig0, zSig1 );\r
997 \r
998 }\r
999 \r
1000 /*----------------------------------------------------------------------------\r
1001 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',\r
1002 | and significand formed by the concatenation of `zSig0' and `zSig1', and\r
1003 | returns the proper quadruple-precision floating-point value corresponding\r
1004 | to the abstract input.  This routine is just like `roundAndPackFloat128'\r
1005 | except that the input significand has fewer bits and does not have to be\r
1006 | normalized.  In all cases, `zExp' must be 1 less than the ``true'' floating-\r
1007 | point exponent.\r
1008 *----------------------------------------------------------------------------*/\r
1009 \r
1010 float128 normalizeRoundAndPackFloat128( softfloat_t* sf, \r
1011      flag zSign, int32_t zExp, bits64 zSig0, bits64 zSig1 )\r
1012 {\r
1013     int8_t shiftCount;\r
1014     bits64 zSig2;\r
1015 \r
1016     if ( zSig0 == 0 ) {\r
1017         zSig0 = zSig1;\r
1018         zSig1 = 0;\r
1019         zExp -= 64;\r
1020     }\r
1021     shiftCount = countLeadingZeros64( zSig0 ) - 15;\r
1022     if ( 0 <= shiftCount ) {\r
1023         zSig2 = 0;\r
1024         shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );\r
1025     }\r
1026     else {\r
1027         shift128ExtraRightJamming(\r
1028             zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );\r
1029     }\r
1030     zExp -= shiftCount;\r
1031     return roundAndPackFloat128( sf,  zSign, zExp, zSig0, zSig1, zSig2 );\r
1032 \r
1033 }\r
1034 \r
1035 #endif\r
1036 \r
1037 /*----------------------------------------------------------------------------\r
1038 | Returns the result of converting the 32-bit two's complement integer `a'\r
1039 | to the single-precision floating-point format.  The conversion is performed\r
1040 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1041 *----------------------------------------------------------------------------*/\r
1042 \r
1043 float32 int32_to_float32( softfloat_t* sf, int32_t a )\r
1044 {\r
1045     flag zSign;\r
1046 \r
1047     if ( a == 0 ) return 0;\r
1048     if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );\r
1049     zSign = ( a < 0 );\r
1050     return normalizeRoundAndPackFloat32( sf, zSign, 0x9C, zSign ? - a : a );\r
1051 \r
1052 }\r
1053 \r
1054 /*----------------------------------------------------------------------------\r
1055 | Returns the result of converting the 32-bit two's complement integer `a'\r
1056 | to the double-precision floating-point format.  The conversion is performed\r
1057 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1058 *----------------------------------------------------------------------------*/\r
1059 \r
1060 float64 int32_to_float64( softfloat_t* sf, int32_t a )\r
1061 {\r
1062     flag zSign;\r
1063     uint32_t absA;\r
1064     int8_t shiftCount;\r
1065     bits64 zSig;\r
1066 \r
1067     if ( a == 0 ) return 0;\r
1068     zSign = ( a < 0 );\r
1069     absA = zSign ? - a : a;\r
1070     shiftCount = countLeadingZeros32( absA ) + 21;\r
1071     zSig = absA;\r
1072     return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );\r
1073 \r
1074 }\r
1075 \r
1076 #ifdef FLOATX80\r
1077 \r
1078 /*----------------------------------------------------------------------------\r
1079 | Returns the result of converting the 32-bit two's complement integer `a'\r
1080 | to the extended double-precision floating-point format.  The conversion\r
1081 | is performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1082 | Arithmetic.\r
1083 *----------------------------------------------------------------------------*/\r
1084 \r
1085 floatx80 int32_to_floatx80( softfloat_t* sf, int32_t a )\r
1086 {\r
1087     flag zSign;\r
1088     uint32_t absA;\r
1089     int8_t shiftCount;\r
1090     bits64 zSig;\r
1091 \r
1092     if ( a == 0 ) return packFloatx80( 0, 0, 0 );\r
1093     zSign = ( a < 0 );\r
1094     absA = zSign ? - a : a;\r
1095     shiftCount = countLeadingZeros32( absA ) + 32;\r
1096     zSig = absA;\r
1097     return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );\r
1098 \r
1099 }\r
1100 \r
1101 #endif\r
1102 \r
1103 #ifdef FLOAT128\r
1104 \r
1105 /*----------------------------------------------------------------------------\r
1106 | Returns the result of converting the 32-bit two's complement integer `a' to\r
1107 | the quadruple-precision floating-point format.  The conversion is performed\r
1108 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1109 *----------------------------------------------------------------------------*/\r
1110 \r
1111 float128 int32_to_float128( softfloat_t* sf, int32_t a )\r
1112 {\r
1113     flag zSign;\r
1114     uint32_t absA;\r
1115     int8_t shiftCount;\r
1116     bits64 zSig0;\r
1117 \r
1118     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );\r
1119     zSign = ( a < 0 );\r
1120     absA = zSign ? - a : a;\r
1121     shiftCount = countLeadingZeros32( absA ) + 17;\r
1122     zSig0 = absA;\r
1123     return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );\r
1124 \r
1125 }\r
1126 \r
1127 #endif\r
1128 \r
1129 /*----------------------------------------------------------------------------\r
1130 | Returns the result of converting the 64-bit two's complement integer `a'\r
1131 | to the single-precision floating-point format.  The conversion is performed\r
1132 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1133 *----------------------------------------------------------------------------*/\r
1134 \r
1135 float32 int64_to_float32( softfloat_t* sf, int64_t a )\r
1136 {\r
1137     flag zSign;\r
1138     uint64_t absA;\r
1139     int8_t shiftCount;\r
1140     bits32 zSig;\r
1141 \r
1142     if ( a == 0 ) return 0;\r
1143     zSign = ( a < 0 );\r
1144     absA = zSign ? - a : a;\r
1145     shiftCount = countLeadingZeros64( absA ) - 40;\r
1146     if ( 0 <= shiftCount ) {\r
1147         return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );\r
1148     }\r
1149     else {\r
1150         shiftCount += 7;\r
1151         if ( shiftCount < 0 ) {\r
1152             shift64RightJamming( absA, - shiftCount, &absA );\r
1153         }\r
1154         else {\r
1155             absA <<= shiftCount;\r
1156         }\r
1157         return roundAndPackFloat32( sf, zSign, 0x9C - shiftCount, absA );\r
1158     }\r
1159 \r
1160 }\r
1161 \r
1162 /*----------------------------------------------------------------------------\r
1163 | Returns the result of converting the 64-bit two's complement integer `a'\r
1164 | to the double-precision floating-point format.  The conversion is performed\r
1165 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1166 *----------------------------------------------------------------------------*/\r
1167 \r
1168 float64 int64_to_float64( softfloat_t* sf, int64_t a )\r
1169 {\r
1170     flag zSign;\r
1171 \r
1172     if ( a == 0 ) return 0;\r
1173     if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {\r
1174         return packFloat64( 1, 0x43E, 0 );\r
1175     }\r
1176     zSign = ( a < 0 );\r
1177     return normalizeRoundAndPackFloat64( sf, zSign, 0x43C, zSign ? - a : a );\r
1178 \r
1179 }\r
1180 \r
1181 #ifdef FLOATX80\r
1182 \r
1183 /*----------------------------------------------------------------------------\r
1184 | Returns the result of converting the 64-bit two's complement integer `a'\r
1185 | to the extended double-precision floating-point format.  The conversion\r
1186 | is performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1187 | Arithmetic.\r
1188 *----------------------------------------------------------------------------*/\r
1189 \r
1190 floatx80 int64_to_floatx80( softfloat_t* sf, int64_t a )\r
1191 {\r
1192     flag zSign;\r
1193     uint64_t absA;\r
1194     int8_t shiftCount;\r
1195 \r
1196     if ( a == 0 ) return packFloatx80( 0, 0, 0 );\r
1197     zSign = ( a < 0 );\r
1198     absA = zSign ? - a : a;\r
1199     shiftCount = countLeadingZeros64( absA );\r
1200     return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );\r
1201 \r
1202 }\r
1203 \r
1204 #endif\r
1205 \r
1206 #ifdef FLOAT128\r
1207 \r
1208 /*----------------------------------------------------------------------------\r
1209 | Returns the result of converting the 64-bit two's complement integer `a' to\r
1210 | the quadruple-precision floating-point format.  The conversion is performed\r
1211 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1212 *----------------------------------------------------------------------------*/\r
1213 \r
1214 float128 int64_to_float128( softfloat_t* sf, int64_t a )\r
1215 {\r
1216     flag zSign;\r
1217     uint64_t absA;\r
1218     int8_t shiftCount;\r
1219     int32_t zExp;\r
1220     bits64 zSig0, zSig1;\r
1221 \r
1222     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );\r
1223     zSign = ( a < 0 );\r
1224     absA = zSign ? - a : a;\r
1225     shiftCount = countLeadingZeros64( absA ) + 49;\r
1226     zExp = 0x406E - shiftCount;\r
1227     if ( 64 <= shiftCount ) {\r
1228         zSig1 = 0;\r
1229         zSig0 = absA;\r
1230         shiftCount -= 64;\r
1231     }\r
1232     else {\r
1233         zSig1 = absA;\r
1234         zSig0 = 0;\r
1235     }\r
1236     shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );\r
1237     return packFloat128( zSign, zExp, zSig0, zSig1 );\r
1238 \r
1239 }\r
1240 \r
1241 #endif\r
1242 \r
1243 /*----------------------------------------------------------------------------\r
1244 | Returns the result of converting the single-precision floating-point value\r
1245 | `a' to the 32-bit two's complement integer format.  The conversion is\r
1246 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1247 | Arithmetic---which means in particular that the conversion is rounded\r
1248 | according to the current rounding mode.  If `a' is a NaN, the largest\r
1249 | positive integer is returned.  Otherwise, if the conversion overflows, the\r
1250 | largest integer with the same sign as `a' is returned.\r
1251 *----------------------------------------------------------------------------*/\r
1252 \r
1253 int32_t float32_to_int32( softfloat_t* sf, float32 a )\r
1254 {\r
1255     flag aSign;\r
1256     int16_t aExp, shiftCount;\r
1257     bits32 aSig;\r
1258     bits64 aSig64;\r
1259 \r
1260     aSig = extractFloat32Frac( a );\r
1261     aExp = extractFloat32Exp( a );\r
1262     aSign = extractFloat32Sign( a );\r
1263     if ( ( aExp == 0xFF ) && aSig ) aSign = 0;\r
1264     if ( aExp ) aSig |= 0x00800000;\r
1265     shiftCount = 0xAF - aExp;\r
1266     aSig64 = aSig;\r
1267     aSig64 <<= 32;\r
1268     if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );\r
1269     return roundAndPackInt32( sf, aSign, aSig64 );\r
1270 \r
1271 }\r
1272 \r
1273 /*----------------------------------------------------------------------------\r
1274 | Returns the result of converting the single-precision floating-point value\r
1275 | `a' to the 32-bit two's complement integer format.  The conversion is\r
1276 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1277 | Arithmetic, except that the conversion is always rounded toward zero.\r
1278 | If `a' is a NaN, the largest positive integer is returned.  Otherwise, if\r
1279 | the conversion overflows, the largest integer with the same sign as `a' is\r
1280 | returned.\r
1281 *----------------------------------------------------------------------------*/\r
1282 \r
1283 int32_t float32_to_int32_round_to_zero( softfloat_t* sf, float32 a )\r
1284 {\r
1285     flag aSign;\r
1286     int16_t aExp, shiftCount;\r
1287     bits32 aSig;\r
1288     int32_t z;\r
1289 \r
1290     aSig = extractFloat32Frac( a );\r
1291     aExp = extractFloat32Exp( a );\r
1292     aSign = extractFloat32Sign( a );\r
1293     shiftCount = aExp - 0x9E;\r
1294     if ( 0 <= shiftCount ) {\r
1295         if ( a != 0xCF000000 ) {\r
1296             float_raise( sf, float_flag_invalid );\r
1297             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;\r
1298         }\r
1299         return (sbits32) 0x80000000;\r
1300     }\r
1301     else if ( aExp <= 0x7E ) {\r
1302         if ( aExp | aSig ) sf->float_exception_flags |= float_flag_inexact;\r
1303         return 0;\r
1304     }\r
1305     aSig = ( aSig | 0x00800000 )<<8;\r
1306     z = aSig>>( - shiftCount );\r
1307     if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {\r
1308         sf->float_exception_flags |= float_flag_inexact;\r
1309     }\r
1310     if ( aSign ) z = - z;\r
1311     return z;\r
1312 \r
1313 }\r
1314 \r
1315 /*----------------------------------------------------------------------------\r
1316 | Returns the result of converting the single-precision floating-point value\r
1317 | `a' to the 64-bit two's complement integer format.  The conversion is\r
1318 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1319 | Arithmetic---which means in particular that the conversion is rounded\r
1320 | according to the current rounding mode.  If `a' is a NaN, the largest\r
1321 | positive integer is returned.  Otherwise, if the conversion overflows, the\r
1322 | largest integer with the same sign as `a' is returned.\r
1323 *----------------------------------------------------------------------------*/\r
1324 \r
1325 int64_t float32_to_int64( softfloat_t* sf, float32 a )\r
1326 {\r
1327     flag aSign;\r
1328     int16_t aExp, shiftCount;\r
1329     bits32 aSig;\r
1330     bits64 aSig64, aSigExtra;\r
1331 \r
1332     aSig = extractFloat32Frac( a );\r
1333     aExp = extractFloat32Exp( a );\r
1334     aSign = extractFloat32Sign( a );\r
1335     shiftCount = 0xBE - aExp;\r
1336     if ( shiftCount < 0 ) {\r
1337         float_raise( sf, float_flag_invalid );\r
1338         if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {\r
1339             return LIT64( 0x7FFFFFFFFFFFFFFF );\r
1340         }\r
1341         return (sbits64) LIT64( 0x8000000000000000 );\r
1342     }\r
1343     if ( aExp ) aSig |= 0x00800000;\r
1344     aSig64 = aSig;\r
1345     aSig64 <<= 40;\r
1346     shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );\r
1347     return roundAndPackInt64( sf, aSign, aSig64, aSigExtra );\r
1348 \r
1349 }\r
1350 \r
1351 /*----------------------------------------------------------------------------\r
1352 | Returns the result of converting the single-precision floating-point value\r
1353 | `a' to the 64-bit two's complement integer format.  The conversion is\r
1354 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1355 | Arithmetic, except that the conversion is always rounded toward zero.  If\r
1356 | `a' is a NaN, the largest positive integer is returned.  Otherwise, if the\r
1357 | conversion overflows, the largest integer with the same sign as `a' is\r
1358 | returned.\r
1359 *----------------------------------------------------------------------------*/\r
1360 \r
1361 int64_t float32_to_int64_round_to_zero( softfloat_t* sf, float32 a )\r
1362 {\r
1363     flag aSign;\r
1364     int16_t aExp, shiftCount;\r
1365     bits32 aSig;\r
1366     bits64 aSig64;\r
1367     int64_t z;\r
1368 \r
1369     aSig = extractFloat32Frac( a );\r
1370     aExp = extractFloat32Exp( a );\r
1371     aSign = extractFloat32Sign( a );\r
1372     shiftCount = aExp - 0xBE;\r
1373     if ( 0 <= shiftCount ) {\r
1374         if ( a != 0xDF000000 ) {\r
1375             float_raise( sf, float_flag_invalid );\r
1376             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {\r
1377                 return LIT64( 0x7FFFFFFFFFFFFFFF );\r
1378             }\r
1379         }\r
1380         return (sbits64) LIT64( 0x8000000000000000 );\r
1381     }\r
1382     else if ( aExp <= 0x7E ) {\r
1383         if ( aExp | aSig ) sf->float_exception_flags |= float_flag_inexact;\r
1384         return 0;\r
1385     }\r
1386     aSig64 = aSig | 0x00800000;\r
1387     aSig64 <<= 40;\r
1388     z = aSig64>>( - shiftCount );\r
1389     if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {\r
1390         sf->float_exception_flags |= float_flag_inexact;\r
1391     }\r
1392     if ( aSign ) z = - z;\r
1393     return z;\r
1394 \r
1395 }\r
1396 \r
1397 /*----------------------------------------------------------------------------\r
1398 | Returns the result of converting the single-precision floating-point value\r
1399 | `a' to the double-precision floating-point format.  The conversion is\r
1400 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1401 | Arithmetic.\r
1402 *----------------------------------------------------------------------------*/\r
1403 \r
1404 float64 float32_to_float64( softfloat_t* sf, float32 a )\r
1405 {\r
1406     flag aSign;\r
1407     int16_t aExp;\r
1408     bits32 aSig;\r
1409 \r
1410     aSig = extractFloat32Frac( a );\r
1411     aExp = extractFloat32Exp( a );\r
1412     aSign = extractFloat32Sign( a );\r
1413     if ( aExp == 0xFF ) {\r
1414         if ( aSig ) return commonNaNToFloat64( sf, float32ToCommonNaN( sf, a ) );\r
1415         return packFloat64( aSign, 0x7FF, 0 );\r
1416     }\r
1417     if ( aExp == 0 ) {\r
1418         if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );\r
1419         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1420         --aExp;\r
1421     }\r
1422     return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );\r
1423 \r
1424 }\r
1425 \r
1426 #ifdef FLOATX80\r
1427 \r
1428 /*----------------------------------------------------------------------------\r
1429 | Returns the result of converting the single-precision floating-point value\r
1430 | `a' to the extended double-precision floating-point format.  The conversion\r
1431 | is performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1432 | Arithmetic.\r
1433 *----------------------------------------------------------------------------*/\r
1434 \r
1435 floatx80 float32_to_floatx80( softfloat_t* sf, float32 a )\r
1436 {\r
1437     flag aSign;\r
1438     int16_t aExp;\r
1439     bits32 aSig;\r
1440 \r
1441     aSig = extractFloat32Frac( a );\r
1442     aExp = extractFloat32Exp( a );\r
1443     aSign = extractFloat32Sign( a );\r
1444     if ( aExp == 0xFF ) {\r
1445         if ( aSig ) return commonNaNToFloatx80( sf, float32ToCommonNaN( sf, a ) );\r
1446         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );\r
1447     }\r
1448     if ( aExp == 0 ) {\r
1449         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );\r
1450         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1451     }\r
1452     aSig |= 0x00800000;\r
1453     return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );\r
1454 \r
1455 }\r
1456 \r
1457 #endif\r
1458 \r
1459 #ifdef FLOAT128\r
1460 \r
1461 /*----------------------------------------------------------------------------\r
1462 | Returns the result of converting the single-precision floating-point value\r
1463 | `a' to the double-precision floating-point format.  The conversion is\r
1464 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
1465 | Arithmetic.\r
1466 *----------------------------------------------------------------------------*/\r
1467 \r
1468 float128 float32_to_float128( softfloat_t* sf, float32 a )\r
1469 {\r
1470     flag aSign;\r
1471     int16_t aExp;\r
1472     bits32 aSig;\r
1473 \r
1474     aSig = extractFloat32Frac( a );\r
1475     aExp = extractFloat32Exp( a );\r
1476     aSign = extractFloat32Sign( a );\r
1477     if ( aExp == 0xFF ) {\r
1478         if ( aSig ) return commonNaNToFloat128( sf, float32ToCommonNaN( sf, a ) );\r
1479         return packFloat128( aSign, 0x7FFF, 0, 0 );\r
1480     }\r
1481     if ( aExp == 0 ) {\r
1482         if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );\r
1483         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1484         --aExp;\r
1485     }\r
1486     return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );\r
1487 \r
1488 }\r
1489 \r
1490 #endif\r
1491 \r
1492 /*----------------------------------------------------------------------------\r
1493 | Rounds the single-precision floating-point value `a' to an integer, and\r
1494 | returns the result as a single-precision floating-point value.  The\r
1495 | operation is performed according to the IEC/IEEE Standard for Binary\r
1496 | Floating-Point Arithmetic.\r
1497 *----------------------------------------------------------------------------*/\r
1498 \r
1499 float32 float32_round_to_int( softfloat_t* sf, float32 a )\r
1500 {\r
1501     flag aSign;\r
1502     int16_t aExp;\r
1503     bits32 lastBitMask, roundBitsMask;\r
1504     int8_t roundingMode;\r
1505     float32 z;\r
1506 \r
1507     aExp = extractFloat32Exp( a );\r
1508     if ( 0x96 <= aExp ) {\r
1509         if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {\r
1510             return propagateFloat32NaN( sf, a, a );\r
1511         }\r
1512         return a;\r
1513     }\r
1514     if ( aExp <= 0x7E ) {\r
1515         if ( (bits32) ( a<<1 ) == 0 ) return a;\r
1516         sf->float_exception_flags |= float_flag_inexact;\r
1517         aSign = extractFloat32Sign( a );\r
1518         switch ( sf->float_rounding_mode ) {\r
1519          case float_round_nearest_even:\r
1520             if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {\r
1521                 return packFloat32( aSign, 0x7F, 0 );\r
1522             }\r
1523             break;\r
1524          case float_round_down:\r
1525             return aSign ? 0xBF800000 : 0;\r
1526          case float_round_up:\r
1527             return aSign ? 0x80000000 : 0x3F800000;\r
1528         }\r
1529         return packFloat32( aSign, 0, 0 );\r
1530     }\r
1531     lastBitMask = 1;\r
1532     lastBitMask <<= 0x96 - aExp;\r
1533     roundBitsMask = lastBitMask - 1;\r
1534     z = a;\r
1535     roundingMode = sf->float_rounding_mode;\r
1536     if ( roundingMode == float_round_nearest_even ) {\r
1537         z += lastBitMask>>1;\r
1538         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;\r
1539     }\r
1540     else if ( roundingMode != float_round_to_zero ) {\r
1541         if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {\r
1542             z += roundBitsMask;\r
1543         }\r
1544     }\r
1545     z &= ~ roundBitsMask;\r
1546     if ( z != a ) sf->float_exception_flags |= float_flag_inexact;\r
1547     return z;\r
1548 \r
1549 }\r
1550 \r
1551 /*----------------------------------------------------------------------------\r
1552 | Returns the result of adding the absolute values of the single-precision\r
1553 | floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated\r
1554 | before being returned.  `zSign' is ignored if the result is a NaN.\r
1555 | The addition is performed according to the IEC/IEEE Standard for Binary\r
1556 | Floating-Point Arithmetic.\r
1557 *----------------------------------------------------------------------------*/\r
1558 \r
1559 float32 addFloat32Sigs( softfloat_t* sf, float32 a, float32 b, flag zSign )\r
1560 {\r
1561     int16_t aExp, bExp, zExp;\r
1562     bits32 aSig, bSig, zSig;\r
1563     int16_t expDiff;\r
1564 \r
1565     aSig = extractFloat32Frac( a );\r
1566     aExp = extractFloat32Exp( a );\r
1567     bSig = extractFloat32Frac( b );\r
1568     bExp = extractFloat32Exp( b );\r
1569     expDiff = aExp - bExp;\r
1570     aSig <<= 6;\r
1571     bSig <<= 6;\r
1572     if ( 0 < expDiff ) {\r
1573         if ( aExp == 0xFF ) {\r
1574             if ( aSig ) return propagateFloat32NaN( sf, a, b );\r
1575             return a;\r
1576         }\r
1577         if ( bExp == 0 ) {\r
1578             --expDiff;\r
1579         }\r
1580         else {\r
1581             bSig |= 0x20000000;\r
1582         }\r
1583         shift32RightJamming( bSig, expDiff, &bSig );\r
1584         zExp = aExp;\r
1585     }\r
1586     else if ( expDiff < 0 ) {\r
1587         if ( bExp == 0xFF ) {\r
1588             if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1589             return packFloat32( zSign, 0xFF, 0 );\r
1590         }\r
1591         if ( aExp == 0 ) {\r
1592             ++expDiff;\r
1593         }\r
1594         else {\r
1595             aSig |= 0x20000000;\r
1596         }\r
1597         shift32RightJamming( aSig, - expDiff, &aSig );\r
1598         zExp = bExp;\r
1599     }\r
1600     else {\r
1601         if ( aExp == 0xFF ) {\r
1602             if ( aSig | bSig ) return propagateFloat32NaN( sf, a, b );\r
1603             return a;\r
1604         }\r
1605         if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );\r
1606         zSig = 0x40000000 + aSig + bSig;\r
1607         zExp = aExp;\r
1608         goto roundAndPack;\r
1609     }\r
1610     aSig |= 0x20000000;\r
1611     zSig = ( aSig + bSig )<<1;\r
1612     --zExp;\r
1613     if ( (sbits32) zSig < 0 ) {\r
1614         zSig = aSig + bSig;\r
1615         ++zExp;\r
1616     }\r
1617  roundAndPack:\r
1618     return roundAndPackFloat32( sf, zSign, zExp, zSig );\r
1619 \r
1620 }\r
1621 \r
1622 /*----------------------------------------------------------------------------\r
1623 | Returns the result of subtracting the absolute values of the single-\r
1624 | precision floating-point values `a' and `b'.  If `zSign' is 1, the\r
1625 | difference is negated before being returned.  `zSign' is ignored if the\r
1626 | result is a NaN.  The subtraction is performed according to the IEC/IEEE\r
1627 | Standard for Binary Floating-Point Arithmetic.\r
1628 *----------------------------------------------------------------------------*/\r
1629 \r
1630 float32 subFloat32Sigs( softfloat_t* sf, float32 a, float32 b, flag zSign )\r
1631 {\r
1632     int16_t aExp, bExp, zExp;\r
1633     bits32 aSig, bSig, zSig;\r
1634     int16_t expDiff;\r
1635 \r
1636     aSig = extractFloat32Frac( a );\r
1637     aExp = extractFloat32Exp( a );\r
1638     bSig = extractFloat32Frac( b );\r
1639     bExp = extractFloat32Exp( b );\r
1640     expDiff = aExp - bExp;\r
1641     aSig <<= 7;\r
1642     bSig <<= 7;\r
1643     if ( 0 < expDiff ) goto aExpBigger;\r
1644     if ( expDiff < 0 ) goto bExpBigger;\r
1645     if ( aExp == 0xFF ) {\r
1646         if ( aSig | bSig ) return propagateFloat32NaN( sf, a, b );\r
1647         float_raise( sf, float_flag_invalid );\r
1648         return float32_default_nan;\r
1649     }\r
1650     if ( aExp == 0 ) {\r
1651         aExp = 1;\r
1652         bExp = 1;\r
1653     }\r
1654     if ( bSig < aSig ) goto aBigger;\r
1655     if ( aSig < bSig ) goto bBigger;\r
1656     return packFloat32( sf->float_rounding_mode == float_round_down, 0, 0 );\r
1657  bExpBigger:\r
1658     if ( bExp == 0xFF ) {\r
1659         if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1660         return packFloat32( zSign ^ 1, 0xFF, 0 );\r
1661     }\r
1662     if ( aExp == 0 ) {\r
1663         ++expDiff;\r
1664     }\r
1665     else {\r
1666         aSig |= 0x40000000;\r
1667     }\r
1668     shift32RightJamming( aSig, - expDiff, &aSig );\r
1669     bSig |= 0x40000000;\r
1670  bBigger:\r
1671     zSig = bSig - aSig;\r
1672     zExp = bExp;\r
1673     zSign ^= 1;\r
1674     goto normalizeRoundAndPack;\r
1675  aExpBigger:\r
1676     if ( aExp == 0xFF ) {\r
1677         if ( aSig ) return propagateFloat32NaN( sf, a, b );\r
1678         return a;\r
1679     }\r
1680     if ( bExp == 0 ) {\r
1681         --expDiff;\r
1682     }\r
1683     else {\r
1684         bSig |= 0x40000000;\r
1685     }\r
1686     shift32RightJamming( bSig, expDiff, &bSig );\r
1687     aSig |= 0x40000000;\r
1688  aBigger:\r
1689     zSig = aSig - bSig;\r
1690     zExp = aExp;\r
1691  normalizeRoundAndPack:\r
1692     --zExp;\r
1693     return normalizeRoundAndPackFloat32( sf, zSign, zExp, zSig );\r
1694 \r
1695 }\r
1696 \r
1697 /*----------------------------------------------------------------------------\r
1698 | Returns the result of adding the single-precision floating-point values `a'\r
1699 | and `b'.  The operation is performed according to the IEC/IEEE Standard for\r
1700 | Binary Floating-Point Arithmetic.\r
1701 *----------------------------------------------------------------------------*/\r
1702 \r
1703 float32 float32_add( softfloat_t* sf, float32 a, float32 b )\r
1704 {\r
1705     flag aSign, bSign;\r
1706 \r
1707     aSign = extractFloat32Sign( a );\r
1708     bSign = extractFloat32Sign( b );\r
1709     if ( aSign == bSign ) {\r
1710         return addFloat32Sigs( sf, a, b, aSign );\r
1711     }\r
1712     else {\r
1713         return subFloat32Sigs( sf, a, b, aSign );\r
1714     }\r
1715 \r
1716 }\r
1717 \r
1718 /*----------------------------------------------------------------------------\r
1719 | Returns the result of subtracting the single-precision floating-point values\r
1720 | `a' and `b'.  The operation is performed according to the IEC/IEEE Standard\r
1721 | for Binary Floating-Point Arithmetic.\r
1722 *----------------------------------------------------------------------------*/\r
1723 \r
1724 float32 float32_sub( softfloat_t* sf, float32 a, float32 b )\r
1725 {\r
1726     flag aSign, bSign;\r
1727 \r
1728     aSign = extractFloat32Sign( a );\r
1729     bSign = extractFloat32Sign( b );\r
1730     if ( aSign == bSign ) {\r
1731         return subFloat32Sigs( sf, a, b, aSign );\r
1732     }\r
1733     else {\r
1734         return addFloat32Sigs( sf, a, b, aSign );\r
1735     }\r
1736 \r
1737 }\r
1738 \r
1739 /*----------------------------------------------------------------------------\r
1740 | Returns the result of multiplying the single-precision floating-point values\r
1741 | `a' and `b'.  The operation is performed according to the IEC/IEEE Standard\r
1742 | for Binary Floating-Point Arithmetic.\r
1743 *----------------------------------------------------------------------------*/\r
1744 \r
1745 float32 float32_mul( softfloat_t* sf, float32 a, float32 b )\r
1746 {\r
1747     flag aSign, bSign, zSign;\r
1748     int16_t aExp, bExp, zExp;\r
1749     bits32 aSig, bSig;\r
1750     bits64 zSig64;\r
1751     bits32 zSig;\r
1752 \r
1753     aSig = extractFloat32Frac( a );\r
1754     aExp = extractFloat32Exp( a );\r
1755     aSign = extractFloat32Sign( a );\r
1756     bSig = extractFloat32Frac( b );\r
1757     bExp = extractFloat32Exp( b );\r
1758     bSign = extractFloat32Sign( b );\r
1759     zSign = aSign ^ bSign;\r
1760     if ( aExp == 0xFF ) {\r
1761         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {\r
1762             return propagateFloat32NaN( sf, a, b );\r
1763         }\r
1764         if ( ( bExp | bSig ) == 0 ) {\r
1765             float_raise( sf, float_flag_invalid );\r
1766             return float32_default_nan;\r
1767         }\r
1768         return packFloat32( zSign, 0xFF, 0 );\r
1769     }\r
1770     if ( bExp == 0xFF ) {\r
1771         if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1772         if ( ( aExp | aSig ) == 0 ) {\r
1773             float_raise( sf, float_flag_invalid );\r
1774             return float32_default_nan;\r
1775         }\r
1776         return packFloat32( zSign, 0xFF, 0 );\r
1777     }\r
1778     if ( aExp == 0 ) {\r
1779         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );\r
1780         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1781     }\r
1782     if ( bExp == 0 ) {\r
1783         if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );\r
1784         normalizeFloat32Subnormal( bSig, &bExp, &bSig );\r
1785     }\r
1786     zExp = aExp + bExp - 0x7F;\r
1787     aSig = ( aSig | 0x00800000 )<<7;\r
1788     bSig = ( bSig | 0x00800000 )<<8;\r
1789     shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );\r
1790     zSig = zSig64;\r
1791     if ( 0 <= (sbits32) ( zSig<<1 ) ) {\r
1792         zSig <<= 1;\r
1793         --zExp;\r
1794     }\r
1795     return roundAndPackFloat32( sf, zSign, zExp, zSig );\r
1796 \r
1797 }\r
1798 \r
1799 /*----------------------------------------------------------------------------\r
1800 | Returns the result of dividing the single-precision floating-point value `a'\r
1801 | by the corresponding value `b'.  The operation is performed according to the\r
1802 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1803 *----------------------------------------------------------------------------*/\r
1804 \r
1805 float32 float32_div( softfloat_t* sf, float32 a, float32 b )\r
1806 {\r
1807     flag aSign, bSign, zSign;\r
1808     int16_t aExp, bExp, zExp;\r
1809     bits32 aSig, bSig, zSig;\r
1810 \r
1811     aSig = extractFloat32Frac( a );\r
1812     aExp = extractFloat32Exp( a );\r
1813     aSign = extractFloat32Sign( a );\r
1814     bSig = extractFloat32Frac( b );\r
1815     bExp = extractFloat32Exp( b );\r
1816     bSign = extractFloat32Sign( b );\r
1817     zSign = aSign ^ bSign;\r
1818     if ( aExp == 0xFF ) {\r
1819         if ( aSig ) return propagateFloat32NaN( sf, a, b );\r
1820         if ( bExp == 0xFF ) {\r
1821             if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1822             float_raise( sf, float_flag_invalid );\r
1823             return float32_default_nan;\r
1824         }\r
1825         return packFloat32( zSign, 0xFF, 0 );\r
1826     }\r
1827     if ( bExp == 0xFF ) {\r
1828         if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1829         return packFloat32( zSign, 0, 0 );\r
1830     }\r
1831     if ( bExp == 0 ) {\r
1832         if ( bSig == 0 ) {\r
1833             if ( ( aExp | aSig ) == 0 ) {\r
1834                 float_raise( sf, float_flag_invalid );\r
1835                 return float32_default_nan;\r
1836             }\r
1837             float_raise( sf, float_flag_divbyzero );\r
1838             return packFloat32( zSign, 0xFF, 0 );\r
1839         }\r
1840         normalizeFloat32Subnormal( bSig, &bExp, &bSig );\r
1841     }\r
1842     if ( aExp == 0 ) {\r
1843         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );\r
1844         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1845     }\r
1846     zExp = aExp - bExp + 0x7D;\r
1847     aSig = ( aSig | 0x00800000 )<<7;\r
1848     bSig = ( bSig | 0x00800000 )<<8;\r
1849     if ( bSig <= ( aSig + aSig ) ) {\r
1850         aSig >>= 1;\r
1851         ++zExp;\r
1852     }\r
1853     zSig = ( ( (bits64) aSig )<<32 ) / bSig;\r
1854     if ( ( zSig & 0x3F ) == 0 ) {\r
1855         zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );\r
1856     }\r
1857     return roundAndPackFloat32( sf, zSign, zExp, zSig );\r
1858 \r
1859 }\r
1860 \r
1861 /*----------------------------------------------------------------------------\r
1862 | Returns the remainder of the single-precision floating-point value `a'\r
1863 | with respect to the corresponding value `b'.  The operation is performed\r
1864 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
1865 *----------------------------------------------------------------------------*/\r
1866 \r
1867 float32 float32_rem( softfloat_t* sf, float32 a, float32 b )\r
1868 {\r
1869     flag aSign, bSign, zSign;\r
1870     int16_t aExp, bExp, expDiff;\r
1871     bits32 aSig, bSig;\r
1872     bits32 q;\r
1873     bits64 aSig64, bSig64, q64;\r
1874     bits32 alternateASig;\r
1875     sbits32 sigMean;\r
1876 \r
1877     aSig = extractFloat32Frac( a );\r
1878     aExp = extractFloat32Exp( a );\r
1879     aSign = extractFloat32Sign( a );\r
1880     bSig = extractFloat32Frac( b );\r
1881     bExp = extractFloat32Exp( b );\r
1882     bSign = extractFloat32Sign( b );\r
1883     if ( aExp == 0xFF ) {\r
1884         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {\r
1885             return propagateFloat32NaN( sf, a, b );\r
1886         }\r
1887         float_raise( sf, float_flag_invalid );\r
1888         return float32_default_nan;\r
1889     }\r
1890     if ( bExp == 0xFF ) {\r
1891         if ( bSig ) return propagateFloat32NaN( sf, a, b );\r
1892         return a;\r
1893     }\r
1894     if ( bExp == 0 ) {\r
1895         if ( bSig == 0 ) {\r
1896             float_raise( sf, float_flag_invalid );\r
1897             return float32_default_nan;\r
1898         }\r
1899         normalizeFloat32Subnormal( bSig, &bExp, &bSig );\r
1900     }\r
1901     if ( aExp == 0 ) {\r
1902         if ( aSig == 0 ) return a;\r
1903         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1904     }\r
1905     expDiff = aExp - bExp;\r
1906     aSig |= 0x00800000;\r
1907     bSig |= 0x00800000;\r
1908     if ( expDiff < 32 ) {\r
1909         aSig <<= 8;\r
1910         bSig <<= 8;\r
1911         if ( expDiff < 0 ) {\r
1912             if ( expDiff < -1 ) return a;\r
1913             aSig >>= 1;\r
1914         }\r
1915         q = ( bSig <= aSig );\r
1916         if ( q ) aSig -= bSig;\r
1917         if ( 0 < expDiff ) {\r
1918             q = ( ( (bits64) aSig )<<32 ) / bSig;\r
1919             q >>= 32 - expDiff;\r
1920             bSig >>= 2;\r
1921             aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;\r
1922         }\r
1923         else {\r
1924             aSig >>= 2;\r
1925             bSig >>= 2;\r
1926         }\r
1927     }\r
1928     else {\r
1929         if ( bSig <= aSig ) aSig -= bSig;\r
1930         aSig64 = ( (bits64) aSig )<<40;\r
1931         bSig64 = ( (bits64) bSig )<<40;\r
1932         expDiff -= 64;\r
1933         while ( 0 < expDiff ) {\r
1934             q64 = estimateDiv128To64( aSig64, 0, bSig64 );\r
1935             q64 = ( 2 < q64 ) ? q64 - 2 : 0;\r
1936             aSig64 = - ( ( bSig * q64 )<<38 );\r
1937             expDiff -= 62;\r
1938         }\r
1939         expDiff += 64;\r
1940         q64 = estimateDiv128To64( aSig64, 0, bSig64 );\r
1941         q64 = ( 2 < q64 ) ? q64 - 2 : 0;\r
1942         q = q64>>( 64 - expDiff );\r
1943         bSig <<= 6;\r
1944         aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;\r
1945     }\r
1946     do {\r
1947         alternateASig = aSig;\r
1948         ++q;\r
1949         aSig -= bSig;\r
1950     } while ( 0 <= (sbits32) aSig );\r
1951     sigMean = aSig + alternateASig;\r
1952     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {\r
1953         aSig = alternateASig;\r
1954     }\r
1955     zSign = ( (sbits32) aSig < 0 );\r
1956     if ( zSign ) aSig = - aSig;\r
1957     return normalizeRoundAndPackFloat32( sf, aSign ^ zSign, bExp, aSig );\r
1958 \r
1959 }\r
1960 \r
1961 /*----------------------------------------------------------------------------\r
1962 | Returns the square root of the single-precision floating-point value `a'.\r
1963 | The operation is performed according to the IEC/IEEE Standard for Binary\r
1964 | Floating-Point Arithmetic.\r
1965 *----------------------------------------------------------------------------*/\r
1966 \r
1967 float32 float32_sqrt( softfloat_t* sf, float32 a )\r
1968 {\r
1969     flag aSign;\r
1970     int16_t aExp, zExp;\r
1971     bits32 aSig, zSig;\r
1972     bits64 rem, term;\r
1973 \r
1974     aSig = extractFloat32Frac( a );\r
1975     aExp = extractFloat32Exp( a );\r
1976     aSign = extractFloat32Sign( a );\r
1977     if ( aExp == 0xFF ) {\r
1978         if ( aSig ) return propagateFloat32NaN( sf, a, 0 );\r
1979         if ( ! aSign ) return a;\r
1980         float_raise( sf, float_flag_invalid );\r
1981         return float32_default_nan;\r
1982     }\r
1983     if ( aSign ) {\r
1984         if ( ( aExp | aSig ) == 0 ) return a;\r
1985         float_raise( sf, float_flag_invalid );\r
1986         return float32_default_nan;\r
1987     }\r
1988     if ( aExp == 0 ) {\r
1989         if ( aSig == 0 ) return 0;\r
1990         normalizeFloat32Subnormal( aSig, &aExp, &aSig );\r
1991     }\r
1992     zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;\r
1993     aSig = ( aSig | 0x00800000 )<<8;\r
1994     zSig = estimateSqrt32( aExp, aSig ) + 2;\r
1995     if ( ( zSig & 0x7F ) <= 5 ) {\r
1996         if ( zSig < 2 ) {\r
1997             zSig = 0x7FFFFFFF;\r
1998             goto roundAndPack;\r
1999         }\r
2000         aSig >>= aExp & 1;\r
2001         term = ( (bits64) zSig ) * zSig;\r
2002         rem = ( ( (bits64) aSig )<<32 ) - term;\r
2003         while ( (sbits64) rem < 0 ) {\r
2004             --zSig;\r
2005             rem += ( ( (bits64) zSig )<<1 ) | 1;\r
2006         }\r
2007         zSig |= ( rem != 0 );\r
2008     }\r
2009     shift32RightJamming( zSig, 1, &zSig );\r
2010  roundAndPack:\r
2011     return roundAndPackFloat32( sf, 0, zExp, zSig );\r
2012 \r
2013 }\r
2014 \r
2015 /*----------------------------------------------------------------------------\r
2016 | Returns 1 if the single-precision floating-point value `a' is equal to\r
2017 | the corresponding value `b', and 0 otherwise.  The comparison is performed\r
2018 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2019 *----------------------------------------------------------------------------*/\r
2020 \r
2021 flag float32_eq( softfloat_t* sf, float32 a, float32 b )\r
2022 {\r
2023 \r
2024     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2025          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2026        ) {\r
2027         if ( float32_is_signaling_nan( sf, a ) || float32_is_signaling_nan( sf, b ) ) {\r
2028             float_raise( sf, float_flag_invalid );\r
2029         }\r
2030         return 0;\r
2031     }\r
2032     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );\r
2033 \r
2034 }\r
2035 \r
2036 /*----------------------------------------------------------------------------\r
2037 | Returns 1 if the single-precision floating-point value `a' is less than\r
2038 | or equal to the corresponding value `b', and 0 otherwise.  The comparison\r
2039 | is performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2040 | Arithmetic.\r
2041 *----------------------------------------------------------------------------*/\r
2042 \r
2043 flag float32_le( softfloat_t* sf, float32 a, float32 b )\r
2044 {\r
2045     flag aSign, bSign;\r
2046 \r
2047     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2048          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2049        ) {\r
2050         float_raise( sf, float_flag_invalid );\r
2051         return 0;\r
2052     }\r
2053     aSign = extractFloat32Sign( a );\r
2054     bSign = extractFloat32Sign( b );\r
2055     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );\r
2056     return ( a == b ) || ( aSign ^ ( a < b ) );\r
2057 \r
2058 }\r
2059 \r
2060 /*----------------------------------------------------------------------------\r
2061 | Returns 1 if the single-precision floating-point value `a' is less than\r
2062 | the corresponding value `b', and 0 otherwise.  The comparison is performed\r
2063 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2064 *----------------------------------------------------------------------------*/\r
2065 \r
2066 flag float32_lt( softfloat_t* sf, float32 a, float32 b )\r
2067 {\r
2068     flag aSign, bSign;\r
2069 \r
2070     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2071          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2072        ) {\r
2073         float_raise( sf, float_flag_invalid );\r
2074         return 0;\r
2075     }\r
2076     aSign = extractFloat32Sign( a );\r
2077     bSign = extractFloat32Sign( b );\r
2078     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );\r
2079     return ( a != b ) && ( aSign ^ ( a < b ) );\r
2080 \r
2081 }\r
2082 \r
2083 /*----------------------------------------------------------------------------\r
2084 | Returns 1 if the single-precision floating-point value `a' is equal to\r
2085 | the corresponding value `b', and 0 otherwise.  The invalid exception is\r
2086 | raised if either operand is a NaN.  Otherwise, the comparison is performed\r
2087 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2088 *----------------------------------------------------------------------------*/\r
2089 \r
2090 flag float32_eq_signaling( softfloat_t* sf, float32 a, float32 b )\r
2091 {\r
2092 \r
2093     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2094          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2095        ) {\r
2096         float_raise( sf, float_flag_invalid );\r
2097         return 0;\r
2098     }\r
2099     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );\r
2100 \r
2101 }\r
2102 \r
2103 /*----------------------------------------------------------------------------\r
2104 | Returns 1 if the single-precision floating-point value `a' is less than or\r
2105 | equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not\r
2106 | cause an exception.  Otherwise, the comparison is performed according to the\r
2107 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2108 *----------------------------------------------------------------------------*/\r
2109 \r
2110 flag float32_le_quiet( softfloat_t* sf, float32 a, float32 b )\r
2111 {\r
2112     flag aSign, bSign;\r
2113     int16_t aExp, bExp;\r
2114 \r
2115     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2116          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2117        ) {\r
2118         if ( float32_is_signaling_nan( sf, a ) || float32_is_signaling_nan( sf, b ) ) {\r
2119             float_raise( sf, float_flag_invalid );\r
2120         }\r
2121         return 0;\r
2122     }\r
2123     aSign = extractFloat32Sign( a );\r
2124     bSign = extractFloat32Sign( b );\r
2125     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );\r
2126     return ( a == b ) || ( aSign ^ ( a < b ) );\r
2127 \r
2128 }\r
2129 \r
2130 /*----------------------------------------------------------------------------\r
2131 | Returns 1 if the single-precision floating-point value `a' is less than\r
2132 | the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an\r
2133 | exception.  Otherwise, the comparison is performed according to the IEC/IEEE\r
2134 | Standard for Binary Floating-Point Arithmetic.\r
2135 *----------------------------------------------------------------------------*/\r
2136 \r
2137 flag float32_lt_quiet( softfloat_t* sf, float32 a, float32 b )\r
2138 {\r
2139     flag aSign, bSign;\r
2140 \r
2141     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )\r
2142          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )\r
2143        ) {\r
2144         if ( float32_is_signaling_nan( sf, a ) || float32_is_signaling_nan( sf, b ) ) {\r
2145             float_raise( sf, float_flag_invalid );\r
2146         }\r
2147         return 0;\r
2148     }\r
2149     aSign = extractFloat32Sign( a );\r
2150     bSign = extractFloat32Sign( b );\r
2151     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );\r
2152     return ( a != b ) && ( aSign ^ ( a < b ) );\r
2153 \r
2154 }\r
2155 \r
2156 /*----------------------------------------------------------------------------\r
2157 | Returns the result of converting the double-precision floating-point value\r
2158 | `a' to the 32-bit two's complement integer format.  The conversion is\r
2159 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2160 | Arithmetic---which means in particular that the conversion is rounded\r
2161 | according to the current rounding mode.  If `a' is a NaN, the largest\r
2162 | positive integer is returned.  Otherwise, if the conversion overflows, the\r
2163 | largest integer with the same sign as `a' is returned.\r
2164 *----------------------------------------------------------------------------*/\r
2165 \r
2166 int32_t float64_to_int32( softfloat_t* sf, float64 a )\r
2167 {\r
2168     flag aSign;\r
2169     int16_t aExp, shiftCount;\r
2170     bits64 aSig;\r
2171 \r
2172     aSig = extractFloat64Frac( a );\r
2173     aExp = extractFloat64Exp( a );\r
2174     aSign = extractFloat64Sign( a );\r
2175     if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;\r
2176     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );\r
2177     shiftCount = 0x42C - aExp;\r
2178     if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );\r
2179     return roundAndPackInt32( sf, aSign, aSig );\r
2180 \r
2181 }\r
2182 \r
2183 /*----------------------------------------------------------------------------\r
2184 | Returns the result of converting the double-precision floating-point value\r
2185 | `a' to the 32-bit two's complement integer format.  The conversion is\r
2186 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2187 | Arithmetic, except that the conversion is always rounded toward zero.\r
2188 | If `a' is a NaN, the largest positive integer is returned.  Otherwise, if\r
2189 | the conversion overflows, the largest integer with the same sign as `a' is\r
2190 | returned.\r
2191 *----------------------------------------------------------------------------*/\r
2192 \r
2193 int32_t float64_to_int32_round_to_zero( softfloat_t* sf, float64 a )\r
2194 {\r
2195     flag aSign;\r
2196     int16_t aExp, shiftCount;\r
2197     bits64 aSig, savedASig;\r
2198     int32_t z;\r
2199 \r
2200     aSig = extractFloat64Frac( a );\r
2201     aExp = extractFloat64Exp( a );\r
2202     aSign = extractFloat64Sign( a );\r
2203     if ( 0x41E < aExp ) {\r
2204         if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;\r
2205         goto invalid;\r
2206     }\r
2207     else if ( aExp < 0x3FF ) {\r
2208         if ( aExp || aSig ) sf->float_exception_flags |= float_flag_inexact;\r
2209         return 0;\r
2210     }\r
2211     aSig |= LIT64( 0x0010000000000000 );\r
2212     shiftCount = 0x433 - aExp;\r
2213     savedASig = aSig;\r
2214     aSig >>= shiftCount;\r
2215     z = aSig;\r
2216     if ( aSign ) z = - z;\r
2217     if ( ( z < 0 ) ^ aSign ) {\r
2218  invalid:\r
2219         float_raise( sf, float_flag_invalid );\r
2220         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;\r
2221     }\r
2222     if ( ( aSig<<shiftCount ) != savedASig ) {\r
2223         sf->float_exception_flags |= float_flag_inexact;\r
2224     }\r
2225     return z;\r
2226 \r
2227 }\r
2228 \r
2229 /*----------------------------------------------------------------------------\r
2230 | Returns the result of converting the double-precision floating-point value\r
2231 | `a' to the 64-bit two's complement integer format.  The conversion is\r
2232 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2233 | Arithmetic---which means in particular that the conversion is rounded\r
2234 | according to the current rounding mode.  If `a' is a NaN, the largest\r
2235 | positive integer is returned.  Otherwise, if the conversion overflows, the\r
2236 | largest integer with the same sign as `a' is returned.\r
2237 *----------------------------------------------------------------------------*/\r
2238 \r
2239 int64_t float64_to_int64( softfloat_t* sf, float64 a )\r
2240 {\r
2241     flag aSign;\r
2242     int16_t aExp, shiftCount;\r
2243     bits64 aSig, aSigExtra;\r
2244 \r
2245     aSig = extractFloat64Frac( a );\r
2246     aExp = extractFloat64Exp( a );\r
2247     aSign = extractFloat64Sign( a );\r
2248     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );\r
2249     shiftCount = 0x433 - aExp;\r
2250     if ( shiftCount <= 0 ) {\r
2251         if ( 0x43E < aExp ) {\r
2252             float_raise( sf, float_flag_invalid );\r
2253             if (    ! aSign\r
2254                  || (    ( aExp == 0x7FF )\r
2255                       && ( aSig != LIT64( 0x0010000000000000 ) ) )\r
2256                ) {\r
2257                 return LIT64( 0x7FFFFFFFFFFFFFFF );\r
2258             }\r
2259             return (sbits64) LIT64( 0x8000000000000000 );\r
2260         }\r
2261         aSigExtra = 0;\r
2262         aSig <<= - shiftCount;\r
2263     }\r
2264     else {\r
2265         shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );\r
2266     }\r
2267     return roundAndPackInt64( sf, aSign, aSig, aSigExtra );\r
2268 \r
2269 }\r
2270 \r
2271 /*----------------------------------------------------------------------------\r
2272 | Returns the result of converting the double-precision floating-point value\r
2273 | `a' to the 64-bit two's complement integer format.  The conversion is\r
2274 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2275 | Arithmetic, except that the conversion is always rounded toward zero.\r
2276 | If `a' is a NaN, the largest positive integer is returned.  Otherwise, if\r
2277 | the conversion overflows, the largest integer with the same sign as `a' is\r
2278 | returned.\r
2279 *----------------------------------------------------------------------------*/\r
2280 \r
2281 int64_t float64_to_int64_round_to_zero( softfloat_t* sf, float64 a )\r
2282 {\r
2283     flag aSign;\r
2284     int16_t aExp, shiftCount;\r
2285     bits64 aSig;\r
2286     int64_t z;\r
2287 \r
2288     aSig = extractFloat64Frac( a );\r
2289     aExp = extractFloat64Exp( a );\r
2290     aSign = extractFloat64Sign( a );\r
2291     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );\r
2292     shiftCount = aExp - 0x433;\r
2293     if ( 0 <= shiftCount ) {\r
2294         if ( 0x43E <= aExp ) {\r
2295             if ( a != LIT64( 0xC3E0000000000000 ) ) {\r
2296                 float_raise( sf, float_flag_invalid );\r
2297                 if (    ! aSign\r
2298                      || (    ( aExp == 0x7FF )\r
2299                           && ( aSig != LIT64( 0x0010000000000000 ) ) )\r
2300                    ) {\r
2301                     return LIT64( 0x7FFFFFFFFFFFFFFF );\r
2302                 }\r
2303             }\r
2304             return (sbits64) LIT64( 0x8000000000000000 );\r
2305         }\r
2306         z = aSig<<shiftCount;\r
2307     }\r
2308     else {\r
2309         if ( aExp < 0x3FE ) {\r
2310             if ( aExp | aSig ) sf->float_exception_flags |= float_flag_inexact;\r
2311             return 0;\r
2312         }\r
2313         z = aSig>>( - shiftCount );\r
2314         if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {\r
2315             sf->float_exception_flags |= float_flag_inexact;\r
2316         }\r
2317     }\r
2318     if ( aSign ) z = - z;\r
2319     return z;\r
2320 \r
2321 }\r
2322 \r
2323 /*----------------------------------------------------------------------------\r
2324 | Returns the result of converting the double-precision floating-point value\r
2325 | `a' to the single-precision floating-point format.  The conversion is\r
2326 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2327 | Arithmetic.\r
2328 *----------------------------------------------------------------------------*/\r
2329 \r
2330 float32 float64_to_float32( softfloat_t* sf, float64 a )\r
2331 {\r
2332     flag aSign;\r
2333     int16_t aExp;\r
2334     bits64 aSig;\r
2335     bits32 zSig;\r
2336 \r
2337     aSig = extractFloat64Frac( a );\r
2338     aExp = extractFloat64Exp( a );\r
2339     aSign = extractFloat64Sign( a );\r
2340     if ( aExp == 0x7FF ) {\r
2341         if ( aSig ) return commonNaNToFloat32( sf, float64ToCommonNaN( sf, a ) );\r
2342         return packFloat32( aSign, 0xFF, 0 );\r
2343     }\r
2344     shift64RightJamming( aSig, 22, &aSig );\r
2345     zSig = aSig;\r
2346     if ( aExp || zSig ) {\r
2347         zSig |= 0x40000000;\r
2348         aExp -= 0x381;\r
2349     }\r
2350     return roundAndPackFloat32( sf, aSign, aExp, zSig );\r
2351 \r
2352 }\r
2353 \r
2354 #ifdef FLOATX80\r
2355 \r
2356 /*----------------------------------------------------------------------------\r
2357 | Returns the result of converting the double-precision floating-point value\r
2358 | `a' to the extended double-precision floating-point format.  The conversion\r
2359 | is performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2360 | Arithmetic.\r
2361 *----------------------------------------------------------------------------*/\r
2362 \r
2363 floatx80 float64_to_floatx80( softfloat_t* sf, float64 a )\r
2364 {\r
2365     flag aSign;\r
2366     int16_t aExp;\r
2367     bits64 aSig;\r
2368 \r
2369     aSig = extractFloat64Frac( a );\r
2370     aExp = extractFloat64Exp( a );\r
2371     aSign = extractFloat64Sign( a );\r
2372     if ( aExp == 0x7FF ) {\r
2373         if ( aSig ) return commonNaNToFloatx80( sf, float64ToCommonNaN( sf, a ) );\r
2374         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );\r
2375     }\r
2376     if ( aExp == 0 ) {\r
2377         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );\r
2378         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2379     }\r
2380     return\r
2381         packFloatx80(\r
2382             aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );\r
2383 \r
2384 }\r
2385 \r
2386 #endif\r
2387 \r
2388 #ifdef FLOAT128\r
2389 \r
2390 /*----------------------------------------------------------------------------\r
2391 | Returns the result of converting the double-precision floating-point value\r
2392 | `a' to the quadruple-precision floating-point format.  The conversion is\r
2393 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2394 | Arithmetic.\r
2395 *----------------------------------------------------------------------------*/\r
2396 \r
2397 float128 float64_to_float128( softfloat_t* sf, float64 a )\r
2398 {\r
2399     flag aSign;\r
2400     int16_t aExp;\r
2401     bits64 aSig, zSig0, zSig1;\r
2402 \r
2403     aSig = extractFloat64Frac( a );\r
2404     aExp = extractFloat64Exp( a );\r
2405     aSign = extractFloat64Sign( a );\r
2406     if ( aExp == 0x7FF ) {\r
2407         if ( aSig ) return commonNaNToFloat128( sf, float64ToCommonNaN( sf, a ) );\r
2408         return packFloat128( aSign, 0x7FFF, 0, 0 );\r
2409     }\r
2410     if ( aExp == 0 ) {\r
2411         if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );\r
2412         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2413         --aExp;\r
2414     }\r
2415     shift128Right( aSig, 0, 4, &zSig0, &zSig1 );\r
2416     return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );\r
2417 \r
2418 }\r
2419 \r
2420 #endif\r
2421 \r
2422 /*----------------------------------------------------------------------------\r
2423 | Rounds the double-precision floating-point value `a' to an integer, and\r
2424 | returns the result as a double-precision floating-point value.  The\r
2425 | operation is performed according to the IEC/IEEE Standard for Binary\r
2426 | Floating-Point Arithmetic.\r
2427 *----------------------------------------------------------------------------*/\r
2428 \r
2429 float64 float64_round_to_int( softfloat_t* sf, float64 a )\r
2430 {\r
2431     flag aSign;\r
2432     int16_t aExp;\r
2433     bits64 lastBitMask, roundBitsMask;\r
2434     int8_t roundingMode;\r
2435     float64 z;\r
2436 \r
2437     aExp = extractFloat64Exp( a );\r
2438     if ( 0x433 <= aExp ) {\r
2439         if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {\r
2440             return propagateFloat64NaN( sf, a, a );\r
2441         }\r
2442         return a;\r
2443     }\r
2444     if ( aExp < 0x3FF ) {\r
2445         if ( (bits64) ( a<<1 ) == 0 ) return a;\r
2446         sf->float_exception_flags |= float_flag_inexact;\r
2447         aSign = extractFloat64Sign( a );\r
2448         switch ( sf->float_rounding_mode ) {\r
2449          case float_round_nearest_even:\r
2450             if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {\r
2451                 return packFloat64( aSign, 0x3FF, 0 );\r
2452             }\r
2453             break;\r
2454          case float_round_down:\r
2455             return aSign ? LIT64( 0xBFF0000000000000 ) : 0;\r
2456          case float_round_up:\r
2457             return\r
2458             aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );\r
2459         }\r
2460         return packFloat64( aSign, 0, 0 );\r
2461     }\r
2462     lastBitMask = 1;\r
2463     lastBitMask <<= 0x433 - aExp;\r
2464     roundBitsMask = lastBitMask - 1;\r
2465     z = a;\r
2466     roundingMode = sf->float_rounding_mode;\r
2467     if ( roundingMode == float_round_nearest_even ) {\r
2468         z += lastBitMask>>1;\r
2469         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;\r
2470     }\r
2471     else if ( roundingMode != float_round_to_zero ) {\r
2472         if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {\r
2473             z += roundBitsMask;\r
2474         }\r
2475     }\r
2476     z &= ~ roundBitsMask;\r
2477     if ( z != a ) sf->float_exception_flags |= float_flag_inexact;\r
2478     return z;\r
2479 \r
2480 }\r
2481 \r
2482 /*----------------------------------------------------------------------------\r
2483 | Returns the result of adding the absolute values of the double-precision\r
2484 | floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated\r
2485 | before being returned.  `zSign' is ignored if the result is a NaN.\r
2486 | The addition is performed according to the IEC/IEEE Standard for Binary\r
2487 | Floating-Point Arithmetic.\r
2488 *----------------------------------------------------------------------------*/\r
2489 \r
2490 float64 addFloat64Sigs( softfloat_t* sf, float64 a, float64 b, flag zSign )\r
2491 {\r
2492     int16_t aExp, bExp, zExp;\r
2493     bits64 aSig, bSig, zSig;\r
2494     int16_t expDiff;\r
2495 \r
2496     aSig = extractFloat64Frac( a );\r
2497     aExp = extractFloat64Exp( a );\r
2498     bSig = extractFloat64Frac( b );\r
2499     bExp = extractFloat64Exp( b );\r
2500     expDiff = aExp - bExp;\r
2501     aSig <<= 9;\r
2502     bSig <<= 9;\r
2503     if ( 0 < expDiff ) {\r
2504         if ( aExp == 0x7FF ) {\r
2505             if ( aSig ) return propagateFloat64NaN( sf, a, b );\r
2506             return a;\r
2507         }\r
2508         if ( bExp == 0 ) {\r
2509             --expDiff;\r
2510         }\r
2511         else {\r
2512             bSig |= LIT64( 0x2000000000000000 );\r
2513         }\r
2514         shift64RightJamming( bSig, expDiff, &bSig );\r
2515         zExp = aExp;\r
2516     }\r
2517     else if ( expDiff < 0 ) {\r
2518         if ( bExp == 0x7FF ) {\r
2519             if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2520             return packFloat64( zSign, 0x7FF, 0 );\r
2521         }\r
2522         if ( aExp == 0 ) {\r
2523             ++expDiff;\r
2524         }\r
2525         else {\r
2526             aSig |= LIT64( 0x2000000000000000 );\r
2527         }\r
2528         shift64RightJamming( aSig, - expDiff, &aSig );\r
2529         zExp = bExp;\r
2530     }\r
2531     else {\r
2532         if ( aExp == 0x7FF ) {\r
2533             if ( aSig | bSig ) return propagateFloat64NaN( sf, a, b );\r
2534             return a;\r
2535         }\r
2536         if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );\r
2537         zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;\r
2538         zExp = aExp;\r
2539         goto roundAndPack;\r
2540     }\r
2541     aSig |= LIT64( 0x2000000000000000 );\r
2542     zSig = ( aSig + bSig )<<1;\r
2543     --zExp;\r
2544     if ( (sbits64) zSig < 0 ) {\r
2545         zSig = aSig + bSig;\r
2546         ++zExp;\r
2547     }\r
2548  roundAndPack:\r
2549     return roundAndPackFloat64( sf, zSign, zExp, zSig );\r
2550 \r
2551 }\r
2552 \r
2553 /*----------------------------------------------------------------------------\r
2554 | Returns the result of subtracting the absolute values of the double-\r
2555 | precision floating-point values `a' and `b'.  If `zSign' is 1, the\r
2556 | difference is negated before being returned.  `zSign' is ignored if the\r
2557 | result is a NaN.  The subtraction is performed according to the IEC/IEEE\r
2558 | Standard for Binary Floating-Point Arithmetic.\r
2559 *----------------------------------------------------------------------------*/\r
2560 \r
2561 float64 subFloat64Sigs( softfloat_t* sf, float64 a, float64 b, flag zSign )\r
2562 {\r
2563     int16_t aExp, bExp, zExp;\r
2564     bits64 aSig, bSig, zSig;\r
2565     int16_t expDiff;\r
2566 \r
2567     aSig = extractFloat64Frac( a );\r
2568     aExp = extractFloat64Exp( a );\r
2569     bSig = extractFloat64Frac( b );\r
2570     bExp = extractFloat64Exp( b );\r
2571     expDiff = aExp - bExp;\r
2572     aSig <<= 10;\r
2573     bSig <<= 10;\r
2574     if ( 0 < expDiff ) goto aExpBigger;\r
2575     if ( expDiff < 0 ) goto bExpBigger;\r
2576     if ( aExp == 0x7FF ) {\r
2577         if ( aSig | bSig ) return propagateFloat64NaN( sf, a, b );\r
2578         float_raise( sf, float_flag_invalid );\r
2579         return float64_default_nan;\r
2580     }\r
2581     if ( aExp == 0 ) {\r
2582         aExp = 1;\r
2583         bExp = 1;\r
2584     }\r
2585     if ( bSig < aSig ) goto aBigger;\r
2586     if ( aSig < bSig ) goto bBigger;\r
2587     return packFloat64( sf->float_rounding_mode == float_round_down, 0, 0 );\r
2588  bExpBigger:\r
2589     if ( bExp == 0x7FF ) {\r
2590         if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2591         return packFloat64( zSign ^ 1, 0x7FF, 0 );\r
2592     }\r
2593     if ( aExp == 0 ) {\r
2594         ++expDiff;\r
2595     }\r
2596     else {\r
2597         aSig |= LIT64( 0x4000000000000000 );\r
2598     }\r
2599     shift64RightJamming( aSig, - expDiff, &aSig );\r
2600     bSig |= LIT64( 0x4000000000000000 );\r
2601  bBigger:\r
2602     zSig = bSig - aSig;\r
2603     zExp = bExp;\r
2604     zSign ^= 1;\r
2605     goto normalizeRoundAndPack;\r
2606  aExpBigger:\r
2607     if ( aExp == 0x7FF ) {\r
2608         if ( aSig ) return propagateFloat64NaN( sf, a, b );\r
2609         return a;\r
2610     }\r
2611     if ( bExp == 0 ) {\r
2612         --expDiff;\r
2613     }\r
2614     else {\r
2615         bSig |= LIT64( 0x4000000000000000 );\r
2616     }\r
2617     shift64RightJamming( bSig, expDiff, &bSig );\r
2618     aSig |= LIT64( 0x4000000000000000 );\r
2619  aBigger:\r
2620     zSig = aSig - bSig;\r
2621     zExp = aExp;\r
2622  normalizeRoundAndPack:\r
2623     --zExp;\r
2624     return normalizeRoundAndPackFloat64( sf, zSign, zExp, zSig );\r
2625 \r
2626 }\r
2627 \r
2628 /*----------------------------------------------------------------------------\r
2629 | Returns the result of adding the double-precision floating-point values `a'\r
2630 | and `b'.  The operation is performed according to the IEC/IEEE Standard for\r
2631 | Binary Floating-Point Arithmetic.\r
2632 *----------------------------------------------------------------------------*/\r
2633 \r
2634 float64 float64_add( softfloat_t* sf, float64 a, float64 b )\r
2635 {\r
2636     flag aSign, bSign;\r
2637 \r
2638     aSign = extractFloat64Sign( a );\r
2639     bSign = extractFloat64Sign( b );\r
2640     if ( aSign == bSign ) {\r
2641         return addFloat64Sigs( sf, a, b, aSign );\r
2642     }\r
2643     else {\r
2644         return subFloat64Sigs( sf, a, b, aSign );\r
2645     }\r
2646 \r
2647 }\r
2648 \r
2649 /*----------------------------------------------------------------------------\r
2650 | Returns the result of subtracting the double-precision floating-point values\r
2651 | `a' and `b'.  The operation is performed according to the IEC/IEEE Standard\r
2652 | for Binary Floating-Point Arithmetic.\r
2653 *----------------------------------------------------------------------------*/\r
2654 \r
2655 float64 float64_sub( softfloat_t* sf, float64 a, float64 b )\r
2656 {\r
2657     flag aSign, bSign;\r
2658 \r
2659     aSign = extractFloat64Sign( a );\r
2660     bSign = extractFloat64Sign( b );\r
2661     if ( aSign == bSign ) {\r
2662         return subFloat64Sigs( sf, a, b, aSign );\r
2663     }\r
2664     else {\r
2665         return addFloat64Sigs( sf, a, b, aSign );\r
2666     }\r
2667 \r
2668 }\r
2669 \r
2670 /*----------------------------------------------------------------------------\r
2671 | Returns the result of multiplying the double-precision floating-point values\r
2672 | `a' and `b'.  The operation is performed according to the IEC/IEEE Standard\r
2673 | for Binary Floating-Point Arithmetic.\r
2674 *----------------------------------------------------------------------------*/\r
2675 \r
2676 float64 float64_mul( softfloat_t* sf, float64 a, float64 b )\r
2677 {\r
2678     flag aSign, bSign, zSign;\r
2679     int16_t aExp, bExp, zExp;\r
2680     bits64 aSig, bSig, zSig0, zSig1;\r
2681 \r
2682     aSig = extractFloat64Frac( a );\r
2683     aExp = extractFloat64Exp( a );\r
2684     aSign = extractFloat64Sign( a );\r
2685     bSig = extractFloat64Frac( b );\r
2686     bExp = extractFloat64Exp( b );\r
2687     bSign = extractFloat64Sign( b );\r
2688     zSign = aSign ^ bSign;\r
2689     if ( aExp == 0x7FF ) {\r
2690         if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {\r
2691             return propagateFloat64NaN( sf, a, b );\r
2692         }\r
2693         if ( ( bExp | bSig ) == 0 ) {\r
2694             float_raise( sf, float_flag_invalid );\r
2695             return float64_default_nan;\r
2696         }\r
2697         return packFloat64( zSign, 0x7FF, 0 );\r
2698     }\r
2699     if ( bExp == 0x7FF ) {\r
2700         if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2701         if ( ( aExp | aSig ) == 0 ) {\r
2702             float_raise( sf, float_flag_invalid );\r
2703             return float64_default_nan;\r
2704         }\r
2705         return packFloat64( zSign, 0x7FF, 0 );\r
2706     }\r
2707     if ( aExp == 0 ) {\r
2708         if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );\r
2709         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2710     }\r
2711     if ( bExp == 0 ) {\r
2712         if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );\r
2713         normalizeFloat64Subnormal( bSig, &bExp, &bSig );\r
2714     }\r
2715     zExp = aExp + bExp - 0x3FF;\r
2716     aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;\r
2717     bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;\r
2718     mul64To128( aSig, bSig, &zSig0, &zSig1 );\r
2719     zSig0 |= ( zSig1 != 0 );\r
2720     if ( 0 <= (sbits64) ( zSig0<<1 ) ) {\r
2721         zSig0 <<= 1;\r
2722         --zExp;\r
2723     }\r
2724 \r
2725     return roundAndPackFloat64( sf, zSign, zExp, zSig0 );\r
2726 \r
2727 }\r
2728 \r
2729 /*----------------------------------------------------------------------------\r
2730 | Returns the result of dividing the double-precision floating-point value `a'\r
2731 | by the corresponding value `b'.  The operation is performed according to\r
2732 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2733 *----------------------------------------------------------------------------*/\r
2734 \r
2735 float64 float64_div( softfloat_t* sf, float64 a, float64 b )\r
2736 {\r
2737     flag aSign, bSign, zSign;\r
2738     int16_t aExp, bExp, zExp;\r
2739     bits64 aSig, bSig, zSig;\r
2740     bits64 rem0, rem1;\r
2741     bits64 term0, term1;\r
2742 \r
2743     aSig = extractFloat64Frac( a );\r
2744     aExp = extractFloat64Exp( a );\r
2745     aSign = extractFloat64Sign( a );\r
2746     bSig = extractFloat64Frac( b );\r
2747     bExp = extractFloat64Exp( b );\r
2748     bSign = extractFloat64Sign( b );\r
2749     zSign = aSign ^ bSign;\r
2750     if ( aExp == 0x7FF ) {\r
2751         if ( aSig ) return propagateFloat64NaN( sf, a, b );\r
2752         if ( bExp == 0x7FF ) {\r
2753             if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2754             float_raise( sf, float_flag_invalid );\r
2755             return float64_default_nan;\r
2756         }\r
2757         return packFloat64( zSign, 0x7FF, 0 );\r
2758     }\r
2759     if ( bExp == 0x7FF ) {\r
2760         if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2761         return packFloat64( zSign, 0, 0 );\r
2762     }\r
2763     if ( bExp == 0 ) {\r
2764         if ( bSig == 0 ) {\r
2765             if ( ( aExp | aSig ) == 0 ) {\r
2766                 float_raise( sf, float_flag_invalid );\r
2767                 return float64_default_nan;\r
2768             }\r
2769             float_raise( sf, float_flag_divbyzero );\r
2770             return packFloat64( zSign, 0x7FF, 0 );\r
2771         }\r
2772         normalizeFloat64Subnormal( bSig, &bExp, &bSig );\r
2773     }\r
2774     if ( aExp == 0 ) {\r
2775         if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );\r
2776         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2777     }\r
2778     zExp = aExp - bExp + 0x3FD;\r
2779     aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;\r
2780     bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;\r
2781     if ( bSig <= ( aSig + aSig ) ) {\r
2782         aSig >>= 1;\r
2783         ++zExp;\r
2784     }\r
2785     zSig = estimateDiv128To64( aSig, 0, bSig );\r
2786     if ( ( zSig & 0x1FF ) <= 2 ) {\r
2787         mul64To128( bSig, zSig, &term0, &term1 );\r
2788         sub128( aSig, 0, term0, term1, &rem0, &rem1 );\r
2789         while ( (sbits64) rem0 < 0 ) {\r
2790             --zSig;\r
2791             add128( rem0, rem1, 0, bSig, &rem0, &rem1 );\r
2792         }\r
2793         zSig |= ( rem1 != 0 );\r
2794     }\r
2795     return roundAndPackFloat64( sf, zSign, zExp, zSig );\r
2796 \r
2797 }\r
2798 \r
2799 /*----------------------------------------------------------------------------\r
2800 | Returns the remainder of the double-precision floating-point value `a'\r
2801 | with respect to the corresponding value `b'.  The operation is performed\r
2802 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2803 *----------------------------------------------------------------------------*/\r
2804 \r
2805 float64 float64_rem( softfloat_t* sf, float64 a, float64 b )\r
2806 {\r
2807     flag aSign, bSign, zSign;\r
2808     int16_t aExp, bExp, expDiff;\r
2809     bits64 aSig, bSig;\r
2810     bits64 q, alternateASig;\r
2811     sbits64 sigMean;\r
2812 \r
2813     aSig = extractFloat64Frac( a );\r
2814     aExp = extractFloat64Exp( a );\r
2815     aSign = extractFloat64Sign( a );\r
2816     bSig = extractFloat64Frac( b );\r
2817     bExp = extractFloat64Exp( b );\r
2818     bSign = extractFloat64Sign( b );\r
2819     if ( aExp == 0x7FF ) {\r
2820         if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {\r
2821             return propagateFloat64NaN( sf, a, b );\r
2822         }\r
2823         float_raise( sf, float_flag_invalid );\r
2824         return float64_default_nan;\r
2825     }\r
2826     if ( bExp == 0x7FF ) {\r
2827         if ( bSig ) return propagateFloat64NaN( sf, a, b );\r
2828         return a;\r
2829     }\r
2830     if ( bExp == 0 ) {\r
2831         if ( bSig == 0 ) {\r
2832             float_raise( sf, float_flag_invalid );\r
2833             return float64_default_nan;\r
2834         }\r
2835         normalizeFloat64Subnormal( bSig, &bExp, &bSig );\r
2836     }\r
2837     if ( aExp == 0 ) {\r
2838         if ( aSig == 0 ) return a;\r
2839         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2840     }\r
2841     expDiff = aExp - bExp;\r
2842     aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;\r
2843     bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;\r
2844     if ( expDiff < 0 ) {\r
2845         if ( expDiff < -1 ) return a;\r
2846         aSig >>= 1;\r
2847     }\r
2848     q = ( bSig <= aSig );\r
2849     if ( q ) aSig -= bSig;\r
2850     expDiff -= 64;\r
2851     while ( 0 < expDiff ) {\r
2852         q = estimateDiv128To64( aSig, 0, bSig );\r
2853         q = ( 2 < q ) ? q - 2 : 0;\r
2854         aSig = - ( ( bSig>>2 ) * q );\r
2855         expDiff -= 62;\r
2856     }\r
2857     expDiff += 64;\r
2858     if ( 0 < expDiff ) {\r
2859         q = estimateDiv128To64( aSig, 0, bSig );\r
2860         q = ( 2 < q ) ? q - 2 : 0;\r
2861         q >>= 64 - expDiff;\r
2862         bSig >>= 2;\r
2863         aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;\r
2864     }\r
2865     else {\r
2866         aSig >>= 2;\r
2867         bSig >>= 2;\r
2868     }\r
2869     do {\r
2870         alternateASig = aSig;\r
2871         ++q;\r
2872         aSig -= bSig;\r
2873     } while ( 0 <= (sbits64) aSig );\r
2874     sigMean = aSig + alternateASig;\r
2875     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {\r
2876         aSig = alternateASig;\r
2877     }\r
2878     zSign = ( (sbits64) aSig < 0 );\r
2879     if ( zSign ) aSig = - aSig;\r
2880     return normalizeRoundAndPackFloat64( sf, aSign ^ zSign, bExp, aSig );\r
2881 \r
2882 }\r
2883 \r
2884 /*----------------------------------------------------------------------------\r
2885 | Returns the square root of the double-precision floating-point value `a'.\r
2886 | The operation is performed according to the IEC/IEEE Standard for Binary\r
2887 | Floating-Point Arithmetic.\r
2888 *----------------------------------------------------------------------------*/\r
2889 \r
2890 float64 float64_sqrt( softfloat_t* sf, float64 a )\r
2891 {\r
2892     flag aSign;\r
2893     int16_t aExp, zExp;\r
2894     bits64 aSig, zSig, doubleZSig;\r
2895     bits64 rem0, rem1, term0, term1;\r
2896     float64 z;\r
2897 \r
2898     aSig = extractFloat64Frac( a );\r
2899     aExp = extractFloat64Exp( a );\r
2900     aSign = extractFloat64Sign( a );\r
2901     if ( aExp == 0x7FF ) {\r
2902         if ( aSig ) return propagateFloat64NaN( sf, a, a );\r
2903         if ( ! aSign ) return a;\r
2904         float_raise( sf, float_flag_invalid );\r
2905         return float64_default_nan;\r
2906     }\r
2907     if ( aSign ) {\r
2908         if ( ( aExp | aSig ) == 0 ) return a;\r
2909         float_raise( sf, float_flag_invalid );\r
2910         return float64_default_nan;\r
2911     }\r
2912     if ( aExp == 0 ) {\r
2913         if ( aSig == 0 ) return 0;\r
2914         normalizeFloat64Subnormal( aSig, &aExp, &aSig );\r
2915     }\r
2916     zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;\r
2917     aSig |= LIT64( 0x0010000000000000 );\r
2918     zSig = estimateSqrt32( aExp, aSig>>21 );\r
2919     aSig <<= 9 - ( aExp & 1 );\r
2920     zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );\r
2921     if ( ( zSig & 0x1FF ) <= 5 ) {\r
2922         doubleZSig = zSig<<1;\r
2923         mul64To128( zSig, zSig, &term0, &term1 );\r
2924         sub128( aSig, 0, term0, term1, &rem0, &rem1 );\r
2925         while ( (sbits64) rem0 < 0 ) {\r
2926             --zSig;\r
2927             doubleZSig -= 2;\r
2928             add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );\r
2929         }\r
2930         zSig |= ( ( rem0 | rem1 ) != 0 );\r
2931     }\r
2932     return roundAndPackFloat64( sf, 0, zExp, zSig );\r
2933 \r
2934 }\r
2935 \r
2936 /*----------------------------------------------------------------------------\r
2937 | Returns 1 if the double-precision floating-point value `a' is equal to the\r
2938 | corresponding value `b', and 0 otherwise.  The comparison is performed\r
2939 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2940 *----------------------------------------------------------------------------*/\r
2941 \r
2942 flag float64_eq( softfloat_t* sf, float64 a, float64 b )\r
2943 {\r
2944 \r
2945     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
2946          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
2947        ) {\r
2948         if ( float64_is_signaling_nan( sf, a ) || float64_is_signaling_nan( sf, b ) ) {\r
2949             float_raise( sf, float_flag_invalid );\r
2950         }\r
2951         return 0;\r
2952     }\r
2953     return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );\r
2954 \r
2955 }\r
2956 \r
2957 /*----------------------------------------------------------------------------\r
2958 | Returns 1 if the double-precision floating-point value `a' is less than or\r
2959 | equal to the corresponding value `b', and 0 otherwise.  The comparison is\r
2960 | performed according to the IEC/IEEE Standard for Binary Floating-Point\r
2961 | Arithmetic.\r
2962 *----------------------------------------------------------------------------*/\r
2963 \r
2964 flag float64_le( softfloat_t* sf, float64 a, float64 b )\r
2965 {\r
2966     flag aSign, bSign;\r
2967 \r
2968     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
2969          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
2970        ) {\r
2971         float_raise( sf, float_flag_invalid );\r
2972         return 0;\r
2973     }\r
2974     aSign = extractFloat64Sign( a );\r
2975     bSign = extractFloat64Sign( b );\r
2976     if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );\r
2977     return ( a == b ) || ( aSign ^ ( a < b ) );\r
2978 \r
2979 }\r
2980 \r
2981 /*----------------------------------------------------------------------------\r
2982 | Returns 1 if the double-precision floating-point value `a' is less than\r
2983 | the corresponding value `b', and 0 otherwise.  The comparison is performed\r
2984 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
2985 *----------------------------------------------------------------------------*/\r
2986 \r
2987 flag float64_lt( softfloat_t* sf, float64 a, float64 b )\r
2988 {\r
2989     flag aSign, bSign;\r
2990 \r
2991     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
2992          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
2993        ) {\r
2994         float_raise( sf, float_flag_invalid );\r
2995         return 0;\r
2996     }\r
2997     aSign = extractFloat64Sign( a );\r
2998     bSign = extractFloat64Sign( b );\r
2999     if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );\r
3000     return ( a != b ) && ( aSign ^ ( a < b ) );\r
3001 \r
3002 }\r
3003 \r
3004 /*----------------------------------------------------------------------------\r
3005 | Returns 1 if the double-precision floating-point value `a' is equal to the\r
3006 | corresponding value `b', and 0 otherwise.  The invalid exception is raised\r
3007 | if either operand is a NaN.  Otherwise, the comparison is performed\r
3008 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
3009 *----------------------------------------------------------------------------*/\r
3010 \r
3011 flag float64_eq_signaling( softfloat_t* sf, float64 a, float64 b )\r
3012 {\r
3013 \r
3014     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
3015          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
3016        ) {\r
3017         float_raise( sf, float_flag_invalid );\r
3018         return 0;\r
3019     }\r
3020     return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );\r
3021 \r
3022 }\r
3023 \r
3024 /*----------------------------------------------------------------------------\r
3025 | Returns 1 if the double-precision floating-point value `a' is less than or\r
3026 | equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not\r
3027 | cause an exception.  Otherwise, the comparison is performed according to the\r
3028 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.\r
3029 *----------------------------------------------------------------------------*/\r
3030 \r
3031 flag float64_le_quiet( softfloat_t* sf, float64 a, float64 b )\r
3032 {\r
3033     flag aSign, bSign;\r
3034     int16_t aExp, bExp;\r
3035 \r
3036     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
3037          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
3038        ) {\r
3039         if ( float64_is_signaling_nan( sf, a ) || float64_is_signaling_nan( sf, b ) ) {\r
3040             float_raise( sf, float_flag_invalid );\r
3041         }\r
3042         return 0;\r
3043     }\r
3044     aSign = extractFloat64Sign( a );\r
3045     bSign = extractFloat64Sign( b );\r
3046     if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );\r
3047     return ( a == b ) || ( aSign ^ ( a < b ) );\r
3048 \r
3049 }\r
3050 \r
3051 /*----------------------------------------------------------------------------\r
3052 | Returns 1 if the double-precision floating-point value `a' is less than\r
3053 | the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an\r
3054 | exception.  Otherwise, the comparison is performed according to the IEC/IEEE\r
3055 | Standard for Binary Floating-Point Arithmetic.\r
3056 *----------------------------------------------------------------------------*/\r
3057 \r
3058 flag float64_lt_quiet( softfloat_t* sf, float64 a, float64 b )\r
3059 {\r
3060     flag aSign, bSign;\r
3061 \r
3062     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )\r
3063          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )\r
3064        ) {\r
3065         if ( float64_is_signaling_nan( sf, a ) || float64_is_signaling_nan( sf, b ) ) {\r
3066             float_raise( sf, float_flag_invalid );\r
3067         }\r
3068         return 0;\r
3069     }\r
3070     aSign = extractFloat64Sign( a );\r
3071     bSign = extractFloat64Sign( b );\r
3072     if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );\r
3073     return ( a != b ) && ( aSign ^ ( a < b ) );\r
3074 \r
3075 }\r
3076 \r
3077 #ifdef FLOATX80\r
3078 \r
3079 /*----------------------------------------------------------------------------\r
3080 | Returns the result of converting the extended double-precision floating-\r
3081 | point value `a' to the 32-bit two's complement integer format.  The\r
3082 | conversion is performed according to the IEC/IEEE Standard for Binary\r
3083 | Floating-Point Arithmetic---which means in particular that the conversion\r
3084 | is rounded according to the current rounding mode.  If `a' is a NaN, the\r
3085 | largest positive integer is returned.  Otherwise, if the conversion\r
3086 | overflows, the largest integer with the same sign as `a' is returned.\r
3087 *----------------------------------------------------------------------------*/\r
3088 \r
3089 int32_t floatx80_to_int32( softfloat_t* sf, floatx80 a )\r
3090 {\r
3091     flag aSign;\r
3092     int32_t aExp, shiftCount;\r
3093     bits64 aSig;\r
3094 \r
3095     aSig = extractFloatx80Frac( a );\r
3096     aExp = extractFloatx80Exp( a );\r
3097     aSign = extractFloatx80Sign( a );\r
3098     if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;\r
3099     shiftCount = 0x4037 - aExp;\r
3100     if ( shiftCount <= 0 ) shiftCount = 1;\r
3101     shift64RightJamming( aSig, shiftCount, &aSig );\r
3102     return roundAndPackInt32( sf, aSign, aSig );\r
3103 \r
3104 }\r
3105 \r
3106 /*----------------------------------------------------------------------------\r
3107 | Returns the result of converting the extended double-precision floating-\r
3108 | point value `a' to the 32-bit two's complement integer format.  The\r
3109 | conversion is performed according to the IEC/IEEE Standard for Binary\r
3110 | Floating-Point Arithmetic, except that the conversion is always rounded\r
3111 | toward zero.  If `a' is a NaN, the largest positive integer is returned.\r
3112 | Otherwise, if the conversion overflows, the largest integer with the same\r
3113 | sign as `a' is returned.\r
3114 *----------------------------------------------------------------------------*/\r
3115 \r
3116 int32_t floatx80_to_int32_round_to_zero( softfloat_t* sf, floatx80 a )\r
3117 {\r
3118     flag aSign;\r
3119     int32_t aExp, shiftCount;\r
3120     bits64 aSig, savedASig;\r
3121     int32_t z;\r
3122 \r
3123     aSig = extractFloatx80Frac( a );\r
3124     aExp = extractFloatx80Exp( a );\r
3125     aSign = extractFloatx80Sign( a );\r
3126     if ( 0x401E < aExp ) {\r
3127         if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;\r
3128         goto invalid;\r
3129     }\r
3130     else if ( aExp < 0x3FFF ) {\r
3131         if ( aExp || aSig ) sf->float_exception_flags |= float_flag_inexact;\r
3132         return 0;\r
3133     }\r
3134     shiftCount = 0x403E - aExp;\r
3135     savedASig = aSig;\r
3136     aSig >>= shiftCount;\r
3137     z = aSig;\r
3138     if ( aSign ) z = - z;\r
3139     if ( ( z < 0 ) ^ aSign ) {\r
3140  invalid:\r
3141         float_raise( sf, float_flag_invalid );\r
3142         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;\r
3143     }\r
3144     if ( ( aSig<<shiftCount ) != savedASig ) {\r
3145         sf->float_exception_flags |= float_flag_inexact;\r
3146     }\r
3147     return z;\r
3148 \r
3149 }\r
3150 \r
3151 /*----------------------------------------------------------------------------\r
3152 | Returns the result of converting the extended double-precision floating-\r
3153 | point value `a' to the 64-bit two's complement integer format.  The\r
3154 | conversion is performed according to the IEC/IEEE Standard for Binary\r
3155 | Floating-Point Arithmetic---which means in particular that the conversion\r
3156 | is rounded according to the current rounding mode.  If `a' is a NaN,\r
3157 | the largest positive integer is returned.  Otherwise, if the conversion\r
3158 | overflows, the largest integer with the same sign as `a' is returned.\r
3159 *----------------------------------------------------------------------------*/\r
3160 \r
3161 int64_t floatx80_to_int64( softfloat_t* sf, floatx80 a )\r
3162 {\r
3163     flag aSign;\r
3164     int32_t aExp, shiftCount;\r
3165     bits64 aSig, aSigExtra;\r
3166 \r
3167     aSig = extractFloatx80Frac( a );\r
3168     aExp = extractFloatx80Exp( a );\r
3169     aSign = extractFloatx80Sign( a );\r
3170     shiftCount = 0x403E - aExp;\r
3171     if ( shiftCount <= 0 ) {\r
3172         if ( shiftCount ) {\r
3173             float_raise( sf, float_flag_invalid );\r
3174             if (    ! aSign\r
3175                  || (    ( aExp == 0x7FFF )\r
3176                       && ( aSig != LIT64( 0x8000000000000000 ) ) )\r
3177                ) {\r
3178                 return LIT64( 0x7FFFFFFFFFFFFFFF );\r
3179             }\r
3180             return (sbits64) LIT64( 0x8000000000000000 );\r
3181         }\r
3182         aSigExtra = 0;\r
3183     }\r
3184     else {\r
3185         shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );\r
3186     }\r
3187     return roundAndPackInt64( sf, aSign, aSig, aSigExtra );\r
3188 \r
3189 }\r
3190 \r
3191 /*----------------------------------------------------------------------------\r
3192 | Returns the result of converting the extended double-precision floating-\r
3193 | point value `a' to the 64-bit two's complement integer format.  The\r
3194 | conversion is performed according to the IEC/IEEE Standard for Binary\r
3195 | Floating-Point Arithmetic, except that the conversion is always rounded\r
3196 | toward zero.  If `a' is a NaN, the largest positive integer is returned.\r
3197 | Otherwise, if the conversion overflows, the largest integer with the same\r
3198 | sign as `a' is returned.\r
3199 *----------------------------------------------------------------------------*/\r
3200 \r
3201 int64_t floatx80_to_int64_round_to_zero( softfloat_t* sf, floatx80 a )\r
3202 {\r
3203     flag aSign;\r
3204     int32_t aExp, shiftCount;\r
3205     bits64 aSig;\r
3206     int64_t z;\r
3207 \r
3208     aSig = extractFloatx80Frac( a );\r
3209     aExp = extractFloatx80Exp( a );\r