4 * Derived from SoftFloat.
7 /*============================================================================
9 This C source file is part of the SoftFloat IEC/IEEE Floating-point Arithmetic
12 Written by John R. Hauser. This work was made possible in part by the
13 International Computer Science Institute, located at Suite 600, 1947 Center
14 Street, Berkeley, California 94704. Funding was partially provided by the
15 National Science Foundation under grant MIP-9311980. The original version
16 of this code was written as part of a project to build a fixed-point vector
17 processor in collaboration with the University of California at Berkeley,
18 overseen by Profs. Nelson Morgan and John Wawrzynek. More information
19 is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
20 arithmetic/SoftFloat.html'.
22 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
23 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
24 RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
25 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
26 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
27 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
28 INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
29 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
31 Derivative works are acceptable, even for commercial purposes, so long as
32 (1) the source code for the derivative work includes prominent notice that
33 the work is derivative, and (2) the source code includes prominent notice with
34 these four paragraphs for those parts of this code that are retained.
36 =============================================================================*/
38 /* softfloat (and in particular the code in softfloat-specialize.h) is
39 * target-dependent and needs the TARGET_* macros.
43 #include "softfloat.h"
45 /*----------------------------------------------------------------------------
46 | Primitive arithmetic functions, including multi-word arithmetic, and
47 | division and square root approximations. (Can be specialized to target if
49 *----------------------------------------------------------------------------*/
50 #include "softfloat-macros.h"
52 /*----------------------------------------------------------------------------
53 | Functions and definitions to determine: (1) whether tininess for underflow
54 | is detected before or after rounding by default, (2) what (if anything)
55 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
56 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
57 | are propagated from function inputs to output. These details are target-
59 *----------------------------------------------------------------------------*/
60 #include "softfloat-specialize.h"
62 void set_float_rounding_mode(int val STATUS_PARAM
)
64 STATUS(float_rounding_mode
) = val
;
67 void set_float_exception_flags(int val STATUS_PARAM
)
69 STATUS(float_exception_flags
) = val
;
72 void set_floatx80_rounding_precision(int val STATUS_PARAM
)
74 STATUS(floatx80_rounding_precision
) = val
;
77 /*----------------------------------------------------------------------------
78 | Returns the fraction bits of the half-precision floating-point value `a'.
79 *----------------------------------------------------------------------------*/
81 INLINE
uint32_t extractFloat16Frac(float16 a
)
83 return float16_val(a
) & 0x3ff;
86 /*----------------------------------------------------------------------------
87 | Returns the exponent bits of the half-precision floating-point value `a'.
88 *----------------------------------------------------------------------------*/
90 INLINE int16
extractFloat16Exp(float16 a
)
92 return (float16_val(a
) >> 10) & 0x1f;
95 /*----------------------------------------------------------------------------
96 | Returns the sign bit of the single-precision floating-point value `a'.
97 *----------------------------------------------------------------------------*/
99 INLINE flag
extractFloat16Sign(float16 a
)
101 return float16_val(a
)>>15;
104 /*----------------------------------------------------------------------------
105 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
106 | and 7, and returns the properly rounded 32-bit integer corresponding to the
107 | input. If `zSign' is 1, the input is negated before being converted to an
108 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
109 | is simply rounded to an integer, with the inexact exception raised if the
110 | input cannot be represented exactly as an integer. However, if the fixed-
111 | point input is too large, the invalid exception is raised and the largest
112 | positive or negative integer is returned.
113 *----------------------------------------------------------------------------*/
115 static int32
roundAndPackInt32( flag zSign
, uint64_t absZ STATUS_PARAM
)
118 flag roundNearestEven
;
119 int8 roundIncrement
, roundBits
;
122 roundingMode
= STATUS(float_rounding_mode
);
123 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
124 roundIncrement
= 0x40;
125 if ( ! roundNearestEven
) {
126 if ( roundingMode
== float_round_to_zero
) {
130 roundIncrement
= 0x7F;
132 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
135 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
139 roundBits
= absZ
& 0x7F;
140 absZ
= ( absZ
+ roundIncrement
)>>7;
141 absZ
&= ~ ( ( ( roundBits
^ 0x40 ) == 0 ) & roundNearestEven
);
143 if ( zSign
) z
= - z
;
144 if ( ( absZ
>>32 ) || ( z
&& ( ( z
< 0 ) ^ zSign
) ) ) {
145 float_raise( float_flag_invalid STATUS_VAR
);
146 return zSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
148 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
153 /*----------------------------------------------------------------------------
154 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
155 | `absZ1', with binary point between bits 63 and 64 (between the input words),
156 | and returns the properly rounded 64-bit integer corresponding to the input.
157 | If `zSign' is 1, the input is negated before being converted to an integer.
158 | Ordinarily, the fixed-point input is simply rounded to an integer, with
159 | the inexact exception raised if the input cannot be represented exactly as
160 | an integer. However, if the fixed-point input is too large, the invalid
161 | exception is raised and the largest positive or negative integer is
163 *----------------------------------------------------------------------------*/
165 static int64
roundAndPackInt64( flag zSign
, uint64_t absZ0
, uint64_t absZ1 STATUS_PARAM
)
168 flag roundNearestEven
, increment
;
171 roundingMode
= STATUS(float_rounding_mode
);
172 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
173 increment
= ( (int64_t) absZ1
< 0 );
174 if ( ! roundNearestEven
) {
175 if ( roundingMode
== float_round_to_zero
) {
180 increment
= ( roundingMode
== float_round_down
) && absZ1
;
183 increment
= ( roundingMode
== float_round_up
) && absZ1
;
189 if ( absZ0
== 0 ) goto overflow
;
190 absZ0
&= ~ ( ( (uint64_t) ( absZ1
<<1 ) == 0 ) & roundNearestEven
);
193 if ( zSign
) z
= - z
;
194 if ( z
&& ( ( z
< 0 ) ^ zSign
) ) {
196 float_raise( float_flag_invalid STATUS_VAR
);
198 zSign
? (int64_t) LIT64( 0x8000000000000000 )
199 : LIT64( 0x7FFFFFFFFFFFFFFF );
201 if ( absZ1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
206 /*----------------------------------------------------------------------------
207 | Returns the fraction bits of the single-precision floating-point value `a'.
208 *----------------------------------------------------------------------------*/
210 INLINE
uint32_t extractFloat32Frac( float32 a
)
213 return float32_val(a
) & 0x007FFFFF;
217 /*----------------------------------------------------------------------------
218 | Returns the exponent bits of the single-precision floating-point value `a'.
219 *----------------------------------------------------------------------------*/
221 INLINE int16
extractFloat32Exp( float32 a
)
224 return ( float32_val(a
)>>23 ) & 0xFF;
228 /*----------------------------------------------------------------------------
229 | Returns the sign bit of the single-precision floating-point value `a'.
230 *----------------------------------------------------------------------------*/
232 INLINE flag
extractFloat32Sign( float32 a
)
235 return float32_val(a
)>>31;
239 /*----------------------------------------------------------------------------
240 | If `a' is denormal and we are in flush-to-zero mode then set the
241 | input-denormal exception and return zero. Otherwise just return the value.
242 *----------------------------------------------------------------------------*/
243 static float32
float32_squash_input_denormal(float32 a STATUS_PARAM
)
245 if (STATUS(flush_inputs_to_zero
)) {
246 if (extractFloat32Exp(a
) == 0 && extractFloat32Frac(a
) != 0) {
247 float_raise(float_flag_input_denormal STATUS_VAR
);
248 return make_float32(float32_val(a
) & 0x80000000);
254 /*----------------------------------------------------------------------------
255 | Normalizes the subnormal single-precision floating-point value represented
256 | by the denormalized significand `aSig'. The normalized exponent and
257 | significand are stored at the locations pointed to by `zExpPtr' and
258 | `zSigPtr', respectively.
259 *----------------------------------------------------------------------------*/
262 normalizeFloat32Subnormal( uint32_t aSig
, int16
*zExpPtr
, uint32_t *zSigPtr
)
266 shiftCount
= countLeadingZeros32( aSig
) - 8;
267 *zSigPtr
= aSig
<<shiftCount
;
268 *zExpPtr
= 1 - shiftCount
;
272 /*----------------------------------------------------------------------------
273 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
274 | single-precision floating-point value, returning the result. After being
275 | shifted into the proper positions, the three fields are simply added
276 | together to form the result. This means that any integer portion of `zSig'
277 | will be added into the exponent. Since a properly normalized significand
278 | will have an integer portion equal to 1, the `zExp' input should be 1 less
279 | than the desired result exponent whenever `zSig' is a complete, normalized
281 *----------------------------------------------------------------------------*/
283 INLINE float32
packFloat32( flag zSign
, int16 zExp
, uint32_t zSig
)
287 ( ( (uint32_t) zSign
)<<31 ) + ( ( (uint32_t) zExp
)<<23 ) + zSig
);
291 /*----------------------------------------------------------------------------
292 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
293 | and significand `zSig', and returns the proper single-precision floating-
294 | point value corresponding to the abstract input. Ordinarily, the abstract
295 | value is simply rounded and packed into the single-precision format, with
296 | the inexact exception raised if the abstract input cannot be represented
297 | exactly. However, if the abstract value is too large, the overflow and
298 | inexact exceptions are raised and an infinity or maximal finite value is
299 | returned. If the abstract value is too small, the input value is rounded to
300 | a subnormal number, and the underflow and inexact exceptions are raised if
301 | the abstract input cannot be represented exactly as a subnormal single-
302 | precision floating-point number.
303 | The input significand `zSig' has its binary point between bits 30
304 | and 29, which is 7 bits to the left of the usual location. This shifted
305 | significand must be normalized or smaller. If `zSig' is not normalized,
306 | `zExp' must be 0; in that case, the result returned is a subnormal number,
307 | and it must not require rounding. In the usual case that `zSig' is
308 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
309 | The handling of underflow and overflow follows the IEC/IEEE Standard for
310 | Binary Floating-Point Arithmetic.
311 *----------------------------------------------------------------------------*/
313 static float32
roundAndPackFloat32( flag zSign
, int16 zExp
, uint32_t zSig STATUS_PARAM
)
316 flag roundNearestEven
;
317 int8 roundIncrement
, roundBits
;
320 roundingMode
= STATUS(float_rounding_mode
);
321 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
322 roundIncrement
= 0x40;
323 if ( ! roundNearestEven
) {
324 if ( roundingMode
== float_round_to_zero
) {
328 roundIncrement
= 0x7F;
330 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
333 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
337 roundBits
= zSig
& 0x7F;
338 if ( 0xFD <= (uint16_t) zExp
) {
340 || ( ( zExp
== 0xFD )
341 && ( (int32_t) ( zSig
+ roundIncrement
) < 0 ) )
343 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
344 return packFloat32( zSign
, 0xFF, - ( roundIncrement
== 0 ));
347 if (STATUS(flush_to_zero
)) {
348 float_raise(float_flag_output_denormal STATUS_VAR
);
349 return packFloat32(zSign
, 0, 0);
352 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
354 || ( zSig
+ roundIncrement
< 0x80000000 );
355 shift32RightJamming( zSig
, - zExp
, &zSig
);
357 roundBits
= zSig
& 0x7F;
358 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
361 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
362 zSig
= ( zSig
+ roundIncrement
)>>7;
363 zSig
&= ~ ( ( ( roundBits
^ 0x40 ) == 0 ) & roundNearestEven
);
364 if ( zSig
== 0 ) zExp
= 0;
365 return packFloat32( zSign
, zExp
, zSig
);
369 /*----------------------------------------------------------------------------
370 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
371 | and significand `zSig', and returns the proper single-precision floating-
372 | point value corresponding to the abstract input. This routine is just like
373 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
374 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
375 | floating-point exponent.
376 *----------------------------------------------------------------------------*/
379 normalizeRoundAndPackFloat32( flag zSign
, int16 zExp
, uint32_t zSig STATUS_PARAM
)
383 shiftCount
= countLeadingZeros32( zSig
) - 1;
384 return roundAndPackFloat32( zSign
, zExp
- shiftCount
, zSig
<<shiftCount STATUS_VAR
);
388 /*----------------------------------------------------------------------------
389 | Returns the fraction bits of the double-precision floating-point value `a'.
390 *----------------------------------------------------------------------------*/
392 INLINE
uint64_t extractFloat64Frac( float64 a
)
395 return float64_val(a
) & LIT64( 0x000FFFFFFFFFFFFF );
399 /*----------------------------------------------------------------------------
400 | Returns the exponent bits of the double-precision floating-point value `a'.
401 *----------------------------------------------------------------------------*/
403 INLINE int16
extractFloat64Exp( float64 a
)
406 return ( float64_val(a
)>>52 ) & 0x7FF;
410 /*----------------------------------------------------------------------------
411 | Returns the sign bit of the double-precision floating-point value `a'.
412 *----------------------------------------------------------------------------*/
414 INLINE flag
extractFloat64Sign( float64 a
)
417 return float64_val(a
)>>63;
421 /*----------------------------------------------------------------------------
422 | If `a' is denormal and we are in flush-to-zero mode then set the
423 | input-denormal exception and return zero. Otherwise just return the value.
424 *----------------------------------------------------------------------------*/
425 static float64
float64_squash_input_denormal(float64 a STATUS_PARAM
)
427 if (STATUS(flush_inputs_to_zero
)) {
428 if (extractFloat64Exp(a
) == 0 && extractFloat64Frac(a
) != 0) {
429 float_raise(float_flag_input_denormal STATUS_VAR
);
430 return make_float64(float64_val(a
) & (1ULL << 63));
436 /*----------------------------------------------------------------------------
437 | Normalizes the subnormal double-precision floating-point value represented
438 | by the denormalized significand `aSig'. The normalized exponent and
439 | significand are stored at the locations pointed to by `zExpPtr' and
440 | `zSigPtr', respectively.
441 *----------------------------------------------------------------------------*/
444 normalizeFloat64Subnormal( uint64_t aSig
, int16
*zExpPtr
, uint64_t *zSigPtr
)
448 shiftCount
= countLeadingZeros64( aSig
) - 11;
449 *zSigPtr
= aSig
<<shiftCount
;
450 *zExpPtr
= 1 - shiftCount
;
454 /*----------------------------------------------------------------------------
455 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
456 | double-precision floating-point value, returning the result. After being
457 | shifted into the proper positions, the three fields are simply added
458 | together to form the result. This means that any integer portion of `zSig'
459 | will be added into the exponent. Since a properly normalized significand
460 | will have an integer portion equal to 1, the `zExp' input should be 1 less
461 | than the desired result exponent whenever `zSig' is a complete, normalized
463 *----------------------------------------------------------------------------*/
465 INLINE float64
packFloat64( flag zSign
, int16 zExp
, uint64_t zSig
)
469 ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<52 ) + zSig
);
473 /*----------------------------------------------------------------------------
474 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
475 | and significand `zSig', and returns the proper double-precision floating-
476 | point value corresponding to the abstract input. Ordinarily, the abstract
477 | value is simply rounded and packed into the double-precision format, with
478 | the inexact exception raised if the abstract input cannot be represented
479 | exactly. However, if the abstract value is too large, the overflow and
480 | inexact exceptions are raised and an infinity or maximal finite value is
481 | returned. If the abstract value is too small, the input value is rounded
482 | to a subnormal number, and the underflow and inexact exceptions are raised
483 | if the abstract input cannot be represented exactly as a subnormal double-
484 | precision floating-point number.
485 | The input significand `zSig' has its binary point between bits 62
486 | and 61, which is 10 bits to the left of the usual location. This shifted
487 | significand must be normalized or smaller. If `zSig' is not normalized,
488 | `zExp' must be 0; in that case, the result returned is a subnormal number,
489 | and it must not require rounding. In the usual case that `zSig' is
490 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
491 | The handling of underflow and overflow follows the IEC/IEEE Standard for
492 | Binary Floating-Point Arithmetic.
493 *----------------------------------------------------------------------------*/
495 static float64
roundAndPackFloat64( flag zSign
, int16 zExp
, uint64_t zSig STATUS_PARAM
)
498 flag roundNearestEven
;
499 int16 roundIncrement
, roundBits
;
502 roundingMode
= STATUS(float_rounding_mode
);
503 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
504 roundIncrement
= 0x200;
505 if ( ! roundNearestEven
) {
506 if ( roundingMode
== float_round_to_zero
) {
510 roundIncrement
= 0x3FF;
512 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
515 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
519 roundBits
= zSig
& 0x3FF;
520 if ( 0x7FD <= (uint16_t) zExp
) {
521 if ( ( 0x7FD < zExp
)
522 || ( ( zExp
== 0x7FD )
523 && ( (int64_t) ( zSig
+ roundIncrement
) < 0 ) )
525 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
526 return packFloat64( zSign
, 0x7FF, - ( roundIncrement
== 0 ));
529 if (STATUS(flush_to_zero
)) {
530 float_raise(float_flag_output_denormal STATUS_VAR
);
531 return packFloat64(zSign
, 0, 0);
534 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
536 || ( zSig
+ roundIncrement
< LIT64( 0x8000000000000000 ) );
537 shift64RightJamming( zSig
, - zExp
, &zSig
);
539 roundBits
= zSig
& 0x3FF;
540 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
543 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
544 zSig
= ( zSig
+ roundIncrement
)>>10;
545 zSig
&= ~ ( ( ( roundBits
^ 0x200 ) == 0 ) & roundNearestEven
);
546 if ( zSig
== 0 ) zExp
= 0;
547 return packFloat64( zSign
, zExp
, zSig
);
551 /*----------------------------------------------------------------------------
552 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
553 | and significand `zSig', and returns the proper double-precision floating-
554 | point value corresponding to the abstract input. This routine is just like
555 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
556 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
557 | floating-point exponent.
558 *----------------------------------------------------------------------------*/
561 normalizeRoundAndPackFloat64( flag zSign
, int16 zExp
, uint64_t zSig STATUS_PARAM
)
565 shiftCount
= countLeadingZeros64( zSig
) - 1;
566 return roundAndPackFloat64( zSign
, zExp
- shiftCount
, zSig
<<shiftCount STATUS_VAR
);
570 /*----------------------------------------------------------------------------
571 | Returns the fraction bits of the extended double-precision floating-point
573 *----------------------------------------------------------------------------*/
575 INLINE
uint64_t extractFloatx80Frac( floatx80 a
)
582 /*----------------------------------------------------------------------------
583 | Returns the exponent bits of the extended double-precision floating-point
585 *----------------------------------------------------------------------------*/
587 INLINE int32
extractFloatx80Exp( floatx80 a
)
590 return a
.high
& 0x7FFF;
594 /*----------------------------------------------------------------------------
595 | Returns the sign bit of the extended double-precision floating-point value
597 *----------------------------------------------------------------------------*/
599 INLINE flag
extractFloatx80Sign( floatx80 a
)
606 /*----------------------------------------------------------------------------
607 | Normalizes the subnormal extended double-precision floating-point value
608 | represented by the denormalized significand `aSig'. The normalized exponent
609 | and significand are stored at the locations pointed to by `zExpPtr' and
610 | `zSigPtr', respectively.
611 *----------------------------------------------------------------------------*/
614 normalizeFloatx80Subnormal( uint64_t aSig
, int32
*zExpPtr
, uint64_t *zSigPtr
)
618 shiftCount
= countLeadingZeros64( aSig
);
619 *zSigPtr
= aSig
<<shiftCount
;
620 *zExpPtr
= 1 - shiftCount
;
624 /*----------------------------------------------------------------------------
625 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
626 | extended double-precision floating-point value, returning the result.
627 *----------------------------------------------------------------------------*/
629 INLINE floatx80
packFloatx80( flag zSign
, int32 zExp
, uint64_t zSig
)
634 z
.high
= ( ( (uint16_t) zSign
)<<15 ) + zExp
;
639 /*----------------------------------------------------------------------------
640 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
641 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
642 | and returns the proper extended double-precision floating-point value
643 | corresponding to the abstract input. Ordinarily, the abstract value is
644 | rounded and packed into the extended double-precision format, with the
645 | inexact exception raised if the abstract input cannot be represented
646 | exactly. However, if the abstract value is too large, the overflow and
647 | inexact exceptions are raised and an infinity or maximal finite value is
648 | returned. If the abstract value is too small, the input value is rounded to
649 | a subnormal number, and the underflow and inexact exceptions are raised if
650 | the abstract input cannot be represented exactly as a subnormal extended
651 | double-precision floating-point number.
652 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
653 | number of bits as single or double precision, respectively. Otherwise, the
654 | result is rounded to the full precision of the extended double-precision
656 | The input significand must be normalized or smaller. If the input
657 | significand is not normalized, `zExp' must be 0; in that case, the result
658 | returned is a subnormal number, and it must not require rounding. The
659 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
660 | Floating-Point Arithmetic.
661 *----------------------------------------------------------------------------*/
664 roundAndPackFloatx80(
665 int8 roundingPrecision
, flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
669 flag roundNearestEven
, increment
, isTiny
;
670 int64 roundIncrement
, roundMask
, roundBits
;
672 roundingMode
= STATUS(float_rounding_mode
);
673 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
674 if ( roundingPrecision
== 80 ) goto precision80
;
675 if ( roundingPrecision
== 64 ) {
676 roundIncrement
= LIT64( 0x0000000000000400 );
677 roundMask
= LIT64( 0x00000000000007FF );
679 else if ( roundingPrecision
== 32 ) {
680 roundIncrement
= LIT64( 0x0000008000000000 );
681 roundMask
= LIT64( 0x000000FFFFFFFFFF );
686 zSig0
|= ( zSig1
!= 0 );
687 if ( ! roundNearestEven
) {
688 if ( roundingMode
== float_round_to_zero
) {
692 roundIncrement
= roundMask
;
694 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
697 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
701 roundBits
= zSig0
& roundMask
;
702 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
703 if ( ( 0x7FFE < zExp
)
704 || ( ( zExp
== 0x7FFE ) && ( zSig0
+ roundIncrement
< zSig0
) )
709 if (STATUS(flush_to_zero
)) {
710 float_raise(float_flag_output_denormal STATUS_VAR
);
711 return packFloatx80(zSign
, 0, 0);
714 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
716 || ( zSig0
<= zSig0
+ roundIncrement
);
717 shift64RightJamming( zSig0
, 1 - zExp
, &zSig0
);
719 roundBits
= zSig0
& roundMask
;
720 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
721 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
722 zSig0
+= roundIncrement
;
723 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
724 roundIncrement
= roundMask
+ 1;
725 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
726 roundMask
|= roundIncrement
;
728 zSig0
&= ~ roundMask
;
729 return packFloatx80( zSign
, zExp
, zSig0
);
732 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
733 zSig0
+= roundIncrement
;
734 if ( zSig0
< roundIncrement
) {
736 zSig0
= LIT64( 0x8000000000000000 );
738 roundIncrement
= roundMask
+ 1;
739 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
740 roundMask
|= roundIncrement
;
742 zSig0
&= ~ roundMask
;
743 if ( zSig0
== 0 ) zExp
= 0;
744 return packFloatx80( zSign
, zExp
, zSig0
);
746 increment
= ( (int64_t) zSig1
< 0 );
747 if ( ! roundNearestEven
) {
748 if ( roundingMode
== float_round_to_zero
) {
753 increment
= ( roundingMode
== float_round_down
) && zSig1
;
756 increment
= ( roundingMode
== float_round_up
) && zSig1
;
760 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
761 if ( ( 0x7FFE < zExp
)
762 || ( ( zExp
== 0x7FFE )
763 && ( zSig0
== LIT64( 0xFFFFFFFFFFFFFFFF ) )
769 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
770 if ( ( roundingMode
== float_round_to_zero
)
771 || ( zSign
&& ( roundingMode
== float_round_up
) )
772 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
774 return packFloatx80( zSign
, 0x7FFE, ~ roundMask
);
776 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
780 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
783 || ( zSig0
< LIT64( 0xFFFFFFFFFFFFFFFF ) );
784 shift64ExtraRightJamming( zSig0
, zSig1
, 1 - zExp
, &zSig0
, &zSig1
);
786 if ( isTiny
&& zSig1
) float_raise( float_flag_underflow STATUS_VAR
);
787 if ( zSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
788 if ( roundNearestEven
) {
789 increment
= ( (int64_t) zSig1
< 0 );
793 increment
= ( roundingMode
== float_round_down
) && zSig1
;
796 increment
= ( roundingMode
== float_round_up
) && zSig1
;
802 ~ ( ( (uint64_t) ( zSig1
<<1 ) == 0 ) & roundNearestEven
);
803 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
805 return packFloatx80( zSign
, zExp
, zSig0
);
808 if ( zSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
813 zSig0
= LIT64( 0x8000000000000000 );
816 zSig0
&= ~ ( ( (uint64_t) ( zSig1
<<1 ) == 0 ) & roundNearestEven
);
820 if ( zSig0
== 0 ) zExp
= 0;
822 return packFloatx80( zSign
, zExp
, zSig0
);
826 /*----------------------------------------------------------------------------
827 | Takes an abstract floating-point value having sign `zSign', exponent
828 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
829 | and returns the proper extended double-precision floating-point value
830 | corresponding to the abstract input. This routine is just like
831 | `roundAndPackFloatx80' except that the input significand does not have to be
833 *----------------------------------------------------------------------------*/
836 normalizeRoundAndPackFloatx80(
837 int8 roundingPrecision
, flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
847 shiftCount
= countLeadingZeros64( zSig0
);
848 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
851 roundAndPackFloatx80( roundingPrecision
, zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
855 /*----------------------------------------------------------------------------
856 | Returns the least-significant 64 fraction bits of the quadruple-precision
857 | floating-point value `a'.
858 *----------------------------------------------------------------------------*/
860 INLINE
uint64_t extractFloat128Frac1( float128 a
)
867 /*----------------------------------------------------------------------------
868 | Returns the most-significant 48 fraction bits of the quadruple-precision
869 | floating-point value `a'.
870 *----------------------------------------------------------------------------*/
872 INLINE
uint64_t extractFloat128Frac0( float128 a
)
875 return a
.high
& LIT64( 0x0000FFFFFFFFFFFF );
879 /*----------------------------------------------------------------------------
880 | Returns the exponent bits of the quadruple-precision floating-point value
882 *----------------------------------------------------------------------------*/
884 INLINE int32
extractFloat128Exp( float128 a
)
887 return ( a
.high
>>48 ) & 0x7FFF;
891 /*----------------------------------------------------------------------------
892 | Returns the sign bit of the quadruple-precision floating-point value `a'.
893 *----------------------------------------------------------------------------*/
895 INLINE flag
extractFloat128Sign( float128 a
)
902 /*----------------------------------------------------------------------------
903 | Normalizes the subnormal quadruple-precision floating-point value
904 | represented by the denormalized significand formed by the concatenation of
905 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
906 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
907 | significand are stored at the location pointed to by `zSig0Ptr', and the
908 | least significant 64 bits of the normalized significand are stored at the
909 | location pointed to by `zSig1Ptr'.
910 *----------------------------------------------------------------------------*/
913 normalizeFloat128Subnormal(
924 shiftCount
= countLeadingZeros64( aSig1
) - 15;
925 if ( shiftCount
< 0 ) {
926 *zSig0Ptr
= aSig1
>>( - shiftCount
);
927 *zSig1Ptr
= aSig1
<<( shiftCount
& 63 );
930 *zSig0Ptr
= aSig1
<<shiftCount
;
933 *zExpPtr
= - shiftCount
- 63;
936 shiftCount
= countLeadingZeros64( aSig0
) - 15;
937 shortShift128Left( aSig0
, aSig1
, shiftCount
, zSig0Ptr
, zSig1Ptr
);
938 *zExpPtr
= 1 - shiftCount
;
943 /*----------------------------------------------------------------------------
944 | Packs the sign `zSign', the exponent `zExp', and the significand formed
945 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
946 | floating-point value, returning the result. After being shifted into the
947 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
948 | added together to form the most significant 32 bits of the result. This
949 | means that any integer portion of `zSig0' will be added into the exponent.
950 | Since a properly normalized significand will have an integer portion equal
951 | to 1, the `zExp' input should be 1 less than the desired result exponent
952 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
954 *----------------------------------------------------------------------------*/
957 packFloat128( flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
)
962 z
.high
= ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<48 ) + zSig0
;
967 /*----------------------------------------------------------------------------
968 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
969 | and extended significand formed by the concatenation of `zSig0', `zSig1',
970 | and `zSig2', and returns the proper quadruple-precision floating-point value
971 | corresponding to the abstract input. Ordinarily, the abstract value is
972 | simply rounded and packed into the quadruple-precision format, with the
973 | inexact exception raised if the abstract input cannot be represented
974 | exactly. However, if the abstract value is too large, the overflow and
975 | inexact exceptions are raised and an infinity or maximal finite value is
976 | returned. If the abstract value is too small, the input value is rounded to
977 | a subnormal number, and the underflow and inexact exceptions are raised if
978 | the abstract input cannot be represented exactly as a subnormal quadruple-
979 | precision floating-point number.
980 | The input significand must be normalized or smaller. If the input
981 | significand is not normalized, `zExp' must be 0; in that case, the result
982 | returned is a subnormal number, and it must not require rounding. In the
983 | usual case that the input significand is normalized, `zExp' must be 1 less
984 | than the ``true'' floating-point exponent. The handling of underflow and
985 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
986 *----------------------------------------------------------------------------*/
989 roundAndPackFloat128(
990 flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
, uint64_t zSig2 STATUS_PARAM
)
993 flag roundNearestEven
, increment
, isTiny
;
995 roundingMode
= STATUS(float_rounding_mode
);
996 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
997 increment
= ( (int64_t) zSig2
< 0 );
998 if ( ! roundNearestEven
) {
999 if ( roundingMode
== float_round_to_zero
) {
1004 increment
= ( roundingMode
== float_round_down
) && zSig2
;
1007 increment
= ( roundingMode
== float_round_up
) && zSig2
;
1011 if ( 0x7FFD <= (uint32_t) zExp
) {
1012 if ( ( 0x7FFD < zExp
)
1013 || ( ( zExp
== 0x7FFD )
1015 LIT64( 0x0001FFFFFFFFFFFF ),
1016 LIT64( 0xFFFFFFFFFFFFFFFF ),
1023 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
1024 if ( ( roundingMode
== float_round_to_zero
)
1025 || ( zSign
&& ( roundingMode
== float_round_up
) )
1026 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
1032 LIT64( 0x0000FFFFFFFFFFFF ),
1033 LIT64( 0xFFFFFFFFFFFFFFFF )
1036 return packFloat128( zSign
, 0x7FFF, 0, 0 );
1039 if (STATUS(flush_to_zero
)) {
1040 float_raise(float_flag_output_denormal STATUS_VAR
);
1041 return packFloat128(zSign
, 0, 0, 0);
1044 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
1050 LIT64( 0x0001FFFFFFFFFFFF ),
1051 LIT64( 0xFFFFFFFFFFFFFFFF )
1053 shift128ExtraRightJamming(
1054 zSig0
, zSig1
, zSig2
, - zExp
, &zSig0
, &zSig1
, &zSig2
);
1056 if ( isTiny
&& zSig2
) float_raise( float_flag_underflow STATUS_VAR
);
1057 if ( roundNearestEven
) {
1058 increment
= ( (int64_t) zSig2
< 0 );
1062 increment
= ( roundingMode
== float_round_down
) && zSig2
;
1065 increment
= ( roundingMode
== float_round_up
) && zSig2
;
1070 if ( zSig2
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1072 add128( zSig0
, zSig1
, 0, 1, &zSig0
, &zSig1
);
1073 zSig1
&= ~ ( ( zSig2
+ zSig2
== 0 ) & roundNearestEven
);
1076 if ( ( zSig0
| zSig1
) == 0 ) zExp
= 0;
1078 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
1082 /*----------------------------------------------------------------------------
1083 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1084 | and significand formed by the concatenation of `zSig0' and `zSig1', and
1085 | returns the proper quadruple-precision floating-point value corresponding
1086 | to the abstract input. This routine is just like `roundAndPackFloat128'
1087 | except that the input significand has fewer bits and does not have to be
1088 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
1090 *----------------------------------------------------------------------------*/
1093 normalizeRoundAndPackFloat128(
1094 flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1 STATUS_PARAM
)
1104 shiftCount
= countLeadingZeros64( zSig0
) - 15;
1105 if ( 0 <= shiftCount
) {
1107 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
1110 shift128ExtraRightJamming(
1111 zSig0
, zSig1
, 0, - shiftCount
, &zSig0
, &zSig1
, &zSig2
);
1114 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
1118 /*----------------------------------------------------------------------------
1119 | Returns the result of converting the 32-bit two's complement integer `a'
1120 | to the single-precision floating-point format. The conversion is performed
1121 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1122 *----------------------------------------------------------------------------*/
1124 float32
int32_to_float32( int32 a STATUS_PARAM
)
1128 if ( a
== 0 ) return float32_zero
;
1129 if ( a
== (int32_t) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
1131 return normalizeRoundAndPackFloat32( zSign
, 0x9C, zSign
? - a
: a STATUS_VAR
);
1135 /*----------------------------------------------------------------------------
1136 | Returns the result of converting the 32-bit two's complement integer `a'
1137 | to the double-precision floating-point format. The conversion is performed
1138 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1139 *----------------------------------------------------------------------------*/
1141 float64
int32_to_float64( int32 a STATUS_PARAM
)
1148 if ( a
== 0 ) return float64_zero
;
1150 absA
= zSign
? - a
: a
;
1151 shiftCount
= countLeadingZeros32( absA
) + 21;
1153 return packFloat64( zSign
, 0x432 - shiftCount
, zSig
<<shiftCount
);
1157 /*----------------------------------------------------------------------------
1158 | Returns the result of converting the 32-bit two's complement integer `a'
1159 | to the extended double-precision floating-point format. The conversion
1160 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1162 *----------------------------------------------------------------------------*/
1164 floatx80
int32_to_floatx80( int32 a STATUS_PARAM
)
1171 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
1173 absA
= zSign
? - a
: a
;
1174 shiftCount
= countLeadingZeros32( absA
) + 32;
1176 return packFloatx80( zSign
, 0x403E - shiftCount
, zSig
<<shiftCount
);
1180 /*----------------------------------------------------------------------------
1181 | Returns the result of converting the 32-bit two's complement integer `a' to
1182 | the quadruple-precision floating-point format. The conversion is performed
1183 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1184 *----------------------------------------------------------------------------*/
1186 float128
int32_to_float128( int32 a STATUS_PARAM
)
1193 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
1195 absA
= zSign
? - a
: a
;
1196 shiftCount
= countLeadingZeros32( absA
) + 17;
1198 return packFloat128( zSign
, 0x402E - shiftCount
, zSig0
<<shiftCount
, 0 );
1202 /*----------------------------------------------------------------------------
1203 | Returns the result of converting the 64-bit two's complement integer `a'
1204 | to the single-precision floating-point format. The conversion is performed
1205 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1206 *----------------------------------------------------------------------------*/
1208 float32
int64_to_float32( int64 a STATUS_PARAM
)
1214 if ( a
== 0 ) return float32_zero
;
1216 absA
= zSign
? - a
: a
;
1217 shiftCount
= countLeadingZeros64( absA
) - 40;
1218 if ( 0 <= shiftCount
) {
1219 return packFloat32( zSign
, 0x95 - shiftCount
, absA
<<shiftCount
);
1223 if ( shiftCount
< 0 ) {
1224 shift64RightJamming( absA
, - shiftCount
, &absA
);
1227 absA
<<= shiftCount
;
1229 return roundAndPackFloat32( zSign
, 0x9C - shiftCount
, absA STATUS_VAR
);
1234 float32
uint64_to_float32( uint64 a STATUS_PARAM
)
1238 if ( a
== 0 ) return float32_zero
;
1239 shiftCount
= countLeadingZeros64( a
) - 40;
1240 if ( 0 <= shiftCount
) {
1241 return packFloat32( 1 > 0, 0x95 - shiftCount
, a
<<shiftCount
);
1245 if ( shiftCount
< 0 ) {
1246 shift64RightJamming( a
, - shiftCount
, &a
);
1251 return roundAndPackFloat32( 1 > 0, 0x9C - shiftCount
, a STATUS_VAR
);
1255 /*----------------------------------------------------------------------------
1256 | Returns the result of converting the 64-bit two's complement integer `a'
1257 | to the double-precision floating-point format. The conversion is performed
1258 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1259 *----------------------------------------------------------------------------*/
1261 float64
int64_to_float64( int64 a STATUS_PARAM
)
1265 if ( a
== 0 ) return float64_zero
;
1266 if ( a
== (int64_t) LIT64( 0x8000000000000000 ) ) {
1267 return packFloat64( 1, 0x43E, 0 );
1270 return normalizeRoundAndPackFloat64( zSign
, 0x43C, zSign
? - a
: a STATUS_VAR
);
1274 float64
uint64_to_float64( uint64 a STATUS_PARAM
)
1276 if ( a
== 0 ) return float64_zero
;
1277 return normalizeRoundAndPackFloat64( 0, 0x43C, a STATUS_VAR
);
1281 /*----------------------------------------------------------------------------
1282 | Returns the result of converting the 64-bit two's complement integer `a'
1283 | to the extended double-precision floating-point format. The conversion
1284 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1286 *----------------------------------------------------------------------------*/
1288 floatx80
int64_to_floatx80( int64 a STATUS_PARAM
)
1294 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
1296 absA
= zSign
? - a
: a
;
1297 shiftCount
= countLeadingZeros64( absA
);
1298 return packFloatx80( zSign
, 0x403E - shiftCount
, absA
<<shiftCount
);
1302 /*----------------------------------------------------------------------------
1303 | Returns the result of converting the 64-bit two's complement integer `a' to
1304 | the quadruple-precision floating-point format. The conversion is performed
1305 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1306 *----------------------------------------------------------------------------*/
1308 float128
int64_to_float128( int64 a STATUS_PARAM
)
1314 uint64_t zSig0
, zSig1
;
1316 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
1318 absA
= zSign
? - a
: a
;
1319 shiftCount
= countLeadingZeros64( absA
) + 49;
1320 zExp
= 0x406E - shiftCount
;
1321 if ( 64 <= shiftCount
) {
1330 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
1331 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
1335 /*----------------------------------------------------------------------------
1336 | Returns the result of converting the single-precision floating-point value
1337 | `a' to the 32-bit two's complement integer format. The conversion is
1338 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1339 | Arithmetic---which means in particular that the conversion is rounded
1340 | according to the current rounding mode. If `a' is a NaN, the largest
1341 | positive integer is returned. Otherwise, if the conversion overflows, the
1342 | largest integer with the same sign as `a' is returned.
1343 *----------------------------------------------------------------------------*/
1345 int32
float32_to_int32( float32 a STATUS_PARAM
)
1348 int16 aExp
, shiftCount
;
1352 a
= float32_squash_input_denormal(a STATUS_VAR
);
1353 aSig
= extractFloat32Frac( a
);
1354 aExp
= extractFloat32Exp( a
);
1355 aSign
= extractFloat32Sign( a
);
1356 if ( ( aExp
== 0xFF ) && aSig
) aSign
= 0;
1357 if ( aExp
) aSig
|= 0x00800000;
1358 shiftCount
= 0xAF - aExp
;
1361 if ( 0 < shiftCount
) shift64RightJamming( aSig64
, shiftCount
, &aSig64
);
1362 return roundAndPackInt32( aSign
, aSig64 STATUS_VAR
);
1366 /*----------------------------------------------------------------------------
1367 | Returns the result of converting the single-precision floating-point value
1368 | `a' to the 32-bit two's complement integer format. The conversion is
1369 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1370 | Arithmetic, except that the conversion is always rounded toward zero.
1371 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1372 | the conversion overflows, the largest integer with the same sign as `a' is
1374 *----------------------------------------------------------------------------*/
1376 int32
float32_to_int32_round_to_zero( float32 a STATUS_PARAM
)
1379 int16 aExp
, shiftCount
;
1382 a
= float32_squash_input_denormal(a STATUS_VAR
);
1384 aSig
= extractFloat32Frac( a
);
1385 aExp
= extractFloat32Exp( a
);
1386 aSign
= extractFloat32Sign( a
);
1387 shiftCount
= aExp
- 0x9E;
1388 if ( 0 <= shiftCount
) {
1389 if ( float32_val(a
) != 0xCF000000 ) {
1390 float_raise( float_flag_invalid STATUS_VAR
);
1391 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) return 0x7FFFFFFF;
1393 return (int32_t) 0x80000000;
1395 else if ( aExp
<= 0x7E ) {
1396 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1399 aSig
= ( aSig
| 0x00800000 )<<8;
1400 z
= aSig
>>( - shiftCount
);
1401 if ( (uint32_t) ( aSig
<<( shiftCount
& 31 ) ) ) {
1402 STATUS(float_exception_flags
) |= float_flag_inexact
;
1404 if ( aSign
) z
= - z
;
1409 /*----------------------------------------------------------------------------
1410 | Returns the result of converting the single-precision floating-point value
1411 | `a' to the 16-bit two's complement integer format. The conversion is
1412 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1413 | Arithmetic, except that the conversion is always rounded toward zero.
1414 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1415 | the conversion overflows, the largest integer with the same sign as `a' is
1417 *----------------------------------------------------------------------------*/
1419 int16
float32_to_int16_round_to_zero( float32 a STATUS_PARAM
)
1422 int16 aExp
, shiftCount
;
1426 aSig
= extractFloat32Frac( a
);
1427 aExp
= extractFloat32Exp( a
);
1428 aSign
= extractFloat32Sign( a
);
1429 shiftCount
= aExp
- 0x8E;
1430 if ( 0 <= shiftCount
) {
1431 if ( float32_val(a
) != 0xC7000000 ) {
1432 float_raise( float_flag_invalid STATUS_VAR
);
1433 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1437 return (int32_t) 0xffff8000;
1439 else if ( aExp
<= 0x7E ) {
1440 if ( aExp
| aSig
) {
1441 STATUS(float_exception_flags
) |= float_flag_inexact
;
1446 aSig
= ( aSig
| 0x00800000 )<<8;
1447 z
= aSig
>>( - shiftCount
);
1448 if ( (uint32_t) ( aSig
<<( shiftCount
& 31 ) ) ) {
1449 STATUS(float_exception_flags
) |= float_flag_inexact
;
1458 /*----------------------------------------------------------------------------
1459 | Returns the result of converting the single-precision floating-point value
1460 | `a' to the 64-bit two's complement integer format. The conversion is
1461 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1462 | Arithmetic---which means in particular that the conversion is rounded
1463 | according to the current rounding mode. If `a' is a NaN, the largest
1464 | positive integer is returned. Otherwise, if the conversion overflows, the
1465 | largest integer with the same sign as `a' is returned.
1466 *----------------------------------------------------------------------------*/
1468 int64
float32_to_int64( float32 a STATUS_PARAM
)
1471 int16 aExp
, shiftCount
;
1473 uint64_t aSig64
, aSigExtra
;
1474 a
= float32_squash_input_denormal(a STATUS_VAR
);
1476 aSig
= extractFloat32Frac( a
);
1477 aExp
= extractFloat32Exp( a
);
1478 aSign
= extractFloat32Sign( a
);
1479 shiftCount
= 0xBE - aExp
;
1480 if ( shiftCount
< 0 ) {
1481 float_raise( float_flag_invalid STATUS_VAR
);
1482 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1483 return LIT64( 0x7FFFFFFFFFFFFFFF );
1485 return (int64_t) LIT64( 0x8000000000000000 );
1487 if ( aExp
) aSig
|= 0x00800000;
1490 shift64ExtraRightJamming( aSig64
, 0, shiftCount
, &aSig64
, &aSigExtra
);
1491 return roundAndPackInt64( aSign
, aSig64
, aSigExtra STATUS_VAR
);
1495 /*----------------------------------------------------------------------------
1496 | Returns the result of converting the single-precision floating-point value
1497 | `a' to the 64-bit two's complement integer format. The conversion is
1498 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1499 | Arithmetic, except that the conversion is always rounded toward zero. If
1500 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
1501 | conversion overflows, the largest integer with the same sign as `a' is
1503 *----------------------------------------------------------------------------*/
1505 int64
float32_to_int64_round_to_zero( float32 a STATUS_PARAM
)
1508 int16 aExp
, shiftCount
;
1512 a
= float32_squash_input_denormal(a STATUS_VAR
);
1514 aSig
= extractFloat32Frac( a
);
1515 aExp
= extractFloat32Exp( a
);
1516 aSign
= extractFloat32Sign( a
);
1517 shiftCount
= aExp
- 0xBE;
1518 if ( 0 <= shiftCount
) {
1519 if ( float32_val(a
) != 0xDF000000 ) {
1520 float_raise( float_flag_invalid STATUS_VAR
);
1521 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1522 return LIT64( 0x7FFFFFFFFFFFFFFF );
1525 return (int64_t) LIT64( 0x8000000000000000 );
1527 else if ( aExp
<= 0x7E ) {
1528 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1531 aSig64
= aSig
| 0x00800000;
1533 z
= aSig64
>>( - shiftCount
);
1534 if ( (uint64_t) ( aSig64
<<( shiftCount
& 63 ) ) ) {
1535 STATUS(float_exception_flags
) |= float_flag_inexact
;
1537 if ( aSign
) z
= - z
;
1542 /*----------------------------------------------------------------------------
1543 | Returns the result of converting the single-precision floating-point value
1544 | `a' to the double-precision floating-point format. The conversion is
1545 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1547 *----------------------------------------------------------------------------*/
1549 float64
float32_to_float64( float32 a STATUS_PARAM
)
1554 a
= float32_squash_input_denormal(a STATUS_VAR
);
1556 aSig
= extractFloat32Frac( a
);
1557 aExp
= extractFloat32Exp( a
);
1558 aSign
= extractFloat32Sign( a
);
1559 if ( aExp
== 0xFF ) {
1560 if ( aSig
) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1561 return packFloat64( aSign
, 0x7FF, 0 );
1564 if ( aSig
== 0 ) return packFloat64( aSign
, 0, 0 );
1565 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1568 return packFloat64( aSign
, aExp
+ 0x380, ( (uint64_t) aSig
)<<29 );
1572 /*----------------------------------------------------------------------------
1573 | Returns the result of converting the single-precision floating-point value
1574 | `a' to the extended double-precision floating-point format. The conversion
1575 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1577 *----------------------------------------------------------------------------*/
1579 floatx80
float32_to_floatx80( float32 a STATUS_PARAM
)
1585 a
= float32_squash_input_denormal(a STATUS_VAR
);
1586 aSig
= extractFloat32Frac( a
);
1587 aExp
= extractFloat32Exp( a
);
1588 aSign
= extractFloat32Sign( a
);
1589 if ( aExp
== 0xFF ) {
1590 if ( aSig
) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1591 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
1594 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
1595 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1598 return packFloatx80( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<40 );
1602 /*----------------------------------------------------------------------------
1603 | Returns the result of converting the single-precision floating-point value
1604 | `a' to the double-precision floating-point format. The conversion is
1605 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1607 *----------------------------------------------------------------------------*/
1609 float128
float32_to_float128( float32 a STATUS_PARAM
)
1615 a
= float32_squash_input_denormal(a STATUS_VAR
);
1616 aSig
= extractFloat32Frac( a
);
1617 aExp
= extractFloat32Exp( a
);
1618 aSign
= extractFloat32Sign( a
);
1619 if ( aExp
== 0xFF ) {
1620 if ( aSig
) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1621 return packFloat128( aSign
, 0x7FFF, 0, 0 );
1624 if ( aSig
== 0 ) return packFloat128( aSign
, 0, 0, 0 );
1625 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1628 return packFloat128( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<25, 0 );
1632 /*----------------------------------------------------------------------------
1633 | Rounds the single-precision floating-point value `a' to an integer, and
1634 | returns the result as a single-precision floating-point value. The
1635 | operation is performed according to the IEC/IEEE Standard for Binary
1636 | Floating-Point Arithmetic.
1637 *----------------------------------------------------------------------------*/
1639 float32
float32_round_to_int( float32 a STATUS_PARAM
)
1643 uint32_t lastBitMask
, roundBitsMask
;
1646 a
= float32_squash_input_denormal(a STATUS_VAR
);
1648 aExp
= extractFloat32Exp( a
);
1649 if ( 0x96 <= aExp
) {
1650 if ( ( aExp
== 0xFF ) && extractFloat32Frac( a
) ) {
1651 return propagateFloat32NaN( a
, a STATUS_VAR
);
1655 if ( aExp
<= 0x7E ) {
1656 if ( (uint32_t) ( float32_val(a
)<<1 ) == 0 ) return a
;
1657 STATUS(float_exception_flags
) |= float_flag_inexact
;
1658 aSign
= extractFloat32Sign( a
);
1659 switch ( STATUS(float_rounding_mode
) ) {
1660 case float_round_nearest_even
:
1661 if ( ( aExp
== 0x7E ) && extractFloat32Frac( a
) ) {
1662 return packFloat32( aSign
, 0x7F, 0 );
1665 case float_round_down
:
1666 return make_float32(aSign
? 0xBF800000 : 0);
1667 case float_round_up
:
1668 return make_float32(aSign
? 0x80000000 : 0x3F800000);
1670 return packFloat32( aSign
, 0, 0 );
1673 lastBitMask
<<= 0x96 - aExp
;
1674 roundBitsMask
= lastBitMask
- 1;
1676 roundingMode
= STATUS(float_rounding_mode
);
1677 if ( roundingMode
== float_round_nearest_even
) {
1678 z
+= lastBitMask
>>1;
1679 if ( ( z
& roundBitsMask
) == 0 ) z
&= ~ lastBitMask
;
1681 else if ( roundingMode
!= float_round_to_zero
) {
1682 if ( extractFloat32Sign( make_float32(z
) ) ^ ( roundingMode
== float_round_up
) ) {
1686 z
&= ~ roundBitsMask
;
1687 if ( z
!= float32_val(a
) ) STATUS(float_exception_flags
) |= float_flag_inexact
;
1688 return make_float32(z
);
1692 /*----------------------------------------------------------------------------
1693 | Returns the result of adding the absolute values of the single-precision
1694 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
1695 | before being returned. `zSign' is ignored if the result is a NaN.
1696 | The addition is performed according to the IEC/IEEE Standard for Binary
1697 | Floating-Point Arithmetic.
1698 *----------------------------------------------------------------------------*/
1700 static float32
addFloat32Sigs( float32 a
, float32 b
, flag zSign STATUS_PARAM
)
1702 int16 aExp
, bExp
, zExp
;
1703 uint32_t aSig
, bSig
, zSig
;
1706 aSig
= extractFloat32Frac( a
);
1707 aExp
= extractFloat32Exp( a
);
1708 bSig
= extractFloat32Frac( b
);
1709 bExp
= extractFloat32Exp( b
);
1710 expDiff
= aExp
- bExp
;
1713 if ( 0 < expDiff
) {
1714 if ( aExp
== 0xFF ) {
1715 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1724 shift32RightJamming( bSig
, expDiff
, &bSig
);
1727 else if ( expDiff
< 0 ) {
1728 if ( bExp
== 0xFF ) {
1729 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1730 return packFloat32( zSign
, 0xFF, 0 );
1738 shift32RightJamming( aSig
, - expDiff
, &aSig
);
1742 if ( aExp
== 0xFF ) {
1743 if ( aSig
| bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1747 if (STATUS(flush_to_zero
)) {
1749 float_raise(float_flag_output_denormal STATUS_VAR
);
1751 return packFloat32(zSign
, 0, 0);
1753 return packFloat32( zSign
, 0, ( aSig
+ bSig
)>>6 );
1755 zSig
= 0x40000000 + aSig
+ bSig
;
1760 zSig
= ( aSig
+ bSig
)<<1;
1762 if ( (int32_t) zSig
< 0 ) {
1767 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1771 /*----------------------------------------------------------------------------
1772 | Returns the result of subtracting the absolute values of the single-
1773 | precision floating-point values `a' and `b'. If `zSign' is 1, the
1774 | difference is negated before being returned. `zSign' is ignored if the
1775 | result is a NaN. The subtraction is performed according to the IEC/IEEE
1776 | Standard for Binary Floating-Point Arithmetic.
1777 *----------------------------------------------------------------------------*/
1779 static float32
subFloat32Sigs( float32 a
, float32 b
, flag zSign STATUS_PARAM
)
1781 int16 aExp
, bExp
, zExp
;
1782 uint32_t aSig
, bSig
, zSig
;
1785 aSig
= extractFloat32Frac( a
);
1786 aExp
= extractFloat32Exp( a
);
1787 bSig
= extractFloat32Frac( b
);
1788 bExp
= extractFloat32Exp( b
);
1789 expDiff
= aExp
- bExp
;
1792 if ( 0 < expDiff
) goto aExpBigger
;
1793 if ( expDiff
< 0 ) goto bExpBigger
;
1794 if ( aExp
== 0xFF ) {
1795 if ( aSig
| bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1796 float_raise( float_flag_invalid STATUS_VAR
);
1797 return float32_default_nan
;
1803 if ( bSig
< aSig
) goto aBigger
;
1804 if ( aSig
< bSig
) goto bBigger
;
1805 return packFloat32( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
1807 if ( bExp
== 0xFF ) {
1808 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1809 return packFloat32( zSign
^ 1, 0xFF, 0 );
1817 shift32RightJamming( aSig
, - expDiff
, &aSig
);
1823 goto normalizeRoundAndPack
;
1825 if ( aExp
== 0xFF ) {
1826 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1835 shift32RightJamming( bSig
, expDiff
, &bSig
);
1840 normalizeRoundAndPack
:
1842 return normalizeRoundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1846 /*----------------------------------------------------------------------------
1847 | Returns the result of adding the single-precision floating-point values `a'
1848 | and `b'. The operation is performed according to the IEC/IEEE Standard for
1849 | Binary Floating-Point Arithmetic.
1850 *----------------------------------------------------------------------------*/
1852 float32
float32_add( float32 a
, float32 b STATUS_PARAM
)
1855 a
= float32_squash_input_denormal(a STATUS_VAR
);
1856 b
= float32_squash_input_denormal(b STATUS_VAR
);
1858 aSign
= extractFloat32Sign( a
);
1859 bSign
= extractFloat32Sign( b
);
1860 if ( aSign
== bSign
) {
1861 return addFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1864 return subFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1869 /*----------------------------------------------------------------------------
1870 | Returns the result of subtracting the single-precision floating-point values
1871 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1872 | for Binary Floating-Point Arithmetic.
1873 *----------------------------------------------------------------------------*/
1875 float32
float32_sub( float32 a
, float32 b STATUS_PARAM
)
1878 a
= float32_squash_input_denormal(a STATUS_VAR
);
1879 b
= float32_squash_input_denormal(b STATUS_VAR
);
1881 aSign
= extractFloat32Sign( a
);
1882 bSign
= extractFloat32Sign( b
);
1883 if ( aSign
== bSign
) {
1884 return subFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1887 return addFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1892 /*----------------------------------------------------------------------------
1893 | Returns the result of multiplying the single-precision floating-point values
1894 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1895 | for Binary Floating-Point Arithmetic.
1896 *----------------------------------------------------------------------------*/
1898 float32
float32_mul( float32 a
, float32 b STATUS_PARAM
)
1900 flag aSign
, bSign
, zSign
;
1901 int16 aExp
, bExp
, zExp
;
1902 uint32_t aSig
, bSig
;
1906 a
= float32_squash_input_denormal(a STATUS_VAR
);
1907 b
= float32_squash_input_denormal(b STATUS_VAR
);
1909 aSig
= extractFloat32Frac( a
);
1910 aExp
= extractFloat32Exp( a
);
1911 aSign
= extractFloat32Sign( a
);
1912 bSig
= extractFloat32Frac( b
);
1913 bExp
= extractFloat32Exp( b
);
1914 bSign
= extractFloat32Sign( b
);
1915 zSign
= aSign
^ bSign
;
1916 if ( aExp
== 0xFF ) {
1917 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
1918 return propagateFloat32NaN( a
, b STATUS_VAR
);
1920 if ( ( bExp
| bSig
) == 0 ) {
1921 float_raise( float_flag_invalid STATUS_VAR
);
1922 return float32_default_nan
;
1924 return packFloat32( zSign
, 0xFF, 0 );
1926 if ( bExp
== 0xFF ) {
1927 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1928 if ( ( aExp
| aSig
) == 0 ) {
1929 float_raise( float_flag_invalid STATUS_VAR
);
1930 return float32_default_nan
;
1932 return packFloat32( zSign
, 0xFF, 0 );
1935 if ( aSig
== 0 ) return packFloat32( zSign
, 0, 0 );
1936 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1939 if ( bSig
== 0 ) return packFloat32( zSign
, 0, 0 );
1940 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
1942 zExp
= aExp
+ bExp
- 0x7F;
1943 aSig
= ( aSig
| 0x00800000 )<<7;
1944 bSig
= ( bSig
| 0x00800000 )<<8;
1945 shift64RightJamming( ( (uint64_t) aSig
) * bSig
, 32, &zSig64
);
1947 if ( 0 <= (int32_t) ( zSig
<<1 ) ) {
1951 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1955 /*----------------------------------------------------------------------------
1956 | Returns the result of dividing the single-precision floating-point value `a'
1957 | by the corresponding value `b'. The operation is performed according to the
1958 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1959 *----------------------------------------------------------------------------*/
1961 float32
float32_div( float32 a
, float32 b STATUS_PARAM
)
1963 flag aSign
, bSign
, zSign
;
1964 int16 aExp
, bExp
, zExp
;
1965 uint32_t aSig
, bSig
, zSig
;
1966 a
= float32_squash_input_denormal(a STATUS_VAR
);
1967 b
= float32_squash_input_denormal(b STATUS_VAR
);
1969 aSig
= extractFloat32Frac( a
);
1970 aExp
= extractFloat32Exp( a
);
1971 aSign
= extractFloat32Sign( a
);
1972 bSig
= extractFloat32Frac( b
);
1973 bExp
= extractFloat32Exp( b
);
1974 bSign
= extractFloat32Sign( b
);
1975 zSign
= aSign
^ bSign
;
1976 if ( aExp
== 0xFF ) {
1977 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1978 if ( bExp
== 0xFF ) {
1979 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1980 float_raise( float_flag_invalid STATUS_VAR
);
1981 return float32_default_nan
;
1983 return packFloat32( zSign
, 0xFF, 0 );
1985 if ( bExp
== 0xFF ) {
1986 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1987 return packFloat32( zSign
, 0, 0 );
1991 if ( ( aExp
| aSig
) == 0 ) {
1992 float_raise( float_flag_invalid STATUS_VAR
);
1993 return float32_default_nan
;
1995 float_raise( float_flag_divbyzero STATUS_VAR
);
1996 return packFloat32( zSign
, 0xFF, 0 );
1998 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
2001 if ( aSig
== 0 ) return packFloat32( zSign
, 0, 0 );
2002 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2004 zExp
= aExp
- bExp
+ 0x7D;
2005 aSig
= ( aSig
| 0x00800000 )<<7;
2006 bSig
= ( bSig
| 0x00800000 )<<8;
2007 if ( bSig
<= ( aSig
+ aSig
) ) {
2011 zSig
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
2012 if ( ( zSig
& 0x3F ) == 0 ) {
2013 zSig
|= ( (uint64_t) bSig
* zSig
!= ( (uint64_t) aSig
)<<32 );
2015 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
2019 /*----------------------------------------------------------------------------
2020 | Returns the remainder of the single-precision floating-point value `a'
2021 | with respect to the corresponding value `b'. The operation is performed
2022 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2023 *----------------------------------------------------------------------------*/
2025 float32
float32_rem( float32 a
, float32 b STATUS_PARAM
)
2028 int16 aExp
, bExp
, expDiff
;
2029 uint32_t aSig
, bSig
;
2031 uint64_t aSig64
, bSig64
, q64
;
2032 uint32_t alternateASig
;
2034 a
= float32_squash_input_denormal(a STATUS_VAR
);
2035 b
= float32_squash_input_denormal(b STATUS_VAR
);
2037 aSig
= extractFloat32Frac( a
);
2038 aExp
= extractFloat32Exp( a
);
2039 aSign
= extractFloat32Sign( a
);
2040 bSig
= extractFloat32Frac( b
);
2041 bExp
= extractFloat32Exp( b
);
2042 if ( aExp
== 0xFF ) {
2043 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
2044 return propagateFloat32NaN( a
, b STATUS_VAR
);
2046 float_raise( float_flag_invalid STATUS_VAR
);
2047 return float32_default_nan
;
2049 if ( bExp
== 0xFF ) {
2050 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
2055 float_raise( float_flag_invalid STATUS_VAR
);
2056 return float32_default_nan
;
2058 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
2061 if ( aSig
== 0 ) return a
;
2062 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2064 expDiff
= aExp
- bExp
;
2067 if ( expDiff
< 32 ) {
2070 if ( expDiff
< 0 ) {
2071 if ( expDiff
< -1 ) return a
;
2074 q
= ( bSig
<= aSig
);
2075 if ( q
) aSig
-= bSig
;
2076 if ( 0 < expDiff
) {
2077 q
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
2080 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
2088 if ( bSig
<= aSig
) aSig
-= bSig
;
2089 aSig64
= ( (uint64_t) aSig
)<<40;
2090 bSig64
= ( (uint64_t) bSig
)<<40;
2092 while ( 0 < expDiff
) {
2093 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
2094 q64
= ( 2 < q64
) ? q64
- 2 : 0;
2095 aSig64
= - ( ( bSig
* q64
)<<38 );
2099 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
2100 q64
= ( 2 < q64
) ? q64
- 2 : 0;
2101 q
= q64
>>( 64 - expDiff
);
2103 aSig
= ( ( aSig64
>>33 )<<( expDiff
- 1 ) ) - bSig
* q
;
2106 alternateASig
= aSig
;
2109 } while ( 0 <= (int32_t) aSig
);
2110 sigMean
= aSig
+ alternateASig
;
2111 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
2112 aSig
= alternateASig
;
2114 zSign
= ( (int32_t) aSig
< 0 );
2115 if ( zSign
) aSig
= - aSig
;
2116 return normalizeRoundAndPackFloat32( aSign
^ zSign
, bExp
, aSig STATUS_VAR
);
2120 /*----------------------------------------------------------------------------
2121 | Returns the result of multiplying the single-precision floating-point values
2122 | `a' and `b' then adding 'c', with no intermediate rounding step after the
2123 | multiplication. The operation is performed according to the IEC/IEEE
2124 | Standard for Binary Floating-Point Arithmetic 754-2008.
2125 | The flags argument allows the caller to select negation of the
2126 | addend, the intermediate product, or the final result. (The difference
2127 | between this and having the caller do a separate negation is that negating
2128 | externally will flip the sign bit on NaNs.)
2129 *----------------------------------------------------------------------------*/
2131 float32
float32_muladd(float32 a
, float32 b
, float32 c
, int flags STATUS_PARAM
)
2133 flag aSign
, bSign
, cSign
, zSign
;
2134 int aExp
, bExp
, cExp
, pExp
, zExp
, expDiff
;
2135 uint32_t aSig
, bSig
, cSig
;
2136 flag pInf
, pZero
, pSign
;
2137 uint64_t pSig64
, cSig64
, zSig64
;
2140 flag signflip
, infzero
;
2142 a
= float32_squash_input_denormal(a STATUS_VAR
);
2143 b
= float32_squash_input_denormal(b STATUS_VAR
);
2144 c
= float32_squash_input_denormal(c STATUS_VAR
);
2145 aSig
= extractFloat32Frac(a
);
2146 aExp
= extractFloat32Exp(a
);
2147 aSign
= extractFloat32Sign(a
);
2148 bSig
= extractFloat32Frac(b
);
2149 bExp
= extractFloat32Exp(b
);
2150 bSign
= extractFloat32Sign(b
);
2151 cSig
= extractFloat32Frac(c
);
2152 cExp
= extractFloat32Exp(c
);
2153 cSign
= extractFloat32Sign(c
);
2155 infzero
= ((aExp
== 0 && aSig
== 0 && bExp
== 0xff && bSig
== 0) ||
2156 (aExp
== 0xff && aSig
== 0 && bExp
== 0 && bSig
== 0));
2158 /* It is implementation-defined whether the cases of (0,inf,qnan)
2159 * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
2160 * they return if they do), so we have to hand this information
2161 * off to the target-specific pick-a-NaN routine.
2163 if (((aExp
== 0xff) && aSig
) ||
2164 ((bExp
== 0xff) && bSig
) ||
2165 ((cExp
== 0xff) && cSig
)) {
2166 return propagateFloat32MulAddNaN(a
, b
, c
, infzero STATUS_VAR
);
2170 float_raise(float_flag_invalid STATUS_VAR
);
2171 return float32_default_nan
;
2174 if (flags
& float_muladd_negate_c
) {
2178 signflip
= (flags
& float_muladd_negate_result
) ? 1 : 0;
2180 /* Work out the sign and type of the product */
2181 pSign
= aSign
^ bSign
;
2182 if (flags
& float_muladd_negate_product
) {
2185 pInf
= (aExp
== 0xff) || (bExp
== 0xff);
2186 pZero
= ((aExp
| aSig
) == 0) || ((bExp
| bSig
) == 0);
2189 if (pInf
&& (pSign
^ cSign
)) {
2190 /* addition of opposite-signed infinities => InvalidOperation */
2191 float_raise(float_flag_invalid STATUS_VAR
);
2192 return float32_default_nan
;
2194 /* Otherwise generate an infinity of the same sign */
2195 return packFloat32(cSign
^ signflip
, 0xff, 0);
2199 return packFloat32(pSign
^ signflip
, 0xff, 0);
2205 /* Adding two exact zeroes */
2206 if (pSign
== cSign
) {
2208 } else if (STATUS(float_rounding_mode
) == float_round_down
) {
2213 return packFloat32(zSign
^ signflip
, 0, 0);
2215 /* Exact zero plus a denorm */
2216 if (STATUS(flush_to_zero
)) {
2217 float_raise(float_flag_output_denormal STATUS_VAR
);
2218 return packFloat32(cSign
^ signflip
, 0, 0);
2221 /* Zero plus something non-zero : just return the something */
2222 return make_float32(float32_val(c
) ^ (signflip
<< 31));
2226 normalizeFloat32Subnormal(aSig
, &aExp
, &aSig
);
2229 normalizeFloat32Subnormal(bSig
, &bExp
, &bSig
);
2232 /* Calculate the actual result a * b + c */
2234 /* Multiply first; this is easy. */
2235 /* NB: we subtract 0x7e where float32_mul() subtracts 0x7f
2236 * because we want the true exponent, not the "one-less-than"
2237 * flavour that roundAndPackFloat32() takes.
2239 pExp
= aExp
+ bExp
- 0x7e;
2240 aSig
= (aSig
| 0x00800000) << 7;
2241 bSig
= (bSig
| 0x00800000) << 8;
2242 pSig64
= (uint64_t)aSig
* bSig
;
2243 if ((int64_t)(pSig64
<< 1) >= 0) {
2248 zSign
= pSign
^ signflip
;
2250 /* Now pSig64 is the significand of the multiply, with the explicit bit in
2255 /* Throw out the special case of c being an exact zero now */
2256 shift64RightJamming(pSig64
, 32, &pSig64
);
2258 return roundAndPackFloat32(zSign
, pExp
- 1,
2261 normalizeFloat32Subnormal(cSig
, &cExp
, &cSig
);
2264 cSig64
= (uint64_t)cSig
<< (62 - 23);
2265 cSig64
|= LIT64(0x4000000000000000);
2266 expDiff
= pExp
- cExp
;
2268 if (pSign
== cSign
) {
2271 /* scale c to match p */
2272 shift64RightJamming(cSig64
, expDiff
, &cSig64
);
2274 } else if (expDiff
< 0) {
2275 /* scale p to match c */
2276 shift64RightJamming(pSig64
, -expDiff
, &pSig64
);
2279 /* no scaling needed */
2282 /* Add significands and make sure explicit bit ends up in posn 62 */
2283 zSig64
= pSig64
+ cSig64
;
2284 if ((int64_t)zSig64
< 0) {
2285 shift64RightJamming(zSig64
, 1, &zSig64
);
2292 shift64RightJamming(cSig64
, expDiff
, &cSig64
);
2293 zSig64
= pSig64
- cSig64
;
2295 } else if (expDiff
< 0) {
2296 shift64RightJamming(pSig64
, -expDiff
, &pSig64
);
2297 zSig64
= cSig64
- pSig64
;
2302 if (cSig64
< pSig64
) {
2303 zSig64
= pSig64
- cSig64
;
2304 } else if (pSig64
< cSig64
) {
2305 zSig64
= cSig64
- pSig64
;
2310 if (STATUS(float_rounding_mode
) == float_round_down
) {
2313 return packFloat32(zSign
, 0, 0);
2317 /* Normalize to put the explicit bit back into bit 62. */
2318 shiftcount
= countLeadingZeros64(zSig64
) - 1;
2319 zSig64
<<= shiftcount
;
2322 shift64RightJamming(zSig64
, 32, &zSig64
);
2323 return roundAndPackFloat32(zSign
, zExp
, zSig64 STATUS_VAR
);
2327 /*----------------------------------------------------------------------------
2328 | Returns the square root of the single-precision floating-point value `a'.
2329 | The operation is performed according to the IEC/IEEE Standard for Binary
2330 | Floating-Point Arithmetic.
2331 *----------------------------------------------------------------------------*/
2333 float32
float32_sqrt( float32 a STATUS_PARAM
)
2337 uint32_t aSig
, zSig
;
2339 a
= float32_squash_input_denormal(a STATUS_VAR
);
2341 aSig
= extractFloat32Frac( a
);
2342 aExp
= extractFloat32Exp( a
);
2343 aSign
= extractFloat32Sign( a
);
2344 if ( aExp
== 0xFF ) {
2345 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2346 if ( ! aSign
) return a
;
2347 float_raise( float_flag_invalid STATUS_VAR
);
2348 return float32_default_nan
;
2351 if ( ( aExp
| aSig
) == 0 ) return a
;
2352 float_raise( float_flag_invalid STATUS_VAR
);
2353 return float32_default_nan
;
2356 if ( aSig
== 0 ) return float32_zero
;
2357 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2359 zExp
= ( ( aExp
- 0x7F )>>1 ) + 0x7E;
2360 aSig
= ( aSig
| 0x00800000 )<<8;
2361 zSig
= estimateSqrt32( aExp
, aSig
) + 2;
2362 if ( ( zSig
& 0x7F ) <= 5 ) {
2368 term
= ( (uint64_t) zSig
) * zSig
;
2369 rem
= ( ( (uint64_t) aSig
)<<32 ) - term
;
2370 while ( (int64_t) rem
< 0 ) {
2372 rem
+= ( ( (uint64_t) zSig
)<<1 ) | 1;
2374 zSig
|= ( rem
!= 0 );
2376 shift32RightJamming( zSig
, 1, &zSig
);
2378 return roundAndPackFloat32( 0, zExp
, zSig STATUS_VAR
);
2382 /*----------------------------------------------------------------------------
2383 | Returns the binary exponential of the single-precision floating-point value
2384 | `a'. The operation is performed according to the IEC/IEEE Standard for
2385 | Binary Floating-Point Arithmetic.
2387 | Uses the following identities:
2389 | 1. -------------------------------------------------------------------------
2393 | 2. -------------------------------------------------------------------------
2396 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
2398 *----------------------------------------------------------------------------*/
2400 static const float64 float32_exp2_coefficients
[15] =
2402 const_float64( 0x3ff0000000000000ll
), /* 1 */
2403 const_float64( 0x3fe0000000000000ll
), /* 2 */
2404 const_float64( 0x3fc5555555555555ll
), /* 3 */
2405 const_float64( 0x3fa5555555555555ll
), /* 4 */
2406 const_float64( 0x3f81111111111111ll
), /* 5 */
2407 const_float64( 0x3f56c16c16c16c17ll
), /* 6 */
2408 const_float64( 0x3f2a01a01a01a01all
), /* 7 */
2409 const_float64( 0x3efa01a01a01a01all
), /* 8 */
2410 const_float64( 0x3ec71de3a556c734ll
), /* 9 */
2411 const_float64( 0x3e927e4fb7789f5cll
), /* 10 */
2412 const_float64( 0x3e5ae64567f544e4ll
), /* 11 */
2413 const_float64( 0x3e21eed8eff8d898ll
), /* 12 */
2414 const_float64( 0x3de6124613a86d09ll
), /* 13 */
2415 const_float64( 0x3da93974a8c07c9dll
), /* 14 */
2416 const_float64( 0x3d6ae7f3e733b81fll
), /* 15 */
2419 float32
float32_exp2( float32 a STATUS_PARAM
)
2426 a
= float32_squash_input_denormal(a STATUS_VAR
);
2428 aSig
= extractFloat32Frac( a
);
2429 aExp
= extractFloat32Exp( a
);
2430 aSign
= extractFloat32Sign( a
);
2432 if ( aExp
== 0xFF) {
2433 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2434 return (aSign
) ? float32_zero
: a
;
2437 if (aSig
== 0) return float32_one
;
2440 float_raise( float_flag_inexact STATUS_VAR
);
2442 /* ******************************* */
2443 /* using float64 for approximation */
2444 /* ******************************* */
2445 x
= float32_to_float64(a STATUS_VAR
);
2446 x
= float64_mul(x
, float64_ln2 STATUS_VAR
);
2450 for (i
= 0 ; i
< 15 ; i
++) {
2453 f
= float64_mul(xn
, float32_exp2_coefficients
[i
] STATUS_VAR
);
2454 r
= float64_add(r
, f STATUS_VAR
);
2456 xn
= float64_mul(xn
, x STATUS_VAR
);
2459 return float64_to_float32(r
, status
);
2462 /*----------------------------------------------------------------------------
2463 | Returns the binary log of the single-precision floating-point value `a'.
2464 | The operation is performed according to the IEC/IEEE Standard for Binary
2465 | Floating-Point Arithmetic.
2466 *----------------------------------------------------------------------------*/
2467 float32
float32_log2( float32 a STATUS_PARAM
)
2471 uint32_t aSig
, zSig
, i
;
2473 a
= float32_squash_input_denormal(a STATUS_VAR
);
2474 aSig
= extractFloat32Frac( a
);
2475 aExp
= extractFloat32Exp( a
);
2476 aSign
= extractFloat32Sign( a
);
2479 if ( aSig
== 0 ) return packFloat32( 1, 0xFF, 0 );
2480 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2483 float_raise( float_flag_invalid STATUS_VAR
);
2484 return float32_default_nan
;
2486 if ( aExp
== 0xFF ) {
2487 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2496 for (i
= 1 << 22; i
> 0; i
>>= 1) {
2497 aSig
= ( (uint64_t)aSig
* aSig
) >> 23;
2498 if ( aSig
& 0x01000000 ) {
2507 return normalizeRoundAndPackFloat32( zSign
, 0x85, zSig STATUS_VAR
);
2510 /*----------------------------------------------------------------------------
2511 | Returns 1 if the single-precision floating-point value `a' is equal to
2512 | the corresponding value `b', and 0 otherwise. The invalid exception is
2513 | raised if either operand is a NaN. Otherwise, the comparison is performed
2514 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2515 *----------------------------------------------------------------------------*/
2517 int float32_eq( float32 a
, float32 b STATUS_PARAM
)
2520 a
= float32_squash_input_denormal(a STATUS_VAR
);
2521 b
= float32_squash_input_denormal(b STATUS_VAR
);
2523 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2524 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2526 float_raise( float_flag_invalid STATUS_VAR
);
2529 av
= float32_val(a
);
2530 bv
= float32_val(b
);
2531 return ( av
== bv
) || ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2534 /*----------------------------------------------------------------------------
2535 | Returns 1 if the single-precision floating-point value `a' is less than
2536 | or equal to the corresponding value `b', and 0 otherwise. The invalid
2537 | exception is raised if either operand is a NaN. The comparison is performed
2538 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2539 *----------------------------------------------------------------------------*/
2541 int float32_le( float32 a
, float32 b STATUS_PARAM
)
2545 a
= float32_squash_input_denormal(a STATUS_VAR
);
2546 b
= float32_squash_input_denormal(b STATUS_VAR
);
2548 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2549 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2551 float_raise( float_flag_invalid STATUS_VAR
);
2554 aSign
= extractFloat32Sign( a
);
2555 bSign
= extractFloat32Sign( b
);
2556 av
= float32_val(a
);
2557 bv
= float32_val(b
);
2558 if ( aSign
!= bSign
) return aSign
|| ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2559 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
2563 /*----------------------------------------------------------------------------
2564 | Returns 1 if the single-precision floating-point value `a' is less than
2565 | the corresponding value `b', and 0 otherwise. The invalid exception is
2566 | raised if either operand is a NaN. The comparison is performed according
2567 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2568 *----------------------------------------------------------------------------*/
2570 int float32_lt( float32 a
, float32 b STATUS_PARAM
)
2574 a
= float32_squash_input_denormal(a STATUS_VAR
);
2575 b
= float32_squash_input_denormal(b STATUS_VAR
);
2577 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2578 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2580 float_raise( float_flag_invalid STATUS_VAR
);
2583 aSign
= extractFloat32Sign( a
);
2584 bSign
= extractFloat32Sign( b
);
2585 av
= float32_val(a
);
2586 bv
= float32_val(b
);
2587 if ( aSign
!= bSign
) return aSign
&& ( (uint32_t) ( ( av
| bv
)<<1 ) != 0 );
2588 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
2592 /*----------------------------------------------------------------------------
2593 | Returns 1 if the single-precision floating-point values `a' and `b' cannot
2594 | be compared, and 0 otherwise. The invalid exception is raised if either
2595 | operand is a NaN. The comparison is performed according to the IEC/IEEE
2596 | Standard for Binary Floating-Point Arithmetic.
2597 *----------------------------------------------------------------------------*/
2599 int float32_unordered( float32 a
, float32 b STATUS_PARAM
)
2601 a
= float32_squash_input_denormal(a STATUS_VAR
);
2602 b
= float32_squash_input_denormal(b STATUS_VAR
);
2604 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2605 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2607 float_raise( float_flag_invalid STATUS_VAR
);
2613 /*----------------------------------------------------------------------------
2614 | Returns 1 if the single-precision floating-point value `a' is equal to
2615 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2616 | exception. The comparison is performed according to the IEC/IEEE Standard
2617 | for Binary Floating-Point Arithmetic.
2618 *----------------------------------------------------------------------------*/
2620 int float32_eq_quiet( float32 a
, float32 b STATUS_PARAM
)
2622 a
= float32_squash_input_denormal(a STATUS_VAR
);
2623 b
= float32_squash_input_denormal(b STATUS_VAR
);
2625 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2626 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2628 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2629 float_raise( float_flag_invalid STATUS_VAR
);
2633 return ( float32_val(a
) == float32_val(b
) ) ||
2634 ( (uint32_t) ( ( float32_val(a
) | float32_val(b
) )<<1 ) == 0 );
2637 /*----------------------------------------------------------------------------
2638 | Returns 1 if the single-precision floating-point value `a' is less than or
2639 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
2640 | cause an exception. Otherwise, the comparison is performed according to the
2641 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2642 *----------------------------------------------------------------------------*/
2644 int float32_le_quiet( float32 a
, float32 b STATUS_PARAM
)
2648 a
= float32_squash_input_denormal(a STATUS_VAR
);
2649 b
= float32_squash_input_denormal(b STATUS_VAR
);
2651 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2652 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2654 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2655 float_raise( float_flag_invalid STATUS_VAR
);
2659 aSign
= extractFloat32Sign( a
);
2660 bSign
= extractFloat32Sign( b
);
2661 av
= float32_val(a
);
2662 bv
= float32_val(b
);
2663 if ( aSign
!= bSign
) return aSign
|| ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2664 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
2668 /*----------------------------------------------------------------------------
2669 | Returns 1 if the single-precision floating-point value `a' is less than
2670 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2671 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
2672 | Standard for Binary Floating-Point Arithmetic.
2673 *----------------------------------------------------------------------------*/
2675 int float32_lt_quiet( float32 a
, float32 b STATUS_PARAM
)
2679 a
= float32_squash_input_denormal(a STATUS_VAR
);
2680 b
= float32_squash_input_denormal(b STATUS_VAR
);
2682 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2683 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2685 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2686 float_raise( float_flag_invalid STATUS_VAR
);
2690 aSign
= extractFloat32Sign( a
);
2691 bSign
= extractFloat32Sign( b
);
2692 av
= float32_val(a
);
2693 bv
= float32_val(b
);
2694 if ( aSign
!= bSign
) return aSign
&& ( (uint32_t) ( ( av
| bv
)<<1 ) != 0 );
2695 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
2699 /*----------------------------------------------------------------------------
2700 | Returns 1 if the single-precision floating-point values `a' and `b' cannot
2701 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
2702 | comparison is performed according to the IEC/IEEE Standard for Binary
2703 | Floating-Point Arithmetic.
2704 *----------------------------------------------------------------------------*/
2706 int float32_unordered_quiet( float32 a
, float32 b STATUS_PARAM
)
2708 a
= float32_squash_input_denormal(a STATUS_VAR
);
2709 b
= float32_squash_input_denormal(b STATUS_VAR
);
2711 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2712 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2714 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2715 float_raise( float_flag_invalid STATUS_VAR
);
2722 /*----------------------------------------------------------------------------
2723 | Returns the result of converting the double-precision floating-point value
2724 | `a' to the 32-bit two's complement integer format. The conversion is
2725 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2726 | Arithmetic---which means in particular that the conversion is rounded
2727 | according to the current rounding mode. If `a' is a NaN, the largest
2728 | positive integer is returned. Otherwise, if the conversion overflows, the
2729 | largest integer with the same sign as `a' is returned.
2730 *----------------------------------------------------------------------------*/
2732 int32
float64_to_int32( float64 a STATUS_PARAM
)
2735 int16 aExp
, shiftCount
;
2737 a
= float64_squash_input_denormal(a STATUS_VAR
);
2739 aSig
= extractFloat64Frac( a
);
2740 aExp
= extractFloat64Exp( a
);
2741 aSign
= extractFloat64Sign( a
);
2742 if ( ( aExp
== 0x7FF ) && aSig
) aSign
= 0;
2743 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2744 shiftCount
= 0x42C - aExp
;
2745 if ( 0 < shiftCount
) shift64RightJamming( aSig
, shiftCount
, &aSig
);
2746 return roundAndPackInt32( aSign
, aSig STATUS_VAR
);
2750 /*----------------------------------------------------------------------------
2751 | Returns the result of converting the double-precision floating-point value
2752 | `a' to the 32-bit two's complement integer format. The conversion is
2753 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2754 | Arithmetic, except that the conversion is always rounded toward zero.
2755 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2756 | the conversion overflows, the largest integer with the same sign as `a' is
2758 *----------------------------------------------------------------------------*/
2760 int32
float64_to_int32_round_to_zero( float64 a STATUS_PARAM
)
2763 int16 aExp
, shiftCount
;
2764 uint64_t aSig
, savedASig
;
2766 a
= float64_squash_input_denormal(a STATUS_VAR
);
2768 aSig
= extractFloat64Frac( a
);
2769 aExp
= extractFloat64Exp( a
);
2770 aSign
= extractFloat64Sign( a
);
2771 if ( 0x41E < aExp
) {
2772 if ( ( aExp
== 0x7FF ) && aSig
) aSign
= 0;
2775 else if ( aExp
< 0x3FF ) {
2776 if ( aExp
|| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
2779 aSig
|= LIT64( 0x0010000000000000 );
2780 shiftCount
= 0x433 - aExp
;
2782 aSig
>>= shiftCount
;
2784 if ( aSign
) z
= - z
;
2785 if ( ( z
< 0 ) ^ aSign
) {
2787 float_raise( float_flag_invalid STATUS_VAR
);
2788 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
2790 if ( ( aSig
<<shiftCount
) != savedASig
) {
2791 STATUS(float_exception_flags
) |= float_flag_inexact
;
2797 /*----------------------------------------------------------------------------
2798 | Returns the result of converting the double-precision floating-point value
2799 | `a' to the 16-bit two's complement integer format. The conversion is
2800 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2801 | Arithmetic, except that the conversion is always rounded toward zero.
2802 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2803 | the conversion overflows, the largest integer with the same sign as `a' is
2805 *----------------------------------------------------------------------------*/
2807 int16
float64_to_int16_round_to_zero( float64 a STATUS_PARAM
)
2810 int16 aExp
, shiftCount
;
2811 uint64_t aSig
, savedASig
;
2814 aSig
= extractFloat64Frac( a
);
2815 aExp
= extractFloat64Exp( a
);
2816 aSign
= extractFloat64Sign( a
);
2817 if ( 0x40E < aExp
) {
2818 if ( ( aExp
== 0x7FF ) && aSig
) {
2823 else if ( aExp
< 0x3FF ) {
2824 if ( aExp
|| aSig
) {
2825 STATUS(float_exception_flags
) |= float_flag_inexact
;
2829 aSig
|= LIT64( 0x0010000000000000 );
2830 shiftCount
= 0x433 - aExp
;
2832 aSig
>>= shiftCount
;
2837 if ( ( (int16_t)z
< 0 ) ^ aSign
) {
2839 float_raise( float_flag_invalid STATUS_VAR
);
2840 return aSign
? (int32_t) 0xffff8000 : 0x7FFF;
2842 if ( ( aSig
<<shiftCount
) != savedASig
) {
2843 STATUS(float_exception_flags
) |= float_flag_inexact
;
2848 /*----------------------------------------------------------------------------
2849 | Returns the result of converting the double-precision floating-point value
2850 | `a' to the 64-bit two's complement integer format. The conversion is
2851 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2852 | Arithmetic---which means in particular that the conversion is rounded
2853 | according to the current rounding mode. If `a' is a NaN, the largest
2854 | positive integer is returned. Otherwise, if the conversion overflows, the
2855 | largest integer with the same sign as `a' is returned.
2856 *----------------------------------------------------------------------------*/
2858 int64
float64_to_int64( float64 a STATUS_PARAM
)
2861 int16 aExp
, shiftCount
;
2862 uint64_t aSig
, aSigExtra
;
2863 a
= float64_squash_input_denormal(a STATUS_VAR
);
2865 aSig
= extractFloat64Frac( a
);
2866 aExp
= extractFloat64Exp( a
);
2867 aSign
= extractFloat64Sign( a
);
2868 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2869 shiftCount
= 0x433 - aExp
;
2870 if ( shiftCount
<= 0 ) {
2871 if ( 0x43E < aExp
) {
2872 float_raise( float_flag_invalid STATUS_VAR
);
2874 || ( ( aExp
== 0x7FF )
2875 && ( aSig
!= LIT64( 0x0010000000000000 ) ) )
2877 return LIT64( 0x7FFFFFFFFFFFFFFF );
2879 return (int64_t) LIT64( 0x8000000000000000 );
2882 aSig
<<= - shiftCount
;
2885 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
2887 return roundAndPackInt64( aSign
, aSig
, aSigExtra STATUS_VAR
);
2891 /*----------------------------------------------------------------------------
2892 | Returns the result of converting the double-precision floating-point value
2893 | `a' to the 64-bit two's complement integer format. The conversion is
2894 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2895 | Arithmetic, except that the conversion is always rounded toward zero.
2896 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2897 | the conversion overflows, the largest integer with the same sign as `a' is
2899 *----------------------------------------------------------------------------*/
2901 int64
float64_to_int64_round_to_zero( float64 a STATUS_PARAM
)
2904 int16 aExp
, shiftCount
;
2907 a
= float64_squash_input_denormal(a STATUS_VAR
);
2909 aSig
= extractFloat64Frac( a
);
2910 aExp
= extractFloat64Exp( a
);
2911 aSign
= extractFloat64Sign( a
);
2912 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2913 shiftCount
= aExp
- 0x433;
2914 if ( 0 <= shiftCount
) {
2915 if ( 0x43E <= aExp
) {
2916 if ( float64_val(a
) != LIT64( 0xC3E0000000000000 ) ) {
2917 float_raise( float_flag_invalid STATUS_VAR
);
2919 || ( ( aExp
== 0x7FF )
2920 && ( aSig
!= LIT64( 0x0010000000000000 ) ) )
2922 return LIT64( 0x7FFFFFFFFFFFFFFF );
2925 return (int64_t) LIT64( 0x8000000000000000 );
2927 z
= aSig
<<shiftCount
;
2930 if ( aExp
< 0x3FE ) {
2931 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
2934 z
= aSig
>>( - shiftCount
);
2935 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
2936 STATUS(float_exception_flags
) |= float_flag_inexact
;
2939 if ( aSign
) z
= - z
;
2944 /*----------------------------------------------------------------------------
2945 | Returns the result of converting the double-precision floating-point value
2946 | `a' to the single-precision floating-point format. The conversion is
2947 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2949 *----------------------------------------------------------------------------*/
2951 float32
float64_to_float32( float64 a STATUS_PARAM
)
2957 a
= float64_squash_input_denormal(a STATUS_VAR
);
2959 aSig
= extractFloat64Frac( a
);
2960 aExp
= extractFloat64Exp( a
);
2961 aSign
= extractFloat64Sign( a
);
2962 if ( aExp
== 0x7FF ) {
2963 if ( aSig
) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
2964 return packFloat32( aSign
, 0xFF, 0 );
2966 shift64RightJamming( aSig
, 22, &aSig
);
2968 if ( aExp
|| zSig
) {
2972 return roundAndPackFloat32( aSign
, aExp
, zSig STATUS_VAR
);
2977 /*----------------------------------------------------------------------------
2978 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
2979 | half-precision floating-point value, returning the result. After being
2980 | shifted into the proper positions, the three fields are simply added
2981 | together to form the result. This means that any integer portion of `zSig'
2982 | will be added into the exponent. Since a properly normalized significand
2983 | will have an integer portion equal to 1, the `zExp' input should be 1 less
2984 | than the desired result exponent whenever `zSig' is a complete, normalized
2986 *----------------------------------------------------------------------------*/
2987 static float16
packFloat16(flag zSign
, int16 zExp
, uint16_t zSig
)
2989 return make_float16(
2990 (((uint32_t)zSign
) << 15) + (((uint32_t)zExp
) << 10) + zSig
);
2993 /* Half precision floats come in two formats: standard IEEE and "ARM" format.
2994 The latter gains extra exponent range by omitting the NaN/Inf encodings. */
2996 float32
float16_to_float32(float16 a
, flag ieee STATUS_PARAM
)
3002 aSign
= extractFloat16Sign(a
);
3003 aExp
= extractFloat16Exp(a
);
3004 aSig
= extractFloat16Frac(a
);
3006 if (aExp
== 0x1f && ieee
) {
3008 return commonNaNToFloat32(float16ToCommonNaN(a STATUS_VAR
) STATUS_VAR
);
3010 return packFloat32(aSign
, 0xff, aSig
<< 13);
3016 return packFloat32(aSign
, 0, 0);
3019 shiftCount
= countLeadingZeros32( aSig
) - 21;
3020 aSig
= aSig
<< shiftCount
;
3023 return packFloat32( aSign
, aExp
+ 0x70, aSig
<< 13);
3026 float16
float32_to_float16(float32 a
, flag ieee STATUS_PARAM
)
3034 a
= float32_squash_input_denormal(a STATUS_VAR
);
3036 aSig
= extractFloat32Frac( a
);
3037 aExp
= extractFloat32Exp( a
);
3038 aSign
= extractFloat32Sign( a
);
3039 if ( aExp
== 0xFF ) {
3041 /* Input is a NaN */
3042 float16 r
= commonNaNToFloat16( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
3044 return packFloat16(aSign
, 0, 0);
3050 float_raise(float_flag_invalid STATUS_VAR
);
3051 return packFloat16(aSign
, 0x1f, 0x3ff);
3053 return packFloat16(aSign
, 0x1f, 0);
3055 if (aExp
== 0 && aSig
== 0) {
3056 return packFloat16(aSign
, 0, 0);
3058 /* Decimal point between bits 22 and 23. */
3070 float_raise( float_flag_underflow STATUS_VAR
);
3071 roundingMode
= STATUS(float_rounding_mode
);
3072 switch (roundingMode
) {
3073 case float_round_nearest_even
:
3074 increment
= (mask
+ 1) >> 1;
3075 if ((aSig
& mask
) == increment
) {
3076 increment
= aSig
& (increment
<< 1);
3079 case float_round_up
:
3080 increment
= aSign
? 0 : mask
;
3082 case float_round_down
:
3083 increment
= aSign
? mask
: 0;
3085 default: /* round_to_zero */
3090 if (aSig
>= 0x01000000) {
3094 } else if (aExp
< -14
3095 && STATUS(float_detect_tininess
) == float_tininess_before_rounding
) {
3096 float_raise( float_flag_underflow STATUS_VAR
);
3101 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
3102 return packFloat16(aSign
, 0x1f, 0);
3106 float_raise(float_flag_invalid
| float_flag_inexact STATUS_VAR
);
3107 return packFloat16(aSign
, 0x1f, 0x3ff);
3111 return packFloat16(aSign
, 0, 0);
3114 aSig
>>= -14 - aExp
;
3117 return packFloat16(aSign
, aExp
+ 14, aSig
>> 13);
3120 /*----------------------------------------------------------------------------
3121 | Returns the result of converting the double-precision floating-point value
3122 | `a' to the extended double-precision floating-point format. The conversion
3123 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
3125 *----------------------------------------------------------------------------*/
3127 floatx80
float64_to_floatx80( float64 a STATUS_PARAM
)
3133 a
= float64_squash_input_denormal(a STATUS_VAR
);
3134 aSig
= extractFloat64Frac( a
);
3135 aExp
= extractFloat64Exp( a
);
3136 aSign
= extractFloat64Sign( a
);
3137 if ( aExp
== 0x7FF ) {
3138 if ( aSig
) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
3139 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
3142 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
3143 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3147 aSign
, aExp
+ 0x3C00, ( aSig
| LIT64( 0x0010000000000000 ) )<<11 );
3151 /*----------------------------------------------------------------------------
3152 | Returns the result of converting the double-precision floating-point value
3153 | `a' to the quadruple-precision floating-point format. The conversion is
3154 | performed according to the IEC/IEEE Standard for Binary Floating-Point
3156 *----------------------------------------------------------------------------*/
3158 float128
float64_to_float128( float64 a STATUS_PARAM
)
3162 uint64_t aSig
, zSig0
, zSig1
;
3164 a
= float64_squash_input_denormal(a STATUS_VAR
);
3165 aSig
= extractFloat64Frac( a
);
3166 aExp
= extractFloat64Exp( a
);
3167 aSign
= extractFloat64Sign( a
);
3168 if ( aExp
== 0x7FF ) {
3169 if ( aSig
) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
3170 return packFloat128( aSign
, 0x7FFF, 0, 0 );
3173 if ( aSig
== 0 ) return packFloat128( aSign
, 0, 0, 0 );
3174 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3177 shift128Right( aSig
, 0, 4, &zSig0
, &zSig1
);
3178 return packFloat128( aSign
, aExp
+ 0x3C00, zSig0
, zSig1
);
3182 /*----------------------------------------------------------------------------
3183 | Rounds the double-precision floating-point value `a' to an integer, and
3184 | returns the result as a double-precision floating-point value. The
3185 | operation is performed according to the IEC/IEEE Standard for Binary
3186 | Floating-Point Arithmetic.
3187 *----------------------------------------------------------------------------*/
3189 float64
float64_round_to_int( float64 a STATUS_PARAM
)
3193 uint64_t lastBitMask
, roundBitsMask
;
3196 a
= float64_squash_input_denormal(a STATUS_VAR
);
3198 aExp
= extractFloat64Exp( a
);
3199 if ( 0x433 <= aExp
) {
3200 if ( ( aExp
== 0x7FF ) && extractFloat64Frac( a
) ) {
3201 return propagateFloat64NaN( a
, a STATUS_VAR
);
3205 if ( aExp
< 0x3FF ) {
3206 if ( (uint64_t) ( float64_val(a
)<<1 ) == 0 ) return a
;
3207 STATUS(float_exception_flags
) |= float_flag_inexact
;
3208 aSign
= extractFloat64Sign( a
);
3209 switch ( STATUS(float_rounding_mode
) ) {
3210 case float_round_nearest_even
:
3211 if ( ( aExp
== 0x3FE ) && extractFloat64Frac( a
) ) {
3212 return packFloat64( aSign
, 0x3FF, 0 );
3215 case float_round_down
:
3216 return make_float64(aSign
? LIT64( 0xBFF0000000000000 ) : 0);
3217 case float_round_up
:
3218 return make_float64(
3219 aSign
? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 ));
3221 return packFloat64( aSign
, 0, 0 );
3224 lastBitMask
<<= 0x433 - aExp
;
3225 roundBitsMask
= lastBitMask
- 1;
3227 roundingMode
= STATUS(float_rounding_mode
);
3228 if ( roundingMode
== float_round_nearest_even
) {
3229 z
+= lastBitMask
>>1;
3230 if ( ( z
& roundBitsMask
) == 0 ) z
&= ~ lastBitMask
;
3232 else if ( roundingMode
!= float_round_to_zero
) {
3233 if ( extractFloat64Sign( make_float64(z
) ) ^ ( roundingMode
== float_round_up
) ) {
3237 z
&= ~ roundBitsMask
;
3238 if ( z
!= float64_val(a
) )
3239 STATUS(float_exception_flags
) |= float_flag_inexact
;
3240 return make_float64(z
);
3244 float64
float64_trunc_to_int( float64 a STATUS_PARAM
)
3248 oldmode
= STATUS(float_rounding_mode
);
3249 STATUS(float_rounding_mode
) = float_round_to_zero
;
3250 res
= float64_round_to_int(a STATUS_VAR
);
3251 STATUS(float_rounding_mode
) = oldmode
;
3255 /*----------------------------------------------------------------------------
3256 | Returns the result of adding the absolute values of the double-precision
3257 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
3258 | before being returned. `zSign' is ignored if the result is a NaN.
3259 | The addition is performed according to the IEC/IEEE Standard for Binary
3260 | Floating-Point Arithmetic.
3261 *----------------------------------------------------------------------------*/
3263 static float64
addFloat64Sigs( float64 a
, float64 b
, flag zSign STATUS_PARAM
)
3265 int16 aExp
, bExp
, zExp
;
3266 uint64_t aSig
, bSig
, zSig
;
3269 aSig
= extractFloat64Frac( a
);
3270 aExp
= extractFloat64Exp( a
);
3271 bSig
= extractFloat64Frac( b
);
3272 bExp
= extractFloat64Exp( b
);
3273 expDiff
= aExp
- bExp
;
3276 if ( 0 < expDiff
) {
3277 if ( aExp
== 0x7FF ) {
3278 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3285 bSig
|= LIT64( 0x2000000000000000 );
3287 shift64RightJamming( bSig
, expDiff
, &bSig
);
3290 else if ( expDiff
< 0 ) {
3291 if ( bExp
== 0x7FF ) {
3292 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3293 return packFloat64( zSign
, 0x7FF, 0 );
3299 aSig
|= LIT64( 0x2000000000000000 );
3301 shift64RightJamming( aSig
, - expDiff
, &aSig
);
3305 if ( aExp
== 0x7FF ) {
3306 if ( aSig
| bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3310 if (STATUS(flush_to_zero
)) {
3312 float_raise(float_flag_output_denormal STATUS_VAR
);
3314 return packFloat64(zSign
, 0, 0);
3316 return packFloat64( zSign
, 0, ( aSig
+ bSig
)>>9 );
3318 zSig
= LIT64( 0x4000000000000000 ) + aSig
+ bSig
;
3322 aSig
|= LIT64( 0x2000000000000000 );
3323 zSig
= ( aSig
+ bSig
)<<1;
3325 if ( (int64_t) zSig
< 0 ) {
3330 return roundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3334 /*----------------------------------------------------------------------------
3335 | Returns the result of subtracting the absolute values of the double-
3336 | precision floating-point values `a' and `b'. If `zSign' is 1, the
3337 | difference is negated before being returned. `zSign' is ignored if the
3338 | result is a NaN. The subtraction is performed according to the IEC/IEEE
3339 | Standard for Binary Floating-Point Arithmetic.
3340 *----------------------------------------------------------------------------*/
3342 static float64
subFloat64Sigs( float64 a
, float64 b
, flag zSign STATUS_PARAM
)
3344 int16 aExp
, bExp
, zExp
;
3345 uint64_t aSig
, bSig
, zSig
;
3348 aSig
= extractFloat64Frac( a
);
3349 aExp
= extractFloat64Exp( a
);
3350 bSig
= extractFloat64Frac( b
);
3351 bExp
= extractFloat64Exp( b
);
3352 expDiff
= aExp
- bExp
;
3355 if ( 0 < expDiff
) goto aExpBigger
;
3356 if ( expDiff
< 0 ) goto bExpBigger
;
3357 if ( aExp
== 0x7FF ) {
3358 if ( aSig
| bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3359 float_raise( float_flag_invalid STATUS_VAR
);
3360 return float64_default_nan
;
3366 if ( bSig
< aSig
) goto aBigger
;
3367 if ( aSig
< bSig
) goto bBigger
;
3368 return packFloat64( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
3370 if ( bExp
== 0x7FF ) {
3371 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3372 return packFloat64( zSign
^ 1, 0x7FF, 0 );
3378 aSig
|= LIT64( 0x4000000000000000 );
3380 shift64RightJamming( aSig
, - expDiff
, &aSig
);
3381 bSig
|= LIT64( 0x4000000000000000 );
3386 goto normalizeRoundAndPack
;
3388 if ( aExp
== 0x7FF ) {
3389 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3396 bSig
|= LIT64( 0x4000000000000000 );
3398 shift64RightJamming( bSig
, expDiff
, &bSig
);
3399 aSig
|= LIT64( 0x4000000000000000 );
3403 normalizeRoundAndPack
:
3405 return normalizeRoundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3409 /*----------------------------------------------------------------------------
3410 | Returns the result of adding the double-precision floating-point values `a'
3411 | and `b'. The operation is performed according to the IEC/IEEE Standard for
3412 | Binary Floating-Point Arithmetic.
3413 *----------------------------------------------------------------------------*/
3415 float64
float64_add( float64 a
, float64 b STATUS_PARAM
)
3418 a
= float64_squash_input_denormal(a STATUS_VAR
);
3419 b
= float64_squash_input_denormal(b STATUS_VAR
);
3421 aSign
= extractFloat64Sign( a
);
3422 bSign
= extractFloat64Sign( b
);
3423 if ( aSign
== bSign
) {
3424 return addFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3427 return subFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3432 /*----------------------------------------------------------------------------
3433 | Returns the result of subtracting the double-precision floating-point values
3434 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
3435 | for Binary Floating-Point Arithmetic.
3436 *----------------------------------------------------------------------------*/
3438 float64
float64_sub( float64 a
, float64 b STATUS_PARAM
)
3441 a
= float64_squash_input_denormal(a STATUS_VAR
);
3442 b
= float64_squash_input_denormal(b STATUS_VAR
);
3444 aSign
= extractFloat64Sign( a
);
3445 bSign
= extractFloat64Sign( b
);
3446 if ( aSign
== bSign
) {
3447 return subFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3450 return addFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3455 /*----------------------------------------------------------------------------
3456 | Returns the result of multiplying the double-precision floating-point values
3457 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
3458 | for Binary Floating-Point Arithmetic.
3459 *----------------------------------------------------------------------------*/
3461 float64
float64_mul( float64 a
, float64 b STATUS_PARAM
)
3463 flag aSign
, bSign
, zSign
;
3464 int16 aExp
, bExp
, zExp
;
3465 uint64_t aSig
, bSig
, zSig0
, zSig1
;
3467 a
= float64_squash_input_denormal(a STATUS_VAR
);
3468 b
= float64_squash_input_denormal(b STATUS_VAR
);
3470 aSig
= extractFloat64Frac( a
);
3471 aExp
= extractFloat64Exp( a
);
3472 aSign
= extractFloat64Sign( a
);
3473 bSig
= extractFloat64Frac( b
);
3474 bExp
= extractFloat64Exp( b
);
3475 bSign
= extractFloat64Sign( b
);
3476 zSign
= aSign
^ bSign
;
3477 if ( aExp
== 0x7FF ) {
3478 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
3479 return propagateFloat64NaN( a
, b STATUS_VAR
);
3481 if ( ( bExp
| bSig
) == 0 ) {
3482 float_raise( float_flag_invalid STATUS_VAR
);
3483 return float64_default_nan
;
3485 return packFloat64( zSign
, 0x7FF, 0 );
3487 if ( bExp
== 0x7FF ) {
3488 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3489 if ( ( aExp
| aSig
) == 0 ) {
3490 float_raise( float_flag_invalid STATUS_VAR
);
3491 return float64_default_nan
;
3493 return packFloat64( zSign
, 0x7FF, 0 );
3496 if ( aSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3497 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3500 if ( bSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3501 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3503 zExp
= aExp
+ bExp
- 0x3FF;
3504 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<10;
3505 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3506 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
3507 zSig0
|= ( zSig1
!= 0 );
3508 if ( 0 <= (int64_t) ( zSig0
<<1 ) ) {
3512 return roundAndPackFloat64( zSign
, zExp
, zSig0 STATUS_VAR
);
3516 /*----------------------------------------------------------------------------
3517 | Returns the result of dividing the double-precision floating-point value `a'
3518 | by the corresponding value `b'. The operation is performed according to
3519 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3520 *----------------------------------------------------------------------------*/
3522 float64
float64_div( float64 a
, float64 b STATUS_PARAM
)
3524 flag aSign
, bSign
, zSign
;
3525 int16 aExp
, bExp
, zExp
;
3526 uint64_t aSig
, bSig
, zSig
;
3527 uint64_t rem0
, rem1
;
3528 uint64_t term0
, term1
;
3529 a
= float64_squash_input_denormal(a STATUS_VAR
);
3530 b
= float64_squash_input_denormal(b STATUS_VAR
);
3532 aSig
= extractFloat64Frac( a
);
3533 aExp
= extractFloat64Exp( a
);
3534 aSign
= extractFloat64Sign( a
);
3535 bSig
= extractFloat64Frac( b
);
3536 bExp
= extractFloat64Exp( b
);
3537 bSign
= extractFloat64Sign( b
);
3538 zSign
= aSign
^ bSign
;
3539 if ( aExp
== 0x7FF ) {
3540 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3541 if ( bExp
== 0x7FF ) {
3542 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3543 float_raise( float_flag_invalid STATUS_VAR
);
3544 return float64_default_nan
;
3546 return packFloat64( zSign
, 0x7FF, 0 );
3548 if ( bExp
== 0x7FF ) {
3549 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3550 return packFloat64( zSign
, 0, 0 );
3554 if ( ( aExp
| aSig
) == 0 ) {
3555 float_raise( float_flag_invalid STATUS_VAR
);
3556 return float64_default_nan
;
3558 float_raise( float_flag_divbyzero STATUS_VAR
);
3559 return packFloat64( zSign
, 0x7FF, 0 );
3561 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3564 if ( aSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3565 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3567 zExp
= aExp
- bExp
+ 0x3FD;
3568 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<10;
3569 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3570 if ( bSig
<= ( aSig
+ aSig
) ) {
3574 zSig
= estimateDiv128To64( aSig
, 0, bSig
);
3575 if ( ( zSig
& 0x1FF ) <= 2 ) {
3576 mul64To128( bSig
, zSig
, &term0
, &term1
);
3577 sub128( aSig
, 0, term0
, term1
, &rem0
, &rem1
);
3578 while ( (int64_t) rem0
< 0 ) {
3580 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
3582 zSig
|= ( rem1
!= 0 );
3584 return roundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3588 /*----------------------------------------------------------------------------
3589 | Returns the remainder of the double-precision floating-point value `a'
3590 | with respect to the corresponding value `b'. The operation is performed
3591 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3592 *----------------------------------------------------------------------------*/
3594 float64
float64_rem( float64 a
, float64 b STATUS_PARAM
)
3597 int16 aExp
, bExp
, expDiff
;
3598 uint64_t aSig
, bSig
;
3599 uint64_t q
, alternateASig
;
3602 a
= float64_squash_input_denormal(a STATUS_VAR
);
3603 b
= float64_squash_input_denormal(b STATUS_VAR
);
3604 aSig
= extractFloat64Frac( a
);
3605 aExp
= extractFloat64Exp( a
);
3606 aSign
= extractFloat64Sign( a
);
3607 bSig
= extractFloat64Frac( b
);
3608 bExp
= extractFloat64Exp( b
);
3609 if ( aExp
== 0x7FF ) {
3610 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
3611 return propagateFloat64NaN( a
, b STATUS_VAR
);
3613 float_raise( float_flag_invalid STATUS_VAR
);
3614 return float64_default_nan
;
3616 if ( bExp
== 0x7FF ) {
3617 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3622 float_raise( float_flag_invalid STATUS_VAR
);
3623 return float64_default_nan
;
3625 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3628 if ( aSig
== 0 ) return a
;
3629 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3631 expDiff
= aExp
- bExp
;
3632 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<11;
3633 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3634 if ( expDiff
< 0 ) {
3635 if ( expDiff
< -1 ) return a
;
3638 q
= ( bSig
<= aSig
);
3639 if ( q
) aSig
-= bSig
;
3641 while ( 0 < expDiff
) {
3642 q
= estimateDiv128To64( aSig
, 0, bSig
);
3643 q
= ( 2 < q
) ? q
- 2 : 0;
3644 aSig
= - ( ( bSig
>>2 ) * q
);
3648 if ( 0 < expDiff
) {
3649 q
= estimateDiv128To64( aSig
, 0, bSig
);
3650 q
= ( 2 < q
) ? q
- 2 : 0;
3653 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
3660 alternateASig
= aSig
;
3663 } while ( 0 <= (int64_t) aSig
);
3664 sigMean
= aSig
+ alternateASig
;
3665 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
3666 aSig
= alternateASig
;
3668 zSign
= ( (int64_t) aSig
< 0 );
3669 if ( zSign
) aSig
= - aSig
;
3670 return normalizeRoundAndPackFloat64( aSign
^ zSign
, bExp
, aSig STATUS_VAR
);
3674 /*----------------------------------------------------------------------------
3675 | Returns the result of multiplying the double-precision floating-point values
3676 | `a' and `b' then adding 'c', with no intermediate rounding step after the
3677 | multiplication. The operation is performed according to the IEC/IEEE
3678 | Standard for Binary Floating-Point Arithmetic 754-2008.
3679 | The flags argument allows the caller to select negation of the
3680 | addend, the intermediate product, or the final result. (The difference
3681 | between this and having the caller do a separate negation is that negating
3682 | externally will flip the sign bit on NaNs.)
3683 *----------------------------------------------------------------------------*/
3685 float64
float64_muladd(float64 a
, float64 b
, float64 c
, int flags STATUS_PARAM
)
3687 flag aSign
, bSign
, cSign
, zSign
;
3688 int aExp
, bExp
, cExp
, pExp
, zExp
, expDiff
;
3689 uint64_t aSig
, bSig
, cSig
;
3690 flag pInf
, pZero
, pSign
;
3691 uint64_t pSig0
, pSig1
, cSig0
, cSig1
, zSig0
, zSig1
;
3693 flag signflip
, infzero
;
3695 a
= float64_squash_input_denormal(a STATUS_VAR
);
3696 b
= float64_squash_input_denormal(b STATUS_VAR
);
3697 c
= float64_squash_input_denormal(c STATUS_VAR
);
3698 aSig
= extractFloat64Frac(a
);
3699 aExp
= extractFloat64Exp(a
);
3700 aSign
= extractFloat64Sign(a
);
3701 bSig
= extractFloat64Frac(b
);
3702 bExp
= extractFloat64Exp(b
);
3703 bSign
= extractFloat64Sign(b
);
3704 cSig
= extractFloat64Frac(c
);
3705 cExp
= extractFloat64Exp(c
);
3706 cSign
= extractFloat64Sign(c
);
3708 infzero
= ((aExp
== 0 && aSig
== 0 && bExp
== 0x7ff && bSig
== 0) ||
3709 (aExp
== 0x7ff && aSig
== 0 && bExp
== 0 && bSig
== 0));
3711 /* It is implementation-defined whether the cases of (0,inf,qnan)
3712 * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
3713 * they return if they do), so we have to hand this information
3714 * off to the target-specific pick-a-NaN routine.
3716 if (((aExp
== 0x7ff) && aSig
) ||
3717 ((bExp
== 0x7ff) && bSig
) ||
3718 ((cExp
== 0x7ff) && cSig
)) {
3719 return propagateFloat64MulAddNaN(a
, b
, c
, infzero STATUS_VAR
);
3723 float_raise(float_flag_invalid STATUS_VAR
);
3724 return float64_default_nan
;
3727 if (flags
& float_muladd_negate_c
) {
3731 signflip
= (flags
& float_muladd_negate_result
) ? 1 : 0;
3733 /* Work out the sign and type of the product */
3734 pSign
= aSign
^ bSign
;
3735 if (flags
& float_muladd_negate_product
) {
3738 pInf
= (aExp
== 0x7ff) || (bExp
== 0x7ff);
3739 pZero
= ((aExp
| aSig
) == 0) || ((bExp
| bSig
) == 0);
3741 if (cExp
== 0x7ff) {
3742 if (pInf
&& (pSign
^ cSign
)) {
3743 /* addition of opposite-signed infinities => InvalidOperation */
3744 float_raise(float_flag_invalid STATUS_VAR
);
3745 return float64_default_nan
;
3747 /* Otherwise generate an infinity of the same sign */
3748 return packFloat64(cSign
^ signflip
, 0x7ff, 0);
3752 return packFloat64(pSign
^ signflip
, 0x7ff, 0);
3758 /* Adding two exact zeroes */
3759 if (pSign
== cSign
) {
3761 } else if (STATUS(float_rounding_mode
) == float_round_down
) {
3766 return packFloat64(zSign
^ signflip
, 0, 0);
3768 /* Exact zero plus a denorm */
3769 if (STATUS(flush_to_zero
)) {
3770 float_raise(float_flag_output_denormal STATUS_VAR
);
3771 return packFloat64(cSign
^ signflip
, 0, 0);
3774 /* Zero plus something non-zero : just return the something */
3775 return make_float64(float64_val(c
) ^ ((uint64_t)signflip
<< 63));
3779 normalizeFloat64Subnormal(aSig
, &aExp
, &aSig
);
3782 normalizeFloat64Subnormal(bSig
, &bExp
, &bSig
);
3785 /* Calculate the actual result a * b + c */
3787 /* Multiply first; this is easy. */
3788 /* NB: we subtract 0x3fe where float64_mul() subtracts 0x3ff
3789 * because we want the true exponent, not the "one-less-than"
3790 * flavour that roundAndPackFloat64() takes.
3792 pExp
= aExp
+ bExp
- 0x3fe;
3793 aSig
= (aSig
| LIT64(0x0010000000000000))<<10;
3794 bSig
= (bSig
| LIT64(0x0010000000000000))<<11;
3795 mul64To128(aSig
, bSig
, &pSig0
, &pSig1
);
3796 if ((int64_t)(pSig0
<< 1) >= 0) {
3797 shortShift128Left(pSig0
, pSig1
, 1, &pSig0
, &pSig1
);
3801 zSign
= pSign
^ signflip
;
3803 /* Now [pSig0:pSig1] is the significand of the multiply, with the explicit
3804 * bit in position 126.
3808 /* Throw out the special case of c being an exact zero now */
3809 shift128RightJamming(pSig0
, pSig1
, 64, &pSig0
, &pSig1
);
3810 return roundAndPackFloat64(zSign
, pExp
- 1,
3813 normalizeFloat64Subnormal(cSig
, &cExp
, &cSig
);
3816 /* Shift cSig and add the explicit bit so [cSig0:cSig1] is the
3817 * significand of the addend, with the explicit bit in position 126.
3819 cSig0
= cSig
<< (126 - 64 - 52);
3821 cSig0
|= LIT64(0x4000000000000000);
3822 expDiff
= pExp
- cExp
;
3824 if (pSign
== cSign
) {
3827 /* scale c to match p */
3828 shift128RightJamming(cSig0
, cSig1
, expDiff
, &cSig0
, &cSig1
);
3830 } else if (expDiff
< 0) {
3831 /* scale p to match c */
3832 shift128RightJamming(pSig0
, pSig1
, -expDiff
, &pSig0
, &pSig1
);
3835 /* no scaling needed */
3838 /* Add significands and make sure explicit bit ends up in posn 126 */
3839 add128(pSig0
, pSig1
, cSig0
, cSig1
, &zSig0
, &zSig1
);
3840 if ((int64_t)zSig0
< 0) {
3841 shift128RightJamming(zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
3845 shift128RightJamming(zSig0
, zSig1
, 64, &zSig0
, &zSig1
);
3846 return roundAndPackFloat64(zSign
, zExp
, zSig1 STATUS_VAR
);
3850 shift128RightJamming(cSig0
, cSig1
, expDiff
, &cSig0
, &cSig1
);
3851 sub128(pSig0
, pSig1
, cSig0
, cSig1
, &zSig0
, &zSig1
);
3853 } else if (expDiff
< 0) {
3854 shift128RightJamming(pSig0
, pSig1
, -expDiff
, &pSig0
, &pSig1
);
3855 sub128(cSig0
, cSig1
, pSig0
, pSig1
, &zSig0
, &zSig1
);
3860 if (lt128(cSig0
, cSig1
, pSig0
, pSig1
)) {
3861 sub128(pSig0
, pSig1
, cSig0
, cSig1
, &zSig0
, &zSig1
);
3862 } else if (lt128(pSig0
, pSig1
, cSig0
, cSig1
)) {
3863 sub128(cSig0
, cSig1
, pSig0
, pSig1
, &zSig0
, &zSig1
);
3868 if (STATUS(float_rounding_mode
) == float_round_down
) {
3871 return packFloat64(zSign
, 0, 0);
3875 /* Do the equivalent of normalizeRoundAndPackFloat64() but
3876 * starting with the significand in a pair of uint64_t.
3879 shiftcount
= countLeadingZeros64(zSig0
) - 1;
3880 shortShift128Left(zSig0
, zSig1
, shiftcount
, &zSig0
, &zSig1
);
3886 shiftcount
= countLeadingZeros64(zSig1
) - 1;
3887 zSig0
= zSig1
<< shiftcount
;
3888 zExp
-= (shiftcount
+ 64);
3890 return roundAndPackFloat64(zSign
, zExp
, zSig0 STATUS_VAR
);
3894 /*----------------------------------------------------------------------------
3895 | Returns the square root of the double-precision floating-point value `a'.
3896 | The operation is performed according to the IEC/IEEE Standard for Binary
3897 | Floating-Point Arithmetic.
3898 *----------------------------------------------------------------------------*/
3900 float64
float64_sqrt( float64 a STATUS_PARAM
)
3904 uint64_t aSig
, zSig
, doubleZSig
;
3905 uint64_t rem0
, rem1
, term0
, term1
;
3906 a
= float64_squash_input_denormal(a STATUS_VAR
);
3908 aSig
= extractFloat64Frac( a
);
3909 aExp
= extractFloat64Exp( a
);
3910 aSign
= extractFloat64Sign( a
);
3911 if ( aExp
== 0x7FF ) {
3912 if ( aSig
) return propagateFloat64NaN( a
, a STATUS_VAR
);
3913 if ( ! aSign
) return a
;
3914 float_raise( float_flag_invalid STATUS_VAR
);
3915 return float64_default_nan
;
3918 if ( ( aExp
| aSig
) == 0 ) return a
;
3919 float_raise( float_flag_invalid STATUS_VAR
);
3920 return float64_default_nan
;
3923 if ( aSig
== 0 ) return float64_zero
;
3924 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3926 zExp
= ( ( aExp
- 0x3FF )>>1 ) + 0x3FE;
3927 aSig
|= LIT64( 0x0010000000000000 );
3928 zSig
= estimateSqrt32( aExp
, aSig
>>21 );
3929 aSig
<<= 9 - ( aExp
& 1 );
3930 zSig
= estimateDiv128To64( aSig
, 0, zSig
<<32 ) + ( zSig
<<30 );
3931 if ( ( zSig
& 0x1FF ) <= 5 ) {
3932 doubleZSig
= zSig
<<1;
3933 mul64To128( zSig
, zSig
, &term0
, &term1
);
3934 sub128( aSig
, 0, term0
, term1
, &rem0
, &rem1
);
3935 while ( (int64_t) rem0
< 0 ) {
3938 add128( rem0
, rem1
, zSig
>>63, doubleZSig
| 1, &rem0
, &rem1
);
3940 zSig
|= ( ( rem0
| rem1
) != 0 );
3942 return roundAndPackFloat64( 0, zExp
, zSig STATUS_VAR
);
3946 /*----------------------------------------------------------------------------
3947 | Returns the binary log of the double-precision floating-point value `a'.
3948 | The operation is performed according to the IEC/IEEE Standard for Binary
3949 | Floating-Point Arithmetic.
3950 *----------------------------------------------------------------------------*/
3951 float64
float64_log2( float64 a STATUS_PARAM
)
3955 uint64_t aSig
, aSig0
, aSig1
, zSig
, i
;
3956 a
= float64_squash_input_denormal(a STATUS_VAR
);
3958 aSig
= extractFloat64Frac( a
);
3959 aExp
= extractFloat64Exp( a
);
3960 aSign
= extractFloat64Sign( a
);
3963 if ( aSig
== 0 ) return packFloat64( 1, 0x7FF, 0 );
3964 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3967 float_raise( float_flag_invalid STATUS_VAR
);
3968 return float64_default_nan
;
3970 if ( aExp
== 0x7FF ) {
3971 if ( aSig
) return propagateFloat64NaN( a
, float64_zero STATUS_VAR
);
3976 aSig
|= LIT64( 0x0010000000000000 );
3978 zSig
= (uint64_t)aExp
<< 52;
3979 for (i
= 1LL << 51; i
> 0; i
>>= 1) {
3980 mul64To128( aSig
, aSig
, &aSig0
, &aSig1
);
3981 aSig
= ( aSig0
<< 12 ) | ( aSig1
>> 52 );
3982 if ( aSig
& LIT64( 0x0020000000000000 ) ) {
3990 return normalizeRoundAndPackFloat64( zSign
, 0x408, zSig STATUS_VAR
);
3993 /*----------------------------------------------------------------------------
3994 | Returns 1 if the double-precision floating-point value `a' is equal to the
3995 | corresponding value `b', and 0 otherwise. The invalid exception is raised
3996 | if either operand is a NaN. Otherwise, the comparison is performed
3997 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3998 *----------------------------------------------------------------------------*/
4000 int float64_eq( float64 a
, float64 b STATUS_PARAM
)
4003 a
= float64_squash_input_denormal(a STATUS_VAR
);
4004 b
= float64_squash_input_denormal(b STATUS_VAR
);
4006 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4007 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4009 float_raise( float_flag_invalid STATUS_VAR
);
4012 av
= float64_val(a
);
4013 bv
= float64_val(b
);
4014 return ( av
== bv
) || ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
4018 /*----------------------------------------------------------------------------
4019 | Returns 1 if the double-precision floating-point value `a' is less than or
4020 | equal to the corresponding value `b', and 0 otherwise. The invalid
4021 | exception is raised if either operand is a NaN. The comparison is performed
4022 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4023 *----------------------------------------------------------------------------*/
4025 int float64_le( float64 a
, float64 b STATUS_PARAM
)
4029 a
= float64_squash_input_denormal(a STATUS_VAR
);
4030 b
= float64_squash_input_denormal(b STATUS_VAR
);
4032 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4033 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4035 float_raise( float_flag_invalid STATUS_VAR
);
4038 aSign
= extractFloat64Sign( a
);
4039 bSign
= extractFloat64Sign( b
);
4040 av
= float64_val(a
);
4041 bv
= float64_val(b
);
4042 if ( aSign
!= bSign
) return aSign
|| ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
4043 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
4047 /*----------------------------------------------------------------------------
4048 | Returns 1 if the double-precision floating-point value `a' is less than
4049 | the corresponding value `b', and 0 otherwise. The invalid exception is
4050 | raised if either operand is a NaN. The comparison is performed according
4051 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4052 *----------------------------------------------------------------------------*/
4054 int float64_lt( float64 a
, float64 b STATUS_PARAM
)
4059 a
= float64_squash_input_denormal(a STATUS_VAR
);
4060 b
= float64_squash_input_denormal(b STATUS_VAR
);
4061 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4062 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4064 float_raise( float_flag_invalid STATUS_VAR
);
4067 aSign
= extractFloat64Sign( a
);
4068 bSign
= extractFloat64Sign( b
);
4069 av
= float64_val(a
);
4070 bv
= float64_val(b
);
4071 if ( aSign
!= bSign
) return aSign
&& ( (uint64_t) ( ( av
| bv
)<<1 ) != 0 );
4072 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
4076 /*----------------------------------------------------------------------------
4077 | Returns 1 if the double-precision floating-point values `a' and `b' cannot
4078 | be compared, and 0 otherwise. The invalid exception is raised if either
4079 | operand is a NaN. The comparison is performed according to the IEC/IEEE
4080 | Standard for Binary Floating-Point Arithmetic.
4081 *----------------------------------------------------------------------------*/
4083 int float64_unordered( float64 a
, float64 b STATUS_PARAM
)
4085 a
= float64_squash_input_denormal(a STATUS_VAR
);
4086 b
= float64_squash_input_denormal(b STATUS_VAR
);
4088 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4089 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4091 float_raise( float_flag_invalid STATUS_VAR
);
4097 /*----------------------------------------------------------------------------
4098 | Returns 1 if the double-precision floating-point value `a' is equal to the
4099 | corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
4100 | exception.The comparison is performed according to the IEC/IEEE Standard
4101 | for Binary Floating-Point Arithmetic.
4102 *----------------------------------------------------------------------------*/
4104 int float64_eq_quiet( float64 a
, float64 b STATUS_PARAM
)
4107 a
= float64_squash_input_denormal(a STATUS_VAR
);
4108 b
= float64_squash_input_denormal(b STATUS_VAR
);
4110 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4111 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4113 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
4114 float_raise( float_flag_invalid STATUS_VAR
);
4118 av
= float64_val(a
);
4119 bv
= float64_val(b
);
4120 return ( av
== bv
) || ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
4124 /*----------------------------------------------------------------------------
4125 | Returns 1 if the double-precision floating-point value `a' is less than or
4126 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
4127 | cause an exception. Otherwise, the comparison is performed according to the
4128 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4129 *----------------------------------------------------------------------------*/
4131 int float64_le_quiet( float64 a
, float64 b STATUS_PARAM
)
4135 a
= float64_squash_input_denormal(a STATUS_VAR
);
4136 b
= float64_squash_input_denormal(b STATUS_VAR
);
4138 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4139 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4141 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
4142 float_raise( float_flag_invalid STATUS_VAR
);
4146 aSign
= extractFloat64Sign( a
);
4147 bSign
= extractFloat64Sign( b
);
4148 av
= float64_val(a
);
4149 bv
= float64_val(b
);
4150 if ( aSign
!= bSign
) return aSign
|| ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
4151 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
4155 /*----------------------------------------------------------------------------
4156 | Returns 1 if the double-precision floating-point value `a' is less than
4157 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
4158 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
4159 | Standard for Binary Floating-Point Arithmetic.
4160 *----------------------------------------------------------------------------*/
4162 int float64_lt_quiet( float64 a
, float64 b STATUS_PARAM
)
4166 a
= float64_squash_input_denormal(a STATUS_VAR
);
4167 b
= float64_squash_input_denormal(b STATUS_VAR
);
4169 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4170 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4172 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
4173 float_raise( float_flag_invalid STATUS_VAR
);
4177 aSign
= extractFloat64Sign( a
);
4178 bSign
= extractFloat64Sign( b
);
4179 av
= float64_val(a
);
4180 bv
= float64_val(b
);
4181 if ( aSign
!= bSign
) return aSign
&& ( (uint64_t) ( ( av
| bv
)<<1 ) != 0 );
4182 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
4186 /*----------------------------------------------------------------------------
4187 | Returns 1 if the double-precision floating-point values `a' and `b' cannot
4188 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
4189 | comparison is performed according to the IEC/IEEE Standard for Binary
4190 | Floating-Point Arithmetic.
4191 *----------------------------------------------------------------------------*/
4193 int float64_unordered_quiet( float64 a
, float64 b STATUS_PARAM
)
4195 a
= float64_squash_input_denormal(a STATUS_VAR
);
4196 b
= float64_squash_input_denormal(b STATUS_VAR
);
4198 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
4199 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
4201 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
4202 float_raise( float_flag_invalid STATUS_VAR
);
4209 /*----------------------------------------------------------------------------
4210 | Returns the result of converting the extended double-precision floating-
4211 | point value `a' to the 32-bit two's complement integer format. The
4212 | conversion is performed according to the IEC/IEEE Standard for Binary
4213 | Floating-Point Arithmetic---which means in particular that the conversion
4214 | is rounded according to the current rounding mode. If `a' is a NaN, the
4215 | largest positive integer is returned. Otherwise, if the conversion
4216 | overflows, the largest integer with the same sign as `a' is returned.
4217 *----------------------------------------------------------------------------*/
4219 int32
floatx80_to_int32( floatx80 a STATUS_PARAM
)
4222 int32 aExp
, shiftCount
;
4225 aSig
= extractFloatx80Frac( a
);
4226 aExp
= extractFloatx80Exp( a
);
4227 aSign
= extractFloatx80Sign( a
);
4228 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
4229 shiftCount
= 0x4037 - aExp
;
4230 if ( shiftCount
<= 0 ) shiftCount
= 1;
4231 shift64RightJamming( aSig
, shiftCount
, &aSig
);
4232 return roundAndPackInt32( aSign
, aSig STATUS_VAR
);
4236 /*----------------------------------------------------------------------------
4237 | Returns the result of converting the extended double-precision floating-
4238 | point value `a' to the 32-bit two's complement integer format. The
4239 | conversion is performed according to the IEC/IEEE Standard for Binary
4240 | Floating-Point Arithmetic, except that the conversion is always rounded
4241 | toward zero. If `a' is a NaN, the largest positive integer is returned.
4242 | Otherwise, if the conversion overflows, the largest integer with the same
4243 | sign as `a' is returned.
4244 *----------------------------------------------------------------------------*/
4246 int32
floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM
)
4249 int32 aExp
, shiftCount
;
4250 uint64_t aSig
, savedASig
;
4253 aSig
= extractFloatx80Frac( a
);
4254 aExp
= extractFloatx80Exp( a
);
4255 aSign
= extractFloatx80Sign( a
);
4256 if ( 0x401E < aExp
) {
4257 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
4260 else if ( aExp
< 0x3FFF ) {
4261 if ( aExp
|| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4264 shiftCount
= 0x403E - aExp
;
4266 aSig
>>= shiftCount
;
4268 if ( aSign
) z
= - z
;
4269 if ( ( z
< 0 ) ^ aSign
) {
4271 float_raise( float_flag_invalid STATUS_VAR
);
4272 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
4274 if ( ( aSig
<<shiftCount
) != savedASig
) {
4275 STATUS(float_exception_flags
) |= float_flag_inexact
;
4281 /*----------------------------------------------------------------------------
4282 | Returns the result of converting the extended double-precision floating-
4283 | point value `a' to the 64-bit two's complement integer format. The
4284 | conversion is performed according to the IEC/IEEE Standard for Binary
4285 | Floating-Point Arithmetic---which means in particular that the conversion
4286 | is rounded according to the current rounding mode. If `a' is a NaN,
4287 | the largest positive integer is returned. Otherwise, if the conversion
4288 | overflows, the largest integer with the same sign as `a' is returned.
4289 *----------------------------------------------------------------------------*/
4291 int64
floatx80_to_int64( floatx80 a STATUS_PARAM
)
4294 int32 aExp
, shiftCount
;
4295 uint64_t aSig
, aSigExtra
;
4297 aSig
= extractFloatx80Frac( a
);
4298 aExp
= extractFloatx80Exp( a
);
4299 aSign
= extractFloatx80Sign( a
);
4300 shiftCount
= 0x403E - aExp
;
4301 if ( shiftCount
<= 0 ) {
4303 float_raise( float_flag_invalid STATUS_VAR
);
4305 || ( ( aExp
== 0x7FFF )
4306 && ( aSig
!= LIT64( 0x8000000000000000 ) ) )
4308 return LIT64( 0x7FFFFFFFFFFFFFFF );
4310 return (int64_t) LIT64( 0x8000000000000000 );
4315 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
4317 return roundAndPackInt64( aSign
, aSig
, aSigExtra STATUS_VAR
);
4321 /*----------------------------------------------------------------------------
4322 | Returns the result of converting the extended double-precision floating-
4323 | point value `a' to the 64-bit two's complement integer format. The
4324 | conversion is performed according to the IEC/IEEE Standard for Binary
4325 | Floating-Point Arithmetic, except that the conversion is always rounded
4326 | toward zero. If `a' is a NaN, the largest positive integer is returned.
4327 | Otherwise, if the conversion overflows, the largest integer with the same
4328 | sign as `a' is returned.
4329 *----------------------------------------------------------------------------*/
4331 int64
floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM
)
4334 int32 aExp
, shiftCount
;
4338 aSig
= extractFloatx80Frac( a
);
4339 aExp
= extractFloatx80Exp( a
);
4340 aSign
= extractFloatx80Sign( a
);
4341 shiftCount
= aExp
- 0x403E;
4342 if ( 0 <= shiftCount
) {
4343 aSig
&= LIT64( 0x7FFFFFFFFFFFFFFF );
4344 if ( ( a
.high
!= 0xC03E ) || aSig
) {
4345 float_raise( float_flag_invalid STATUS_VAR
);
4346 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && aSig
) ) {
4347 return LIT64( 0x7FFFFFFFFFFFFFFF );
4350 return (int64_t) LIT64( 0x8000000000000000 );
4352 else if ( aExp
< 0x3FFF ) {
4353 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4356 z
= aSig
>>( - shiftCount
);
4357 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
4358 STATUS(float_exception_flags
) |= float_flag_inexact
;
4360 if ( aSign
) z
= - z
;
4365 /*----------------------------------------------------------------------------
4366 | Returns the result of converting the extended double-precision floating-
4367 | point value `a' to the single-precision floating-point format. The
4368 | conversion is performed according to the IEC/IEEE Standard for Binary
4369 | Floating-Point Arithmetic.
4370 *----------------------------------------------------------------------------*/
4372 float32
floatx80_to_float32( floatx80 a STATUS_PARAM
)
4378 aSig
= extractFloatx80Frac( a
);
4379 aExp
= extractFloatx80Exp( a
);
4380 aSign
= extractFloatx80Sign( a
);
4381 if ( aExp
== 0x7FFF ) {
4382 if ( (uint64_t) ( aSig
<<1 ) ) {
4383 return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
4385 return packFloat32( aSign
, 0xFF, 0 );
4387 shift64RightJamming( aSig
, 33, &aSig
);
4388 if ( aExp
|| aSig
) aExp
-= 0x3F81;
4389 return roundAndPackFloat32( aSign
, aExp
, aSig STATUS_VAR
);
4393 /*----------------------------------------------------------------------------
4394 | Returns the result of converting the extended double-precision floating-
4395 | point value `a' to the double-precision floating-point format. The
4396 | conversion is performed according to the IEC/IEEE Standard for Binary
4397 | Floating-Point Arithmetic.
4398 *----------------------------------------------------------------------------*/
4400 float64
floatx80_to_float64( floatx80 a STATUS_PARAM
)
4404 uint64_t aSig
, zSig
;
4406 aSig
= extractFloatx80Frac( a
);
4407 aExp
= extractFloatx80Exp( a
);
4408 aSign
= extractFloatx80Sign( a
);
4409 if ( aExp
== 0x7FFF ) {
4410 if ( (uint64_t) ( aSig
<<1 ) ) {
4411 return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
4413 return packFloat64( aSign
, 0x7FF, 0 );
4415 shift64RightJamming( aSig
, 1, &zSig
);
4416 if ( aExp
|| aSig
) aExp
-= 0x3C01;
4417 return roundAndPackFloat64( aSign
, aExp
, zSig STATUS_VAR
);
4421 /*----------------------------------------------------------------------------
4422 | Returns the result of converting the extended double-precision floating-
4423 | point value `a' to the quadruple-precision floating-point format. The
4424 | conversion is performed according to the IEC/IEEE Standard for Binary
4425 | Floating-Point Arithmetic.
4426 *----------------------------------------------------------------------------*/
4428 float128
floatx80_to_float128( floatx80 a STATUS_PARAM
)
4432 uint64_t aSig
, zSig0
, zSig1
;
4434 aSig
= extractFloatx80Frac( a
);
4435 aExp
= extractFloatx80Exp( a
);
4436 aSign
= extractFloatx80Sign( a
);
4437 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) {
4438 return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
4440 shift128Right( aSig
<<1, 0, 16, &zSig0
, &zSig1
);
4441 return packFloat128( aSign
, aExp
, zSig0
, zSig1
);
4445 /*----------------------------------------------------------------------------
4446 | Rounds the extended double-precision floating-point value `a' to an integer,
4447 | and returns the result as an extended quadruple-precision floating-point
4448 | value. The operation is performed according to the IEC/IEEE Standard for
4449 | Binary Floating-Point Arithmetic.
4450 *----------------------------------------------------------------------------*/
4452 floatx80
floatx80_round_to_int( floatx80 a STATUS_PARAM
)
4456 uint64_t lastBitMask
, roundBitsMask
;
4460 aExp
= extractFloatx80Exp( a
);
4461 if ( 0x403E <= aExp
) {
4462 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) ) {
4463 return propagateFloatx80NaN( a
, a STATUS_VAR
);
4467 if ( aExp
< 0x3FFF ) {
4469 && ( (uint64_t) ( extractFloatx80Frac( a
)<<1 ) == 0 ) ) {
4472 STATUS(float_exception_flags
) |= float_flag_inexact
;
4473 aSign
= extractFloatx80Sign( a
);
4474 switch ( STATUS(float_rounding_mode
) ) {
4475 case float_round_nearest_even
:
4476 if ( ( aExp
== 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 )
4479 packFloatx80( aSign
, 0x3FFF, LIT64( 0x8000000000000000 ) );
4482 case float_round_down
:
4485 packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
4486 : packFloatx80( 0, 0, 0 );
4487 case float_round_up
:
4489 aSign
? packFloatx80( 1, 0, 0 )
4490 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
4492 return packFloatx80( aSign
, 0, 0 );
4495 lastBitMask
<<= 0x403E - aExp
;
4496 roundBitsMask
= lastBitMask
- 1;
4498 roundingMode
= STATUS(float_rounding_mode
);
4499 if ( roundingMode
== float_round_nearest_even
) {
4500 z
.low
+= lastBitMask
>>1;
4501 if ( ( z
.low
& roundBitsMask
) == 0 ) z
.low
&= ~ lastBitMask
;
4503 else if ( roundingMode
!= float_round_to_zero
) {
4504 if ( extractFloatx80Sign( z
) ^ ( roundingMode
== float_round_up
) ) {
4505 z
.low
+= roundBitsMask
;
4508 z
.low
&= ~ roundBitsMask
;
4511 z
.low
= LIT64( 0x8000000000000000 );
4513 if ( z
.low
!= a
.low
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4518 /*----------------------------------------------------------------------------
4519 | Returns the result of adding the absolute values of the extended double-
4520 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
4521 | negated before being returned. `zSign' is ignored if the result is a NaN.
4522 | The addition is performed according to the IEC/IEEE Standard for Binary
4523 | Floating-Point Arithmetic.
4524 *----------------------------------------------------------------------------*/
4526 static floatx80
addFloatx80Sigs( floatx80 a
, floatx80 b
, flag zSign STATUS_PARAM
)
4528 int32 aExp
, bExp
, zExp
;
4529 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4532 aSig
= extractFloatx80Frac( a
);
4533 aExp
= extractFloatx80Exp( a
);
4534 bSig
= extractFloatx80Frac( b
);
4535 bExp
= extractFloatx80Exp( b
);
4536 expDiff
= aExp
- bExp
;
4537 if ( 0 < expDiff
) {
4538 if ( aExp
== 0x7FFF ) {
4539 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4542 if ( bExp
== 0 ) --expDiff
;
4543 shift64ExtraRightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
4546 else if ( expDiff
< 0 ) {
4547 if ( bExp
== 0x7FFF ) {
4548 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4549 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4551 if ( aExp
== 0 ) ++expDiff
;
4552 shift64ExtraRightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
4556 if ( aExp
== 0x7FFF ) {
4557 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
4558 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4563 zSig0
= aSig
+ bSig
;
4565 normalizeFloatx80Subnormal( zSig0
, &zExp
, &zSig0
);
4571 zSig0
= aSig
+ bSig
;
4572 if ( (int64_t) zSig0
< 0 ) goto roundAndPack
;
4574 shift64ExtraRightJamming( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
4575 zSig0
|= LIT64( 0x8000000000000000 );
4579 roundAndPackFloatx80(
4580 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4584 /*----------------------------------------------------------------------------
4585 | Returns the result of subtracting the absolute values of the extended
4586 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
4587 | difference is negated before being returned. `zSign' is ignored if the
4588 | result is a NaN. The subtraction is performed according to the IEC/IEEE
4589 | Standard for Binary Floating-Point Arithmetic.
4590 *----------------------------------------------------------------------------*/
4592 static floatx80
subFloatx80Sigs( floatx80 a
, floatx80 b
, flag zSign STATUS_PARAM
)
4594 int32 aExp
, bExp
, zExp
;
4595 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4599 aSig
= extractFloatx80Frac( a
);
4600 aExp
= extractFloatx80Exp( a
);
4601 bSig
= extractFloatx80Frac( b
);
4602 bExp
= extractFloatx80Exp( b
);
4603 expDiff
= aExp
- bExp
;
4604 if ( 0 < expDiff
) goto aExpBigger
;
4605 if ( expDiff
< 0 ) goto bExpBigger
;
4606 if ( aExp
== 0x7FFF ) {
4607 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
4608 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4610 float_raise( float_flag_invalid STATUS_VAR
);
4611 z
.low
= floatx80_default_nan_low
;
4612 z
.high
= floatx80_default_nan_high
;
4620 if ( bSig
< aSig
) goto aBigger
;
4621 if ( aSig
< bSig
) goto bBigger
;
4622 return packFloatx80( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
4624 if ( bExp
== 0x7FFF ) {
4625 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4626 return packFloatx80( zSign
^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
4628 if ( aExp
== 0 ) ++expDiff
;
4629 shift128RightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
4631 sub128( bSig
, 0, aSig
, zSig1
, &zSig0
, &zSig1
);
4634 goto normalizeRoundAndPack
;
4636 if ( aExp
== 0x7FFF ) {
4637 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4640 if ( bExp
== 0 ) --expDiff
;
4641 shift128RightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
4643 sub128( aSig
, 0, bSig
, zSig1
, &zSig0
, &zSig1
);
4645 normalizeRoundAndPack
:
4647 normalizeRoundAndPackFloatx80(
4648 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4652 /*----------------------------------------------------------------------------
4653 | Returns the result of adding the extended double-precision floating-point
4654 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4655 | Standard for Binary Floating-Point Arithmetic.
4656 *----------------------------------------------------------------------------*/
4658 floatx80
floatx80_add( floatx80 a
, floatx80 b STATUS_PARAM
)
4662 aSign
= extractFloatx80Sign( a
);
4663 bSign
= extractFloatx80Sign( b
);
4664 if ( aSign
== bSign
) {
4665 return addFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4668 return subFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4673 /*----------------------------------------------------------------------------
4674 | Returns the result of subtracting the extended double-precision floating-
4675 | point values `a' and `b'. The operation is performed according to the
4676 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4677 *----------------------------------------------------------------------------*/
4679 floatx80
floatx80_sub( floatx80 a
, floatx80 b STATUS_PARAM
)
4683 aSign
= extractFloatx80Sign( a
);
4684 bSign
= extractFloatx80Sign( b
);
4685 if ( aSign
== bSign
) {
4686 return subFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4689 return addFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4694 /*----------------------------------------------------------------------------
4695 | Returns the result of multiplying the extended double-precision floating-
4696 | point values `a' and `b'. The operation is performed according to the
4697 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4698 *----------------------------------------------------------------------------*/
4700 floatx80
floatx80_mul( floatx80 a
, floatx80 b STATUS_PARAM
)
4702 flag aSign
, bSign
, zSign
;
4703 int32 aExp
, bExp
, zExp
;
4704 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4707 aSig
= extractFloatx80Frac( a
);
4708 aExp
= extractFloatx80Exp( a
);
4709 aSign
= extractFloatx80Sign( a
);
4710 bSig
= extractFloatx80Frac( b
);
4711 bExp
= extractFloatx80Exp( b
);
4712 bSign
= extractFloatx80Sign( b
);
4713 zSign
= aSign
^ bSign
;
4714 if ( aExp
== 0x7FFF ) {
4715 if ( (uint64_t) ( aSig
<<1 )
4716 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
4717 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4719 if ( ( bExp
| bSig
) == 0 ) goto invalid
;
4720 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4722 if ( bExp
== 0x7FFF ) {
4723 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4724 if ( ( aExp
| aSig
) == 0 ) {
4726 float_raise( float_flag_invalid STATUS_VAR
);
4727 z
.low
= floatx80_default_nan_low
;
4728 z
.high
= floatx80_default_nan_high
;
4731 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4734 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4735 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
4738 if ( bSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4739 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4741 zExp
= aExp
+ bExp
- 0x3FFE;
4742 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
4743 if ( 0 < (int64_t) zSig0
) {
4744 shortShift128Left( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
4748 roundAndPackFloatx80(
4749 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4753 /*----------------------------------------------------------------------------
4754 | Returns the result of dividing the extended double-precision floating-point
4755 | value `a' by the corresponding value `b'. The operation is performed
4756 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4757 *----------------------------------------------------------------------------*/
4759 floatx80
floatx80_div( floatx80 a
, floatx80 b STATUS_PARAM
)
4761 flag aSign
, bSign
, zSign
;
4762 int32 aExp
, bExp
, zExp
;
4763 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4764 uint64_t rem0
, rem1
, rem2
, term0
, term1
, term2
;
4767 aSig
= extractFloatx80Frac( a
);
4768 aExp
= extractFloatx80Exp( a
);
4769 aSign
= extractFloatx80Sign( a
);
4770 bSig
= extractFloatx80Frac( b
);
4771 bExp
= extractFloatx80Exp( b
);
4772 bSign
= extractFloatx80Sign( b
);
4773 zSign
= aSign
^ bSign
;
4774 if ( aExp
== 0x7FFF ) {
4775 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4776 if ( bExp
== 0x7FFF ) {
4777 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4780 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4782 if ( bExp
== 0x7FFF ) {
4783 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4784 return packFloatx80( zSign
, 0, 0 );
4788 if ( ( aExp
| aSig
) == 0 ) {
4790 float_raise( float_flag_invalid STATUS_VAR
);
4791 z
.low
= floatx80_default_nan_low
;
4792 z
.high
= floatx80_default_nan_high
;
4795 float_raise( float_flag_divbyzero STATUS_VAR
);
4796 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4798 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4801 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4802 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
4804 zExp
= aExp
- bExp
+ 0x3FFE;
4806 if ( bSig
<= aSig
) {
4807 shift128Right( aSig
, 0, 1, &aSig
, &rem1
);
4810 zSig0
= estimateDiv128To64( aSig
, rem1
, bSig
);
4811 mul64To128( bSig
, zSig0
, &term0
, &term1
);
4812 sub128( aSig
, rem1
, term0
, term1
, &rem0
, &rem1
);
4813 while ( (int64_t) rem0
< 0 ) {
4815 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
4817 zSig1
= estimateDiv128To64( rem1
, 0, bSig
);
4818 if ( (uint64_t) ( zSig1
<<1 ) <= 8 ) {
4819 mul64To128( bSig
, zSig1
, &term1
, &term2
);
4820 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
4821 while ( (int64_t) rem1
< 0 ) {
4823 add128( rem1
, rem2
, 0, bSig
, &rem1
, &rem2
);
4825 zSig1
|= ( ( rem1
| rem2
) != 0 );
4828 roundAndPackFloatx80(
4829 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4833 /*----------------------------------------------------------------------------
4834 | Returns the remainder of the extended double-precision floating-point value
4835 | `a' with respect to the corresponding value `b'. The operation is performed
4836 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4837 *----------------------------------------------------------------------------*/
4839 floatx80
floatx80_rem( floatx80 a
, floatx80 b STATUS_PARAM
)
4842 int32 aExp
, bExp
, expDiff
;
4843 uint64_t aSig0
, aSig1
, bSig
;
4844 uint64_t q
, term0
, term1
, alternateASig0
, alternateASig1
;
4847 aSig0
= extractFloatx80Frac( a
);
4848 aExp
= extractFloatx80Exp( a
);
4849 aSign
= extractFloatx80Sign( a
);
4850 bSig
= extractFloatx80Frac( b
);
4851 bExp
= extractFloatx80Exp( b
);
4852 if ( aExp
== 0x7FFF ) {
4853 if ( (uint64_t) ( aSig0
<<1 )
4854 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
4855 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4859 if ( bExp
== 0x7FFF ) {
4860 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4866 float_raise( float_flag_invalid STATUS_VAR
);
4867 z
.low
= floatx80_default_nan_low
;
4868 z
.high
= floatx80_default_nan_high
;
4871 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4874 if ( (uint64_t) ( aSig0
<<1 ) == 0 ) return a
;
4875 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
4877 bSig
|= LIT64( 0x8000000000000000 );
4879 expDiff
= aExp
- bExp
;
4881 if ( expDiff
< 0 ) {
4882 if ( expDiff
< -1 ) return a
;
4883 shift128Right( aSig0
, 0, 1, &aSig0
, &aSig1
);
4886 q
= ( bSig
<= aSig0
);
4887 if ( q
) aSig0
-= bSig
;
4889 while ( 0 < expDiff
) {
4890 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
4891 q
= ( 2 < q
) ? q
- 2 : 0;
4892 mul64To128( bSig
, q
, &term0
, &term1
);
4893 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4894 shortShift128Left( aSig0
, aSig1
, 62, &aSig0
, &aSig1
);
4898 if ( 0 < expDiff
) {
4899 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
4900 q
= ( 2 < q
) ? q
- 2 : 0;
4902 mul64To128( bSig
, q
<<( 64 - expDiff
), &term0
, &term1
);
4903 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4904 shortShift128Left( 0, bSig
, 64 - expDiff
, &term0
, &term1
);
4905 while ( le128( term0
, term1
, aSig0
, aSig1
) ) {
4907 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4914 sub128( term0
, term1
, aSig0
, aSig1
, &alternateASig0
, &alternateASig1
);
4915 if ( lt128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
4916 || ( eq128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
4919 aSig0
= alternateASig0
;
4920 aSig1
= alternateASig1
;
4924 normalizeRoundAndPackFloatx80(
4925 80, zSign
, bExp
+ expDiff
, aSig0
, aSig1 STATUS_VAR
);
4929 /*----------------------------------------------------------------------------
4930 | Returns the square root of the extended double-precision floating-point
4931 | value `a'. The operation is performed according to the IEC/IEEE Standard
4932 | for Binary Floating-Point Arithmetic.
4933 *----------------------------------------------------------------------------*/
4935 floatx80
floatx80_sqrt( floatx80 a STATUS_PARAM
)
4939 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, doubleZSig0
;
4940 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
4943 aSig0
= extractFloatx80Frac( a
);
4944 aExp
= extractFloatx80Exp( a
);
4945 aSign
= extractFloatx80Sign( a
);
4946 if ( aExp
== 0x7FFF ) {
4947 if ( (uint64_t) ( aSig0
<<1 ) ) return propagateFloatx80NaN( a
, a STATUS_VAR
);
4948 if ( ! aSign
) return a
;
4952 if ( ( aExp
| aSig0
) == 0 ) return a
;
4954 float_raise( float_flag_invalid STATUS_VAR
);
4955 z
.low
= floatx80_default_nan_low
;
4956 z
.high
= floatx80_default_nan_high
;
4960 if ( aSig0
== 0 ) return packFloatx80( 0, 0, 0 );
4961 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
4963 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFF;
4964 zSig0
= estimateSqrt32( aExp
, aSig0
>>32 );
4965 shift128Right( aSig0
, 0, 2 + ( aExp
& 1 ), &aSig0
, &aSig1
);
4966 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
4967 doubleZSig0
= zSig0
<<1;
4968 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
4969 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
4970 while ( (int64_t) rem0
< 0 ) {
4973 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
4975 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
4976 if ( ( zSig1
& LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
4977 if ( zSig1
== 0 ) zSig1
= 1;
4978 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
4979 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
4980 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
4981 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
4982 while ( (int64_t) rem1
< 0 ) {
4984 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
4986 term2
|= doubleZSig0
;
4987 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
4989 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
4991 shortShift128Left( 0, zSig1
, 1, &zSig0
, &zSig1
);
4992 zSig0
|= doubleZSig0
;
4994 roundAndPackFloatx80(
4995 STATUS(floatx80_rounding_precision
), 0, zExp
, zSig0
, zSig1 STATUS_VAR
);
4999 /*----------------------------------------------------------------------------
5000 | Returns 1 if the extended double-precision floating-point value `a' is equal
5001 | to the corresponding value `b', and 0 otherwise. The invalid exception is
5002 | raised if either operand is a NaN. Otherwise, the comparison is performed
5003 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5004 *----------------------------------------------------------------------------*/
5006 int floatx80_eq( floatx80 a
, floatx80 b STATUS_PARAM
)
5009 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5010 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5011 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5012 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5014 float_raise( float_flag_invalid STATUS_VAR
);
5019 && ( ( a
.high
== b
.high
)
5021 && ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
5026 /*----------------------------------------------------------------------------
5027 | Returns 1 if the extended double-precision floating-point value `a' is
5028 | less than or equal to the corresponding value `b', and 0 otherwise. The
5029 | invalid exception is raised if either operand is a NaN. The comparison is
5030 | performed according to the IEC/IEEE Standard for Binary Floating-Point
5032 *----------------------------------------------------------------------------*/
5034 int floatx80_le( floatx80 a
, floatx80 b STATUS_PARAM
)
5038 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5039 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5040 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5041 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5043 float_raise( float_flag_invalid STATUS_VAR
);
5046 aSign
= extractFloatx80Sign( a
);
5047 bSign
= extractFloatx80Sign( b
);
5048 if ( aSign
!= bSign
) {
5051 || ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5055 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
5056 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
5060 /*----------------------------------------------------------------------------
5061 | Returns 1 if the extended double-precision floating-point value `a' is
5062 | less than the corresponding value `b', and 0 otherwise. The invalid
5063 | exception is raised if either operand is a NaN. The comparison is performed
5064 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5065 *----------------------------------------------------------------------------*/
5067 int floatx80_lt( floatx80 a
, floatx80 b STATUS_PARAM
)
5071 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5072 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5073 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5074 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5076 float_raise( float_flag_invalid STATUS_VAR
);
5079 aSign
= extractFloatx80Sign( a
);
5080 bSign
= extractFloatx80Sign( b
);
5081 if ( aSign
!= bSign
) {
5084 && ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5088 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
5089 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
5093 /*----------------------------------------------------------------------------
5094 | Returns 1 if the extended double-precision floating-point values `a' and `b'
5095 | cannot be compared, and 0 otherwise. The invalid exception is raised if
5096 | either operand is a NaN. The comparison is performed according to the
5097 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5098 *----------------------------------------------------------------------------*/
5099 int floatx80_unordered( floatx80 a
, floatx80 b STATUS_PARAM
)
5101 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5102 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5103 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5104 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5106 float_raise( float_flag_invalid STATUS_VAR
);
5112 /*----------------------------------------------------------------------------
5113 | Returns 1 if the extended double-precision floating-point value `a' is
5114 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
5115 | cause an exception. The comparison is performed according to the IEC/IEEE
5116 | Standard for Binary Floating-Point Arithmetic.
5117 *----------------------------------------------------------------------------*/
5119 int floatx80_eq_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
5122 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5123 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5124 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5125 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5127 if ( floatx80_is_signaling_nan( a
)
5128 || floatx80_is_signaling_nan( b
) ) {
5129 float_raise( float_flag_invalid STATUS_VAR
);
5135 && ( ( a
.high
== b
.high
)
5137 && ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
5142 /*----------------------------------------------------------------------------
5143 | Returns 1 if the extended double-precision floating-point value `a' is less
5144 | than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
5145 | do not cause an exception. Otherwise, the comparison is performed according
5146 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5147 *----------------------------------------------------------------------------*/
5149 int floatx80_le_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
5153 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5154 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5155 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5156 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5158 if ( floatx80_is_signaling_nan( a
)
5159 || floatx80_is_signaling_nan( b
) ) {
5160 float_raise( float_flag_invalid STATUS_VAR
);
5164 aSign
= extractFloatx80Sign( a
);
5165 bSign
= extractFloatx80Sign( b
);
5166 if ( aSign
!= bSign
) {
5169 || ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5173 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
5174 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
5178 /*----------------------------------------------------------------------------
5179 | Returns 1 if the extended double-precision floating-point value `a' is less
5180 | than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
5181 | an exception. Otherwise, the comparison is performed according to the
5182 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5183 *----------------------------------------------------------------------------*/
5185 int floatx80_lt_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
5189 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5190 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5191 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5192 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5194 if ( floatx80_is_signaling_nan( a
)
5195 || floatx80_is_signaling_nan( b
) ) {
5196 float_raise( float_flag_invalid STATUS_VAR
);
5200 aSign
= extractFloatx80Sign( a
);
5201 bSign
= extractFloatx80Sign( b
);
5202 if ( aSign
!= bSign
) {
5205 && ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5209 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
5210 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
5214 /*----------------------------------------------------------------------------
5215 | Returns 1 if the extended double-precision floating-point values `a' and `b'
5216 | cannot be compared, and 0 otherwise. Quiet NaNs do not cause an exception.
5217 | The comparison is performed according to the IEC/IEEE Standard for Binary
5218 | Floating-Point Arithmetic.
5219 *----------------------------------------------------------------------------*/
5220 int floatx80_unordered_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
5222 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
5223 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
5224 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
5225 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
5227 if ( floatx80_is_signaling_nan( a
)
5228 || floatx80_is_signaling_nan( b
) ) {
5229 float_raise( float_flag_invalid STATUS_VAR
);
5236 /*----------------------------------------------------------------------------
5237 | Returns the result of converting the quadruple-precision floating-point
5238 | value `a' to the 32-bit two's complement integer format. The conversion
5239 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5240 | Arithmetic---which means in particular that the conversion is rounded
5241 | according to the current rounding mode. If `a' is a NaN, the largest
5242 | positive integer is returned. Otherwise, if the conversion overflows, the
5243 | largest integer with the same sign as `a' is returned.
5244 *----------------------------------------------------------------------------*/
5246 int32
float128_to_int32( float128 a STATUS_PARAM
)
5249 int32 aExp
, shiftCount
;
5250 uint64_t aSig0
, aSig1
;
5252 aSig1
= extractFloat128Frac1( a
);
5253 aSig0
= extractFloat128Frac0( a
);
5254 aExp
= extractFloat128Exp( a
);
5255 aSign
= extractFloat128Sign( a
);
5256 if ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) aSign
= 0;
5257 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
5258 aSig0
|= ( aSig1
!= 0 );
5259 shiftCount
= 0x4028 - aExp
;
5260 if ( 0 < shiftCount
) shift64RightJamming( aSig0
, shiftCount
, &aSig0
);
5261 return roundAndPackInt32( aSign
, aSig0 STATUS_VAR
);
5265 /*----------------------------------------------------------------------------
5266 | Returns the result of converting the quadruple-precision floating-point
5267 | value `a' to the 32-bit two's complement integer format. The conversion
5268 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5269 | Arithmetic, except that the conversion is always rounded toward zero. If
5270 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
5271 | conversion overflows, the largest integer with the same sign as `a' is
5273 *----------------------------------------------------------------------------*/
5275 int32
float128_to_int32_round_to_zero( float128 a STATUS_PARAM
)
5278 int32 aExp
, shiftCount
;
5279 uint64_t aSig0
, aSig1
, savedASig
;
5282 aSig1
= extractFloat128Frac1( a
);
5283 aSig0
= extractFloat128Frac0( a
);
5284 aExp
= extractFloat128Exp( a
);
5285 aSign
= extractFloat128Sign( a
);
5286 aSig0
|= ( aSig1
!= 0 );
5287 if ( 0x401E < aExp
) {
5288 if ( ( aExp
== 0x7FFF ) && aSig0
) aSign
= 0;
5291 else if ( aExp
< 0x3FFF ) {
5292 if ( aExp
|| aSig0
) STATUS(float_exception_flags
) |= float_flag_inexact
;
5295 aSig0
|= LIT64( 0x0001000000000000 );
5296 shiftCount
= 0x402F - aExp
;
5298 aSig0
>>= shiftCount
;
5300 if ( aSign
) z
= - z
;
5301 if ( ( z
< 0 ) ^ aSign
) {
5303 float_raise( float_flag_invalid STATUS_VAR
);
5304 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
5306 if ( ( aSig0
<<shiftCount
) != savedASig
) {
5307 STATUS(float_exception_flags
) |= float_flag_inexact
;
5313 /*----------------------------------------------------------------------------
5314 | Returns the result of converting the quadruple-precision floating-point
5315 | value `a' to the 64-bit two's complement integer format. The conversion
5316 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5317 | Arithmetic---which means in particular that the conversion is rounded
5318 | according to the current rounding mode. If `a' is a NaN, the largest
5319 | positive integer is returned. Otherwise, if the conversion overflows, the
5320 | largest integer with the same sign as `a' is returned.
5321 *----------------------------------------------------------------------------*/
5323 int64
float128_to_int64( float128 a STATUS_PARAM
)
5326 int32 aExp
, shiftCount
;
5327 uint64_t aSig0
, aSig1
;
5329 aSig1
= extractFloat128Frac1( a
);
5330 aSig0
= extractFloat128Frac0( a
);
5331 aExp
= extractFloat128Exp( a
);
5332 aSign
= extractFloat128Sign( a
);
5333 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
5334 shiftCount
= 0x402F - aExp
;
5335 if ( shiftCount
<= 0 ) {
5336 if ( 0x403E < aExp
) {
5337 float_raise( float_flag_invalid STATUS_VAR
);
5339 || ( ( aExp
== 0x7FFF )
5340 && ( aSig1
|| ( aSig0
!= LIT64( 0x0001000000000000 ) ) )
5343 return LIT64( 0x7FFFFFFFFFFFFFFF );
5345 return (int64_t) LIT64( 0x8000000000000000 );
5347 shortShift128Left( aSig0
, aSig1
, - shiftCount
, &aSig0
, &aSig1
);
5350 shift64ExtraRightJamming( aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
5352 return roundAndPackInt64( aSign
, aSig0
, aSig1 STATUS_VAR
);
5356 /*----------------------------------------------------------------------------
5357 | Returns the result of converting the quadruple-precision floating-point
5358 | value `a' to the 64-bit two's complement integer format. The conversion
5359 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5360 | Arithmetic, except that the conversion is always rounded toward zero.
5361 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
5362 | the conversion overflows, the largest integer with the same sign as `a' is
5364 *----------------------------------------------------------------------------*/
5366 int64
float128_to_int64_round_to_zero( float128 a STATUS_PARAM
)
5369 int32 aExp
, shiftCount
;
5370 uint64_t aSig0
, aSig1
;
5373 aSig1
= extractFloat128Frac1( a
);
5374 aSig0
= extractFloat128Frac0( a
);
5375 aExp
= extractFloat128Exp( a
);
5376 aSign
= extractFloat128Sign( a
);
5377 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
5378 shiftCount
= aExp
- 0x402F;
5379 if ( 0 < shiftCount
) {
5380 if ( 0x403E <= aExp
) {
5381 aSig0
&= LIT64( 0x0000FFFFFFFFFFFF );
5382 if ( ( a
.high
== LIT64( 0xC03E000000000000 ) )
5383 && ( aSig1
< LIT64( 0x0002000000000000 ) ) ) {
5384 if ( aSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
5387 float_raise( float_flag_invalid STATUS_VAR
);
5388 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) ) {
5389 return LIT64( 0x7FFFFFFFFFFFFFFF );
5392 return (int64_t) LIT64( 0x8000000000000000 );
5394 z
= ( aSig0
<<shiftCount
) | ( aSig1
>>( ( - shiftCount
) & 63 ) );
5395 if ( (uint64_t) ( aSig1
<<shiftCount
) ) {
5396 STATUS(float_exception_flags
) |= float_flag_inexact
;
5400 if ( aExp
< 0x3FFF ) {
5401 if ( aExp
| aSig0
| aSig1
) {
5402 STATUS(float_exception_flags
) |= float_flag_inexact
;
5406 z
= aSig0
>>( - shiftCount
);
5408 || ( shiftCount
&& (uint64_t) ( aSig0
<<( shiftCount
& 63 ) ) ) ) {
5409 STATUS(float_exception_flags
) |= float_flag_inexact
;
5412 if ( aSign
) z
= - z
;
5417 /*----------------------------------------------------------------------------
5418 | Returns the result of converting the quadruple-precision floating-point
5419 | value `a' to the single-precision floating-point format. The conversion
5420 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5422 *----------------------------------------------------------------------------*/
5424 float32
float128_to_float32( float128 a STATUS_PARAM
)
5428 uint64_t aSig0
, aSig1
;
5431 aSig1
= extractFloat128Frac1( a
);
5432 aSig0
= extractFloat128Frac0( a
);
5433 aExp
= extractFloat128Exp( a
);
5434 aSign
= extractFloat128Sign( a
);
5435 if ( aExp
== 0x7FFF ) {
5436 if ( aSig0
| aSig1
) {
5437 return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5439 return packFloat32( aSign
, 0xFF, 0 );
5441 aSig0
|= ( aSig1
!= 0 );
5442 shift64RightJamming( aSig0
, 18, &aSig0
);
5444 if ( aExp
|| zSig
) {
5448 return roundAndPackFloat32( aSign
, aExp
, zSig STATUS_VAR
);
5452 /*----------------------------------------------------------------------------
5453 | Returns the result of converting the quadruple-precision floating-point
5454 | value `a' to the double-precision floating-point format. The conversion
5455 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5457 *----------------------------------------------------------------------------*/
5459 float64
float128_to_float64( float128 a STATUS_PARAM
)
5463 uint64_t aSig0
, aSig1
;
5465 aSig1
= extractFloat128Frac1( a
);
5466 aSig0
= extractFloat128Frac0( a
);
5467 aExp
= extractFloat128Exp( a
);
5468 aSign
= extractFloat128Sign( a
);
5469 if ( aExp
== 0x7FFF ) {
5470 if ( aSig0
| aSig1
) {
5471 return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5473 return packFloat64( aSign
, 0x7FF, 0 );
5475 shortShift128Left( aSig0
, aSig1
, 14, &aSig0
, &aSig1
);
5476 aSig0
|= ( aSig1
!= 0 );
5477 if ( aExp
|| aSig0
) {
5478 aSig0
|= LIT64( 0x4000000000000000 );
5481 return roundAndPackFloat64( aSign
, aExp
, aSig0 STATUS_VAR
);
5485 /*----------------------------------------------------------------------------
5486 | Returns the result of converting the quadruple-precision floating-point
5487 | value `a' to the extended double-precision floating-point format. The
5488 | conversion is performed according to the IEC/IEEE Standard for Binary
5489 | Floating-Point Arithmetic.
5490 *----------------------------------------------------------------------------*/
5492 floatx80
float128_to_floatx80( float128 a STATUS_PARAM
)
5496 uint64_t aSig0
, aSig1
;
5498 aSig1
= extractFloat128Frac1( a
);
5499 aSig0
= extractFloat128Frac0( a
);
5500 aExp
= extractFloat128Exp( a
);
5501 aSign
= extractFloat128Sign( a
);
5502 if ( aExp
== 0x7FFF ) {
5503 if ( aSig0
| aSig1
) {
5504 return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5506 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
5509 if ( ( aSig0
| aSig1
) == 0 ) return packFloatx80( aSign
, 0, 0 );
5510 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5513 aSig0
|= LIT64( 0x0001000000000000 );
5515 shortShift128Left( aSig0
, aSig1
, 15, &aSig0
, &aSig1
);
5516 return roundAndPackFloatx80( 80, aSign
, aExp
, aSig0
, aSig1 STATUS_VAR
);
5520 /*----------------------------------------------------------------------------
5521 | Rounds the quadruple-precision floating-point value `a' to an integer, and
5522 | returns the result as a quadruple-precision floating-point value. The
5523 | operation is performed according to the IEC/IEEE Standard for Binary
5524 | Floating-Point Arithmetic.
5525 *----------------------------------------------------------------------------*/
5527 float128
float128_round_to_int( float128 a STATUS_PARAM
)
5531 uint64_t lastBitMask
, roundBitsMask
;
5535 aExp
= extractFloat128Exp( a
);
5536 if ( 0x402F <= aExp
) {
5537 if ( 0x406F <= aExp
) {
5538 if ( ( aExp
== 0x7FFF )
5539 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) )
5541 return propagateFloat128NaN( a
, a STATUS_VAR
);
5546 lastBitMask
= ( lastBitMask
<<( 0x406E - aExp
) )<<1;
5547 roundBitsMask
= lastBitMask
- 1;
5549 roundingMode
= STATUS(float_rounding_mode
);
5550 if ( roundingMode
== float_round_nearest_even
) {
5551 if ( lastBitMask
) {
5552 add128( z
.high
, z
.low
, 0, lastBitMask
>>1, &z
.high
, &z
.low
);
5553 if ( ( z
.low
& roundBitsMask
) == 0 ) z
.low
&= ~ lastBitMask
;
5556 if ( (int64_t) z
.low
< 0 ) {
5558 if ( (uint64_t) ( z
.low
<<1 ) == 0 ) z
.high
&= ~1;
5562 else if ( roundingMode
!= float_round_to_zero
) {
5563 if ( extractFloat128Sign( z
)
5564 ^ ( roundingMode
== float_round_up
) ) {
5565 add128( z
.high
, z
.low
, 0, roundBitsMask
, &z
.high
, &z
.low
);
5568 z
.low
&= ~ roundBitsMask
;
5571 if ( aExp
< 0x3FFF ) {
5572 if ( ( ( (uint64_t) ( a
.high
<<1 ) ) | a
.low
) == 0 ) return a
;
5573 STATUS(float_exception_flags
) |= float_flag_inexact
;
5574 aSign
= extractFloat128Sign( a
);
5575 switch ( STATUS(float_rounding_mode
) ) {
5576 case float_round_nearest_even
:
5577 if ( ( aExp
== 0x3FFE )
5578 && ( extractFloat128Frac0( a
)
5579 | extractFloat128Frac1( a
) )
5581 return packFloat128( aSign
, 0x3FFF, 0, 0 );
5584 case float_round_down
:
5586 aSign
? packFloat128( 1, 0x3FFF, 0, 0 )
5587 : packFloat128( 0, 0, 0, 0 );
5588 case float_round_up
:
5590 aSign
? packFloat128( 1, 0, 0, 0 )
5591 : packFloat128( 0, 0x3FFF, 0, 0 );
5593 return packFloat128( aSign
, 0, 0, 0 );
5596 lastBitMask
<<= 0x402F - aExp
;
5597 roundBitsMask
= lastBitMask
- 1;
5600 roundingMode
= STATUS(float_rounding_mode
);
5601 if ( roundingMode
== float_round_nearest_even
) {
5602 z
.high
+= lastBitMask
>>1;
5603 if ( ( ( z
.high
& roundBitsMask
) | a
.low
) == 0 ) {
5604 z
.high
&= ~ lastBitMask
;
5607 else if ( roundingMode
!= float_round_to_zero
) {
5608 if ( extractFloat128Sign( z
)
5609 ^ ( roundingMode
== float_round_up
) ) {
5610 z
.high
|= ( a
.low
!= 0 );
5611 z
.high
+= roundBitsMask
;
5614 z
.high
&= ~ roundBitsMask
;
5616 if ( ( z
.low
!= a
.low
) || ( z
.high
!= a
.high
) ) {
5617 STATUS(float_exception_flags
) |= float_flag_inexact
;
5623 /*----------------------------------------------------------------------------
5624 | Returns the result of adding the absolute values of the quadruple-precision
5625 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
5626 | before being returned. `zSign' is ignored if the result is a NaN.
5627 | The addition is performed according to the IEC/IEEE Standard for Binary
5628 | Floating-Point Arithmetic.
5629 *----------------------------------------------------------------------------*/
5631 static float128
addFloat128Sigs( float128 a
, float128 b
, flag zSign STATUS_PARAM
)
5633 int32 aExp
, bExp
, zExp
;
5634 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
;
5637 aSig1
= extractFloat128Frac1( a
);
5638 aSig0
= extractFloat128Frac0( a
);
5639 aExp
= extractFloat128Exp( a
);
5640 bSig1
= extractFloat128Frac1( b
);
5641 bSig0
= extractFloat128Frac0( b
);
5642 bExp
= extractFloat128Exp( b
);
5643 expDiff
= aExp
- bExp
;
5644 if ( 0 < expDiff
) {
5645 if ( aExp
== 0x7FFF ) {
5646 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5653 bSig0
|= LIT64( 0x0001000000000000 );
5655 shift128ExtraRightJamming(
5656 bSig0
, bSig1
, 0, expDiff
, &bSig0
, &bSig1
, &zSig2
);
5659 else if ( expDiff
< 0 ) {
5660 if ( bExp
== 0x7FFF ) {
5661 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5662 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5668 aSig0
|= LIT64( 0x0001000000000000 );
5670 shift128ExtraRightJamming(
5671 aSig0
, aSig1
, 0, - expDiff
, &aSig0
, &aSig1
, &zSig2
);
5675 if ( aExp
== 0x7FFF ) {
5676 if ( aSig0
| aSig1
| bSig0
| bSig1
) {
5677 return propagateFloat128NaN( a
, b STATUS_VAR
);
5681 add128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5683 if (STATUS(flush_to_zero
)) {
5684 if (zSig0
| zSig1
) {
5685 float_raise(float_flag_output_denormal STATUS_VAR
);
5687 return packFloat128(zSign
, 0, 0, 0);
5689 return packFloat128( zSign
, 0, zSig0
, zSig1
);
5692 zSig0
|= LIT64( 0x0002000000000000 );
5696 aSig0
|= LIT64( 0x0001000000000000 );
5697 add128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5699 if ( zSig0
< LIT64( 0x0002000000000000 ) ) goto roundAndPack
;
5702 shift128ExtraRightJamming(
5703 zSig0
, zSig1
, zSig2
, 1, &zSig0
, &zSig1
, &zSig2
);
5705 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5709 /*----------------------------------------------------------------------------
5710 | Returns the result of subtracting the absolute values of the quadruple-
5711 | precision floating-point values `a' and `b'. If `zSign' is 1, the
5712 | difference is negated before being returned. `zSign' is ignored if the
5713 | result is a NaN. The subtraction is performed according to the IEC/IEEE
5714 | Standard for Binary Floating-Point Arithmetic.
5715 *----------------------------------------------------------------------------*/
5717 static float128
subFloat128Sigs( float128 a
, float128 b
, flag zSign STATUS_PARAM
)
5719 int32 aExp
, bExp
, zExp
;
5720 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
;
5724 aSig1
= extractFloat128Frac1( a
);
5725 aSig0
= extractFloat128Frac0( a
);
5726 aExp
= extractFloat128Exp( a
);
5727 bSig1
= extractFloat128Frac1( b
);
5728 bSig0
= extractFloat128Frac0( b
);
5729 bExp
= extractFloat128Exp( b
);
5730 expDiff
= aExp
- bExp
;
5731 shortShift128Left( aSig0
, aSig1
, 14, &aSig0
, &aSig1
);
5732 shortShift128Left( bSig0
, bSig1
, 14, &bSig0
, &bSig1
);
5733 if ( 0 < expDiff
) goto aExpBigger
;
5734 if ( expDiff
< 0 ) goto bExpBigger
;
5735 if ( aExp
== 0x7FFF ) {
5736 if ( aSig0
| aSig1
| bSig0
| bSig1
) {
5737 return propagateFloat128NaN( a
, b STATUS_VAR
);
5739 float_raise( float_flag_invalid STATUS_VAR
);
5740 z
.low
= float128_default_nan_low
;
5741 z
.high
= float128_default_nan_high
;
5748 if ( bSig0
< aSig0
) goto aBigger
;
5749 if ( aSig0
< bSig0
) goto bBigger
;
5750 if ( bSig1
< aSig1
) goto aBigger
;
5751 if ( aSig1
< bSig1
) goto bBigger
;
5752 return packFloat128( STATUS(float_rounding_mode
) == float_round_down
, 0, 0, 0 );
5754 if ( bExp
== 0x7FFF ) {
5755 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5756 return packFloat128( zSign
^ 1, 0x7FFF, 0, 0 );
5762 aSig0
|= LIT64( 0x4000000000000000 );
5764 shift128RightJamming( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
5765 bSig0
|= LIT64( 0x4000000000000000 );
5767 sub128( bSig0
, bSig1
, aSig0
, aSig1
, &zSig0
, &zSig1
);
5770 goto normalizeRoundAndPack
;
5772 if ( aExp
== 0x7FFF ) {
5773 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5780 bSig0
|= LIT64( 0x4000000000000000 );
5782 shift128RightJamming( bSig0
, bSig1
, expDiff
, &bSig0
, &bSig1
);
5783 aSig0
|= LIT64( 0x4000000000000000 );
5785 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5787 normalizeRoundAndPack
:
5789 return normalizeRoundAndPackFloat128( zSign
, zExp
- 14, zSig0
, zSig1 STATUS_VAR
);
5793 /*----------------------------------------------------------------------------
5794 | Returns the result of adding the quadruple-precision floating-point values
5795 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
5796 | for Binary Floating-Point Arithmetic.
5797 *----------------------------------------------------------------------------*/
5799 float128
float128_add( float128 a
, float128 b STATUS_PARAM
)
5803 aSign
= extractFloat128Sign( a
);
5804 bSign
= extractFloat128Sign( b
);
5805 if ( aSign
== bSign
) {
5806 return addFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5809 return subFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5814 /*----------------------------------------------------------------------------
5815 | Returns the result of subtracting the quadruple-precision floating-point
5816 | values `a' and `b'. The operation is performed according to the IEC/IEEE
5817 | Standard for Binary Floating-Point Arithmetic.
5818 *----------------------------------------------------------------------------*/
5820 float128
float128_sub( float128 a
, float128 b STATUS_PARAM
)
5824 aSign
= extractFloat128Sign( a
);
5825 bSign
= extractFloat128Sign( b
);
5826 if ( aSign
== bSign
) {
5827 return subFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5830 return addFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5835 /*----------------------------------------------------------------------------
5836 | Returns the result of multiplying the quadruple-precision floating-point
5837 | values `a' and `b'. The operation is performed according to the IEC/IEEE
5838 | Standard for Binary Floating-Point Arithmetic.
5839 *----------------------------------------------------------------------------*/
5841 float128
float128_mul( float128 a
, float128 b STATUS_PARAM
)
5843 flag aSign
, bSign
, zSign
;
5844 int32 aExp
, bExp
, zExp
;
5845 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
, zSig3
;
5848 aSig1
= extractFloat128Frac1( a
);
5849 aSig0
= extractFloat128Frac0( a
);
5850 aExp
= extractFloat128Exp( a
);
5851 aSign
= extractFloat128Sign( a
);
5852 bSig1
= extractFloat128Frac1( b
);
5853 bSig0
= extractFloat128Frac0( b
);
5854 bExp
= extractFloat128Exp( b
);
5855 bSign
= extractFloat128Sign( b
);
5856 zSign
= aSign
^ bSign
;
5857 if ( aExp
== 0x7FFF ) {
5858 if ( ( aSig0
| aSig1
)
5859 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
5860 return propagateFloat128NaN( a
, b STATUS_VAR
);
5862 if ( ( bExp
| bSig0
| bSig1
) == 0 ) goto invalid
;
5863 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5865 if ( bExp
== 0x7FFF ) {
5866 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5867 if ( ( aExp
| aSig0
| aSig1
) == 0 ) {
5869 float_raise( float_flag_invalid STATUS_VAR
);
5870 z
.low
= float128_default_nan_low
;
5871 z
.high
= float128_default_nan_high
;
5874 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5877 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5878 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5881 if ( ( bSig0
| bSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5882 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5884 zExp
= aExp
+ bExp
- 0x4000;
5885 aSig0
|= LIT64( 0x0001000000000000 );
5886 shortShift128Left( bSig0
, bSig1
, 16, &bSig0
, &bSig1
);
5887 mul128To256( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
, &zSig2
, &zSig3
);
5888 add128( zSig0
, zSig1
, aSig0
, aSig1
, &zSig0
, &zSig1
);
5889 zSig2
|= ( zSig3
!= 0 );
5890 if ( LIT64( 0x0002000000000000 ) <= zSig0
) {
5891 shift128ExtraRightJamming(
5892 zSig0
, zSig1
, zSig2
, 1, &zSig0
, &zSig1
, &zSig2
);
5895 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5899 /*----------------------------------------------------------------------------
5900 | Returns the result of dividing the quadruple-precision floating-point value
5901 | `a' by the corresponding value `b'. The operation is performed according to
5902 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5903 *----------------------------------------------------------------------------*/
5905 float128
float128_div( float128 a
, float128 b STATUS_PARAM
)
5907 flag aSign
, bSign
, zSign
;
5908 int32 aExp
, bExp
, zExp
;
5909 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
;
5910 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
5913 aSig1
= extractFloat128Frac1( a
);
5914 aSig0
= extractFloat128Frac0( a
);
5915 aExp
= extractFloat128Exp( a
);
5916 aSign
= extractFloat128Sign( a
);
5917 bSig1
= extractFloat128Frac1( b
);
5918 bSig0
= extractFloat128Frac0( b
);
5919 bExp
= extractFloat128Exp( b
);
5920 bSign
= extractFloat128Sign( b
);
5921 zSign
= aSign
^ bSign
;
5922 if ( aExp
== 0x7FFF ) {
5923 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5924 if ( bExp
== 0x7FFF ) {
5925 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5928 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5930 if ( bExp
== 0x7FFF ) {
5931 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5932 return packFloat128( zSign
, 0, 0, 0 );
5935 if ( ( bSig0
| bSig1
) == 0 ) {
5936 if ( ( aExp
| aSig0
| aSig1
) == 0 ) {
5938 float_raise( float_flag_invalid STATUS_VAR
);
5939 z
.low
= float128_default_nan_low
;
5940 z
.high
= float128_default_nan_high
;
5943 float_raise( float_flag_divbyzero STATUS_VAR
);
5944 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5946 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5949 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5950 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5952 zExp
= aExp
- bExp
+ 0x3FFD;
5954 aSig0
| LIT64( 0x0001000000000000 ), aSig1
, 15, &aSig0
, &aSig1
);
5956 bSig0
| LIT64( 0x0001000000000000 ), bSig1
, 15, &bSig0
, &bSig1
);
5957 if ( le128( bSig0
, bSig1
, aSig0
, aSig1
) ) {
5958 shift128Right( aSig0
, aSig1
, 1, &aSig0
, &aSig1
);
5961 zSig0
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5962 mul128By64To192( bSig0
, bSig1
, zSig0
, &term0
, &term1
, &term2
);
5963 sub192( aSig0
, aSig1
, 0, term0
, term1
, term2
, &rem0
, &rem1
, &rem2
);
5964 while ( (int64_t) rem0
< 0 ) {
5966 add192( rem0
, rem1
, rem2
, 0, bSig0
, bSig1
, &rem0
, &rem1
, &rem2
);
5968 zSig1
= estimateDiv128To64( rem1
, rem2
, bSig0
);
5969 if ( ( zSig1
& 0x3FFF ) <= 4 ) {
5970 mul128By64To192( bSig0
, bSig1
, zSig1
, &term1
, &term2
, &term3
);
5971 sub192( rem1
, rem2
, 0, term1
, term2
, term3
, &rem1
, &rem2
, &rem3
);
5972 while ( (int64_t) rem1
< 0 ) {
5974 add192( rem1
, rem2
, rem3
, 0, bSig0
, bSig1
, &rem1
, &rem2
, &rem3
);
5976 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
5978 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 15, &zSig0
, &zSig1
, &zSig2
);
5979 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5983 /*----------------------------------------------------------------------------
5984 | Returns the remainder of the quadruple-precision floating-point value `a'
5985 | with respect to the corresponding value `b'. The operation is performed
5986 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5987 *----------------------------------------------------------------------------*/
5989 float128
float128_rem( float128 a
, float128 b STATUS_PARAM
)
5992 int32 aExp
, bExp
, expDiff
;
5993 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, q
, term0
, term1
, term2
;
5994 uint64_t allZero
, alternateASig0
, alternateASig1
, sigMean1
;
5998 aSig1
= extractFloat128Frac1( a
);
5999 aSig0
= extractFloat128Frac0( a
);
6000 aExp
= extractFloat128Exp( a
);
6001 aSign
= extractFloat128Sign( a
);
6002 bSig1
= extractFloat128Frac1( b
);
6003 bSig0
= extractFloat128Frac0( b
);
6004 bExp
= extractFloat128Exp( b
);
6005 if ( aExp
== 0x7FFF ) {
6006 if ( ( aSig0
| aSig1
)
6007 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
6008 return propagateFloat128NaN( a
, b STATUS_VAR
);
6012 if ( bExp
== 0x7FFF ) {
6013 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
6017 if ( ( bSig0
| bSig1
) == 0 ) {
6019 float_raise( float_flag_invalid STATUS_VAR
);
6020 z
.low
= float128_default_nan_low
;
6021 z
.high
= float128_default_nan_high
;
6024 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
6027 if ( ( aSig0
| aSig1
) == 0 ) return a
;
6028 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
6030 expDiff
= aExp
- bExp
;
6031 if ( expDiff
< -1 ) return a
;
6033 aSig0
| LIT64( 0x0001000000000000 ),
6035 15 - ( expDiff
< 0 ),
6040 bSig0
| LIT64( 0x0001000000000000 ), bSig1
, 15, &bSig0
, &bSig1
);
6041 q
= le128( bSig0
, bSig1
, aSig0
, aSig1
);
6042 if ( q
) sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
6044 while ( 0 < expDiff
) {
6045 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
6046 q
= ( 4 < q
) ? q
- 4 : 0;
6047 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
6048 shortShift192Left( term0
, term1
, term2
, 61, &term1
, &term2
, &allZero
);
6049 shortShift128Left( aSig0
, aSig1
, 61, &aSig0
, &allZero
);
6050 sub128( aSig0
, 0, term1
, term2
, &aSig0
, &aSig1
);
6053 if ( -64 < expDiff
) {
6054 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
6055 q
= ( 4 < q
) ? q
- 4 : 0;
6057 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
6059 if ( expDiff
< 0 ) {
6060 shift128Right( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
6063 shortShift128Left( aSig0
, aSig1
, expDiff
, &aSig0
, &aSig1
);
6065 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
6066 sub128( aSig0
, aSig1
, term1
, term2
, &aSig0
, &aSig1
);
6069 shift128Right( aSig0
, aSig1
, 12, &aSig0
, &aSig1
);
6070 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
6073 alternateASig0
= aSig0
;
6074 alternateASig1
= aSig1
;
6076 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
6077 } while ( 0 <= (int64_t) aSig0
);
6079 aSig0
, aSig1
, alternateASig0
, alternateASig1
, (uint64_t *)&sigMean0
, &sigMean1
);
6080 if ( ( sigMean0
< 0 )
6081 || ( ( ( sigMean0
| sigMean1
) == 0 ) && ( q
& 1 ) ) ) {
6082 aSig0
= alternateASig0
;
6083 aSig1
= alternateASig1
;
6085 zSign
= ( (int64_t) aSig0
< 0 );
6086 if ( zSign
) sub128( 0, 0, aSig0
, aSig1
, &aSig0
, &aSig1
);
6088 normalizeRoundAndPackFloat128( aSign
^ zSign
, bExp
- 4, aSig0
, aSig1 STATUS_VAR
);
6092 /*----------------------------------------------------------------------------
6093 | Returns the square root of the quadruple-precision floating-point value `a'.
6094 | The operation is performed according to the IEC/IEEE Standard for Binary
6095 | Floating-Point Arithmetic.
6096 *----------------------------------------------------------------------------*/
6098 float128
float128_sqrt( float128 a STATUS_PARAM
)
6102 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, zSig2
, doubleZSig0
;
6103 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
6106 aSig1
= extractFloat128Frac1( a
);
6107 aSig0
= extractFloat128Frac0( a
);
6108 aExp
= extractFloat128Exp( a
);
6109 aSign
= extractFloat128Sign( a
);
6110 if ( aExp
== 0x7FFF ) {
6111 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, a STATUS_VAR
);
6112 if ( ! aSign
) return a
;
6116 if ( ( aExp
| aSig0
| aSig1
) == 0 ) return a
;
6118 float_raise( float_flag_invalid STATUS_VAR
);
6119 z
.low
= float128_default_nan_low
;
6120 z
.high
= float128_default_nan_high
;
6124 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( 0, 0, 0, 0 );
6125 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
6127 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFE;
6128 aSig0
|= LIT64( 0x0001000000000000 );
6129 zSig0
= estimateSqrt32( aExp
, aSig0
>>17 );
6130 shortShift128Left( aSig0
, aSig1
, 13 - ( aExp
& 1 ), &aSig0
, &aSig1
);
6131 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
6132 doubleZSig0
= zSig0
<<1;
6133 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
6134 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
6135 while ( (int64_t) rem0
< 0 ) {
6138 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
6140 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
6141 if ( ( zSig1
& 0x1FFF ) <= 5 ) {
6142 if ( zSig1
== 0 ) zSig1
= 1;
6143 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
6144 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
6145 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
6146 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6147 while ( (int64_t) rem1
< 0 ) {
6149 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
6151 term2
|= doubleZSig0
;
6152 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
6154 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
6156 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 14, &zSig0
, &zSig1
, &zSig2
);
6157 return roundAndPackFloat128( 0, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
6161 /*----------------------------------------------------------------------------
6162 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
6163 | the corresponding value `b', and 0 otherwise. The invalid exception is
6164 | raised if either operand is a NaN. Otherwise, the comparison is performed
6165 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6166 *----------------------------------------------------------------------------*/
6168 int float128_eq( float128 a
, float128 b STATUS_PARAM
)
6171 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6172 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6173 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6174 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6176 float_raise( float_flag_invalid STATUS_VAR
);
6181 && ( ( a
.high
== b
.high
)
6183 && ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
6188 /*----------------------------------------------------------------------------
6189 | Returns 1 if the quadruple-precision floating-point value `a' is less than
6190 | or equal to the corresponding value `b', and 0 otherwise. The invalid
6191 | exception is raised if either operand is a NaN. The comparison is performed
6192 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6193 *----------------------------------------------------------------------------*/
6195 int float128_le( float128 a
, float128 b STATUS_PARAM
)
6199 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6200 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6201 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6202 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6204 float_raise( float_flag_invalid STATUS_VAR
);
6207 aSign
= extractFloat128Sign( a
);
6208 bSign
= extractFloat128Sign( b
);
6209 if ( aSign
!= bSign
) {
6212 || ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
6216 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
6217 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
6221 /*----------------------------------------------------------------------------
6222 | Returns 1 if the quadruple-precision floating-point value `a' is less than
6223 | the corresponding value `b', and 0 otherwise. The invalid exception is
6224 | raised if either operand is a NaN. The comparison is performed according
6225 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6226 *----------------------------------------------------------------------------*/
6228 int float128_lt( float128 a
, float128 b STATUS_PARAM
)
6232 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6233 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6234 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6235 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6237 float_raise( float_flag_invalid STATUS_VAR
);
6240 aSign
= extractFloat128Sign( a
);
6241 bSign
= extractFloat128Sign( b
);
6242 if ( aSign
!= bSign
) {
6245 && ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
6249 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
6250 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
6254 /*----------------------------------------------------------------------------
6255 | Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
6256 | be compared, and 0 otherwise. The invalid exception is raised if either
6257 | operand is a NaN. The comparison is performed according to the IEC/IEEE
6258 | Standard for Binary Floating-Point Arithmetic.
6259 *----------------------------------------------------------------------------*/
6261 int float128_unordered( float128 a
, float128 b STATUS_PARAM
)
6263 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6264 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6265 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6266 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6268 float_raise( float_flag_invalid STATUS_VAR
);
6274 /*----------------------------------------------------------------------------
6275 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
6276 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
6277 | exception. The comparison is performed according to the IEC/IEEE Standard
6278 | for Binary Floating-Point Arithmetic.
6279 *----------------------------------------------------------------------------*/
6281 int float128_eq_quiet( float128 a
, float128 b STATUS_PARAM
)
6284 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6285 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6286 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6287 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6289 if ( float128_is_signaling_nan( a
)
6290 || float128_is_signaling_nan( b
) ) {
6291 float_raise( float_flag_invalid STATUS_VAR
);
6297 && ( ( a
.high
== b
.high
)
6299 && ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
6304 /*----------------------------------------------------------------------------
6305 | Returns 1 if the quadruple-precision floating-point value `a' is less than
6306 | or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
6307 | cause an exception. Otherwise, the comparison is performed according to the
6308 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
6309 *----------------------------------------------------------------------------*/
6311 int float128_le_quiet( float128 a
, float128 b STATUS_PARAM
)
6315 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6316 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6317 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6318 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6320 if ( float128_is_signaling_nan( a
)
6321 || float128_is_signaling_nan( b
) ) {
6322 float_raise( float_flag_invalid STATUS_VAR
);
6326 aSign
= extractFloat128Sign( a
);
6327 bSign
= extractFloat128Sign( b
);
6328 if ( aSign
!= bSign
) {
6331 || ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
6335 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
6336 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
6340 /*----------------------------------------------------------------------------
6341 | Returns 1 if the quadruple-precision floating-point value `a' is less than
6342 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
6343 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
6344 | Standard for Binary Floating-Point Arithmetic.
6345 *----------------------------------------------------------------------------*/
6347 int float128_lt_quiet( float128 a
, float128 b STATUS_PARAM
)
6351 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6352 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6353 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6354 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6356 if ( float128_is_signaling_nan( a
)
6357 || float128_is_signaling_nan( b
) ) {
6358 float_raise( float_flag_invalid STATUS_VAR
);
6362 aSign
= extractFloat128Sign( a
);
6363 bSign
= extractFloat128Sign( b
);
6364 if ( aSign
!= bSign
) {
6367 && ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
6371 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
6372 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
6376 /*----------------------------------------------------------------------------
6377 | Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
6378 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
6379 | comparison is performed according to the IEC/IEEE Standard for Binary
6380 | Floating-Point Arithmetic.
6381 *----------------------------------------------------------------------------*/
6383 int float128_unordered_quiet( float128 a
, float128 b STATUS_PARAM
)
6385 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
6386 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
6387 || ( ( extractFloat128Exp( b
) == 0x7FFF )
6388 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
6390 if ( float128_is_signaling_nan( a
)
6391 || float128_is_signaling_nan( b
) ) {
6392 float_raise( float_flag_invalid STATUS_VAR
);
6399 /* misc functions */
6400 float32
uint32_to_float32( uint32 a STATUS_PARAM
)
6402 return int64_to_float32(a STATUS_VAR
);
6405 float64
uint32_to_float64( uint32 a STATUS_PARAM
)
6407 return int64_to_float64(a STATUS_VAR
);
6410 uint32
float32_to_uint32( float32 a STATUS_PARAM
)
6415 v
= float32_to_int64(a STATUS_VAR
);
6418 float_raise( float_flag_invalid STATUS_VAR
);
6419 } else if (v
> 0xffffffff) {
6421 float_raise( float_flag_invalid STATUS_VAR
);
6428 uint32
float32_to_uint32_round_to_zero( float32 a STATUS_PARAM
)
6433 v
= float32_to_int64_round_to_zero(a STATUS_VAR
);
6436 float_raise( float_flag_invalid STATUS_VAR
);
6437 } else if (v
> 0xffffffff) {
6439 float_raise( float_flag_invalid STATUS_VAR
);
6446 uint16
float32_to_uint16_round_to_zero( float32 a STATUS_PARAM
)
6451 v
= float32_to_int64_round_to_zero(a STATUS_VAR
);
6454 float_raise( float_flag_invalid STATUS_VAR
);
6455 } else if (v
> 0xffff) {
6457 float_raise( float_flag_invalid STATUS_VAR
);
6464 uint32
float64_to_uint32( float64 a STATUS_PARAM
)
6469 v
= float64_to_int64(a STATUS_VAR
);
6472 float_raise( float_flag_invalid STATUS_VAR
);
6473 } else if (v
> 0xffffffff) {
6475 float_raise( float_flag_invalid STATUS_VAR
);
6482 uint32
float64_to_uint32_round_to_zero( float64 a STATUS_PARAM
)
6487 v
= float64_to_int64_round_to_zero(a STATUS_VAR
);
6490 float_raise( float_flag_invalid STATUS_VAR
);
6491 } else if (v
> 0xffffffff) {
6493 float_raise( float_flag_invalid STATUS_VAR
);
6500 uint16
float64_to_uint16_round_to_zero( float64 a STATUS_PARAM
)
6505 v
= float64_to_int64_round_to_zero(a STATUS_VAR
);
6508 float_raise( float_flag_invalid STATUS_VAR
);
6509 } else if (v
> 0xffff) {
6511 float_raise( float_flag_invalid STATUS_VAR
);
6518 /* FIXME: This looks broken. */
6519 uint64_t float64_to_uint64 (float64 a STATUS_PARAM
)
6523 v
= float64_val(int64_to_float64(INT64_MIN STATUS_VAR
));
6524 v
+= float64_val(a
);
6525 v
= float64_to_int64(make_float64(v
) STATUS_VAR
);
6527 return v
- INT64_MIN
;
6530 uint64_t float64_to_uint64_round_to_zero (float64 a STATUS_PARAM
)
6534 v
= float64_val(int64_to_float64(INT64_MIN STATUS_VAR
));
6535 v
+= float64_val(a
);
6536 v
= float64_to_int64_round_to_zero(make_float64(v
) STATUS_VAR
);
6538 return v
- INT64_MIN
;
6541 #define COMPARE(s, nan_exp) \
6542 INLINE int float ## s ## _compare_internal( float ## s a, float ## s b, \
6543 int is_quiet STATUS_PARAM ) \
6545 flag aSign, bSign; \
6546 uint ## s ## _t av, bv; \
6547 a = float ## s ## _squash_input_denormal(a STATUS_VAR); \
6548 b = float ## s ## _squash_input_denormal(b STATUS_VAR); \
6550 if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) && \
6551 extractFloat ## s ## Frac( a ) ) || \
6552 ( ( extractFloat ## s ## Exp( b ) == nan_exp ) && \
6553 extractFloat ## s ## Frac( b ) )) { \
6555 float ## s ## _is_signaling_nan( a ) || \
6556 float ## s ## _is_signaling_nan( b ) ) { \
6557 float_raise( float_flag_invalid STATUS_VAR); \
6559 return float_relation_unordered; \
6561 aSign = extractFloat ## s ## Sign( a ); \
6562 bSign = extractFloat ## s ## Sign( b ); \
6563 av = float ## s ## _val(a); \
6564 bv = float ## s ## _val(b); \
6565 if ( aSign != bSign ) { \
6566 if ( (uint ## s ## _t) ( ( av | bv )<<1 ) == 0 ) { \
6568 return float_relation_equal; \
6570 return 1 - (2 * aSign); \
6574 return float_relation_equal; \
6576 return 1 - 2 * (aSign ^ ( av < bv )); \
6581 int float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM ) \
6583 return float ## s ## _compare_internal(a, b, 0 STATUS_VAR); \
6586 int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
6588 return float ## s ## _compare_internal(a, b, 1 STATUS_VAR); \
6594 INLINE
int floatx80_compare_internal( floatx80 a
, floatx80 b
,
6595 int is_quiet STATUS_PARAM
)
6599 if (( ( extractFloatx80Exp( a
) == 0x7fff ) &&
6600 ( extractFloatx80Frac( a
)<<1 ) ) ||
6601 ( ( extractFloatx80Exp( b
) == 0x7fff ) &&
6602 ( extractFloatx80Frac( b
)<<1 ) )) {
6604 floatx80_is_signaling_nan( a
) ||
6605 floatx80_is_signaling_nan( b
) ) {
6606 float_raise( float_flag_invalid STATUS_VAR
);
6608 return float_relation_unordered
;
6610 aSign
= extractFloatx80Sign( a
);
6611 bSign
= extractFloatx80Sign( b
);
6612 if ( aSign
!= bSign
) {
6614 if ( ( ( (uint16_t) ( ( a
.high
| b
.high
) << 1 ) ) == 0) &&
6615 ( ( a
.low
| b
.low
) == 0 ) ) {
6617 return float_relation_equal
;
6619 return 1 - (2 * aSign
);
6622 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
6623 return float_relation_equal
;
6625 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
6630 int floatx80_compare( floatx80 a
, floatx80 b STATUS_PARAM
)
6632 return floatx80_compare_internal(a
, b
, 0 STATUS_VAR
);
6635 int floatx80_compare_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
6637 return floatx80_compare_internal(a
, b
, 1 STATUS_VAR
);
6640 INLINE
int float128_compare_internal( float128 a
, float128 b
,
6641 int is_quiet STATUS_PARAM
)
6645 if (( ( extractFloat128Exp( a
) == 0x7fff ) &&
6646 ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) ) ||
6647 ( ( extractFloat128Exp( b
) == 0x7fff ) &&
6648 ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )) {
6650 float128_is_signaling_nan( a
) ||
6651 float128_is_signaling_nan( b
) ) {
6652 float_raise( float_flag_invalid STATUS_VAR
);
6654 return float_relation_unordered
;
6656 aSign
= extractFloat128Sign( a
);
6657 bSign
= extractFloat128Sign( b
);
6658 if ( aSign
!= bSign
) {
6659 if ( ( ( ( a
.high
| b
.high
)<<1 ) | a
.low
| b
.low
) == 0 ) {
6661 return float_relation_equal
;
6663 return 1 - (2 * aSign
);
6666 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
6667 return float_relation_equal
;
6669 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
6674 int float128_compare( float128 a
, float128 b STATUS_PARAM
)
6676 return float128_compare_internal(a
, b
, 0 STATUS_VAR
);
6679 int float128_compare_quiet( float128 a
, float128 b STATUS_PARAM
)
6681 return float128_compare_internal(a
, b
, 1 STATUS_VAR
);
6684 /* min() and max() functions. These can't be implemented as
6685 * 'compare and pick one input' because that would mishandle
6686 * NaNs and +0 vs -0.
6688 #define MINMAX(s, nan_exp) \
6689 INLINE float ## s float ## s ## _minmax(float ## s a, float ## s b, \
6690 int ismin STATUS_PARAM ) \
6692 flag aSign, bSign; \
6693 uint ## s ## _t av, bv; \
6694 a = float ## s ## _squash_input_denormal(a STATUS_VAR); \
6695 b = float ## s ## _squash_input_denormal(b STATUS_VAR); \
6696 if (float ## s ## _is_any_nan(a) || \
6697 float ## s ## _is_any_nan(b)) { \
6698 return propagateFloat ## s ## NaN(a, b STATUS_VAR); \
6700 aSign = extractFloat ## s ## Sign(a); \
6701 bSign = extractFloat ## s ## Sign(b); \
6702 av = float ## s ## _val(a); \
6703 bv = float ## s ## _val(b); \
6704 if (aSign != bSign) { \
6706 return aSign ? a : b; \
6708 return aSign ? b : a; \
6712 return (aSign ^ (av < bv)) ? a : b; \
6714 return (aSign ^ (av < bv)) ? b : a; \
6719 float ## s float ## s ## _min(float ## s a, float ## s b STATUS_PARAM) \
6721 return float ## s ## _minmax(a, b, 1 STATUS_VAR); \
6724 float ## s float ## s ## _max(float ## s a, float ## s b STATUS_PARAM) \
6726 return float ## s ## _minmax(a, b, 0 STATUS_VAR); \
6733 /* Multiply A by 2 raised to the power N. */
6734 float32
float32_scalbn( float32 a
, int n STATUS_PARAM
)
6740 a
= float32_squash_input_denormal(a STATUS_VAR
);
6741 aSig
= extractFloat32Frac( a
);
6742 aExp
= extractFloat32Exp( a
);
6743 aSign
= extractFloat32Sign( a
);
6745 if ( aExp
== 0xFF ) {
6747 return propagateFloat32NaN( a
, a STATUS_VAR
);
6753 else if ( aSig
== 0 )
6758 } else if (n
< -0x200) {
6764 return normalizeRoundAndPackFloat32( aSign
, aExp
, aSig STATUS_VAR
);
6767 float64
float64_scalbn( float64 a
, int n STATUS_PARAM
)
6773 a
= float64_squash_input_denormal(a STATUS_VAR
);
6774 aSig
= extractFloat64Frac( a
);
6775 aExp
= extractFloat64Exp( a
);
6776 aSign
= extractFloat64Sign( a
);
6778 if ( aExp
== 0x7FF ) {
6780 return propagateFloat64NaN( a
, a STATUS_VAR
);
6785 aSig
|= LIT64( 0x0010000000000000 );
6786 else if ( aSig
== 0 )
6791 } else if (n
< -0x1000) {
6797 return normalizeRoundAndPackFloat64( aSign
, aExp
, aSig STATUS_VAR
);
6800 floatx80
floatx80_scalbn( floatx80 a
, int n STATUS_PARAM
)
6806 aSig
= extractFloatx80Frac( a
);
6807 aExp
= extractFloatx80Exp( a
);
6808 aSign
= extractFloatx80Sign( a
);
6810 if ( aExp
== 0x7FFF ) {
6812 return propagateFloatx80NaN( a
, a STATUS_VAR
);
6817 if (aExp
== 0 && aSig
== 0)
6822 } else if (n
< -0x10000) {
6827 return normalizeRoundAndPackFloatx80( STATUS(floatx80_rounding_precision
),
6828 aSign
, aExp
, aSig
, 0 STATUS_VAR
);
6831 float128
float128_scalbn( float128 a
, int n STATUS_PARAM
)
6835 uint64_t aSig0
, aSig1
;
6837 aSig1
= extractFloat128Frac1( a
);
6838 aSig0
= extractFloat128Frac0( a
);
6839 aExp
= extractFloat128Exp( a
);
6840 aSign
= extractFloat128Sign( a
);
6841 if ( aExp
== 0x7FFF ) {
6842 if ( aSig0
| aSig1
) {
6843 return propagateFloat128NaN( a
, a STATUS_VAR
);
6848 aSig0
|= LIT64( 0x0001000000000000 );
6849 else if ( aSig0
== 0 && aSig1
== 0 )
6854 } else if (n
< -0x10000) {
6859 return normalizeRoundAndPackFloat128( aSign
, aExp
, aSig0
, aSig1