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 #include "softfloat.h"
40 /*----------------------------------------------------------------------------
41 | Primitive arithmetic functions, including multi-word arithmetic, and
42 | division and square root approximations. (Can be specialized to target if
44 *----------------------------------------------------------------------------*/
45 #include "softfloat-macros.h"
47 /*----------------------------------------------------------------------------
48 | Functions and definitions to determine: (1) whether tininess for underflow
49 | is detected before or after rounding by default, (2) what (if anything)
50 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
51 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
52 | are propagated from function inputs to output. These details are target-
54 *----------------------------------------------------------------------------*/
55 #include "softfloat-specialize.h"
57 void set_float_rounding_mode(int val STATUS_PARAM
)
59 STATUS(float_rounding_mode
) = val
;
62 void set_float_exception_flags(int val STATUS_PARAM
)
64 STATUS(float_exception_flags
) = val
;
68 void set_floatx80_rounding_precision(int val STATUS_PARAM
)
70 STATUS(floatx80_rounding_precision
) = val
;
74 /*----------------------------------------------------------------------------
75 | Returns the fraction bits of the half-precision floating-point value `a'.
76 *----------------------------------------------------------------------------*/
78 INLINE
uint32_t extractFloat16Frac(float16 a
)
80 return float16_val(a
) & 0x3ff;
83 /*----------------------------------------------------------------------------
84 | Returns the exponent bits of the half-precision floating-point value `a'.
85 *----------------------------------------------------------------------------*/
87 INLINE int16
extractFloat16Exp(float16 a
)
89 return (float16_val(a
) >> 10) & 0x1f;
92 /*----------------------------------------------------------------------------
93 | Returns the sign bit of the single-precision floating-point value `a'.
94 *----------------------------------------------------------------------------*/
96 INLINE flag
extractFloat16Sign(float16 a
)
98 return float16_val(a
)>>15;
101 /*----------------------------------------------------------------------------
102 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
103 | and 7, and returns the properly rounded 32-bit integer corresponding to the
104 | input. If `zSign' is 1, the input is negated before being converted to an
105 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
106 | is simply rounded to an integer, with the inexact exception raised if the
107 | input cannot be represented exactly as an integer. However, if the fixed-
108 | point input is too large, the invalid exception is raised and the largest
109 | positive or negative integer is returned.
110 *----------------------------------------------------------------------------*/
112 static int32
roundAndPackInt32( flag zSign
, uint64_t absZ STATUS_PARAM
)
115 flag roundNearestEven
;
116 int8 roundIncrement
, roundBits
;
119 roundingMode
= STATUS(float_rounding_mode
);
120 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
121 roundIncrement
= 0x40;
122 if ( ! roundNearestEven
) {
123 if ( roundingMode
== float_round_to_zero
) {
127 roundIncrement
= 0x7F;
129 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
132 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
136 roundBits
= absZ
& 0x7F;
137 absZ
= ( absZ
+ roundIncrement
)>>7;
138 absZ
&= ~ ( ( ( roundBits
^ 0x40 ) == 0 ) & roundNearestEven
);
140 if ( zSign
) z
= - z
;
141 if ( ( absZ
>>32 ) || ( z
&& ( ( z
< 0 ) ^ zSign
) ) ) {
142 float_raise( float_flag_invalid STATUS_VAR
);
143 return zSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
145 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
150 /*----------------------------------------------------------------------------
151 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
152 | `absZ1', with binary point between bits 63 and 64 (between the input words),
153 | and returns the properly rounded 64-bit integer corresponding to the input.
154 | If `zSign' is 1, the input is negated before being converted to an integer.
155 | Ordinarily, the fixed-point input is simply rounded to an integer, with
156 | the inexact exception raised if the input cannot be represented exactly as
157 | an integer. However, if the fixed-point input is too large, the invalid
158 | exception is raised and the largest positive or negative integer is
160 *----------------------------------------------------------------------------*/
162 static int64
roundAndPackInt64( flag zSign
, uint64_t absZ0
, uint64_t absZ1 STATUS_PARAM
)
165 flag roundNearestEven
, increment
;
168 roundingMode
= STATUS(float_rounding_mode
);
169 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
170 increment
= ( (int64_t) absZ1
< 0 );
171 if ( ! roundNearestEven
) {
172 if ( roundingMode
== float_round_to_zero
) {
177 increment
= ( roundingMode
== float_round_down
) && absZ1
;
180 increment
= ( roundingMode
== float_round_up
) && absZ1
;
186 if ( absZ0
== 0 ) goto overflow
;
187 absZ0
&= ~ ( ( (uint64_t) ( absZ1
<<1 ) == 0 ) & roundNearestEven
);
190 if ( zSign
) z
= - z
;
191 if ( z
&& ( ( z
< 0 ) ^ zSign
) ) {
193 float_raise( float_flag_invalid STATUS_VAR
);
195 zSign
? (int64_t) LIT64( 0x8000000000000000 )
196 : LIT64( 0x7FFFFFFFFFFFFFFF );
198 if ( absZ1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
203 /*----------------------------------------------------------------------------
204 | Returns the fraction bits of the single-precision floating-point value `a'.
205 *----------------------------------------------------------------------------*/
207 INLINE
uint32_t extractFloat32Frac( float32 a
)
210 return float32_val(a
) & 0x007FFFFF;
214 /*----------------------------------------------------------------------------
215 | Returns the exponent bits of the single-precision floating-point value `a'.
216 *----------------------------------------------------------------------------*/
218 INLINE int16
extractFloat32Exp( float32 a
)
221 return ( float32_val(a
)>>23 ) & 0xFF;
225 /*----------------------------------------------------------------------------
226 | Returns the sign bit of the single-precision floating-point value `a'.
227 *----------------------------------------------------------------------------*/
229 INLINE flag
extractFloat32Sign( float32 a
)
232 return float32_val(a
)>>31;
236 /*----------------------------------------------------------------------------
237 | If `a' is denormal and we are in flush-to-zero mode then set the
238 | input-denormal exception and return zero. Otherwise just return the value.
239 *----------------------------------------------------------------------------*/
240 static float32
float32_squash_input_denormal(float32 a STATUS_PARAM
)
242 if (STATUS(flush_inputs_to_zero
)) {
243 if (extractFloat32Exp(a
) == 0 && extractFloat32Frac(a
) != 0) {
244 float_raise(float_flag_input_denormal STATUS_VAR
);
245 return make_float32(float32_val(a
) & 0x80000000);
251 /*----------------------------------------------------------------------------
252 | Normalizes the subnormal single-precision floating-point value represented
253 | by the denormalized significand `aSig'. The normalized exponent and
254 | significand are stored at the locations pointed to by `zExpPtr' and
255 | `zSigPtr', respectively.
256 *----------------------------------------------------------------------------*/
259 normalizeFloat32Subnormal( uint32_t aSig
, int16
*zExpPtr
, uint32_t *zSigPtr
)
263 shiftCount
= countLeadingZeros32( aSig
) - 8;
264 *zSigPtr
= aSig
<<shiftCount
;
265 *zExpPtr
= 1 - shiftCount
;
269 /*----------------------------------------------------------------------------
270 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
271 | single-precision floating-point value, returning the result. After being
272 | shifted into the proper positions, the three fields are simply added
273 | together to form the result. This means that any integer portion of `zSig'
274 | will be added into the exponent. Since a properly normalized significand
275 | will have an integer portion equal to 1, the `zExp' input should be 1 less
276 | than the desired result exponent whenever `zSig' is a complete, normalized
278 *----------------------------------------------------------------------------*/
280 INLINE float32
packFloat32( flag zSign
, int16 zExp
, uint32_t zSig
)
284 ( ( (uint32_t) zSign
)<<31 ) + ( ( (uint32_t) zExp
)<<23 ) + zSig
);
288 /*----------------------------------------------------------------------------
289 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
290 | and significand `zSig', and returns the proper single-precision floating-
291 | point value corresponding to the abstract input. Ordinarily, the abstract
292 | value is simply rounded and packed into the single-precision format, with
293 | the inexact exception raised if the abstract input cannot be represented
294 | exactly. However, if the abstract value is too large, the overflow and
295 | inexact exceptions are raised and an infinity or maximal finite value is
296 | returned. If the abstract value is too small, the input value is rounded to
297 | a subnormal number, and the underflow and inexact exceptions are raised if
298 | the abstract input cannot be represented exactly as a subnormal single-
299 | precision floating-point number.
300 | The input significand `zSig' has its binary point between bits 30
301 | and 29, which is 7 bits to the left of the usual location. This shifted
302 | significand must be normalized or smaller. If `zSig' is not normalized,
303 | `zExp' must be 0; in that case, the result returned is a subnormal number,
304 | and it must not require rounding. In the usual case that `zSig' is
305 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
306 | The handling of underflow and overflow follows the IEC/IEEE Standard for
307 | Binary Floating-Point Arithmetic.
308 *----------------------------------------------------------------------------*/
310 static float32
roundAndPackFloat32( flag zSign
, int16 zExp
, uint32_t zSig STATUS_PARAM
)
313 flag roundNearestEven
;
314 int8 roundIncrement
, roundBits
;
317 roundingMode
= STATUS(float_rounding_mode
);
318 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
319 roundIncrement
= 0x40;
320 if ( ! roundNearestEven
) {
321 if ( roundingMode
== float_round_to_zero
) {
325 roundIncrement
= 0x7F;
327 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
330 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
334 roundBits
= zSig
& 0x7F;
335 if ( 0xFD <= (uint16_t) zExp
) {
337 || ( ( zExp
== 0xFD )
338 && ( (int32_t) ( zSig
+ roundIncrement
) < 0 ) )
340 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
341 return packFloat32( zSign
, 0xFF, - ( roundIncrement
== 0 ));
344 if ( STATUS(flush_to_zero
) ) return packFloat32( zSign
, 0, 0 );
346 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
348 || ( zSig
+ roundIncrement
< 0x80000000 );
349 shift32RightJamming( zSig
, - zExp
, &zSig
);
351 roundBits
= zSig
& 0x7F;
352 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
355 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
356 zSig
= ( zSig
+ roundIncrement
)>>7;
357 zSig
&= ~ ( ( ( roundBits
^ 0x40 ) == 0 ) & roundNearestEven
);
358 if ( zSig
== 0 ) zExp
= 0;
359 return packFloat32( zSign
, zExp
, zSig
);
363 /*----------------------------------------------------------------------------
364 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
365 | and significand `zSig', and returns the proper single-precision floating-
366 | point value corresponding to the abstract input. This routine is just like
367 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
368 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
369 | floating-point exponent.
370 *----------------------------------------------------------------------------*/
373 normalizeRoundAndPackFloat32( flag zSign
, int16 zExp
, uint32_t zSig STATUS_PARAM
)
377 shiftCount
= countLeadingZeros32( zSig
) - 1;
378 return roundAndPackFloat32( zSign
, zExp
- shiftCount
, zSig
<<shiftCount STATUS_VAR
);
382 /*----------------------------------------------------------------------------
383 | Returns the fraction bits of the double-precision floating-point value `a'.
384 *----------------------------------------------------------------------------*/
386 INLINE
uint64_t extractFloat64Frac( float64 a
)
389 return float64_val(a
) & LIT64( 0x000FFFFFFFFFFFFF );
393 /*----------------------------------------------------------------------------
394 | Returns the exponent bits of the double-precision floating-point value `a'.
395 *----------------------------------------------------------------------------*/
397 INLINE int16
extractFloat64Exp( float64 a
)
400 return ( float64_val(a
)>>52 ) & 0x7FF;
404 /*----------------------------------------------------------------------------
405 | Returns the sign bit of the double-precision floating-point value `a'.
406 *----------------------------------------------------------------------------*/
408 INLINE flag
extractFloat64Sign( float64 a
)
411 return float64_val(a
)>>63;
415 /*----------------------------------------------------------------------------
416 | If `a' is denormal and we are in flush-to-zero mode then set the
417 | input-denormal exception and return zero. Otherwise just return the value.
418 *----------------------------------------------------------------------------*/
419 static float64
float64_squash_input_denormal(float64 a STATUS_PARAM
)
421 if (STATUS(flush_inputs_to_zero
)) {
422 if (extractFloat64Exp(a
) == 0 && extractFloat64Frac(a
) != 0) {
423 float_raise(float_flag_input_denormal STATUS_VAR
);
424 return make_float64(float64_val(a
) & (1ULL << 63));
430 /*----------------------------------------------------------------------------
431 | Normalizes the subnormal double-precision floating-point value represented
432 | by the denormalized significand `aSig'. The normalized exponent and
433 | significand are stored at the locations pointed to by `zExpPtr' and
434 | `zSigPtr', respectively.
435 *----------------------------------------------------------------------------*/
438 normalizeFloat64Subnormal( uint64_t aSig
, int16
*zExpPtr
, uint64_t *zSigPtr
)
442 shiftCount
= countLeadingZeros64( aSig
) - 11;
443 *zSigPtr
= aSig
<<shiftCount
;
444 *zExpPtr
= 1 - shiftCount
;
448 /*----------------------------------------------------------------------------
449 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
450 | double-precision floating-point value, returning the result. After being
451 | shifted into the proper positions, the three fields are simply added
452 | together to form the result. This means that any integer portion of `zSig'
453 | will be added into the exponent. Since a properly normalized significand
454 | will have an integer portion equal to 1, the `zExp' input should be 1 less
455 | than the desired result exponent whenever `zSig' is a complete, normalized
457 *----------------------------------------------------------------------------*/
459 INLINE float64
packFloat64( flag zSign
, int16 zExp
, uint64_t zSig
)
463 ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<52 ) + zSig
);
467 /*----------------------------------------------------------------------------
468 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
469 | and significand `zSig', and returns the proper double-precision floating-
470 | point value corresponding to the abstract input. Ordinarily, the abstract
471 | value is simply rounded and packed into the double-precision format, with
472 | the inexact exception raised if the abstract input cannot be represented
473 | exactly. However, if the abstract value is too large, the overflow and
474 | inexact exceptions are raised and an infinity or maximal finite value is
475 | returned. If the abstract value is too small, the input value is rounded
476 | to a subnormal number, and the underflow and inexact exceptions are raised
477 | if the abstract input cannot be represented exactly as a subnormal double-
478 | precision floating-point number.
479 | The input significand `zSig' has its binary point between bits 62
480 | and 61, which is 10 bits to the left of the usual location. This shifted
481 | significand must be normalized or smaller. If `zSig' is not normalized,
482 | `zExp' must be 0; in that case, the result returned is a subnormal number,
483 | and it must not require rounding. In the usual case that `zSig' is
484 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
485 | The handling of underflow and overflow follows the IEC/IEEE Standard for
486 | Binary Floating-Point Arithmetic.
487 *----------------------------------------------------------------------------*/
489 static float64
roundAndPackFloat64( flag zSign
, int16 zExp
, uint64_t zSig STATUS_PARAM
)
492 flag roundNearestEven
;
493 int16 roundIncrement
, roundBits
;
496 roundingMode
= STATUS(float_rounding_mode
);
497 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
498 roundIncrement
= 0x200;
499 if ( ! roundNearestEven
) {
500 if ( roundingMode
== float_round_to_zero
) {
504 roundIncrement
= 0x3FF;
506 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
509 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
513 roundBits
= zSig
& 0x3FF;
514 if ( 0x7FD <= (uint16_t) zExp
) {
515 if ( ( 0x7FD < zExp
)
516 || ( ( zExp
== 0x7FD )
517 && ( (int64_t) ( zSig
+ roundIncrement
) < 0 ) )
519 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
520 return packFloat64( zSign
, 0x7FF, - ( roundIncrement
== 0 ));
523 if ( STATUS(flush_to_zero
) ) return packFloat64( zSign
, 0, 0 );
525 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
527 || ( zSig
+ roundIncrement
< LIT64( 0x8000000000000000 ) );
528 shift64RightJamming( zSig
, - zExp
, &zSig
);
530 roundBits
= zSig
& 0x3FF;
531 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
534 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
535 zSig
= ( zSig
+ roundIncrement
)>>10;
536 zSig
&= ~ ( ( ( roundBits
^ 0x200 ) == 0 ) & roundNearestEven
);
537 if ( zSig
== 0 ) zExp
= 0;
538 return packFloat64( zSign
, zExp
, zSig
);
542 /*----------------------------------------------------------------------------
543 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
544 | and significand `zSig', and returns the proper double-precision floating-
545 | point value corresponding to the abstract input. This routine is just like
546 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
547 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
548 | floating-point exponent.
549 *----------------------------------------------------------------------------*/
552 normalizeRoundAndPackFloat64( flag zSign
, int16 zExp
, uint64_t zSig STATUS_PARAM
)
556 shiftCount
= countLeadingZeros64( zSig
) - 1;
557 return roundAndPackFloat64( zSign
, zExp
- shiftCount
, zSig
<<shiftCount STATUS_VAR
);
563 /*----------------------------------------------------------------------------
564 | Returns the fraction bits of the extended double-precision floating-point
566 *----------------------------------------------------------------------------*/
568 INLINE
uint64_t extractFloatx80Frac( floatx80 a
)
575 /*----------------------------------------------------------------------------
576 | Returns the exponent bits of the extended double-precision floating-point
578 *----------------------------------------------------------------------------*/
580 INLINE int32
extractFloatx80Exp( floatx80 a
)
583 return a
.high
& 0x7FFF;
587 /*----------------------------------------------------------------------------
588 | Returns the sign bit of the extended double-precision floating-point value
590 *----------------------------------------------------------------------------*/
592 INLINE flag
extractFloatx80Sign( floatx80 a
)
599 /*----------------------------------------------------------------------------
600 | Normalizes the subnormal extended double-precision floating-point value
601 | represented by the denormalized significand `aSig'. The normalized exponent
602 | and significand are stored at the locations pointed to by `zExpPtr' and
603 | `zSigPtr', respectively.
604 *----------------------------------------------------------------------------*/
607 normalizeFloatx80Subnormal( uint64_t aSig
, int32
*zExpPtr
, uint64_t *zSigPtr
)
611 shiftCount
= countLeadingZeros64( aSig
);
612 *zSigPtr
= aSig
<<shiftCount
;
613 *zExpPtr
= 1 - shiftCount
;
617 /*----------------------------------------------------------------------------
618 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
619 | extended double-precision floating-point value, returning the result.
620 *----------------------------------------------------------------------------*/
622 INLINE floatx80
packFloatx80( flag zSign
, int32 zExp
, uint64_t zSig
)
627 z
.high
= ( ( (uint16_t) zSign
)<<15 ) + zExp
;
632 /*----------------------------------------------------------------------------
633 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
634 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
635 | and returns the proper extended double-precision floating-point value
636 | corresponding to the abstract input. Ordinarily, the abstract value is
637 | rounded and packed into the extended double-precision format, with the
638 | inexact exception raised if the abstract input cannot be represented
639 | exactly. However, if the abstract value is too large, the overflow and
640 | inexact exceptions are raised and an infinity or maximal finite value is
641 | returned. If the abstract value is too small, the input value is rounded to
642 | a subnormal number, and the underflow and inexact exceptions are raised if
643 | the abstract input cannot be represented exactly as a subnormal extended
644 | double-precision floating-point number.
645 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
646 | number of bits as single or double precision, respectively. Otherwise, the
647 | result is rounded to the full precision of the extended double-precision
649 | The input significand must be normalized or smaller. If the input
650 | significand is not normalized, `zExp' must be 0; in that case, the result
651 | returned is a subnormal number, and it must not require rounding. The
652 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
653 | Floating-Point Arithmetic.
654 *----------------------------------------------------------------------------*/
657 roundAndPackFloatx80(
658 int8 roundingPrecision
, flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
662 flag roundNearestEven
, increment
, isTiny
;
663 int64 roundIncrement
, roundMask
, roundBits
;
665 roundingMode
= STATUS(float_rounding_mode
);
666 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
667 if ( roundingPrecision
== 80 ) goto precision80
;
668 if ( roundingPrecision
== 64 ) {
669 roundIncrement
= LIT64( 0x0000000000000400 );
670 roundMask
= LIT64( 0x00000000000007FF );
672 else if ( roundingPrecision
== 32 ) {
673 roundIncrement
= LIT64( 0x0000008000000000 );
674 roundMask
= LIT64( 0x000000FFFFFFFFFF );
679 zSig0
|= ( zSig1
!= 0 );
680 if ( ! roundNearestEven
) {
681 if ( roundingMode
== float_round_to_zero
) {
685 roundIncrement
= roundMask
;
687 if ( roundingMode
== float_round_up
) roundIncrement
= 0;
690 if ( roundingMode
== float_round_down
) roundIncrement
= 0;
694 roundBits
= zSig0
& roundMask
;
695 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
696 if ( ( 0x7FFE < zExp
)
697 || ( ( zExp
== 0x7FFE ) && ( zSig0
+ roundIncrement
< zSig0
) )
702 if ( STATUS(flush_to_zero
) ) return packFloatx80( zSign
, 0, 0 );
704 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
706 || ( zSig0
<= zSig0
+ roundIncrement
);
707 shift64RightJamming( zSig0
, 1 - zExp
, &zSig0
);
709 roundBits
= zSig0
& roundMask
;
710 if ( isTiny
&& roundBits
) float_raise( float_flag_underflow STATUS_VAR
);
711 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
712 zSig0
+= roundIncrement
;
713 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
714 roundIncrement
= roundMask
+ 1;
715 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
716 roundMask
|= roundIncrement
;
718 zSig0
&= ~ roundMask
;
719 return packFloatx80( zSign
, zExp
, zSig0
);
722 if ( roundBits
) STATUS(float_exception_flags
) |= float_flag_inexact
;
723 zSig0
+= roundIncrement
;
724 if ( zSig0
< roundIncrement
) {
726 zSig0
= LIT64( 0x8000000000000000 );
728 roundIncrement
= roundMask
+ 1;
729 if ( roundNearestEven
&& ( roundBits
<<1 == roundIncrement
) ) {
730 roundMask
|= roundIncrement
;
732 zSig0
&= ~ roundMask
;
733 if ( zSig0
== 0 ) zExp
= 0;
734 return packFloatx80( zSign
, zExp
, zSig0
);
736 increment
= ( (int64_t) zSig1
< 0 );
737 if ( ! roundNearestEven
) {
738 if ( roundingMode
== float_round_to_zero
) {
743 increment
= ( roundingMode
== float_round_down
) && zSig1
;
746 increment
= ( roundingMode
== float_round_up
) && zSig1
;
750 if ( 0x7FFD <= (uint32_t) ( zExp
- 1 ) ) {
751 if ( ( 0x7FFE < zExp
)
752 || ( ( zExp
== 0x7FFE )
753 && ( zSig0
== LIT64( 0xFFFFFFFFFFFFFFFF ) )
759 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
760 if ( ( roundingMode
== float_round_to_zero
)
761 || ( zSign
&& ( roundingMode
== float_round_up
) )
762 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
764 return packFloatx80( zSign
, 0x7FFE, ~ roundMask
);
766 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
770 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
773 || ( zSig0
< LIT64( 0xFFFFFFFFFFFFFFFF ) );
774 shift64ExtraRightJamming( zSig0
, zSig1
, 1 - zExp
, &zSig0
, &zSig1
);
776 if ( isTiny
&& zSig1
) float_raise( float_flag_underflow STATUS_VAR
);
777 if ( zSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
778 if ( roundNearestEven
) {
779 increment
= ( (int64_t) zSig1
< 0 );
783 increment
= ( roundingMode
== float_round_down
) && zSig1
;
786 increment
= ( roundingMode
== float_round_up
) && zSig1
;
792 ~ ( ( (uint64_t) ( zSig1
<<1 ) == 0 ) & roundNearestEven
);
793 if ( (int64_t) zSig0
< 0 ) zExp
= 1;
795 return packFloatx80( zSign
, zExp
, zSig0
);
798 if ( zSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
803 zSig0
= LIT64( 0x8000000000000000 );
806 zSig0
&= ~ ( ( (uint64_t) ( zSig1
<<1 ) == 0 ) & roundNearestEven
);
810 if ( zSig0
== 0 ) zExp
= 0;
812 return packFloatx80( zSign
, zExp
, zSig0
);
816 /*----------------------------------------------------------------------------
817 | Takes an abstract floating-point value having sign `zSign', exponent
818 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
819 | and returns the proper extended double-precision floating-point value
820 | corresponding to the abstract input. This routine is just like
821 | `roundAndPackFloatx80' except that the input significand does not have to be
823 *----------------------------------------------------------------------------*/
826 normalizeRoundAndPackFloatx80(
827 int8 roundingPrecision
, flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
837 shiftCount
= countLeadingZeros64( zSig0
);
838 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
841 roundAndPackFloatx80( roundingPrecision
, zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
849 /*----------------------------------------------------------------------------
850 | Returns the least-significant 64 fraction bits of the quadruple-precision
851 | floating-point value `a'.
852 *----------------------------------------------------------------------------*/
854 INLINE
uint64_t extractFloat128Frac1( float128 a
)
861 /*----------------------------------------------------------------------------
862 | Returns the most-significant 48 fraction bits of the quadruple-precision
863 | floating-point value `a'.
864 *----------------------------------------------------------------------------*/
866 INLINE
uint64_t extractFloat128Frac0( float128 a
)
869 return a
.high
& LIT64( 0x0000FFFFFFFFFFFF );
873 /*----------------------------------------------------------------------------
874 | Returns the exponent bits of the quadruple-precision floating-point value
876 *----------------------------------------------------------------------------*/
878 INLINE int32
extractFloat128Exp( float128 a
)
881 return ( a
.high
>>48 ) & 0x7FFF;
885 /*----------------------------------------------------------------------------
886 | Returns the sign bit of the quadruple-precision floating-point value `a'.
887 *----------------------------------------------------------------------------*/
889 INLINE flag
extractFloat128Sign( float128 a
)
896 /*----------------------------------------------------------------------------
897 | Normalizes the subnormal quadruple-precision floating-point value
898 | represented by the denormalized significand formed by the concatenation of
899 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
900 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
901 | significand are stored at the location pointed to by `zSig0Ptr', and the
902 | least significant 64 bits of the normalized significand are stored at the
903 | location pointed to by `zSig1Ptr'.
904 *----------------------------------------------------------------------------*/
907 normalizeFloat128Subnormal(
918 shiftCount
= countLeadingZeros64( aSig1
) - 15;
919 if ( shiftCount
< 0 ) {
920 *zSig0Ptr
= aSig1
>>( - shiftCount
);
921 *zSig1Ptr
= aSig1
<<( shiftCount
& 63 );
924 *zSig0Ptr
= aSig1
<<shiftCount
;
927 *zExpPtr
= - shiftCount
- 63;
930 shiftCount
= countLeadingZeros64( aSig0
) - 15;
931 shortShift128Left( aSig0
, aSig1
, shiftCount
, zSig0Ptr
, zSig1Ptr
);
932 *zExpPtr
= 1 - shiftCount
;
937 /*----------------------------------------------------------------------------
938 | Packs the sign `zSign', the exponent `zExp', and the significand formed
939 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
940 | floating-point value, returning the result. After being shifted into the
941 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
942 | added together to form the most significant 32 bits of the result. This
943 | means that any integer portion of `zSig0' will be added into the exponent.
944 | Since a properly normalized significand will have an integer portion equal
945 | to 1, the `zExp' input should be 1 less than the desired result exponent
946 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
948 *----------------------------------------------------------------------------*/
951 packFloat128( flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
)
956 z
.high
= ( ( (uint64_t) zSign
)<<63 ) + ( ( (uint64_t) zExp
)<<48 ) + zSig0
;
961 /*----------------------------------------------------------------------------
962 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
963 | and extended significand formed by the concatenation of `zSig0', `zSig1',
964 | and `zSig2', and returns the proper quadruple-precision floating-point value
965 | corresponding to the abstract input. Ordinarily, the abstract value is
966 | simply rounded and packed into the quadruple-precision format, with the
967 | inexact exception raised if the abstract input cannot be represented
968 | exactly. However, if the abstract value is too large, the overflow and
969 | inexact exceptions are raised and an infinity or maximal finite value is
970 | returned. If the abstract value is too small, the input value is rounded to
971 | a subnormal number, and the underflow and inexact exceptions are raised if
972 | the abstract input cannot be represented exactly as a subnormal quadruple-
973 | precision floating-point number.
974 | The input significand must be normalized or smaller. If the input
975 | significand is not normalized, `zExp' must be 0; in that case, the result
976 | returned is a subnormal number, and it must not require rounding. In the
977 | usual case that the input significand is normalized, `zExp' must be 1 less
978 | than the ``true'' floating-point exponent. The handling of underflow and
979 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
980 *----------------------------------------------------------------------------*/
983 roundAndPackFloat128(
984 flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1
, uint64_t zSig2 STATUS_PARAM
)
987 flag roundNearestEven
, increment
, isTiny
;
989 roundingMode
= STATUS(float_rounding_mode
);
990 roundNearestEven
= ( roundingMode
== float_round_nearest_even
);
991 increment
= ( (int64_t) zSig2
< 0 );
992 if ( ! roundNearestEven
) {
993 if ( roundingMode
== float_round_to_zero
) {
998 increment
= ( roundingMode
== float_round_down
) && zSig2
;
1001 increment
= ( roundingMode
== float_round_up
) && zSig2
;
1005 if ( 0x7FFD <= (uint32_t) zExp
) {
1006 if ( ( 0x7FFD < zExp
)
1007 || ( ( zExp
== 0x7FFD )
1009 LIT64( 0x0001FFFFFFFFFFFF ),
1010 LIT64( 0xFFFFFFFFFFFFFFFF ),
1017 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
1018 if ( ( roundingMode
== float_round_to_zero
)
1019 || ( zSign
&& ( roundingMode
== float_round_up
) )
1020 || ( ! zSign
&& ( roundingMode
== float_round_down
) )
1026 LIT64( 0x0000FFFFFFFFFFFF ),
1027 LIT64( 0xFFFFFFFFFFFFFFFF )
1030 return packFloat128( zSign
, 0x7FFF, 0, 0 );
1033 if ( STATUS(flush_to_zero
) ) return packFloat128( zSign
, 0, 0, 0 );
1035 ( STATUS(float_detect_tininess
) == float_tininess_before_rounding
)
1041 LIT64( 0x0001FFFFFFFFFFFF ),
1042 LIT64( 0xFFFFFFFFFFFFFFFF )
1044 shift128ExtraRightJamming(
1045 zSig0
, zSig1
, zSig2
, - zExp
, &zSig0
, &zSig1
, &zSig2
);
1047 if ( isTiny
&& zSig2
) float_raise( float_flag_underflow STATUS_VAR
);
1048 if ( roundNearestEven
) {
1049 increment
= ( (int64_t) zSig2
< 0 );
1053 increment
= ( roundingMode
== float_round_down
) && zSig2
;
1056 increment
= ( roundingMode
== float_round_up
) && zSig2
;
1061 if ( zSig2
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1063 add128( zSig0
, zSig1
, 0, 1, &zSig0
, &zSig1
);
1064 zSig1
&= ~ ( ( zSig2
+ zSig2
== 0 ) & roundNearestEven
);
1067 if ( ( zSig0
| zSig1
) == 0 ) zExp
= 0;
1069 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
1073 /*----------------------------------------------------------------------------
1074 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1075 | and significand formed by the concatenation of `zSig0' and `zSig1', and
1076 | returns the proper quadruple-precision floating-point value corresponding
1077 | to the abstract input. This routine is just like `roundAndPackFloat128'
1078 | except that the input significand has fewer bits and does not have to be
1079 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
1081 *----------------------------------------------------------------------------*/
1084 normalizeRoundAndPackFloat128(
1085 flag zSign
, int32 zExp
, uint64_t zSig0
, uint64_t zSig1 STATUS_PARAM
)
1095 shiftCount
= countLeadingZeros64( zSig0
) - 15;
1096 if ( 0 <= shiftCount
) {
1098 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
1101 shift128ExtraRightJamming(
1102 zSig0
, zSig1
, 0, - shiftCount
, &zSig0
, &zSig1
, &zSig2
);
1105 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
1111 /*----------------------------------------------------------------------------
1112 | Returns the result of converting the 32-bit two's complement integer `a'
1113 | to the single-precision floating-point format. The conversion is performed
1114 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1115 *----------------------------------------------------------------------------*/
1117 float32
int32_to_float32( int32 a STATUS_PARAM
)
1121 if ( a
== 0 ) return float32_zero
;
1122 if ( a
== (int32_t) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
1124 return normalizeRoundAndPackFloat32( zSign
, 0x9C, zSign
? - a
: a STATUS_VAR
);
1128 /*----------------------------------------------------------------------------
1129 | Returns the result of converting the 32-bit two's complement integer `a'
1130 | to the double-precision floating-point format. The conversion is performed
1131 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1132 *----------------------------------------------------------------------------*/
1134 float64
int32_to_float64( int32 a STATUS_PARAM
)
1141 if ( a
== 0 ) return float64_zero
;
1143 absA
= zSign
? - a
: a
;
1144 shiftCount
= countLeadingZeros32( absA
) + 21;
1146 return packFloat64( zSign
, 0x432 - shiftCount
, zSig
<<shiftCount
);
1152 /*----------------------------------------------------------------------------
1153 | Returns the result of converting the 32-bit two's complement integer `a'
1154 | to the extended double-precision floating-point format. The conversion
1155 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1157 *----------------------------------------------------------------------------*/
1159 floatx80
int32_to_floatx80( int32 a STATUS_PARAM
)
1166 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
1168 absA
= zSign
? - a
: a
;
1169 shiftCount
= countLeadingZeros32( absA
) + 32;
1171 return packFloatx80( zSign
, 0x403E - shiftCount
, zSig
<<shiftCount
);
1179 /*----------------------------------------------------------------------------
1180 | Returns the result of converting the 32-bit two's complement integer `a' to
1181 | the quadruple-precision floating-point format. The conversion is performed
1182 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1183 *----------------------------------------------------------------------------*/
1185 float128
int32_to_float128( int32 a STATUS_PARAM
)
1192 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
1194 absA
= zSign
? - a
: a
;
1195 shiftCount
= countLeadingZeros32( absA
) + 17;
1197 return packFloat128( zSign
, 0x402E - shiftCount
, zSig0
<<shiftCount
, 0 );
1203 /*----------------------------------------------------------------------------
1204 | Returns the result of converting the 64-bit two's complement integer `a'
1205 | to the single-precision floating-point format. The conversion is performed
1206 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1207 *----------------------------------------------------------------------------*/
1209 float32
int64_to_float32( int64 a STATUS_PARAM
)
1215 if ( a
== 0 ) return float32_zero
;
1217 absA
= zSign
? - a
: a
;
1218 shiftCount
= countLeadingZeros64( absA
) - 40;
1219 if ( 0 <= shiftCount
) {
1220 return packFloat32( zSign
, 0x95 - shiftCount
, absA
<<shiftCount
);
1224 if ( shiftCount
< 0 ) {
1225 shift64RightJamming( absA
, - shiftCount
, &absA
);
1228 absA
<<= shiftCount
;
1230 return roundAndPackFloat32( zSign
, 0x9C - shiftCount
, absA STATUS_VAR
);
1235 float32
uint64_to_float32( uint64 a STATUS_PARAM
)
1239 if ( a
== 0 ) return float32_zero
;
1240 shiftCount
= countLeadingZeros64( a
) - 40;
1241 if ( 0 <= shiftCount
) {
1242 return packFloat32( 1 > 0, 0x95 - shiftCount
, a
<<shiftCount
);
1246 if ( shiftCount
< 0 ) {
1247 shift64RightJamming( a
, - shiftCount
, &a
);
1252 return roundAndPackFloat32( 1 > 0, 0x9C - shiftCount
, a STATUS_VAR
);
1256 /*----------------------------------------------------------------------------
1257 | Returns the result of converting the 64-bit two's complement integer `a'
1258 | to the double-precision floating-point format. The conversion is performed
1259 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1260 *----------------------------------------------------------------------------*/
1262 float64
int64_to_float64( int64 a STATUS_PARAM
)
1266 if ( a
== 0 ) return float64_zero
;
1267 if ( a
== (int64_t) LIT64( 0x8000000000000000 ) ) {
1268 return packFloat64( 1, 0x43E, 0 );
1271 return normalizeRoundAndPackFloat64( zSign
, 0x43C, zSign
? - a
: a STATUS_VAR
);
1275 float64
uint64_to_float64( uint64 a STATUS_PARAM
)
1277 if ( a
== 0 ) return float64_zero
;
1278 return normalizeRoundAndPackFloat64( 0, 0x43C, a STATUS_VAR
);
1284 /*----------------------------------------------------------------------------
1285 | Returns the result of converting the 64-bit two's complement integer `a'
1286 | to the extended double-precision floating-point format. The conversion
1287 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1289 *----------------------------------------------------------------------------*/
1291 floatx80
int64_to_floatx80( int64 a STATUS_PARAM
)
1297 if ( a
== 0 ) return packFloatx80( 0, 0, 0 );
1299 absA
= zSign
? - a
: a
;
1300 shiftCount
= countLeadingZeros64( absA
);
1301 return packFloatx80( zSign
, 0x403E - shiftCount
, absA
<<shiftCount
);
1309 /*----------------------------------------------------------------------------
1310 | Returns the result of converting the 64-bit two's complement integer `a' to
1311 | the quadruple-precision floating-point format. The conversion is performed
1312 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1313 *----------------------------------------------------------------------------*/
1315 float128
int64_to_float128( int64 a STATUS_PARAM
)
1321 uint64_t zSig0
, zSig1
;
1323 if ( a
== 0 ) return packFloat128( 0, 0, 0, 0 );
1325 absA
= zSign
? - a
: a
;
1326 shiftCount
= countLeadingZeros64( absA
) + 49;
1327 zExp
= 0x406E - shiftCount
;
1328 if ( 64 <= shiftCount
) {
1337 shortShift128Left( zSig0
, zSig1
, shiftCount
, &zSig0
, &zSig1
);
1338 return packFloat128( zSign
, zExp
, zSig0
, zSig1
);
1344 /*----------------------------------------------------------------------------
1345 | Returns the result of converting the single-precision floating-point value
1346 | `a' to the 32-bit two's complement integer format. The conversion is
1347 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1348 | Arithmetic---which means in particular that the conversion is rounded
1349 | according to the current rounding mode. If `a' is a NaN, the largest
1350 | positive integer is returned. Otherwise, if the conversion overflows, the
1351 | largest integer with the same sign as `a' is returned.
1352 *----------------------------------------------------------------------------*/
1354 int32
float32_to_int32( float32 a STATUS_PARAM
)
1357 int16 aExp
, shiftCount
;
1361 a
= float32_squash_input_denormal(a STATUS_VAR
);
1362 aSig
= extractFloat32Frac( a
);
1363 aExp
= extractFloat32Exp( a
);
1364 aSign
= extractFloat32Sign( a
);
1365 if ( ( aExp
== 0xFF ) && aSig
) aSign
= 0;
1366 if ( aExp
) aSig
|= 0x00800000;
1367 shiftCount
= 0xAF - aExp
;
1370 if ( 0 < shiftCount
) shift64RightJamming( aSig64
, shiftCount
, &aSig64
);
1371 return roundAndPackInt32( aSign
, aSig64 STATUS_VAR
);
1375 /*----------------------------------------------------------------------------
1376 | Returns the result of converting the single-precision floating-point value
1377 | `a' to the 32-bit two's complement integer format. The conversion is
1378 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1379 | Arithmetic, except that the conversion is always rounded toward zero.
1380 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1381 | the conversion overflows, the largest integer with the same sign as `a' is
1383 *----------------------------------------------------------------------------*/
1385 int32
float32_to_int32_round_to_zero( float32 a STATUS_PARAM
)
1388 int16 aExp
, shiftCount
;
1391 a
= float32_squash_input_denormal(a STATUS_VAR
);
1393 aSig
= extractFloat32Frac( a
);
1394 aExp
= extractFloat32Exp( a
);
1395 aSign
= extractFloat32Sign( a
);
1396 shiftCount
= aExp
- 0x9E;
1397 if ( 0 <= shiftCount
) {
1398 if ( float32_val(a
) != 0xCF000000 ) {
1399 float_raise( float_flag_invalid STATUS_VAR
);
1400 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) return 0x7FFFFFFF;
1402 return (int32_t) 0x80000000;
1404 else if ( aExp
<= 0x7E ) {
1405 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1408 aSig
= ( aSig
| 0x00800000 )<<8;
1409 z
= aSig
>>( - shiftCount
);
1410 if ( (uint32_t) ( aSig
<<( shiftCount
& 31 ) ) ) {
1411 STATUS(float_exception_flags
) |= float_flag_inexact
;
1413 if ( aSign
) z
= - z
;
1418 /*----------------------------------------------------------------------------
1419 | Returns the result of converting the single-precision floating-point value
1420 | `a' to the 16-bit two's complement integer format. The conversion is
1421 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1422 | Arithmetic, except that the conversion is always rounded toward zero.
1423 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1424 | the conversion overflows, the largest integer with the same sign as `a' is
1426 *----------------------------------------------------------------------------*/
1428 int16
float32_to_int16_round_to_zero( float32 a STATUS_PARAM
)
1431 int16 aExp
, shiftCount
;
1435 aSig
= extractFloat32Frac( a
);
1436 aExp
= extractFloat32Exp( a
);
1437 aSign
= extractFloat32Sign( a
);
1438 shiftCount
= aExp
- 0x8E;
1439 if ( 0 <= shiftCount
) {
1440 if ( float32_val(a
) != 0xC7000000 ) {
1441 float_raise( float_flag_invalid STATUS_VAR
);
1442 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1446 return (int32_t) 0xffff8000;
1448 else if ( aExp
<= 0x7E ) {
1449 if ( aExp
| aSig
) {
1450 STATUS(float_exception_flags
) |= float_flag_inexact
;
1455 aSig
= ( aSig
| 0x00800000 )<<8;
1456 z
= aSig
>>( - shiftCount
);
1457 if ( (uint32_t) ( aSig
<<( shiftCount
& 31 ) ) ) {
1458 STATUS(float_exception_flags
) |= float_flag_inexact
;
1467 /*----------------------------------------------------------------------------
1468 | Returns the result of converting the single-precision floating-point value
1469 | `a' to the 64-bit two's complement integer format. The conversion is
1470 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1471 | Arithmetic---which means in particular that the conversion is rounded
1472 | according to the current rounding mode. If `a' is a NaN, the largest
1473 | positive integer is returned. Otherwise, if the conversion overflows, the
1474 | largest integer with the same sign as `a' is returned.
1475 *----------------------------------------------------------------------------*/
1477 int64
float32_to_int64( float32 a STATUS_PARAM
)
1480 int16 aExp
, shiftCount
;
1482 uint64_t aSig64
, aSigExtra
;
1483 a
= float32_squash_input_denormal(a STATUS_VAR
);
1485 aSig
= extractFloat32Frac( a
);
1486 aExp
= extractFloat32Exp( a
);
1487 aSign
= extractFloat32Sign( a
);
1488 shiftCount
= 0xBE - aExp
;
1489 if ( shiftCount
< 0 ) {
1490 float_raise( float_flag_invalid STATUS_VAR
);
1491 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1492 return LIT64( 0x7FFFFFFFFFFFFFFF );
1494 return (int64_t) LIT64( 0x8000000000000000 );
1496 if ( aExp
) aSig
|= 0x00800000;
1499 shift64ExtraRightJamming( aSig64
, 0, shiftCount
, &aSig64
, &aSigExtra
);
1500 return roundAndPackInt64( aSign
, aSig64
, aSigExtra STATUS_VAR
);
1504 /*----------------------------------------------------------------------------
1505 | Returns the result of converting the single-precision floating-point value
1506 | `a' to the 64-bit two's complement integer format. The conversion is
1507 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1508 | Arithmetic, except that the conversion is always rounded toward zero. If
1509 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
1510 | conversion overflows, the largest integer with the same sign as `a' is
1512 *----------------------------------------------------------------------------*/
1514 int64
float32_to_int64_round_to_zero( float32 a STATUS_PARAM
)
1517 int16 aExp
, shiftCount
;
1521 a
= float32_squash_input_denormal(a STATUS_VAR
);
1523 aSig
= extractFloat32Frac( a
);
1524 aExp
= extractFloat32Exp( a
);
1525 aSign
= extractFloat32Sign( a
);
1526 shiftCount
= aExp
- 0xBE;
1527 if ( 0 <= shiftCount
) {
1528 if ( float32_val(a
) != 0xDF000000 ) {
1529 float_raise( float_flag_invalid STATUS_VAR
);
1530 if ( ! aSign
|| ( ( aExp
== 0xFF ) && aSig
) ) {
1531 return LIT64( 0x7FFFFFFFFFFFFFFF );
1534 return (int64_t) LIT64( 0x8000000000000000 );
1536 else if ( aExp
<= 0x7E ) {
1537 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
1540 aSig64
= aSig
| 0x00800000;
1542 z
= aSig64
>>( - shiftCount
);
1543 if ( (uint64_t) ( aSig64
<<( shiftCount
& 63 ) ) ) {
1544 STATUS(float_exception_flags
) |= float_flag_inexact
;
1546 if ( aSign
) z
= - z
;
1551 /*----------------------------------------------------------------------------
1552 | Returns the result of converting the single-precision floating-point value
1553 | `a' to the double-precision floating-point format. The conversion is
1554 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1556 *----------------------------------------------------------------------------*/
1558 float64
float32_to_float64( float32 a STATUS_PARAM
)
1563 a
= float32_squash_input_denormal(a STATUS_VAR
);
1565 aSig
= extractFloat32Frac( a
);
1566 aExp
= extractFloat32Exp( a
);
1567 aSign
= extractFloat32Sign( a
);
1568 if ( aExp
== 0xFF ) {
1569 if ( aSig
) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1570 return packFloat64( aSign
, 0x7FF, 0 );
1573 if ( aSig
== 0 ) return packFloat64( aSign
, 0, 0 );
1574 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1577 return packFloat64( aSign
, aExp
+ 0x380, ( (uint64_t) aSig
)<<29 );
1583 /*----------------------------------------------------------------------------
1584 | Returns the result of converting the single-precision floating-point value
1585 | `a' to the extended double-precision floating-point format. The conversion
1586 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1588 *----------------------------------------------------------------------------*/
1590 floatx80
float32_to_floatx80( float32 a STATUS_PARAM
)
1596 a
= float32_squash_input_denormal(a STATUS_VAR
);
1597 aSig
= extractFloat32Frac( a
);
1598 aExp
= extractFloat32Exp( a
);
1599 aSign
= extractFloat32Sign( a
);
1600 if ( aExp
== 0xFF ) {
1601 if ( aSig
) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1602 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
1605 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
1606 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1609 return packFloatx80( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<40 );
1617 /*----------------------------------------------------------------------------
1618 | Returns the result of converting the single-precision floating-point value
1619 | `a' to the double-precision floating-point format. The conversion is
1620 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1622 *----------------------------------------------------------------------------*/
1624 float128
float32_to_float128( float32 a STATUS_PARAM
)
1630 a
= float32_squash_input_denormal(a STATUS_VAR
);
1631 aSig
= extractFloat32Frac( a
);
1632 aExp
= extractFloat32Exp( a
);
1633 aSign
= extractFloat32Sign( a
);
1634 if ( aExp
== 0xFF ) {
1635 if ( aSig
) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
1636 return packFloat128( aSign
, 0x7FFF, 0, 0 );
1639 if ( aSig
== 0 ) return packFloat128( aSign
, 0, 0, 0 );
1640 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1643 return packFloat128( aSign
, aExp
+ 0x3F80, ( (uint64_t) aSig
)<<25, 0 );
1649 /*----------------------------------------------------------------------------
1650 | Rounds the single-precision floating-point value `a' to an integer, and
1651 | returns the result as a single-precision floating-point value. The
1652 | operation is performed according to the IEC/IEEE Standard for Binary
1653 | Floating-Point Arithmetic.
1654 *----------------------------------------------------------------------------*/
1656 float32
float32_round_to_int( float32 a STATUS_PARAM
)
1660 uint32_t lastBitMask
, roundBitsMask
;
1663 a
= float32_squash_input_denormal(a STATUS_VAR
);
1665 aExp
= extractFloat32Exp( a
);
1666 if ( 0x96 <= aExp
) {
1667 if ( ( aExp
== 0xFF ) && extractFloat32Frac( a
) ) {
1668 return propagateFloat32NaN( a
, a STATUS_VAR
);
1672 if ( aExp
<= 0x7E ) {
1673 if ( (uint32_t) ( float32_val(a
)<<1 ) == 0 ) return a
;
1674 STATUS(float_exception_flags
) |= float_flag_inexact
;
1675 aSign
= extractFloat32Sign( a
);
1676 switch ( STATUS(float_rounding_mode
) ) {
1677 case float_round_nearest_even
:
1678 if ( ( aExp
== 0x7E ) && extractFloat32Frac( a
) ) {
1679 return packFloat32( aSign
, 0x7F, 0 );
1682 case float_round_down
:
1683 return make_float32(aSign
? 0xBF800000 : 0);
1684 case float_round_up
:
1685 return make_float32(aSign
? 0x80000000 : 0x3F800000);
1687 return packFloat32( aSign
, 0, 0 );
1690 lastBitMask
<<= 0x96 - aExp
;
1691 roundBitsMask
= lastBitMask
- 1;
1693 roundingMode
= STATUS(float_rounding_mode
);
1694 if ( roundingMode
== float_round_nearest_even
) {
1695 z
+= lastBitMask
>>1;
1696 if ( ( z
& roundBitsMask
) == 0 ) z
&= ~ lastBitMask
;
1698 else if ( roundingMode
!= float_round_to_zero
) {
1699 if ( extractFloat32Sign( make_float32(z
) ) ^ ( roundingMode
== float_round_up
) ) {
1703 z
&= ~ roundBitsMask
;
1704 if ( z
!= float32_val(a
) ) STATUS(float_exception_flags
) |= float_flag_inexact
;
1705 return make_float32(z
);
1709 /*----------------------------------------------------------------------------
1710 | Returns the result of adding the absolute values of the single-precision
1711 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
1712 | before being returned. `zSign' is ignored if the result is a NaN.
1713 | The addition is performed according to the IEC/IEEE Standard for Binary
1714 | Floating-Point Arithmetic.
1715 *----------------------------------------------------------------------------*/
1717 static float32
addFloat32Sigs( float32 a
, float32 b
, flag zSign STATUS_PARAM
)
1719 int16 aExp
, bExp
, zExp
;
1720 uint32_t aSig
, bSig
, zSig
;
1723 aSig
= extractFloat32Frac( a
);
1724 aExp
= extractFloat32Exp( a
);
1725 bSig
= extractFloat32Frac( b
);
1726 bExp
= extractFloat32Exp( b
);
1727 expDiff
= aExp
- bExp
;
1730 if ( 0 < expDiff
) {
1731 if ( aExp
== 0xFF ) {
1732 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1741 shift32RightJamming( bSig
, expDiff
, &bSig
);
1744 else if ( expDiff
< 0 ) {
1745 if ( bExp
== 0xFF ) {
1746 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1747 return packFloat32( zSign
, 0xFF, 0 );
1755 shift32RightJamming( aSig
, - expDiff
, &aSig
);
1759 if ( aExp
== 0xFF ) {
1760 if ( aSig
| bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1764 if ( STATUS(flush_to_zero
) ) return packFloat32( zSign
, 0, 0 );
1765 return packFloat32( zSign
, 0, ( aSig
+ bSig
)>>6 );
1767 zSig
= 0x40000000 + aSig
+ bSig
;
1772 zSig
= ( aSig
+ bSig
)<<1;
1774 if ( (int32_t) zSig
< 0 ) {
1779 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1783 /*----------------------------------------------------------------------------
1784 | Returns the result of subtracting the absolute values of the single-
1785 | precision floating-point values `a' and `b'. If `zSign' is 1, the
1786 | difference is negated before being returned. `zSign' is ignored if the
1787 | result is a NaN. The subtraction is performed according to the IEC/IEEE
1788 | Standard for Binary Floating-Point Arithmetic.
1789 *----------------------------------------------------------------------------*/
1791 static float32
subFloat32Sigs( float32 a
, float32 b
, flag zSign STATUS_PARAM
)
1793 int16 aExp
, bExp
, zExp
;
1794 uint32_t aSig
, bSig
, zSig
;
1797 aSig
= extractFloat32Frac( a
);
1798 aExp
= extractFloat32Exp( a
);
1799 bSig
= extractFloat32Frac( b
);
1800 bExp
= extractFloat32Exp( b
);
1801 expDiff
= aExp
- bExp
;
1804 if ( 0 < expDiff
) goto aExpBigger
;
1805 if ( expDiff
< 0 ) goto bExpBigger
;
1806 if ( aExp
== 0xFF ) {
1807 if ( aSig
| bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1808 float_raise( float_flag_invalid STATUS_VAR
);
1809 return float32_default_nan
;
1815 if ( bSig
< aSig
) goto aBigger
;
1816 if ( aSig
< bSig
) goto bBigger
;
1817 return packFloat32( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
1819 if ( bExp
== 0xFF ) {
1820 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1821 return packFloat32( zSign
^ 1, 0xFF, 0 );
1829 shift32RightJamming( aSig
, - expDiff
, &aSig
);
1835 goto normalizeRoundAndPack
;
1837 if ( aExp
== 0xFF ) {
1838 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1847 shift32RightJamming( bSig
, expDiff
, &bSig
);
1852 normalizeRoundAndPack
:
1854 return normalizeRoundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1858 /*----------------------------------------------------------------------------
1859 | Returns the result of adding the single-precision floating-point values `a'
1860 | and `b'. The operation is performed according to the IEC/IEEE Standard for
1861 | Binary Floating-Point Arithmetic.
1862 *----------------------------------------------------------------------------*/
1864 float32
float32_add( float32 a
, float32 b STATUS_PARAM
)
1867 a
= float32_squash_input_denormal(a STATUS_VAR
);
1868 b
= float32_squash_input_denormal(b STATUS_VAR
);
1870 aSign
= extractFloat32Sign( a
);
1871 bSign
= extractFloat32Sign( b
);
1872 if ( aSign
== bSign
) {
1873 return addFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1876 return subFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1881 /*----------------------------------------------------------------------------
1882 | Returns the result of subtracting the single-precision floating-point values
1883 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1884 | for Binary Floating-Point Arithmetic.
1885 *----------------------------------------------------------------------------*/
1887 float32
float32_sub( float32 a
, float32 b STATUS_PARAM
)
1890 a
= float32_squash_input_denormal(a STATUS_VAR
);
1891 b
= float32_squash_input_denormal(b STATUS_VAR
);
1893 aSign
= extractFloat32Sign( a
);
1894 bSign
= extractFloat32Sign( b
);
1895 if ( aSign
== bSign
) {
1896 return subFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1899 return addFloat32Sigs( a
, b
, aSign STATUS_VAR
);
1904 /*----------------------------------------------------------------------------
1905 | Returns the result of multiplying the single-precision floating-point values
1906 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1907 | for Binary Floating-Point Arithmetic.
1908 *----------------------------------------------------------------------------*/
1910 float32
float32_mul( float32 a
, float32 b STATUS_PARAM
)
1912 flag aSign
, bSign
, zSign
;
1913 int16 aExp
, bExp
, zExp
;
1914 uint32_t aSig
, bSig
;
1918 a
= float32_squash_input_denormal(a STATUS_VAR
);
1919 b
= float32_squash_input_denormal(b STATUS_VAR
);
1921 aSig
= extractFloat32Frac( a
);
1922 aExp
= extractFloat32Exp( a
);
1923 aSign
= extractFloat32Sign( a
);
1924 bSig
= extractFloat32Frac( b
);
1925 bExp
= extractFloat32Exp( b
);
1926 bSign
= extractFloat32Sign( b
);
1927 zSign
= aSign
^ bSign
;
1928 if ( aExp
== 0xFF ) {
1929 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
1930 return propagateFloat32NaN( a
, b STATUS_VAR
);
1932 if ( ( bExp
| bSig
) == 0 ) {
1933 float_raise( float_flag_invalid STATUS_VAR
);
1934 return float32_default_nan
;
1936 return packFloat32( zSign
, 0xFF, 0 );
1938 if ( bExp
== 0xFF ) {
1939 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1940 if ( ( aExp
| aSig
) == 0 ) {
1941 float_raise( float_flag_invalid STATUS_VAR
);
1942 return float32_default_nan
;
1944 return packFloat32( zSign
, 0xFF, 0 );
1947 if ( aSig
== 0 ) return packFloat32( zSign
, 0, 0 );
1948 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
1951 if ( bSig
== 0 ) return packFloat32( zSign
, 0, 0 );
1952 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
1954 zExp
= aExp
+ bExp
- 0x7F;
1955 aSig
= ( aSig
| 0x00800000 )<<7;
1956 bSig
= ( bSig
| 0x00800000 )<<8;
1957 shift64RightJamming( ( (uint64_t) aSig
) * bSig
, 32, &zSig64
);
1959 if ( 0 <= (int32_t) ( zSig
<<1 ) ) {
1963 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
1967 /*----------------------------------------------------------------------------
1968 | Returns the result of dividing the single-precision floating-point value `a'
1969 | by the corresponding value `b'. The operation is performed according to the
1970 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1971 *----------------------------------------------------------------------------*/
1973 float32
float32_div( float32 a
, float32 b STATUS_PARAM
)
1975 flag aSign
, bSign
, zSign
;
1976 int16 aExp
, bExp
, zExp
;
1977 uint32_t aSig
, bSig
, zSig
;
1978 a
= float32_squash_input_denormal(a STATUS_VAR
);
1979 b
= float32_squash_input_denormal(b STATUS_VAR
);
1981 aSig
= extractFloat32Frac( a
);
1982 aExp
= extractFloat32Exp( a
);
1983 aSign
= extractFloat32Sign( a
);
1984 bSig
= extractFloat32Frac( b
);
1985 bExp
= extractFloat32Exp( b
);
1986 bSign
= extractFloat32Sign( b
);
1987 zSign
= aSign
^ bSign
;
1988 if ( aExp
== 0xFF ) {
1989 if ( aSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1990 if ( bExp
== 0xFF ) {
1991 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1992 float_raise( float_flag_invalid STATUS_VAR
);
1993 return float32_default_nan
;
1995 return packFloat32( zSign
, 0xFF, 0 );
1997 if ( bExp
== 0xFF ) {
1998 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
1999 return packFloat32( zSign
, 0, 0 );
2003 if ( ( aExp
| aSig
) == 0 ) {
2004 float_raise( float_flag_invalid STATUS_VAR
);
2005 return float32_default_nan
;
2007 float_raise( float_flag_divbyzero STATUS_VAR
);
2008 return packFloat32( zSign
, 0xFF, 0 );
2010 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
2013 if ( aSig
== 0 ) return packFloat32( zSign
, 0, 0 );
2014 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2016 zExp
= aExp
- bExp
+ 0x7D;
2017 aSig
= ( aSig
| 0x00800000 )<<7;
2018 bSig
= ( bSig
| 0x00800000 )<<8;
2019 if ( bSig
<= ( aSig
+ aSig
) ) {
2023 zSig
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
2024 if ( ( zSig
& 0x3F ) == 0 ) {
2025 zSig
|= ( (uint64_t) bSig
* zSig
!= ( (uint64_t) aSig
)<<32 );
2027 return roundAndPackFloat32( zSign
, zExp
, zSig STATUS_VAR
);
2031 /*----------------------------------------------------------------------------
2032 | Returns the remainder of the single-precision floating-point value `a'
2033 | with respect to the corresponding value `b'. The operation is performed
2034 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2035 *----------------------------------------------------------------------------*/
2037 float32
float32_rem( float32 a
, float32 b STATUS_PARAM
)
2040 int16 aExp
, bExp
, expDiff
;
2041 uint32_t aSig
, bSig
;
2043 uint64_t aSig64
, bSig64
, q64
;
2044 uint32_t alternateASig
;
2046 a
= float32_squash_input_denormal(a STATUS_VAR
);
2047 b
= float32_squash_input_denormal(b STATUS_VAR
);
2049 aSig
= extractFloat32Frac( a
);
2050 aExp
= extractFloat32Exp( a
);
2051 aSign
= extractFloat32Sign( a
);
2052 bSig
= extractFloat32Frac( b
);
2053 bExp
= extractFloat32Exp( b
);
2054 if ( aExp
== 0xFF ) {
2055 if ( aSig
|| ( ( bExp
== 0xFF ) && bSig
) ) {
2056 return propagateFloat32NaN( a
, b STATUS_VAR
);
2058 float_raise( float_flag_invalid STATUS_VAR
);
2059 return float32_default_nan
;
2061 if ( bExp
== 0xFF ) {
2062 if ( bSig
) return propagateFloat32NaN( a
, b STATUS_VAR
);
2067 float_raise( float_flag_invalid STATUS_VAR
);
2068 return float32_default_nan
;
2070 normalizeFloat32Subnormal( bSig
, &bExp
, &bSig
);
2073 if ( aSig
== 0 ) return a
;
2074 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2076 expDiff
= aExp
- bExp
;
2079 if ( expDiff
< 32 ) {
2082 if ( expDiff
< 0 ) {
2083 if ( expDiff
< -1 ) return a
;
2086 q
= ( bSig
<= aSig
);
2087 if ( q
) aSig
-= bSig
;
2088 if ( 0 < expDiff
) {
2089 q
= ( ( (uint64_t) aSig
)<<32 ) / bSig
;
2092 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
2100 if ( bSig
<= aSig
) aSig
-= bSig
;
2101 aSig64
= ( (uint64_t) aSig
)<<40;
2102 bSig64
= ( (uint64_t) bSig
)<<40;
2104 while ( 0 < expDiff
) {
2105 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
2106 q64
= ( 2 < q64
) ? q64
- 2 : 0;
2107 aSig64
= - ( ( bSig
* q64
)<<38 );
2111 q64
= estimateDiv128To64( aSig64
, 0, bSig64
);
2112 q64
= ( 2 < q64
) ? q64
- 2 : 0;
2113 q
= q64
>>( 64 - expDiff
);
2115 aSig
= ( ( aSig64
>>33 )<<( expDiff
- 1 ) ) - bSig
* q
;
2118 alternateASig
= aSig
;
2121 } while ( 0 <= (int32_t) aSig
);
2122 sigMean
= aSig
+ alternateASig
;
2123 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
2124 aSig
= alternateASig
;
2126 zSign
= ( (int32_t) aSig
< 0 );
2127 if ( zSign
) aSig
= - aSig
;
2128 return normalizeRoundAndPackFloat32( aSign
^ zSign
, bExp
, aSig STATUS_VAR
);
2132 /*----------------------------------------------------------------------------
2133 | Returns the square root of the single-precision floating-point value `a'.
2134 | The operation is performed according to the IEC/IEEE Standard for Binary
2135 | Floating-Point Arithmetic.
2136 *----------------------------------------------------------------------------*/
2138 float32
float32_sqrt( float32 a STATUS_PARAM
)
2142 uint32_t aSig
, zSig
;
2144 a
= float32_squash_input_denormal(a STATUS_VAR
);
2146 aSig
= extractFloat32Frac( a
);
2147 aExp
= extractFloat32Exp( a
);
2148 aSign
= extractFloat32Sign( a
);
2149 if ( aExp
== 0xFF ) {
2150 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2151 if ( ! aSign
) return a
;
2152 float_raise( float_flag_invalid STATUS_VAR
);
2153 return float32_default_nan
;
2156 if ( ( aExp
| aSig
) == 0 ) return a
;
2157 float_raise( float_flag_invalid STATUS_VAR
);
2158 return float32_default_nan
;
2161 if ( aSig
== 0 ) return float32_zero
;
2162 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2164 zExp
= ( ( aExp
- 0x7F )>>1 ) + 0x7E;
2165 aSig
= ( aSig
| 0x00800000 )<<8;
2166 zSig
= estimateSqrt32( aExp
, aSig
) + 2;
2167 if ( ( zSig
& 0x7F ) <= 5 ) {
2173 term
= ( (uint64_t) zSig
) * zSig
;
2174 rem
= ( ( (uint64_t) aSig
)<<32 ) - term
;
2175 while ( (int64_t) rem
< 0 ) {
2177 rem
+= ( ( (uint64_t) zSig
)<<1 ) | 1;
2179 zSig
|= ( rem
!= 0 );
2181 shift32RightJamming( zSig
, 1, &zSig
);
2183 return roundAndPackFloat32( 0, zExp
, zSig STATUS_VAR
);
2187 /*----------------------------------------------------------------------------
2188 | Returns the binary exponential of the single-precision floating-point value
2189 | `a'. The operation is performed according to the IEC/IEEE Standard for
2190 | Binary Floating-Point Arithmetic.
2192 | Uses the following identities:
2194 | 1. -------------------------------------------------------------------------
2198 | 2. -------------------------------------------------------------------------
2201 | e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
2203 *----------------------------------------------------------------------------*/
2205 static const float64 float32_exp2_coefficients
[15] =
2207 const_float64( 0x3ff0000000000000ll
), /* 1 */
2208 const_float64( 0x3fe0000000000000ll
), /* 2 */
2209 const_float64( 0x3fc5555555555555ll
), /* 3 */
2210 const_float64( 0x3fa5555555555555ll
), /* 4 */
2211 const_float64( 0x3f81111111111111ll
), /* 5 */
2212 const_float64( 0x3f56c16c16c16c17ll
), /* 6 */
2213 const_float64( 0x3f2a01a01a01a01all
), /* 7 */
2214 const_float64( 0x3efa01a01a01a01all
), /* 8 */
2215 const_float64( 0x3ec71de3a556c734ll
), /* 9 */
2216 const_float64( 0x3e927e4fb7789f5cll
), /* 10 */
2217 const_float64( 0x3e5ae64567f544e4ll
), /* 11 */
2218 const_float64( 0x3e21eed8eff8d898ll
), /* 12 */
2219 const_float64( 0x3de6124613a86d09ll
), /* 13 */
2220 const_float64( 0x3da93974a8c07c9dll
), /* 14 */
2221 const_float64( 0x3d6ae7f3e733b81fll
), /* 15 */
2224 float32
float32_exp2( float32 a STATUS_PARAM
)
2231 a
= float32_squash_input_denormal(a STATUS_VAR
);
2233 aSig
= extractFloat32Frac( a
);
2234 aExp
= extractFloat32Exp( a
);
2235 aSign
= extractFloat32Sign( a
);
2237 if ( aExp
== 0xFF) {
2238 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2239 return (aSign
) ? float32_zero
: a
;
2242 if (aSig
== 0) return float32_one
;
2245 float_raise( float_flag_inexact STATUS_VAR
);
2247 /* ******************************* */
2248 /* using float64 for approximation */
2249 /* ******************************* */
2250 x
= float32_to_float64(a STATUS_VAR
);
2251 x
= float64_mul(x
, float64_ln2 STATUS_VAR
);
2255 for (i
= 0 ; i
< 15 ; i
++) {
2258 f
= float64_mul(xn
, float32_exp2_coefficients
[i
] STATUS_VAR
);
2259 r
= float64_add(r
, f STATUS_VAR
);
2261 xn
= float64_mul(xn
, x STATUS_VAR
);
2264 return float64_to_float32(r
, status
);
2267 /*----------------------------------------------------------------------------
2268 | Returns the binary log of the single-precision floating-point value `a'.
2269 | The operation is performed according to the IEC/IEEE Standard for Binary
2270 | Floating-Point Arithmetic.
2271 *----------------------------------------------------------------------------*/
2272 float32
float32_log2( float32 a STATUS_PARAM
)
2276 uint32_t aSig
, zSig
, i
;
2278 a
= float32_squash_input_denormal(a STATUS_VAR
);
2279 aSig
= extractFloat32Frac( a
);
2280 aExp
= extractFloat32Exp( a
);
2281 aSign
= extractFloat32Sign( a
);
2284 if ( aSig
== 0 ) return packFloat32( 1, 0xFF, 0 );
2285 normalizeFloat32Subnormal( aSig
, &aExp
, &aSig
);
2288 float_raise( float_flag_invalid STATUS_VAR
);
2289 return float32_default_nan
;
2291 if ( aExp
== 0xFF ) {
2292 if ( aSig
) return propagateFloat32NaN( a
, float32_zero STATUS_VAR
);
2301 for (i
= 1 << 22; i
> 0; i
>>= 1) {
2302 aSig
= ( (uint64_t)aSig
* aSig
) >> 23;
2303 if ( aSig
& 0x01000000 ) {
2312 return normalizeRoundAndPackFloat32( zSign
, 0x85, zSig STATUS_VAR
);
2315 /*----------------------------------------------------------------------------
2316 | Returns 1 if the single-precision floating-point value `a' is equal to
2317 | the corresponding value `b', and 0 otherwise. The invalid exception is
2318 | raised if either operand is a NaN. Otherwise, the comparison is performed
2319 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2320 *----------------------------------------------------------------------------*/
2322 int float32_eq( float32 a
, float32 b STATUS_PARAM
)
2325 a
= float32_squash_input_denormal(a STATUS_VAR
);
2326 b
= float32_squash_input_denormal(b STATUS_VAR
);
2328 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2329 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2331 float_raise( float_flag_invalid STATUS_VAR
);
2334 av
= float32_val(a
);
2335 bv
= float32_val(b
);
2336 return ( av
== bv
) || ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2339 /*----------------------------------------------------------------------------
2340 | Returns 1 if the single-precision floating-point value `a' is less than
2341 | or equal to the corresponding value `b', and 0 otherwise. The comparison
2342 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2344 *----------------------------------------------------------------------------*/
2346 int float32_le( float32 a
, float32 b STATUS_PARAM
)
2350 a
= float32_squash_input_denormal(a STATUS_VAR
);
2351 b
= float32_squash_input_denormal(b STATUS_VAR
);
2353 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2354 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2356 float_raise( float_flag_invalid STATUS_VAR
);
2359 aSign
= extractFloat32Sign( a
);
2360 bSign
= extractFloat32Sign( b
);
2361 av
= float32_val(a
);
2362 bv
= float32_val(b
);
2363 if ( aSign
!= bSign
) return aSign
|| ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2364 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
2368 /*----------------------------------------------------------------------------
2369 | Returns 1 if the single-precision floating-point value `a' is less than
2370 | the corresponding value `b', and 0 otherwise. The comparison is performed
2371 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2372 *----------------------------------------------------------------------------*/
2374 int float32_lt( float32 a
, float32 b STATUS_PARAM
)
2378 a
= float32_squash_input_denormal(a STATUS_VAR
);
2379 b
= float32_squash_input_denormal(b STATUS_VAR
);
2381 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2382 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2384 float_raise( float_flag_invalid STATUS_VAR
);
2387 aSign
= extractFloat32Sign( a
);
2388 bSign
= extractFloat32Sign( b
);
2389 av
= float32_val(a
);
2390 bv
= float32_val(b
);
2391 if ( aSign
!= bSign
) return aSign
&& ( (uint32_t) ( ( av
| bv
)<<1 ) != 0 );
2392 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
2396 /*----------------------------------------------------------------------------
2397 | Returns 1 if the single-precision floating-point values `a' and `b' cannot
2398 | be compared, and 0 otherwise. The comparison is performed according to the
2399 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2400 *----------------------------------------------------------------------------*/
2402 int float32_unordered( float32 a
, float32 b STATUS_PARAM
)
2404 a
= float32_squash_input_denormal(a STATUS_VAR
);
2405 b
= float32_squash_input_denormal(b STATUS_VAR
);
2407 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2408 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2410 float_raise( float_flag_invalid STATUS_VAR
);
2416 /*----------------------------------------------------------------------------
2417 | Returns 1 if the single-precision floating-point value `a' is equal to
2418 | the corresponding value `b', and 0 otherwise. The comparison is performed
2419 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2420 *----------------------------------------------------------------------------*/
2422 int float32_eq_quiet( float32 a
, float32 b STATUS_PARAM
)
2424 a
= float32_squash_input_denormal(a STATUS_VAR
);
2425 b
= float32_squash_input_denormal(b STATUS_VAR
);
2427 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2428 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2430 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2431 float_raise( float_flag_invalid STATUS_VAR
);
2435 return ( float32_val(a
) == float32_val(b
) ) ||
2436 ( (uint32_t) ( ( float32_val(a
) | float32_val(b
) )<<1 ) == 0 );
2439 /*----------------------------------------------------------------------------
2440 | Returns 1 if the single-precision floating-point value `a' is less than or
2441 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
2442 | cause an exception. Otherwise, the comparison is performed according to the
2443 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2444 *----------------------------------------------------------------------------*/
2446 int float32_le_quiet( float32 a
, float32 b STATUS_PARAM
)
2450 a
= float32_squash_input_denormal(a STATUS_VAR
);
2451 b
= float32_squash_input_denormal(b STATUS_VAR
);
2453 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2454 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2456 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2457 float_raise( float_flag_invalid STATUS_VAR
);
2461 aSign
= extractFloat32Sign( a
);
2462 bSign
= extractFloat32Sign( b
);
2463 av
= float32_val(a
);
2464 bv
= float32_val(b
);
2465 if ( aSign
!= bSign
) return aSign
|| ( (uint32_t) ( ( av
| bv
)<<1 ) == 0 );
2466 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
2470 /*----------------------------------------------------------------------------
2471 | Returns 1 if the single-precision floating-point value `a' is less than
2472 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2473 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
2474 | Standard for Binary Floating-Point Arithmetic.
2475 *----------------------------------------------------------------------------*/
2477 int float32_lt_quiet( float32 a
, float32 b STATUS_PARAM
)
2481 a
= float32_squash_input_denormal(a STATUS_VAR
);
2482 b
= float32_squash_input_denormal(b STATUS_VAR
);
2484 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2485 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2487 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2488 float_raise( float_flag_invalid STATUS_VAR
);
2492 aSign
= extractFloat32Sign( a
);
2493 bSign
= extractFloat32Sign( b
);
2494 av
= float32_val(a
);
2495 bv
= float32_val(b
);
2496 if ( aSign
!= bSign
) return aSign
&& ( (uint32_t) ( ( av
| bv
)<<1 ) != 0 );
2497 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
2501 /*----------------------------------------------------------------------------
2502 | Returns 1 if the single-precision floating-point values `a' and `b' cannot
2503 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
2504 | comparison is performed according to the IEC/IEEE Standard for Binary
2505 | Floating-Point Arithmetic.
2506 *----------------------------------------------------------------------------*/
2508 int float32_unordered_quiet( float32 a
, float32 b STATUS_PARAM
)
2510 a
= float32_squash_input_denormal(a STATUS_VAR
);
2511 b
= float32_squash_input_denormal(b STATUS_VAR
);
2513 if ( ( ( extractFloat32Exp( a
) == 0xFF ) && extractFloat32Frac( a
) )
2514 || ( ( extractFloat32Exp( b
) == 0xFF ) && extractFloat32Frac( b
) )
2516 if ( float32_is_signaling_nan( a
) || float32_is_signaling_nan( b
) ) {
2517 float_raise( float_flag_invalid STATUS_VAR
);
2524 /*----------------------------------------------------------------------------
2525 | Returns the result of converting the double-precision floating-point value
2526 | `a' to the 32-bit two's complement integer format. The conversion is
2527 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2528 | Arithmetic---which means in particular that the conversion is rounded
2529 | according to the current rounding mode. If `a' is a NaN, the largest
2530 | positive integer is returned. Otherwise, if the conversion overflows, the
2531 | largest integer with the same sign as `a' is returned.
2532 *----------------------------------------------------------------------------*/
2534 int32
float64_to_int32( float64 a STATUS_PARAM
)
2537 int16 aExp
, shiftCount
;
2539 a
= float64_squash_input_denormal(a STATUS_VAR
);
2541 aSig
= extractFloat64Frac( a
);
2542 aExp
= extractFloat64Exp( a
);
2543 aSign
= extractFloat64Sign( a
);
2544 if ( ( aExp
== 0x7FF ) && aSig
) aSign
= 0;
2545 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2546 shiftCount
= 0x42C - aExp
;
2547 if ( 0 < shiftCount
) shift64RightJamming( aSig
, shiftCount
, &aSig
);
2548 return roundAndPackInt32( aSign
, aSig STATUS_VAR
);
2552 /*----------------------------------------------------------------------------
2553 | Returns the result of converting the double-precision floating-point value
2554 | `a' to the 32-bit two's complement integer format. The conversion is
2555 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2556 | Arithmetic, except that the conversion is always rounded toward zero.
2557 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2558 | the conversion overflows, the largest integer with the same sign as `a' is
2560 *----------------------------------------------------------------------------*/
2562 int32
float64_to_int32_round_to_zero( float64 a STATUS_PARAM
)
2565 int16 aExp
, shiftCount
;
2566 uint64_t aSig
, savedASig
;
2568 a
= float64_squash_input_denormal(a STATUS_VAR
);
2570 aSig
= extractFloat64Frac( a
);
2571 aExp
= extractFloat64Exp( a
);
2572 aSign
= extractFloat64Sign( a
);
2573 if ( 0x41E < aExp
) {
2574 if ( ( aExp
== 0x7FF ) && aSig
) aSign
= 0;
2577 else if ( aExp
< 0x3FF ) {
2578 if ( aExp
|| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
2581 aSig
|= LIT64( 0x0010000000000000 );
2582 shiftCount
= 0x433 - aExp
;
2584 aSig
>>= shiftCount
;
2586 if ( aSign
) z
= - z
;
2587 if ( ( z
< 0 ) ^ aSign
) {
2589 float_raise( float_flag_invalid STATUS_VAR
);
2590 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
2592 if ( ( aSig
<<shiftCount
) != savedASig
) {
2593 STATUS(float_exception_flags
) |= float_flag_inexact
;
2599 /*----------------------------------------------------------------------------
2600 | Returns the result of converting the double-precision floating-point value
2601 | `a' to the 16-bit two's complement integer format. The conversion is
2602 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2603 | Arithmetic, except that the conversion is always rounded toward zero.
2604 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2605 | the conversion overflows, the largest integer with the same sign as `a' is
2607 *----------------------------------------------------------------------------*/
2609 int16
float64_to_int16_round_to_zero( float64 a STATUS_PARAM
)
2612 int16 aExp
, shiftCount
;
2613 uint64_t aSig
, savedASig
;
2616 aSig
= extractFloat64Frac( a
);
2617 aExp
= extractFloat64Exp( a
);
2618 aSign
= extractFloat64Sign( a
);
2619 if ( 0x40E < aExp
) {
2620 if ( ( aExp
== 0x7FF ) && aSig
) {
2625 else if ( aExp
< 0x3FF ) {
2626 if ( aExp
|| aSig
) {
2627 STATUS(float_exception_flags
) |= float_flag_inexact
;
2631 aSig
|= LIT64( 0x0010000000000000 );
2632 shiftCount
= 0x433 - aExp
;
2634 aSig
>>= shiftCount
;
2639 if ( ( (int16_t)z
< 0 ) ^ aSign
) {
2641 float_raise( float_flag_invalid STATUS_VAR
);
2642 return aSign
? (int32_t) 0xffff8000 : 0x7FFF;
2644 if ( ( aSig
<<shiftCount
) != savedASig
) {
2645 STATUS(float_exception_flags
) |= float_flag_inexact
;
2650 /*----------------------------------------------------------------------------
2651 | Returns the result of converting the double-precision floating-point value
2652 | `a' to the 64-bit two's complement integer format. The conversion is
2653 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2654 | Arithmetic---which means in particular that the conversion is rounded
2655 | according to the current rounding mode. If `a' is a NaN, the largest
2656 | positive integer is returned. Otherwise, if the conversion overflows, the
2657 | largest integer with the same sign as `a' is returned.
2658 *----------------------------------------------------------------------------*/
2660 int64
float64_to_int64( float64 a STATUS_PARAM
)
2663 int16 aExp
, shiftCount
;
2664 uint64_t aSig
, aSigExtra
;
2665 a
= float64_squash_input_denormal(a STATUS_VAR
);
2667 aSig
= extractFloat64Frac( a
);
2668 aExp
= extractFloat64Exp( a
);
2669 aSign
= extractFloat64Sign( a
);
2670 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2671 shiftCount
= 0x433 - aExp
;
2672 if ( shiftCount
<= 0 ) {
2673 if ( 0x43E < aExp
) {
2674 float_raise( float_flag_invalid STATUS_VAR
);
2676 || ( ( aExp
== 0x7FF )
2677 && ( aSig
!= LIT64( 0x0010000000000000 ) ) )
2679 return LIT64( 0x7FFFFFFFFFFFFFFF );
2681 return (int64_t) LIT64( 0x8000000000000000 );
2684 aSig
<<= - shiftCount
;
2687 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
2689 return roundAndPackInt64( aSign
, aSig
, aSigExtra STATUS_VAR
);
2693 /*----------------------------------------------------------------------------
2694 | Returns the result of converting the double-precision floating-point value
2695 | `a' to the 64-bit two's complement integer format. The conversion is
2696 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2697 | Arithmetic, except that the conversion is always rounded toward zero.
2698 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2699 | the conversion overflows, the largest integer with the same sign as `a' is
2701 *----------------------------------------------------------------------------*/
2703 int64
float64_to_int64_round_to_zero( float64 a STATUS_PARAM
)
2706 int16 aExp
, shiftCount
;
2709 a
= float64_squash_input_denormal(a STATUS_VAR
);
2711 aSig
= extractFloat64Frac( a
);
2712 aExp
= extractFloat64Exp( a
);
2713 aSign
= extractFloat64Sign( a
);
2714 if ( aExp
) aSig
|= LIT64( 0x0010000000000000 );
2715 shiftCount
= aExp
- 0x433;
2716 if ( 0 <= shiftCount
) {
2717 if ( 0x43E <= aExp
) {
2718 if ( float64_val(a
) != LIT64( 0xC3E0000000000000 ) ) {
2719 float_raise( float_flag_invalid STATUS_VAR
);
2721 || ( ( aExp
== 0x7FF )
2722 && ( aSig
!= LIT64( 0x0010000000000000 ) ) )
2724 return LIT64( 0x7FFFFFFFFFFFFFFF );
2727 return (int64_t) LIT64( 0x8000000000000000 );
2729 z
= aSig
<<shiftCount
;
2732 if ( aExp
< 0x3FE ) {
2733 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
2736 z
= aSig
>>( - shiftCount
);
2737 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
2738 STATUS(float_exception_flags
) |= float_flag_inexact
;
2741 if ( aSign
) z
= - z
;
2746 /*----------------------------------------------------------------------------
2747 | Returns the result of converting the double-precision floating-point value
2748 | `a' to the single-precision floating-point format. The conversion is
2749 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2751 *----------------------------------------------------------------------------*/
2753 float32
float64_to_float32( float64 a STATUS_PARAM
)
2759 a
= float64_squash_input_denormal(a STATUS_VAR
);
2761 aSig
= extractFloat64Frac( a
);
2762 aExp
= extractFloat64Exp( a
);
2763 aSign
= extractFloat64Sign( a
);
2764 if ( aExp
== 0x7FF ) {
2765 if ( aSig
) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
2766 return packFloat32( aSign
, 0xFF, 0 );
2768 shift64RightJamming( aSig
, 22, &aSig
);
2770 if ( aExp
|| zSig
) {
2774 return roundAndPackFloat32( aSign
, aExp
, zSig STATUS_VAR
);
2779 /*----------------------------------------------------------------------------
2780 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
2781 | half-precision floating-point value, returning the result. After being
2782 | shifted into the proper positions, the three fields are simply added
2783 | together to form the result. This means that any integer portion of `zSig'
2784 | will be added into the exponent. Since a properly normalized significand
2785 | will have an integer portion equal to 1, the `zExp' input should be 1 less
2786 | than the desired result exponent whenever `zSig' is a complete, normalized
2788 *----------------------------------------------------------------------------*/
2789 static float16
packFloat16(flag zSign
, int16 zExp
, uint16_t zSig
)
2791 return make_float16(
2792 (((uint32_t)zSign
) << 15) + (((uint32_t)zExp
) << 10) + zSig
);
2795 /* Half precision floats come in two formats: standard IEEE and "ARM" format.
2796 The latter gains extra exponent range by omitting the NaN/Inf encodings. */
2798 float32
float16_to_float32(float16 a
, flag ieee STATUS_PARAM
)
2804 aSign
= extractFloat16Sign(a
);
2805 aExp
= extractFloat16Exp(a
);
2806 aSig
= extractFloat16Frac(a
);
2808 if (aExp
== 0x1f && ieee
) {
2810 return commonNaNToFloat32(float16ToCommonNaN(a STATUS_VAR
) STATUS_VAR
);
2812 return packFloat32(aSign
, 0xff, aSig
<< 13);
2818 return packFloat32(aSign
, 0, 0);
2821 shiftCount
= countLeadingZeros32( aSig
) - 21;
2822 aSig
= aSig
<< shiftCount
;
2825 return packFloat32( aSign
, aExp
+ 0x70, aSig
<< 13);
2828 float16
float32_to_float16(float32 a
, flag ieee STATUS_PARAM
)
2836 a
= float32_squash_input_denormal(a STATUS_VAR
);
2838 aSig
= extractFloat32Frac( a
);
2839 aExp
= extractFloat32Exp( a
);
2840 aSign
= extractFloat32Sign( a
);
2841 if ( aExp
== 0xFF ) {
2843 /* Input is a NaN */
2844 float16 r
= commonNaNToFloat16( float32ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
2846 return packFloat16(aSign
, 0, 0);
2852 float_raise(float_flag_invalid STATUS_VAR
);
2853 return packFloat16(aSign
, 0x1f, 0x3ff);
2855 return packFloat16(aSign
, 0x1f, 0);
2857 if (aExp
== 0 && aSig
== 0) {
2858 return packFloat16(aSign
, 0, 0);
2860 /* Decimal point between bits 22 and 23. */
2872 float_raise( float_flag_underflow STATUS_VAR
);
2873 roundingMode
= STATUS(float_rounding_mode
);
2874 switch (roundingMode
) {
2875 case float_round_nearest_even
:
2876 increment
= (mask
+ 1) >> 1;
2877 if ((aSig
& mask
) == increment
) {
2878 increment
= aSig
& (increment
<< 1);
2881 case float_round_up
:
2882 increment
= aSign
? 0 : mask
;
2884 case float_round_down
:
2885 increment
= aSign
? mask
: 0;
2887 default: /* round_to_zero */
2892 if (aSig
>= 0x01000000) {
2896 } else if (aExp
< -14
2897 && STATUS(float_detect_tininess
) == float_tininess_before_rounding
) {
2898 float_raise( float_flag_underflow STATUS_VAR
);
2903 float_raise( float_flag_overflow
| float_flag_inexact STATUS_VAR
);
2904 return packFloat16(aSign
, 0x1f, 0);
2908 float_raise(float_flag_invalid
| float_flag_inexact STATUS_VAR
);
2909 return packFloat16(aSign
, 0x1f, 0x3ff);
2913 return packFloat16(aSign
, 0, 0);
2916 aSig
>>= -14 - aExp
;
2919 return packFloat16(aSign
, aExp
+ 14, aSig
>> 13);
2924 /*----------------------------------------------------------------------------
2925 | Returns the result of converting the double-precision floating-point value
2926 | `a' to the extended double-precision floating-point format. The conversion
2927 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2929 *----------------------------------------------------------------------------*/
2931 floatx80
float64_to_floatx80( float64 a STATUS_PARAM
)
2937 a
= float64_squash_input_denormal(a STATUS_VAR
);
2938 aSig
= extractFloat64Frac( a
);
2939 aExp
= extractFloat64Exp( a
);
2940 aSign
= extractFloat64Sign( a
);
2941 if ( aExp
== 0x7FF ) {
2942 if ( aSig
) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
2943 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
2946 if ( aSig
== 0 ) return packFloatx80( aSign
, 0, 0 );
2947 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
2951 aSign
, aExp
+ 0x3C00, ( aSig
| LIT64( 0x0010000000000000 ) )<<11 );
2959 /*----------------------------------------------------------------------------
2960 | Returns the result of converting the double-precision floating-point value
2961 | `a' to the quadruple-precision floating-point format. The conversion is
2962 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2964 *----------------------------------------------------------------------------*/
2966 float128
float64_to_float128( float64 a STATUS_PARAM
)
2970 uint64_t aSig
, zSig0
, zSig1
;
2972 a
= float64_squash_input_denormal(a STATUS_VAR
);
2973 aSig
= extractFloat64Frac( a
);
2974 aExp
= extractFloat64Exp( a
);
2975 aSign
= extractFloat64Sign( a
);
2976 if ( aExp
== 0x7FF ) {
2977 if ( aSig
) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
2978 return packFloat128( aSign
, 0x7FFF, 0, 0 );
2981 if ( aSig
== 0 ) return packFloat128( aSign
, 0, 0, 0 );
2982 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
2985 shift128Right( aSig
, 0, 4, &zSig0
, &zSig1
);
2986 return packFloat128( aSign
, aExp
+ 0x3C00, zSig0
, zSig1
);
2992 /*----------------------------------------------------------------------------
2993 | Rounds the double-precision floating-point value `a' to an integer, and
2994 | returns the result as a double-precision floating-point value. The
2995 | operation is performed according to the IEC/IEEE Standard for Binary
2996 | Floating-Point Arithmetic.
2997 *----------------------------------------------------------------------------*/
2999 float64
float64_round_to_int( float64 a STATUS_PARAM
)
3003 uint64_t lastBitMask
, roundBitsMask
;
3006 a
= float64_squash_input_denormal(a STATUS_VAR
);
3008 aExp
= extractFloat64Exp( a
);
3009 if ( 0x433 <= aExp
) {
3010 if ( ( aExp
== 0x7FF ) && extractFloat64Frac( a
) ) {
3011 return propagateFloat64NaN( a
, a STATUS_VAR
);
3015 if ( aExp
< 0x3FF ) {
3016 if ( (uint64_t) ( float64_val(a
)<<1 ) == 0 ) return a
;
3017 STATUS(float_exception_flags
) |= float_flag_inexact
;
3018 aSign
= extractFloat64Sign( a
);
3019 switch ( STATUS(float_rounding_mode
) ) {
3020 case float_round_nearest_even
:
3021 if ( ( aExp
== 0x3FE ) && extractFloat64Frac( a
) ) {
3022 return packFloat64( aSign
, 0x3FF, 0 );
3025 case float_round_down
:
3026 return make_float64(aSign
? LIT64( 0xBFF0000000000000 ) : 0);
3027 case float_round_up
:
3028 return make_float64(
3029 aSign
? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 ));
3031 return packFloat64( aSign
, 0, 0 );
3034 lastBitMask
<<= 0x433 - aExp
;
3035 roundBitsMask
= lastBitMask
- 1;
3037 roundingMode
= STATUS(float_rounding_mode
);
3038 if ( roundingMode
== float_round_nearest_even
) {
3039 z
+= lastBitMask
>>1;
3040 if ( ( z
& roundBitsMask
) == 0 ) z
&= ~ lastBitMask
;
3042 else if ( roundingMode
!= float_round_to_zero
) {
3043 if ( extractFloat64Sign( make_float64(z
) ) ^ ( roundingMode
== float_round_up
) ) {
3047 z
&= ~ roundBitsMask
;
3048 if ( z
!= float64_val(a
) )
3049 STATUS(float_exception_flags
) |= float_flag_inexact
;
3050 return make_float64(z
);
3054 float64
float64_trunc_to_int( float64 a STATUS_PARAM
)
3058 oldmode
= STATUS(float_rounding_mode
);
3059 STATUS(float_rounding_mode
) = float_round_to_zero
;
3060 res
= float64_round_to_int(a STATUS_VAR
);
3061 STATUS(float_rounding_mode
) = oldmode
;
3065 /*----------------------------------------------------------------------------
3066 | Returns the result of adding the absolute values of the double-precision
3067 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
3068 | before being returned. `zSign' is ignored if the result is a NaN.
3069 | The addition is performed according to the IEC/IEEE Standard for Binary
3070 | Floating-Point Arithmetic.
3071 *----------------------------------------------------------------------------*/
3073 static float64
addFloat64Sigs( float64 a
, float64 b
, flag zSign STATUS_PARAM
)
3075 int16 aExp
, bExp
, zExp
;
3076 uint64_t aSig
, bSig
, zSig
;
3079 aSig
= extractFloat64Frac( a
);
3080 aExp
= extractFloat64Exp( a
);
3081 bSig
= extractFloat64Frac( b
);
3082 bExp
= extractFloat64Exp( b
);
3083 expDiff
= aExp
- bExp
;
3086 if ( 0 < expDiff
) {
3087 if ( aExp
== 0x7FF ) {
3088 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3095 bSig
|= LIT64( 0x2000000000000000 );
3097 shift64RightJamming( bSig
, expDiff
, &bSig
);
3100 else if ( expDiff
< 0 ) {
3101 if ( bExp
== 0x7FF ) {
3102 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3103 return packFloat64( zSign
, 0x7FF, 0 );
3109 aSig
|= LIT64( 0x2000000000000000 );
3111 shift64RightJamming( aSig
, - expDiff
, &aSig
);
3115 if ( aExp
== 0x7FF ) {
3116 if ( aSig
| bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3120 if ( STATUS(flush_to_zero
) ) return packFloat64( zSign
, 0, 0 );
3121 return packFloat64( zSign
, 0, ( aSig
+ bSig
)>>9 );
3123 zSig
= LIT64( 0x4000000000000000 ) + aSig
+ bSig
;
3127 aSig
|= LIT64( 0x2000000000000000 );
3128 zSig
= ( aSig
+ bSig
)<<1;
3130 if ( (int64_t) zSig
< 0 ) {
3135 return roundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3139 /*----------------------------------------------------------------------------
3140 | Returns the result of subtracting the absolute values of the double-
3141 | precision floating-point values `a' and `b'. If `zSign' is 1, the
3142 | difference is negated before being returned. `zSign' is ignored if the
3143 | result is a NaN. The subtraction is performed according to the IEC/IEEE
3144 | Standard for Binary Floating-Point Arithmetic.
3145 *----------------------------------------------------------------------------*/
3147 static float64
subFloat64Sigs( float64 a
, float64 b
, flag zSign STATUS_PARAM
)
3149 int16 aExp
, bExp
, zExp
;
3150 uint64_t aSig
, bSig
, zSig
;
3153 aSig
= extractFloat64Frac( a
);
3154 aExp
= extractFloat64Exp( a
);
3155 bSig
= extractFloat64Frac( b
);
3156 bExp
= extractFloat64Exp( b
);
3157 expDiff
= aExp
- bExp
;
3160 if ( 0 < expDiff
) goto aExpBigger
;
3161 if ( expDiff
< 0 ) goto bExpBigger
;
3162 if ( aExp
== 0x7FF ) {
3163 if ( aSig
| bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3164 float_raise( float_flag_invalid STATUS_VAR
);
3165 return float64_default_nan
;
3171 if ( bSig
< aSig
) goto aBigger
;
3172 if ( aSig
< bSig
) goto bBigger
;
3173 return packFloat64( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
3175 if ( bExp
== 0x7FF ) {
3176 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3177 return packFloat64( zSign
^ 1, 0x7FF, 0 );
3183 aSig
|= LIT64( 0x4000000000000000 );
3185 shift64RightJamming( aSig
, - expDiff
, &aSig
);
3186 bSig
|= LIT64( 0x4000000000000000 );
3191 goto normalizeRoundAndPack
;
3193 if ( aExp
== 0x7FF ) {
3194 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3201 bSig
|= LIT64( 0x4000000000000000 );
3203 shift64RightJamming( bSig
, expDiff
, &bSig
);
3204 aSig
|= LIT64( 0x4000000000000000 );
3208 normalizeRoundAndPack
:
3210 return normalizeRoundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3214 /*----------------------------------------------------------------------------
3215 | Returns the result of adding the double-precision floating-point values `a'
3216 | and `b'. The operation is performed according to the IEC/IEEE Standard for
3217 | Binary Floating-Point Arithmetic.
3218 *----------------------------------------------------------------------------*/
3220 float64
float64_add( float64 a
, float64 b STATUS_PARAM
)
3223 a
= float64_squash_input_denormal(a STATUS_VAR
);
3224 b
= float64_squash_input_denormal(b STATUS_VAR
);
3226 aSign
= extractFloat64Sign( a
);
3227 bSign
= extractFloat64Sign( b
);
3228 if ( aSign
== bSign
) {
3229 return addFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3232 return subFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3237 /*----------------------------------------------------------------------------
3238 | Returns the result of subtracting the double-precision floating-point values
3239 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
3240 | for Binary Floating-Point Arithmetic.
3241 *----------------------------------------------------------------------------*/
3243 float64
float64_sub( float64 a
, float64 b STATUS_PARAM
)
3246 a
= float64_squash_input_denormal(a STATUS_VAR
);
3247 b
= float64_squash_input_denormal(b STATUS_VAR
);
3249 aSign
= extractFloat64Sign( a
);
3250 bSign
= extractFloat64Sign( b
);
3251 if ( aSign
== bSign
) {
3252 return subFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3255 return addFloat64Sigs( a
, b
, aSign STATUS_VAR
);
3260 /*----------------------------------------------------------------------------
3261 | Returns the result of multiplying the double-precision floating-point values
3262 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
3263 | for Binary Floating-Point Arithmetic.
3264 *----------------------------------------------------------------------------*/
3266 float64
float64_mul( float64 a
, float64 b STATUS_PARAM
)
3268 flag aSign
, bSign
, zSign
;
3269 int16 aExp
, bExp
, zExp
;
3270 uint64_t aSig
, bSig
, zSig0
, zSig1
;
3272 a
= float64_squash_input_denormal(a STATUS_VAR
);
3273 b
= float64_squash_input_denormal(b STATUS_VAR
);
3275 aSig
= extractFloat64Frac( a
);
3276 aExp
= extractFloat64Exp( a
);
3277 aSign
= extractFloat64Sign( a
);
3278 bSig
= extractFloat64Frac( b
);
3279 bExp
= extractFloat64Exp( b
);
3280 bSign
= extractFloat64Sign( b
);
3281 zSign
= aSign
^ bSign
;
3282 if ( aExp
== 0x7FF ) {
3283 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
3284 return propagateFloat64NaN( a
, b STATUS_VAR
);
3286 if ( ( bExp
| bSig
) == 0 ) {
3287 float_raise( float_flag_invalid STATUS_VAR
);
3288 return float64_default_nan
;
3290 return packFloat64( zSign
, 0x7FF, 0 );
3292 if ( bExp
== 0x7FF ) {
3293 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3294 if ( ( aExp
| aSig
) == 0 ) {
3295 float_raise( float_flag_invalid STATUS_VAR
);
3296 return float64_default_nan
;
3298 return packFloat64( zSign
, 0x7FF, 0 );
3301 if ( aSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3302 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3305 if ( bSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3306 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3308 zExp
= aExp
+ bExp
- 0x3FF;
3309 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<10;
3310 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3311 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
3312 zSig0
|= ( zSig1
!= 0 );
3313 if ( 0 <= (int64_t) ( zSig0
<<1 ) ) {
3317 return roundAndPackFloat64( zSign
, zExp
, zSig0 STATUS_VAR
);
3321 /*----------------------------------------------------------------------------
3322 | Returns the result of dividing the double-precision floating-point value `a'
3323 | by the corresponding value `b'. The operation is performed according to
3324 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3325 *----------------------------------------------------------------------------*/
3327 float64
float64_div( float64 a
, float64 b STATUS_PARAM
)
3329 flag aSign
, bSign
, zSign
;
3330 int16 aExp
, bExp
, zExp
;
3331 uint64_t aSig
, bSig
, zSig
;
3332 uint64_t rem0
, rem1
;
3333 uint64_t term0
, term1
;
3334 a
= float64_squash_input_denormal(a STATUS_VAR
);
3335 b
= float64_squash_input_denormal(b STATUS_VAR
);
3337 aSig
= extractFloat64Frac( a
);
3338 aExp
= extractFloat64Exp( a
);
3339 aSign
= extractFloat64Sign( a
);
3340 bSig
= extractFloat64Frac( b
);
3341 bExp
= extractFloat64Exp( b
);
3342 bSign
= extractFloat64Sign( b
);
3343 zSign
= aSign
^ bSign
;
3344 if ( aExp
== 0x7FF ) {
3345 if ( aSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3346 if ( bExp
== 0x7FF ) {
3347 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3348 float_raise( float_flag_invalid STATUS_VAR
);
3349 return float64_default_nan
;
3351 return packFloat64( zSign
, 0x7FF, 0 );
3353 if ( bExp
== 0x7FF ) {
3354 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3355 return packFloat64( zSign
, 0, 0 );
3359 if ( ( aExp
| aSig
) == 0 ) {
3360 float_raise( float_flag_invalid STATUS_VAR
);
3361 return float64_default_nan
;
3363 float_raise( float_flag_divbyzero STATUS_VAR
);
3364 return packFloat64( zSign
, 0x7FF, 0 );
3366 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3369 if ( aSig
== 0 ) return packFloat64( zSign
, 0, 0 );
3370 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3372 zExp
= aExp
- bExp
+ 0x3FD;
3373 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<10;
3374 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3375 if ( bSig
<= ( aSig
+ aSig
) ) {
3379 zSig
= estimateDiv128To64( aSig
, 0, bSig
);
3380 if ( ( zSig
& 0x1FF ) <= 2 ) {
3381 mul64To128( bSig
, zSig
, &term0
, &term1
);
3382 sub128( aSig
, 0, term0
, term1
, &rem0
, &rem1
);
3383 while ( (int64_t) rem0
< 0 ) {
3385 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
3387 zSig
|= ( rem1
!= 0 );
3389 return roundAndPackFloat64( zSign
, zExp
, zSig STATUS_VAR
);
3393 /*----------------------------------------------------------------------------
3394 | Returns the remainder of the double-precision floating-point value `a'
3395 | with respect to the corresponding value `b'. The operation is performed
3396 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3397 *----------------------------------------------------------------------------*/
3399 float64
float64_rem( float64 a
, float64 b STATUS_PARAM
)
3402 int16 aExp
, bExp
, expDiff
;
3403 uint64_t aSig
, bSig
;
3404 uint64_t q
, alternateASig
;
3407 a
= float64_squash_input_denormal(a STATUS_VAR
);
3408 b
= float64_squash_input_denormal(b STATUS_VAR
);
3409 aSig
= extractFloat64Frac( a
);
3410 aExp
= extractFloat64Exp( a
);
3411 aSign
= extractFloat64Sign( a
);
3412 bSig
= extractFloat64Frac( b
);
3413 bExp
= extractFloat64Exp( b
);
3414 if ( aExp
== 0x7FF ) {
3415 if ( aSig
|| ( ( bExp
== 0x7FF ) && bSig
) ) {
3416 return propagateFloat64NaN( a
, b STATUS_VAR
);
3418 float_raise( float_flag_invalid STATUS_VAR
);
3419 return float64_default_nan
;
3421 if ( bExp
== 0x7FF ) {
3422 if ( bSig
) return propagateFloat64NaN( a
, b STATUS_VAR
);
3427 float_raise( float_flag_invalid STATUS_VAR
);
3428 return float64_default_nan
;
3430 normalizeFloat64Subnormal( bSig
, &bExp
, &bSig
);
3433 if ( aSig
== 0 ) return a
;
3434 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3436 expDiff
= aExp
- bExp
;
3437 aSig
= ( aSig
| LIT64( 0x0010000000000000 ) )<<11;
3438 bSig
= ( bSig
| LIT64( 0x0010000000000000 ) )<<11;
3439 if ( expDiff
< 0 ) {
3440 if ( expDiff
< -1 ) return a
;
3443 q
= ( bSig
<= aSig
);
3444 if ( q
) aSig
-= bSig
;
3446 while ( 0 < expDiff
) {
3447 q
= estimateDiv128To64( aSig
, 0, bSig
);
3448 q
= ( 2 < q
) ? q
- 2 : 0;
3449 aSig
= - ( ( bSig
>>2 ) * q
);
3453 if ( 0 < expDiff
) {
3454 q
= estimateDiv128To64( aSig
, 0, bSig
);
3455 q
= ( 2 < q
) ? q
- 2 : 0;
3458 aSig
= ( ( aSig
>>1 )<<( expDiff
- 1 ) ) - bSig
* q
;
3465 alternateASig
= aSig
;
3468 } while ( 0 <= (int64_t) aSig
);
3469 sigMean
= aSig
+ alternateASig
;
3470 if ( ( sigMean
< 0 ) || ( ( sigMean
== 0 ) && ( q
& 1 ) ) ) {
3471 aSig
= alternateASig
;
3473 zSign
= ( (int64_t) aSig
< 0 );
3474 if ( zSign
) aSig
= - aSig
;
3475 return normalizeRoundAndPackFloat64( aSign
^ zSign
, bExp
, aSig STATUS_VAR
);
3479 /*----------------------------------------------------------------------------
3480 | Returns the square root of the double-precision floating-point value `a'.
3481 | The operation is performed according to the IEC/IEEE Standard for Binary
3482 | Floating-Point Arithmetic.
3483 *----------------------------------------------------------------------------*/
3485 float64
float64_sqrt( float64 a STATUS_PARAM
)
3489 uint64_t aSig
, zSig
, doubleZSig
;
3490 uint64_t rem0
, rem1
, term0
, term1
;
3491 a
= float64_squash_input_denormal(a STATUS_VAR
);
3493 aSig
= extractFloat64Frac( a
);
3494 aExp
= extractFloat64Exp( a
);
3495 aSign
= extractFloat64Sign( a
);
3496 if ( aExp
== 0x7FF ) {
3497 if ( aSig
) return propagateFloat64NaN( a
, a STATUS_VAR
);
3498 if ( ! aSign
) return a
;
3499 float_raise( float_flag_invalid STATUS_VAR
);
3500 return float64_default_nan
;
3503 if ( ( aExp
| aSig
) == 0 ) return a
;
3504 float_raise( float_flag_invalid STATUS_VAR
);
3505 return float64_default_nan
;
3508 if ( aSig
== 0 ) return float64_zero
;
3509 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3511 zExp
= ( ( aExp
- 0x3FF )>>1 ) + 0x3FE;
3512 aSig
|= LIT64( 0x0010000000000000 );
3513 zSig
= estimateSqrt32( aExp
, aSig
>>21 );
3514 aSig
<<= 9 - ( aExp
& 1 );
3515 zSig
= estimateDiv128To64( aSig
, 0, zSig
<<32 ) + ( zSig
<<30 );
3516 if ( ( zSig
& 0x1FF ) <= 5 ) {
3517 doubleZSig
= zSig
<<1;
3518 mul64To128( zSig
, zSig
, &term0
, &term1
);
3519 sub128( aSig
, 0, term0
, term1
, &rem0
, &rem1
);
3520 while ( (int64_t) rem0
< 0 ) {
3523 add128( rem0
, rem1
, zSig
>>63, doubleZSig
| 1, &rem0
, &rem1
);
3525 zSig
|= ( ( rem0
| rem1
) != 0 );
3527 return roundAndPackFloat64( 0, zExp
, zSig STATUS_VAR
);
3531 /*----------------------------------------------------------------------------
3532 | Returns the binary log of the double-precision floating-point value `a'.
3533 | The operation is performed according to the IEC/IEEE Standard for Binary
3534 | Floating-Point Arithmetic.
3535 *----------------------------------------------------------------------------*/
3536 float64
float64_log2( float64 a STATUS_PARAM
)
3540 uint64_t aSig
, aSig0
, aSig1
, zSig
, i
;
3541 a
= float64_squash_input_denormal(a STATUS_VAR
);
3543 aSig
= extractFloat64Frac( a
);
3544 aExp
= extractFloat64Exp( a
);
3545 aSign
= extractFloat64Sign( a
);
3548 if ( aSig
== 0 ) return packFloat64( 1, 0x7FF, 0 );
3549 normalizeFloat64Subnormal( aSig
, &aExp
, &aSig
);
3552 float_raise( float_flag_invalid STATUS_VAR
);
3553 return float64_default_nan
;
3555 if ( aExp
== 0x7FF ) {
3556 if ( aSig
) return propagateFloat64NaN( a
, float64_zero STATUS_VAR
);
3561 aSig
|= LIT64( 0x0010000000000000 );
3563 zSig
= (uint64_t)aExp
<< 52;
3564 for (i
= 1LL << 51; i
> 0; i
>>= 1) {
3565 mul64To128( aSig
, aSig
, &aSig0
, &aSig1
);
3566 aSig
= ( aSig0
<< 12 ) | ( aSig1
>> 52 );
3567 if ( aSig
& LIT64( 0x0020000000000000 ) ) {
3575 return normalizeRoundAndPackFloat64( zSign
, 0x408, zSig STATUS_VAR
);
3578 /*----------------------------------------------------------------------------
3579 | Returns 1 if the double-precision floating-point value `a' is equal to the
3580 | corresponding value `b', and 0 otherwise. The invalid exception is raised
3581 | if either operand is a NaN. Otherwise, the comparison is performed
3582 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3583 *----------------------------------------------------------------------------*/
3585 int float64_eq( float64 a
, float64 b STATUS_PARAM
)
3588 a
= float64_squash_input_denormal(a STATUS_VAR
);
3589 b
= float64_squash_input_denormal(b STATUS_VAR
);
3591 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3592 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3594 float_raise( float_flag_invalid STATUS_VAR
);
3597 av
= float64_val(a
);
3598 bv
= float64_val(b
);
3599 return ( av
== bv
) || ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
3603 /*----------------------------------------------------------------------------
3604 | Returns 1 if the double-precision floating-point value `a' is less than or
3605 | equal to the corresponding value `b', and 0 otherwise. The comparison is
3606 | performed according to the IEC/IEEE Standard for Binary Floating-Point
3608 *----------------------------------------------------------------------------*/
3610 int float64_le( float64 a
, float64 b STATUS_PARAM
)
3614 a
= float64_squash_input_denormal(a STATUS_VAR
);
3615 b
= float64_squash_input_denormal(b STATUS_VAR
);
3617 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3618 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3620 float_raise( float_flag_invalid STATUS_VAR
);
3623 aSign
= extractFloat64Sign( a
);
3624 bSign
= extractFloat64Sign( b
);
3625 av
= float64_val(a
);
3626 bv
= float64_val(b
);
3627 if ( aSign
!= bSign
) return aSign
|| ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
3628 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
3632 /*----------------------------------------------------------------------------
3633 | Returns 1 if the double-precision floating-point value `a' is less than
3634 | the corresponding value `b', and 0 otherwise. The comparison is performed
3635 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3636 *----------------------------------------------------------------------------*/
3638 int float64_lt( float64 a
, float64 b STATUS_PARAM
)
3643 a
= float64_squash_input_denormal(a STATUS_VAR
);
3644 b
= float64_squash_input_denormal(b STATUS_VAR
);
3645 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3646 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3648 float_raise( float_flag_invalid STATUS_VAR
);
3651 aSign
= extractFloat64Sign( a
);
3652 bSign
= extractFloat64Sign( b
);
3653 av
= float64_val(a
);
3654 bv
= float64_val(b
);
3655 if ( aSign
!= bSign
) return aSign
&& ( (uint64_t) ( ( av
| bv
)<<1 ) != 0 );
3656 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
3660 /*----------------------------------------------------------------------------
3661 | Returns 1 if the double-precision floating-point values `a' and `b' cannot
3662 | be compared, and 0 otherwise. The comparison is performed according to the
3663 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3664 *----------------------------------------------------------------------------*/
3666 int float64_unordered( float64 a
, float64 b STATUS_PARAM
)
3668 a
= float64_squash_input_denormal(a STATUS_VAR
);
3669 b
= float64_squash_input_denormal(b STATUS_VAR
);
3671 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3672 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3674 float_raise( float_flag_invalid STATUS_VAR
);
3680 /*----------------------------------------------------------------------------
3681 | Returns 1 if the double-precision floating-point value `a' is equal to the
3682 | corresponding value `b', and 0 otherwise. The comparison is performed
3683 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3684 *----------------------------------------------------------------------------*/
3686 int float64_eq_quiet( float64 a
, float64 b STATUS_PARAM
)
3689 a
= float64_squash_input_denormal(a STATUS_VAR
);
3690 b
= float64_squash_input_denormal(b STATUS_VAR
);
3692 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3693 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3695 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
3696 float_raise( float_flag_invalid STATUS_VAR
);
3700 av
= float64_val(a
);
3701 bv
= float64_val(b
);
3702 return ( av
== bv
) || ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
3706 /*----------------------------------------------------------------------------
3707 | Returns 1 if the double-precision floating-point value `a' is less than or
3708 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
3709 | cause an exception. Otherwise, the comparison is performed according to the
3710 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3711 *----------------------------------------------------------------------------*/
3713 int float64_le_quiet( float64 a
, float64 b STATUS_PARAM
)
3717 a
= float64_squash_input_denormal(a STATUS_VAR
);
3718 b
= float64_squash_input_denormal(b STATUS_VAR
);
3720 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3721 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3723 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
3724 float_raise( float_flag_invalid STATUS_VAR
);
3728 aSign
= extractFloat64Sign( a
);
3729 bSign
= extractFloat64Sign( b
);
3730 av
= float64_val(a
);
3731 bv
= float64_val(b
);
3732 if ( aSign
!= bSign
) return aSign
|| ( (uint64_t) ( ( av
| bv
)<<1 ) == 0 );
3733 return ( av
== bv
) || ( aSign
^ ( av
< bv
) );
3737 /*----------------------------------------------------------------------------
3738 | Returns 1 if the double-precision floating-point value `a' is less than
3739 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
3740 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
3741 | Standard for Binary Floating-Point Arithmetic.
3742 *----------------------------------------------------------------------------*/
3744 int float64_lt_quiet( float64 a
, float64 b STATUS_PARAM
)
3748 a
= float64_squash_input_denormal(a STATUS_VAR
);
3749 b
= float64_squash_input_denormal(b STATUS_VAR
);
3751 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3752 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3754 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
3755 float_raise( float_flag_invalid STATUS_VAR
);
3759 aSign
= extractFloat64Sign( a
);
3760 bSign
= extractFloat64Sign( b
);
3761 av
= float64_val(a
);
3762 bv
= float64_val(b
);
3763 if ( aSign
!= bSign
) return aSign
&& ( (uint64_t) ( ( av
| bv
)<<1 ) != 0 );
3764 return ( av
!= bv
) && ( aSign
^ ( av
< bv
) );
3768 /*----------------------------------------------------------------------------
3769 | Returns 1 if the double-precision floating-point values `a' and `b' cannot
3770 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
3771 | comparison is performed according to the IEC/IEEE Standard for Binary
3772 | Floating-Point Arithmetic.
3773 *----------------------------------------------------------------------------*/
3775 int float64_unordered_quiet( float64 a
, float64 b STATUS_PARAM
)
3777 a
= float64_squash_input_denormal(a STATUS_VAR
);
3778 b
= float64_squash_input_denormal(b STATUS_VAR
);
3780 if ( ( ( extractFloat64Exp( a
) == 0x7FF ) && extractFloat64Frac( a
) )
3781 || ( ( extractFloat64Exp( b
) == 0x7FF ) && extractFloat64Frac( b
) )
3783 if ( float64_is_signaling_nan( a
) || float64_is_signaling_nan( b
) ) {
3784 float_raise( float_flag_invalid STATUS_VAR
);
3793 /*----------------------------------------------------------------------------
3794 | Returns the result of converting the extended double-precision floating-
3795 | point value `a' to the 32-bit two's complement integer format. The
3796 | conversion is performed according to the IEC/IEEE Standard for Binary
3797 | Floating-Point Arithmetic---which means in particular that the conversion
3798 | is rounded according to the current rounding mode. If `a' is a NaN, the
3799 | largest positive integer is returned. Otherwise, if the conversion
3800 | overflows, the largest integer with the same sign as `a' is returned.
3801 *----------------------------------------------------------------------------*/
3803 int32
floatx80_to_int32( floatx80 a STATUS_PARAM
)
3806 int32 aExp
, shiftCount
;
3809 aSig
= extractFloatx80Frac( a
);
3810 aExp
= extractFloatx80Exp( a
);
3811 aSign
= extractFloatx80Sign( a
);
3812 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
3813 shiftCount
= 0x4037 - aExp
;
3814 if ( shiftCount
<= 0 ) shiftCount
= 1;
3815 shift64RightJamming( aSig
, shiftCount
, &aSig
);
3816 return roundAndPackInt32( aSign
, aSig STATUS_VAR
);
3820 /*----------------------------------------------------------------------------
3821 | Returns the result of converting the extended double-precision floating-
3822 | point value `a' to the 32-bit two's complement integer format. The
3823 | conversion is performed according to the IEC/IEEE Standard for Binary
3824 | Floating-Point Arithmetic, except that the conversion is always rounded
3825 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3826 | Otherwise, if the conversion overflows, the largest integer with the same
3827 | sign as `a' is returned.
3828 *----------------------------------------------------------------------------*/
3830 int32
floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM
)
3833 int32 aExp
, shiftCount
;
3834 uint64_t aSig
, savedASig
;
3837 aSig
= extractFloatx80Frac( a
);
3838 aExp
= extractFloatx80Exp( a
);
3839 aSign
= extractFloatx80Sign( a
);
3840 if ( 0x401E < aExp
) {
3841 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) aSign
= 0;
3844 else if ( aExp
< 0x3FFF ) {
3845 if ( aExp
|| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
3848 shiftCount
= 0x403E - aExp
;
3850 aSig
>>= shiftCount
;
3852 if ( aSign
) z
= - z
;
3853 if ( ( z
< 0 ) ^ aSign
) {
3855 float_raise( float_flag_invalid STATUS_VAR
);
3856 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
3858 if ( ( aSig
<<shiftCount
) != savedASig
) {
3859 STATUS(float_exception_flags
) |= float_flag_inexact
;
3865 /*----------------------------------------------------------------------------
3866 | Returns the result of converting the extended double-precision floating-
3867 | point value `a' to the 64-bit two's complement integer format. The
3868 | conversion is performed according to the IEC/IEEE Standard for Binary
3869 | Floating-Point Arithmetic---which means in particular that the conversion
3870 | is rounded according to the current rounding mode. If `a' is a NaN,
3871 | the largest positive integer is returned. Otherwise, if the conversion
3872 | overflows, the largest integer with the same sign as `a' is returned.
3873 *----------------------------------------------------------------------------*/
3875 int64
floatx80_to_int64( floatx80 a STATUS_PARAM
)
3878 int32 aExp
, shiftCount
;
3879 uint64_t aSig
, aSigExtra
;
3881 aSig
= extractFloatx80Frac( a
);
3882 aExp
= extractFloatx80Exp( a
);
3883 aSign
= extractFloatx80Sign( a
);
3884 shiftCount
= 0x403E - aExp
;
3885 if ( shiftCount
<= 0 ) {
3887 float_raise( float_flag_invalid STATUS_VAR
);
3889 || ( ( aExp
== 0x7FFF )
3890 && ( aSig
!= LIT64( 0x8000000000000000 ) ) )
3892 return LIT64( 0x7FFFFFFFFFFFFFFF );
3894 return (int64_t) LIT64( 0x8000000000000000 );
3899 shift64ExtraRightJamming( aSig
, 0, shiftCount
, &aSig
, &aSigExtra
);
3901 return roundAndPackInt64( aSign
, aSig
, aSigExtra STATUS_VAR
);
3905 /*----------------------------------------------------------------------------
3906 | Returns the result of converting the extended double-precision floating-
3907 | point value `a' to the 64-bit two's complement integer format. The
3908 | conversion is performed according to the IEC/IEEE Standard for Binary
3909 | Floating-Point Arithmetic, except that the conversion is always rounded
3910 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3911 | Otherwise, if the conversion overflows, the largest integer with the same
3912 | sign as `a' is returned.
3913 *----------------------------------------------------------------------------*/
3915 int64
floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM
)
3918 int32 aExp
, shiftCount
;
3922 aSig
= extractFloatx80Frac( a
);
3923 aExp
= extractFloatx80Exp( a
);
3924 aSign
= extractFloatx80Sign( a
);
3925 shiftCount
= aExp
- 0x403E;
3926 if ( 0 <= shiftCount
) {
3927 aSig
&= LIT64( 0x7FFFFFFFFFFFFFFF );
3928 if ( ( a
.high
!= 0xC03E ) || aSig
) {
3929 float_raise( float_flag_invalid STATUS_VAR
);
3930 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && aSig
) ) {
3931 return LIT64( 0x7FFFFFFFFFFFFFFF );
3934 return (int64_t) LIT64( 0x8000000000000000 );
3936 else if ( aExp
< 0x3FFF ) {
3937 if ( aExp
| aSig
) STATUS(float_exception_flags
) |= float_flag_inexact
;
3940 z
= aSig
>>( - shiftCount
);
3941 if ( (uint64_t) ( aSig
<<( shiftCount
& 63 ) ) ) {
3942 STATUS(float_exception_flags
) |= float_flag_inexact
;
3944 if ( aSign
) z
= - z
;
3949 /*----------------------------------------------------------------------------
3950 | Returns the result of converting the extended double-precision floating-
3951 | point value `a' to the single-precision floating-point format. The
3952 | conversion is performed according to the IEC/IEEE Standard for Binary
3953 | Floating-Point Arithmetic.
3954 *----------------------------------------------------------------------------*/
3956 float32
floatx80_to_float32( floatx80 a STATUS_PARAM
)
3962 aSig
= extractFloatx80Frac( a
);
3963 aExp
= extractFloatx80Exp( a
);
3964 aSign
= extractFloatx80Sign( a
);
3965 if ( aExp
== 0x7FFF ) {
3966 if ( (uint64_t) ( aSig
<<1 ) ) {
3967 return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
3969 return packFloat32( aSign
, 0xFF, 0 );
3971 shift64RightJamming( aSig
, 33, &aSig
);
3972 if ( aExp
|| aSig
) aExp
-= 0x3F81;
3973 return roundAndPackFloat32( aSign
, aExp
, aSig STATUS_VAR
);
3977 /*----------------------------------------------------------------------------
3978 | Returns the result of converting the extended double-precision floating-
3979 | point value `a' to the double-precision floating-point format. The
3980 | conversion is performed according to the IEC/IEEE Standard for Binary
3981 | Floating-Point Arithmetic.
3982 *----------------------------------------------------------------------------*/
3984 float64
floatx80_to_float64( floatx80 a STATUS_PARAM
)
3988 uint64_t aSig
, zSig
;
3990 aSig
= extractFloatx80Frac( a
);
3991 aExp
= extractFloatx80Exp( a
);
3992 aSign
= extractFloatx80Sign( a
);
3993 if ( aExp
== 0x7FFF ) {
3994 if ( (uint64_t) ( aSig
<<1 ) ) {
3995 return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
3997 return packFloat64( aSign
, 0x7FF, 0 );
3999 shift64RightJamming( aSig
, 1, &zSig
);
4000 if ( aExp
|| aSig
) aExp
-= 0x3C01;
4001 return roundAndPackFloat64( aSign
, aExp
, zSig STATUS_VAR
);
4007 /*----------------------------------------------------------------------------
4008 | Returns the result of converting the extended double-precision floating-
4009 | point value `a' to the quadruple-precision floating-point format. The
4010 | conversion is performed according to the IEC/IEEE Standard for Binary
4011 | Floating-Point Arithmetic.
4012 *----------------------------------------------------------------------------*/
4014 float128
floatx80_to_float128( floatx80 a STATUS_PARAM
)
4018 uint64_t aSig
, zSig0
, zSig1
;
4020 aSig
= extractFloatx80Frac( a
);
4021 aExp
= extractFloatx80Exp( a
);
4022 aSign
= extractFloatx80Sign( a
);
4023 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( aSig
<<1 ) ) {
4024 return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
4026 shift128Right( aSig
<<1, 0, 16, &zSig0
, &zSig1
);
4027 return packFloat128( aSign
, aExp
, zSig0
, zSig1
);
4033 /*----------------------------------------------------------------------------
4034 | Rounds the extended double-precision floating-point value `a' to an integer,
4035 | and returns the result as an extended quadruple-precision floating-point
4036 | value. The operation is performed according to the IEC/IEEE Standard for
4037 | Binary Floating-Point Arithmetic.
4038 *----------------------------------------------------------------------------*/
4040 floatx80
floatx80_round_to_int( floatx80 a STATUS_PARAM
)
4044 uint64_t lastBitMask
, roundBitsMask
;
4048 aExp
= extractFloatx80Exp( a
);
4049 if ( 0x403E <= aExp
) {
4050 if ( ( aExp
== 0x7FFF ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) ) {
4051 return propagateFloatx80NaN( a
, a STATUS_VAR
);
4055 if ( aExp
< 0x3FFF ) {
4057 && ( (uint64_t) ( extractFloatx80Frac( a
)<<1 ) == 0 ) ) {
4060 STATUS(float_exception_flags
) |= float_flag_inexact
;
4061 aSign
= extractFloatx80Sign( a
);
4062 switch ( STATUS(float_rounding_mode
) ) {
4063 case float_round_nearest_even
:
4064 if ( ( aExp
== 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a
)<<1 )
4067 packFloatx80( aSign
, 0x3FFF, LIT64( 0x8000000000000000 ) );
4070 case float_round_down
:
4073 packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
4074 : packFloatx80( 0, 0, 0 );
4075 case float_round_up
:
4077 aSign
? packFloatx80( 1, 0, 0 )
4078 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
4080 return packFloatx80( aSign
, 0, 0 );
4083 lastBitMask
<<= 0x403E - aExp
;
4084 roundBitsMask
= lastBitMask
- 1;
4086 roundingMode
= STATUS(float_rounding_mode
);
4087 if ( roundingMode
== float_round_nearest_even
) {
4088 z
.low
+= lastBitMask
>>1;
4089 if ( ( z
.low
& roundBitsMask
) == 0 ) z
.low
&= ~ lastBitMask
;
4091 else if ( roundingMode
!= float_round_to_zero
) {
4092 if ( extractFloatx80Sign( z
) ^ ( roundingMode
== float_round_up
) ) {
4093 z
.low
+= roundBitsMask
;
4096 z
.low
&= ~ roundBitsMask
;
4099 z
.low
= LIT64( 0x8000000000000000 );
4101 if ( z
.low
!= a
.low
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4106 /*----------------------------------------------------------------------------
4107 | Returns the result of adding the absolute values of the extended double-
4108 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
4109 | negated before being returned. `zSign' is ignored if the result is a NaN.
4110 | The addition is performed according to the IEC/IEEE Standard for Binary
4111 | Floating-Point Arithmetic.
4112 *----------------------------------------------------------------------------*/
4114 static floatx80
addFloatx80Sigs( floatx80 a
, floatx80 b
, flag zSign STATUS_PARAM
)
4116 int32 aExp
, bExp
, zExp
;
4117 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4120 aSig
= extractFloatx80Frac( a
);
4121 aExp
= extractFloatx80Exp( a
);
4122 bSig
= extractFloatx80Frac( b
);
4123 bExp
= extractFloatx80Exp( b
);
4124 expDiff
= aExp
- bExp
;
4125 if ( 0 < expDiff
) {
4126 if ( aExp
== 0x7FFF ) {
4127 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4130 if ( bExp
== 0 ) --expDiff
;
4131 shift64ExtraRightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
4134 else if ( expDiff
< 0 ) {
4135 if ( bExp
== 0x7FFF ) {
4136 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4137 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4139 if ( aExp
== 0 ) ++expDiff
;
4140 shift64ExtraRightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
4144 if ( aExp
== 0x7FFF ) {
4145 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
4146 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4151 zSig0
= aSig
+ bSig
;
4153 normalizeFloatx80Subnormal( zSig0
, &zExp
, &zSig0
);
4159 zSig0
= aSig
+ bSig
;
4160 if ( (int64_t) zSig0
< 0 ) goto roundAndPack
;
4162 shift64ExtraRightJamming( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
4163 zSig0
|= LIT64( 0x8000000000000000 );
4167 roundAndPackFloatx80(
4168 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4172 /*----------------------------------------------------------------------------
4173 | Returns the result of subtracting the absolute values of the extended
4174 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
4175 | difference is negated before being returned. `zSign' is ignored if the
4176 | result is a NaN. The subtraction is performed according to the IEC/IEEE
4177 | Standard for Binary Floating-Point Arithmetic.
4178 *----------------------------------------------------------------------------*/
4180 static floatx80
subFloatx80Sigs( floatx80 a
, floatx80 b
, flag zSign STATUS_PARAM
)
4182 int32 aExp
, bExp
, zExp
;
4183 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4187 aSig
= extractFloatx80Frac( a
);
4188 aExp
= extractFloatx80Exp( a
);
4189 bSig
= extractFloatx80Frac( b
);
4190 bExp
= extractFloatx80Exp( b
);
4191 expDiff
= aExp
- bExp
;
4192 if ( 0 < expDiff
) goto aExpBigger
;
4193 if ( expDiff
< 0 ) goto bExpBigger
;
4194 if ( aExp
== 0x7FFF ) {
4195 if ( (uint64_t) ( ( aSig
| bSig
)<<1 ) ) {
4196 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4198 float_raise( float_flag_invalid STATUS_VAR
);
4199 z
.low
= floatx80_default_nan_low
;
4200 z
.high
= floatx80_default_nan_high
;
4208 if ( bSig
< aSig
) goto aBigger
;
4209 if ( aSig
< bSig
) goto bBigger
;
4210 return packFloatx80( STATUS(float_rounding_mode
) == float_round_down
, 0, 0 );
4212 if ( bExp
== 0x7FFF ) {
4213 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4214 return packFloatx80( zSign
^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
4216 if ( aExp
== 0 ) ++expDiff
;
4217 shift128RightJamming( aSig
, 0, - expDiff
, &aSig
, &zSig1
);
4219 sub128( bSig
, 0, aSig
, zSig1
, &zSig0
, &zSig1
);
4222 goto normalizeRoundAndPack
;
4224 if ( aExp
== 0x7FFF ) {
4225 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4228 if ( bExp
== 0 ) --expDiff
;
4229 shift128RightJamming( bSig
, 0, expDiff
, &bSig
, &zSig1
);
4231 sub128( aSig
, 0, bSig
, zSig1
, &zSig0
, &zSig1
);
4233 normalizeRoundAndPack
:
4235 normalizeRoundAndPackFloatx80(
4236 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4240 /*----------------------------------------------------------------------------
4241 | Returns the result of adding the extended double-precision floating-point
4242 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4243 | Standard for Binary Floating-Point Arithmetic.
4244 *----------------------------------------------------------------------------*/
4246 floatx80
floatx80_add( floatx80 a
, floatx80 b STATUS_PARAM
)
4250 aSign
= extractFloatx80Sign( a
);
4251 bSign
= extractFloatx80Sign( b
);
4252 if ( aSign
== bSign
) {
4253 return addFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4256 return subFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4261 /*----------------------------------------------------------------------------
4262 | Returns the result of subtracting the extended double-precision floating-
4263 | point values `a' and `b'. The operation is performed according to the
4264 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4265 *----------------------------------------------------------------------------*/
4267 floatx80
floatx80_sub( floatx80 a
, floatx80 b STATUS_PARAM
)
4271 aSign
= extractFloatx80Sign( a
);
4272 bSign
= extractFloatx80Sign( b
);
4273 if ( aSign
== bSign
) {
4274 return subFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4277 return addFloatx80Sigs( a
, b
, aSign STATUS_VAR
);
4282 /*----------------------------------------------------------------------------
4283 | Returns the result of multiplying the extended double-precision floating-
4284 | point values `a' and `b'. The operation is performed according to the
4285 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4286 *----------------------------------------------------------------------------*/
4288 floatx80
floatx80_mul( floatx80 a
, floatx80 b STATUS_PARAM
)
4290 flag aSign
, bSign
, zSign
;
4291 int32 aExp
, bExp
, zExp
;
4292 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4295 aSig
= extractFloatx80Frac( a
);
4296 aExp
= extractFloatx80Exp( a
);
4297 aSign
= extractFloatx80Sign( a
);
4298 bSig
= extractFloatx80Frac( b
);
4299 bExp
= extractFloatx80Exp( b
);
4300 bSign
= extractFloatx80Sign( b
);
4301 zSign
= aSign
^ bSign
;
4302 if ( aExp
== 0x7FFF ) {
4303 if ( (uint64_t) ( aSig
<<1 )
4304 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
4305 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4307 if ( ( bExp
| bSig
) == 0 ) goto invalid
;
4308 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4310 if ( bExp
== 0x7FFF ) {
4311 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4312 if ( ( aExp
| aSig
) == 0 ) {
4314 float_raise( float_flag_invalid STATUS_VAR
);
4315 z
.low
= floatx80_default_nan_low
;
4316 z
.high
= floatx80_default_nan_high
;
4319 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4322 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4323 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
4326 if ( bSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4327 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4329 zExp
= aExp
+ bExp
- 0x3FFE;
4330 mul64To128( aSig
, bSig
, &zSig0
, &zSig1
);
4331 if ( 0 < (int64_t) zSig0
) {
4332 shortShift128Left( zSig0
, zSig1
, 1, &zSig0
, &zSig1
);
4336 roundAndPackFloatx80(
4337 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4341 /*----------------------------------------------------------------------------
4342 | Returns the result of dividing the extended double-precision floating-point
4343 | value `a' by the corresponding value `b'. The operation is performed
4344 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4345 *----------------------------------------------------------------------------*/
4347 floatx80
floatx80_div( floatx80 a
, floatx80 b STATUS_PARAM
)
4349 flag aSign
, bSign
, zSign
;
4350 int32 aExp
, bExp
, zExp
;
4351 uint64_t aSig
, bSig
, zSig0
, zSig1
;
4352 uint64_t rem0
, rem1
, rem2
, term0
, term1
, term2
;
4355 aSig
= extractFloatx80Frac( a
);
4356 aExp
= extractFloatx80Exp( a
);
4357 aSign
= extractFloatx80Sign( a
);
4358 bSig
= extractFloatx80Frac( b
);
4359 bExp
= extractFloatx80Exp( b
);
4360 bSign
= extractFloatx80Sign( b
);
4361 zSign
= aSign
^ bSign
;
4362 if ( aExp
== 0x7FFF ) {
4363 if ( (uint64_t) ( aSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4364 if ( bExp
== 0x7FFF ) {
4365 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4368 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4370 if ( bExp
== 0x7FFF ) {
4371 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4372 return packFloatx80( zSign
, 0, 0 );
4376 if ( ( aExp
| aSig
) == 0 ) {
4378 float_raise( float_flag_invalid STATUS_VAR
);
4379 z
.low
= floatx80_default_nan_low
;
4380 z
.high
= floatx80_default_nan_high
;
4383 float_raise( float_flag_divbyzero STATUS_VAR
);
4384 return packFloatx80( zSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
4386 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4389 if ( aSig
== 0 ) return packFloatx80( zSign
, 0, 0 );
4390 normalizeFloatx80Subnormal( aSig
, &aExp
, &aSig
);
4392 zExp
= aExp
- bExp
+ 0x3FFE;
4394 if ( bSig
<= aSig
) {
4395 shift128Right( aSig
, 0, 1, &aSig
, &rem1
);
4398 zSig0
= estimateDiv128To64( aSig
, rem1
, bSig
);
4399 mul64To128( bSig
, zSig0
, &term0
, &term1
);
4400 sub128( aSig
, rem1
, term0
, term1
, &rem0
, &rem1
);
4401 while ( (int64_t) rem0
< 0 ) {
4403 add128( rem0
, rem1
, 0, bSig
, &rem0
, &rem1
);
4405 zSig1
= estimateDiv128To64( rem1
, 0, bSig
);
4406 if ( (uint64_t) ( zSig1
<<1 ) <= 8 ) {
4407 mul64To128( bSig
, zSig1
, &term1
, &term2
);
4408 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
4409 while ( (int64_t) rem1
< 0 ) {
4411 add128( rem1
, rem2
, 0, bSig
, &rem1
, &rem2
);
4413 zSig1
|= ( ( rem1
| rem2
) != 0 );
4416 roundAndPackFloatx80(
4417 STATUS(floatx80_rounding_precision
), zSign
, zExp
, zSig0
, zSig1 STATUS_VAR
);
4421 /*----------------------------------------------------------------------------
4422 | Returns the remainder of the extended double-precision floating-point value
4423 | `a' with respect to the corresponding value `b'. The operation is performed
4424 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4425 *----------------------------------------------------------------------------*/
4427 floatx80
floatx80_rem( floatx80 a
, floatx80 b STATUS_PARAM
)
4430 int32 aExp
, bExp
, expDiff
;
4431 uint64_t aSig0
, aSig1
, bSig
;
4432 uint64_t q
, term0
, term1
, alternateASig0
, alternateASig1
;
4435 aSig0
= extractFloatx80Frac( a
);
4436 aExp
= extractFloatx80Exp( a
);
4437 aSign
= extractFloatx80Sign( a
);
4438 bSig
= extractFloatx80Frac( b
);
4439 bExp
= extractFloatx80Exp( b
);
4440 if ( aExp
== 0x7FFF ) {
4441 if ( (uint64_t) ( aSig0
<<1 )
4442 || ( ( bExp
== 0x7FFF ) && (uint64_t) ( bSig
<<1 ) ) ) {
4443 return propagateFloatx80NaN( a
, b STATUS_VAR
);
4447 if ( bExp
== 0x7FFF ) {
4448 if ( (uint64_t) ( bSig
<<1 ) ) return propagateFloatx80NaN( a
, b STATUS_VAR
);
4454 float_raise( float_flag_invalid STATUS_VAR
);
4455 z
.low
= floatx80_default_nan_low
;
4456 z
.high
= floatx80_default_nan_high
;
4459 normalizeFloatx80Subnormal( bSig
, &bExp
, &bSig
);
4462 if ( (uint64_t) ( aSig0
<<1 ) == 0 ) return a
;
4463 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
4465 bSig
|= LIT64( 0x8000000000000000 );
4467 expDiff
= aExp
- bExp
;
4469 if ( expDiff
< 0 ) {
4470 if ( expDiff
< -1 ) return a
;
4471 shift128Right( aSig0
, 0, 1, &aSig0
, &aSig1
);
4474 q
= ( bSig
<= aSig0
);
4475 if ( q
) aSig0
-= bSig
;
4477 while ( 0 < expDiff
) {
4478 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
4479 q
= ( 2 < q
) ? q
- 2 : 0;
4480 mul64To128( bSig
, q
, &term0
, &term1
);
4481 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4482 shortShift128Left( aSig0
, aSig1
, 62, &aSig0
, &aSig1
);
4486 if ( 0 < expDiff
) {
4487 q
= estimateDiv128To64( aSig0
, aSig1
, bSig
);
4488 q
= ( 2 < q
) ? q
- 2 : 0;
4490 mul64To128( bSig
, q
<<( 64 - expDiff
), &term0
, &term1
);
4491 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4492 shortShift128Left( 0, bSig
, 64 - expDiff
, &term0
, &term1
);
4493 while ( le128( term0
, term1
, aSig0
, aSig1
) ) {
4495 sub128( aSig0
, aSig1
, term0
, term1
, &aSig0
, &aSig1
);
4502 sub128( term0
, term1
, aSig0
, aSig1
, &alternateASig0
, &alternateASig1
);
4503 if ( lt128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
4504 || ( eq128( alternateASig0
, alternateASig1
, aSig0
, aSig1
)
4507 aSig0
= alternateASig0
;
4508 aSig1
= alternateASig1
;
4512 normalizeRoundAndPackFloatx80(
4513 80, zSign
, bExp
+ expDiff
, aSig0
, aSig1 STATUS_VAR
);
4517 /*----------------------------------------------------------------------------
4518 | Returns the square root of the extended double-precision floating-point
4519 | value `a'. The operation is performed according to the IEC/IEEE Standard
4520 | for Binary Floating-Point Arithmetic.
4521 *----------------------------------------------------------------------------*/
4523 floatx80
floatx80_sqrt( floatx80 a STATUS_PARAM
)
4527 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, doubleZSig0
;
4528 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
4531 aSig0
= extractFloatx80Frac( a
);
4532 aExp
= extractFloatx80Exp( a
);
4533 aSign
= extractFloatx80Sign( a
);
4534 if ( aExp
== 0x7FFF ) {
4535 if ( (uint64_t) ( aSig0
<<1 ) ) return propagateFloatx80NaN( a
, a STATUS_VAR
);
4536 if ( ! aSign
) return a
;
4540 if ( ( aExp
| aSig0
) == 0 ) return a
;
4542 float_raise( float_flag_invalid STATUS_VAR
);
4543 z
.low
= floatx80_default_nan_low
;
4544 z
.high
= floatx80_default_nan_high
;
4548 if ( aSig0
== 0 ) return packFloatx80( 0, 0, 0 );
4549 normalizeFloatx80Subnormal( aSig0
, &aExp
, &aSig0
);
4551 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFF;
4552 zSig0
= estimateSqrt32( aExp
, aSig0
>>32 );
4553 shift128Right( aSig0
, 0, 2 + ( aExp
& 1 ), &aSig0
, &aSig1
);
4554 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
4555 doubleZSig0
= zSig0
<<1;
4556 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
4557 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
4558 while ( (int64_t) rem0
< 0 ) {
4561 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
4563 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
4564 if ( ( zSig1
& LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
4565 if ( zSig1
== 0 ) zSig1
= 1;
4566 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
4567 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
4568 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
4569 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
4570 while ( (int64_t) rem1
< 0 ) {
4572 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
4574 term2
|= doubleZSig0
;
4575 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
4577 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
4579 shortShift128Left( 0, zSig1
, 1, &zSig0
, &zSig1
);
4580 zSig0
|= doubleZSig0
;
4582 roundAndPackFloatx80(
4583 STATUS(floatx80_rounding_precision
), 0, zExp
, zSig0
, zSig1 STATUS_VAR
);
4587 /*----------------------------------------------------------------------------
4588 | Returns 1 if the extended double-precision floating-point value `a' is equal
4589 | to the corresponding value `b', and 0 otherwise. The invalid exception is
4590 | raised if either operand is a NaN. Otherwise, the comparison is performed
4591 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4592 *----------------------------------------------------------------------------*/
4594 int floatx80_eq( floatx80 a
, floatx80 b STATUS_PARAM
)
4597 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4598 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4599 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4600 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4602 float_raise( float_flag_invalid STATUS_VAR
);
4607 && ( ( a
.high
== b
.high
)
4609 && ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
4614 /*----------------------------------------------------------------------------
4615 | Returns 1 if the extended double-precision floating-point value `a' is
4616 | less than or equal to the corresponding value `b', and 0 otherwise. The
4617 | comparison is performed according to the IEC/IEEE Standard for Binary
4618 | Floating-Point Arithmetic.
4619 *----------------------------------------------------------------------------*/
4621 int floatx80_le( floatx80 a
, floatx80 b STATUS_PARAM
)
4625 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4626 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4627 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4628 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4630 float_raise( float_flag_invalid STATUS_VAR
);
4633 aSign
= extractFloatx80Sign( a
);
4634 bSign
= extractFloatx80Sign( b
);
4635 if ( aSign
!= bSign
) {
4638 || ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
4642 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
4643 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
4647 /*----------------------------------------------------------------------------
4648 | Returns 1 if the extended double-precision floating-point value `a' is
4649 | less than the corresponding value `b', and 0 otherwise. The comparison
4650 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4652 *----------------------------------------------------------------------------*/
4654 int floatx80_lt( floatx80 a
, floatx80 b STATUS_PARAM
)
4658 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4659 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4660 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4661 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4663 float_raise( float_flag_invalid STATUS_VAR
);
4666 aSign
= extractFloatx80Sign( a
);
4667 bSign
= extractFloatx80Sign( b
);
4668 if ( aSign
!= bSign
) {
4671 && ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
4675 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
4676 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
4680 /*----------------------------------------------------------------------------
4681 | Returns 1 if the extended double-precision floating-point values `a' and `b'
4682 | cannot be compared, and 0 otherwise. The comparison is performed according
4683 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4684 *----------------------------------------------------------------------------*/
4685 int floatx80_unordered( floatx80 a
, floatx80 b STATUS_PARAM
)
4687 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4688 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4689 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4690 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4692 float_raise( float_flag_invalid STATUS_VAR
);
4698 /*----------------------------------------------------------------------------
4699 | Returns 1 if the extended double-precision floating-point value `a' is
4700 | equal to the corresponding value `b', and 0 otherwise. The comparison is
4701 | performed according to the IEC/IEEE Standard for Binary Floating-Point
4703 *----------------------------------------------------------------------------*/
4705 int floatx80_eq_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
4708 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4709 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4710 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4711 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4713 if ( floatx80_is_signaling_nan( a
)
4714 || floatx80_is_signaling_nan( b
) ) {
4715 float_raise( float_flag_invalid STATUS_VAR
);
4721 && ( ( a
.high
== b
.high
)
4723 && ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
4728 /*----------------------------------------------------------------------------
4729 | Returns 1 if the extended double-precision floating-point value `a' is less
4730 | than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
4731 | do not cause an exception. Otherwise, the comparison is performed according
4732 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4733 *----------------------------------------------------------------------------*/
4735 int floatx80_le_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
4739 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4740 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4741 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4742 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4744 if ( floatx80_is_signaling_nan( a
)
4745 || floatx80_is_signaling_nan( b
) ) {
4746 float_raise( float_flag_invalid STATUS_VAR
);
4750 aSign
= extractFloatx80Sign( a
);
4751 bSign
= extractFloatx80Sign( b
);
4752 if ( aSign
!= bSign
) {
4755 || ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
4759 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
4760 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
4764 /*----------------------------------------------------------------------------
4765 | Returns 1 if the extended double-precision floating-point value `a' is less
4766 | than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
4767 | an exception. Otherwise, the comparison is performed according to the
4768 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4769 *----------------------------------------------------------------------------*/
4771 int floatx80_lt_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
4775 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4776 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4777 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4778 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4780 if ( floatx80_is_signaling_nan( a
)
4781 || floatx80_is_signaling_nan( b
) ) {
4782 float_raise( float_flag_invalid STATUS_VAR
);
4786 aSign
= extractFloatx80Sign( a
);
4787 bSign
= extractFloatx80Sign( b
);
4788 if ( aSign
!= bSign
) {
4791 && ( ( ( (uint16_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
4795 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
4796 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
4800 /*----------------------------------------------------------------------------
4801 | Returns 1 if the extended double-precision floating-point values `a' and `b'
4802 | cannot be compared, and 0 otherwise. Quiet NaNs do not cause an exception.
4803 | The comparison is performed according to the IEC/IEEE Standard for Binary
4804 | Floating-Point Arithmetic.
4805 *----------------------------------------------------------------------------*/
4806 int floatx80_unordered_quiet( floatx80 a
, floatx80 b STATUS_PARAM
)
4808 if ( ( ( extractFloatx80Exp( a
) == 0x7FFF )
4809 && (uint64_t) ( extractFloatx80Frac( a
)<<1 ) )
4810 || ( ( extractFloatx80Exp( b
) == 0x7FFF )
4811 && (uint64_t) ( extractFloatx80Frac( b
)<<1 ) )
4813 if ( floatx80_is_signaling_nan( a
)
4814 || floatx80_is_signaling_nan( b
) ) {
4815 float_raise( float_flag_invalid STATUS_VAR
);
4826 /*----------------------------------------------------------------------------
4827 | Returns the result of converting the quadruple-precision floating-point
4828 | value `a' to the 32-bit two's complement integer format. The conversion
4829 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4830 | Arithmetic---which means in particular that the conversion is rounded
4831 | according to the current rounding mode. If `a' is a NaN, the largest
4832 | positive integer is returned. Otherwise, if the conversion overflows, the
4833 | largest integer with the same sign as `a' is returned.
4834 *----------------------------------------------------------------------------*/
4836 int32
float128_to_int32( float128 a STATUS_PARAM
)
4839 int32 aExp
, shiftCount
;
4840 uint64_t aSig0
, aSig1
;
4842 aSig1
= extractFloat128Frac1( a
);
4843 aSig0
= extractFloat128Frac0( a
);
4844 aExp
= extractFloat128Exp( a
);
4845 aSign
= extractFloat128Sign( a
);
4846 if ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) aSign
= 0;
4847 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
4848 aSig0
|= ( aSig1
!= 0 );
4849 shiftCount
= 0x4028 - aExp
;
4850 if ( 0 < shiftCount
) shift64RightJamming( aSig0
, shiftCount
, &aSig0
);
4851 return roundAndPackInt32( aSign
, aSig0 STATUS_VAR
);
4855 /*----------------------------------------------------------------------------
4856 | Returns the result of converting the quadruple-precision floating-point
4857 | value `a' to the 32-bit two's complement integer format. The conversion
4858 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4859 | Arithmetic, except that the conversion is always rounded toward zero. If
4860 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
4861 | conversion overflows, the largest integer with the same sign as `a' is
4863 *----------------------------------------------------------------------------*/
4865 int32
float128_to_int32_round_to_zero( float128 a STATUS_PARAM
)
4868 int32 aExp
, shiftCount
;
4869 uint64_t aSig0
, aSig1
, savedASig
;
4872 aSig1
= extractFloat128Frac1( a
);
4873 aSig0
= extractFloat128Frac0( a
);
4874 aExp
= extractFloat128Exp( a
);
4875 aSign
= extractFloat128Sign( a
);
4876 aSig0
|= ( aSig1
!= 0 );
4877 if ( 0x401E < aExp
) {
4878 if ( ( aExp
== 0x7FFF ) && aSig0
) aSign
= 0;
4881 else if ( aExp
< 0x3FFF ) {
4882 if ( aExp
|| aSig0
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4885 aSig0
|= LIT64( 0x0001000000000000 );
4886 shiftCount
= 0x402F - aExp
;
4888 aSig0
>>= shiftCount
;
4890 if ( aSign
) z
= - z
;
4891 if ( ( z
< 0 ) ^ aSign
) {
4893 float_raise( float_flag_invalid STATUS_VAR
);
4894 return aSign
? (int32_t) 0x80000000 : 0x7FFFFFFF;
4896 if ( ( aSig0
<<shiftCount
) != savedASig
) {
4897 STATUS(float_exception_flags
) |= float_flag_inexact
;
4903 /*----------------------------------------------------------------------------
4904 | Returns the result of converting the quadruple-precision floating-point
4905 | value `a' to the 64-bit two's complement integer format. The conversion
4906 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4907 | Arithmetic---which means in particular that the conversion is rounded
4908 | according to the current rounding mode. If `a' is a NaN, the largest
4909 | positive integer is returned. Otherwise, if the conversion overflows, the
4910 | largest integer with the same sign as `a' is returned.
4911 *----------------------------------------------------------------------------*/
4913 int64
float128_to_int64( float128 a STATUS_PARAM
)
4916 int32 aExp
, shiftCount
;
4917 uint64_t aSig0
, aSig1
;
4919 aSig1
= extractFloat128Frac1( a
);
4920 aSig0
= extractFloat128Frac0( a
);
4921 aExp
= extractFloat128Exp( a
);
4922 aSign
= extractFloat128Sign( a
);
4923 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
4924 shiftCount
= 0x402F - aExp
;
4925 if ( shiftCount
<= 0 ) {
4926 if ( 0x403E < aExp
) {
4927 float_raise( float_flag_invalid STATUS_VAR
);
4929 || ( ( aExp
== 0x7FFF )
4930 && ( aSig1
|| ( aSig0
!= LIT64( 0x0001000000000000 ) ) )
4933 return LIT64( 0x7FFFFFFFFFFFFFFF );
4935 return (int64_t) LIT64( 0x8000000000000000 );
4937 shortShift128Left( aSig0
, aSig1
, - shiftCount
, &aSig0
, &aSig1
);
4940 shift64ExtraRightJamming( aSig0
, aSig1
, shiftCount
, &aSig0
, &aSig1
);
4942 return roundAndPackInt64( aSign
, aSig0
, aSig1 STATUS_VAR
);
4946 /*----------------------------------------------------------------------------
4947 | Returns the result of converting the quadruple-precision floating-point
4948 | value `a' to the 64-bit two's complement integer format. The conversion
4949 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4950 | Arithmetic, except that the conversion is always rounded toward zero.
4951 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
4952 | the conversion overflows, the largest integer with the same sign as `a' is
4954 *----------------------------------------------------------------------------*/
4956 int64
float128_to_int64_round_to_zero( float128 a STATUS_PARAM
)
4959 int32 aExp
, shiftCount
;
4960 uint64_t aSig0
, aSig1
;
4963 aSig1
= extractFloat128Frac1( a
);
4964 aSig0
= extractFloat128Frac0( a
);
4965 aExp
= extractFloat128Exp( a
);
4966 aSign
= extractFloat128Sign( a
);
4967 if ( aExp
) aSig0
|= LIT64( 0x0001000000000000 );
4968 shiftCount
= aExp
- 0x402F;
4969 if ( 0 < shiftCount
) {
4970 if ( 0x403E <= aExp
) {
4971 aSig0
&= LIT64( 0x0000FFFFFFFFFFFF );
4972 if ( ( a
.high
== LIT64( 0xC03E000000000000 ) )
4973 && ( aSig1
< LIT64( 0x0002000000000000 ) ) ) {
4974 if ( aSig1
) STATUS(float_exception_flags
) |= float_flag_inexact
;
4977 float_raise( float_flag_invalid STATUS_VAR
);
4978 if ( ! aSign
|| ( ( aExp
== 0x7FFF ) && ( aSig0
| aSig1
) ) ) {
4979 return LIT64( 0x7FFFFFFFFFFFFFFF );
4982 return (int64_t) LIT64( 0x8000000000000000 );
4984 z
= ( aSig0
<<shiftCount
) | ( aSig1
>>( ( - shiftCount
) & 63 ) );
4985 if ( (uint64_t) ( aSig1
<<shiftCount
) ) {
4986 STATUS(float_exception_flags
) |= float_flag_inexact
;
4990 if ( aExp
< 0x3FFF ) {
4991 if ( aExp
| aSig0
| aSig1
) {
4992 STATUS(float_exception_flags
) |= float_flag_inexact
;
4996 z
= aSig0
>>( - shiftCount
);
4998 || ( shiftCount
&& (uint64_t) ( aSig0
<<( shiftCount
& 63 ) ) ) ) {
4999 STATUS(float_exception_flags
) |= float_flag_inexact
;
5002 if ( aSign
) z
= - z
;
5007 /*----------------------------------------------------------------------------
5008 | Returns the result of converting the quadruple-precision floating-point
5009 | value `a' to the single-precision floating-point format. The conversion
5010 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5012 *----------------------------------------------------------------------------*/
5014 float32
float128_to_float32( float128 a STATUS_PARAM
)
5018 uint64_t aSig0
, aSig1
;
5021 aSig1
= extractFloat128Frac1( a
);
5022 aSig0
= extractFloat128Frac0( a
);
5023 aExp
= extractFloat128Exp( a
);
5024 aSign
= extractFloat128Sign( a
);
5025 if ( aExp
== 0x7FFF ) {
5026 if ( aSig0
| aSig1
) {
5027 return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5029 return packFloat32( aSign
, 0xFF, 0 );
5031 aSig0
|= ( aSig1
!= 0 );
5032 shift64RightJamming( aSig0
, 18, &aSig0
);
5034 if ( aExp
|| zSig
) {
5038 return roundAndPackFloat32( aSign
, aExp
, zSig STATUS_VAR
);
5042 /*----------------------------------------------------------------------------
5043 | Returns the result of converting the quadruple-precision floating-point
5044 | value `a' to the double-precision floating-point format. The conversion
5045 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5047 *----------------------------------------------------------------------------*/
5049 float64
float128_to_float64( float128 a STATUS_PARAM
)
5053 uint64_t aSig0
, aSig1
;
5055 aSig1
= extractFloat128Frac1( a
);
5056 aSig0
= extractFloat128Frac0( a
);
5057 aExp
= extractFloat128Exp( a
);
5058 aSign
= extractFloat128Sign( a
);
5059 if ( aExp
== 0x7FFF ) {
5060 if ( aSig0
| aSig1
) {
5061 return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5063 return packFloat64( aSign
, 0x7FF, 0 );
5065 shortShift128Left( aSig0
, aSig1
, 14, &aSig0
, &aSig1
);
5066 aSig0
|= ( aSig1
!= 0 );
5067 if ( aExp
|| aSig0
) {
5068 aSig0
|= LIT64( 0x4000000000000000 );
5071 return roundAndPackFloat64( aSign
, aExp
, aSig0 STATUS_VAR
);
5077 /*----------------------------------------------------------------------------
5078 | Returns the result of converting the quadruple-precision floating-point
5079 | value `a' to the extended double-precision floating-point format. The
5080 | conversion is performed according to the IEC/IEEE Standard for Binary
5081 | Floating-Point Arithmetic.
5082 *----------------------------------------------------------------------------*/
5084 floatx80
float128_to_floatx80( float128 a STATUS_PARAM
)
5088 uint64_t aSig0
, aSig1
;
5090 aSig1
= extractFloat128Frac1( a
);
5091 aSig0
= extractFloat128Frac0( a
);
5092 aExp
= extractFloat128Exp( a
);
5093 aSign
= extractFloat128Sign( a
);
5094 if ( aExp
== 0x7FFF ) {
5095 if ( aSig0
| aSig1
) {
5096 return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR
) STATUS_VAR
);
5098 return packFloatx80( aSign
, 0x7FFF, LIT64( 0x8000000000000000 ) );
5101 if ( ( aSig0
| aSig1
) == 0 ) return packFloatx80( aSign
, 0, 0 );
5102 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5105 aSig0
|= LIT64( 0x0001000000000000 );
5107 shortShift128Left( aSig0
, aSig1
, 15, &aSig0
, &aSig1
);
5108 return roundAndPackFloatx80( 80, aSign
, aExp
, aSig0
, aSig1 STATUS_VAR
);
5114 /*----------------------------------------------------------------------------
5115 | Rounds the quadruple-precision floating-point value `a' to an integer, and
5116 | returns the result as a quadruple-precision floating-point value. The
5117 | operation is performed according to the IEC/IEEE Standard for Binary
5118 | Floating-Point Arithmetic.
5119 *----------------------------------------------------------------------------*/
5121 float128
float128_round_to_int( float128 a STATUS_PARAM
)
5125 uint64_t lastBitMask
, roundBitsMask
;
5129 aExp
= extractFloat128Exp( a
);
5130 if ( 0x402F <= aExp
) {
5131 if ( 0x406F <= aExp
) {
5132 if ( ( aExp
== 0x7FFF )
5133 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) )
5135 return propagateFloat128NaN( a
, a STATUS_VAR
);
5140 lastBitMask
= ( lastBitMask
<<( 0x406E - aExp
) )<<1;
5141 roundBitsMask
= lastBitMask
- 1;
5143 roundingMode
= STATUS(float_rounding_mode
);
5144 if ( roundingMode
== float_round_nearest_even
) {
5145 if ( lastBitMask
) {
5146 add128( z
.high
, z
.low
, 0, lastBitMask
>>1, &z
.high
, &z
.low
);
5147 if ( ( z
.low
& roundBitsMask
) == 0 ) z
.low
&= ~ lastBitMask
;
5150 if ( (int64_t) z
.low
< 0 ) {
5152 if ( (uint64_t) ( z
.low
<<1 ) == 0 ) z
.high
&= ~1;
5156 else if ( roundingMode
!= float_round_to_zero
) {
5157 if ( extractFloat128Sign( z
)
5158 ^ ( roundingMode
== float_round_up
) ) {
5159 add128( z
.high
, z
.low
, 0, roundBitsMask
, &z
.high
, &z
.low
);
5162 z
.low
&= ~ roundBitsMask
;
5165 if ( aExp
< 0x3FFF ) {
5166 if ( ( ( (uint64_t) ( a
.high
<<1 ) ) | a
.low
) == 0 ) return a
;
5167 STATUS(float_exception_flags
) |= float_flag_inexact
;
5168 aSign
= extractFloat128Sign( a
);
5169 switch ( STATUS(float_rounding_mode
) ) {
5170 case float_round_nearest_even
:
5171 if ( ( aExp
== 0x3FFE )
5172 && ( extractFloat128Frac0( a
)
5173 | extractFloat128Frac1( a
) )
5175 return packFloat128( aSign
, 0x3FFF, 0, 0 );
5178 case float_round_down
:
5180 aSign
? packFloat128( 1, 0x3FFF, 0, 0 )
5181 : packFloat128( 0, 0, 0, 0 );
5182 case float_round_up
:
5184 aSign
? packFloat128( 1, 0, 0, 0 )
5185 : packFloat128( 0, 0x3FFF, 0, 0 );
5187 return packFloat128( aSign
, 0, 0, 0 );
5190 lastBitMask
<<= 0x402F - aExp
;
5191 roundBitsMask
= lastBitMask
- 1;
5194 roundingMode
= STATUS(float_rounding_mode
);
5195 if ( roundingMode
== float_round_nearest_even
) {
5196 z
.high
+= lastBitMask
>>1;
5197 if ( ( ( z
.high
& roundBitsMask
) | a
.low
) == 0 ) {
5198 z
.high
&= ~ lastBitMask
;
5201 else if ( roundingMode
!= float_round_to_zero
) {
5202 if ( extractFloat128Sign( z
)
5203 ^ ( roundingMode
== float_round_up
) ) {
5204 z
.high
|= ( a
.low
!= 0 );
5205 z
.high
+= roundBitsMask
;
5208 z
.high
&= ~ roundBitsMask
;
5210 if ( ( z
.low
!= a
.low
) || ( z
.high
!= a
.high
) ) {
5211 STATUS(float_exception_flags
) |= float_flag_inexact
;
5217 /*----------------------------------------------------------------------------
5218 | Returns the result of adding the absolute values of the quadruple-precision
5219 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
5220 | before being returned. `zSign' is ignored if the result is a NaN.
5221 | The addition is performed according to the IEC/IEEE Standard for Binary
5222 | Floating-Point Arithmetic.
5223 *----------------------------------------------------------------------------*/
5225 static float128
addFloat128Sigs( float128 a
, float128 b
, flag zSign STATUS_PARAM
)
5227 int32 aExp
, bExp
, zExp
;
5228 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
;
5231 aSig1
= extractFloat128Frac1( a
);
5232 aSig0
= extractFloat128Frac0( a
);
5233 aExp
= extractFloat128Exp( a
);
5234 bSig1
= extractFloat128Frac1( b
);
5235 bSig0
= extractFloat128Frac0( b
);
5236 bExp
= extractFloat128Exp( b
);
5237 expDiff
= aExp
- bExp
;
5238 if ( 0 < expDiff
) {
5239 if ( aExp
== 0x7FFF ) {
5240 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5247 bSig0
|= LIT64( 0x0001000000000000 );
5249 shift128ExtraRightJamming(
5250 bSig0
, bSig1
, 0, expDiff
, &bSig0
, &bSig1
, &zSig2
);
5253 else if ( expDiff
< 0 ) {
5254 if ( bExp
== 0x7FFF ) {
5255 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5256 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5262 aSig0
|= LIT64( 0x0001000000000000 );
5264 shift128ExtraRightJamming(
5265 aSig0
, aSig1
, 0, - expDiff
, &aSig0
, &aSig1
, &zSig2
);
5269 if ( aExp
== 0x7FFF ) {
5270 if ( aSig0
| aSig1
| bSig0
| bSig1
) {
5271 return propagateFloat128NaN( a
, b STATUS_VAR
);
5275 add128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5277 if ( STATUS(flush_to_zero
) ) return packFloat128( zSign
, 0, 0, 0 );
5278 return packFloat128( zSign
, 0, zSig0
, zSig1
);
5281 zSig0
|= LIT64( 0x0002000000000000 );
5285 aSig0
|= LIT64( 0x0001000000000000 );
5286 add128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5288 if ( zSig0
< LIT64( 0x0002000000000000 ) ) goto roundAndPack
;
5291 shift128ExtraRightJamming(
5292 zSig0
, zSig1
, zSig2
, 1, &zSig0
, &zSig1
, &zSig2
);
5294 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5298 /*----------------------------------------------------------------------------
5299 | Returns the result of subtracting the absolute values of the quadruple-
5300 | precision floating-point values `a' and `b'. If `zSign' is 1, the
5301 | difference is negated before being returned. `zSign' is ignored if the
5302 | result is a NaN. The subtraction is performed according to the IEC/IEEE
5303 | Standard for Binary Floating-Point Arithmetic.
5304 *----------------------------------------------------------------------------*/
5306 static float128
subFloat128Sigs( float128 a
, float128 b
, flag zSign STATUS_PARAM
)
5308 int32 aExp
, bExp
, zExp
;
5309 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
;
5313 aSig1
= extractFloat128Frac1( a
);
5314 aSig0
= extractFloat128Frac0( a
);
5315 aExp
= extractFloat128Exp( a
);
5316 bSig1
= extractFloat128Frac1( b
);
5317 bSig0
= extractFloat128Frac0( b
);
5318 bExp
= extractFloat128Exp( b
);
5319 expDiff
= aExp
- bExp
;
5320 shortShift128Left( aSig0
, aSig1
, 14, &aSig0
, &aSig1
);
5321 shortShift128Left( bSig0
, bSig1
, 14, &bSig0
, &bSig1
);
5322 if ( 0 < expDiff
) goto aExpBigger
;
5323 if ( expDiff
< 0 ) goto bExpBigger
;
5324 if ( aExp
== 0x7FFF ) {
5325 if ( aSig0
| aSig1
| bSig0
| bSig1
) {
5326 return propagateFloat128NaN( a
, b STATUS_VAR
);
5328 float_raise( float_flag_invalid STATUS_VAR
);
5329 z
.low
= float128_default_nan_low
;
5330 z
.high
= float128_default_nan_high
;
5337 if ( bSig0
< aSig0
) goto aBigger
;
5338 if ( aSig0
< bSig0
) goto bBigger
;
5339 if ( bSig1
< aSig1
) goto aBigger
;
5340 if ( aSig1
< bSig1
) goto bBigger
;
5341 return packFloat128( STATUS(float_rounding_mode
) == float_round_down
, 0, 0, 0 );
5343 if ( bExp
== 0x7FFF ) {
5344 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5345 return packFloat128( zSign
^ 1, 0x7FFF, 0, 0 );
5351 aSig0
|= LIT64( 0x4000000000000000 );
5353 shift128RightJamming( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
5354 bSig0
|= LIT64( 0x4000000000000000 );
5356 sub128( bSig0
, bSig1
, aSig0
, aSig1
, &zSig0
, &zSig1
);
5359 goto normalizeRoundAndPack
;
5361 if ( aExp
== 0x7FFF ) {
5362 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5369 bSig0
|= LIT64( 0x4000000000000000 );
5371 shift128RightJamming( bSig0
, bSig1
, expDiff
, &bSig0
, &bSig1
);
5372 aSig0
|= LIT64( 0x4000000000000000 );
5374 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
);
5376 normalizeRoundAndPack
:
5378 return normalizeRoundAndPackFloat128( zSign
, zExp
- 14, zSig0
, zSig1 STATUS_VAR
);
5382 /*----------------------------------------------------------------------------
5383 | Returns the result of adding the quadruple-precision floating-point values
5384 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
5385 | for Binary Floating-Point Arithmetic.
5386 *----------------------------------------------------------------------------*/
5388 float128
float128_add( float128 a
, float128 b STATUS_PARAM
)
5392 aSign
= extractFloat128Sign( a
);
5393 bSign
= extractFloat128Sign( b
);
5394 if ( aSign
== bSign
) {
5395 return addFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5398 return subFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5403 /*----------------------------------------------------------------------------
5404 | Returns the result of subtracting the quadruple-precision floating-point
5405 | values `a' and `b'. The operation is performed according to the IEC/IEEE
5406 | Standard for Binary Floating-Point Arithmetic.
5407 *----------------------------------------------------------------------------*/
5409 float128
float128_sub( float128 a
, float128 b STATUS_PARAM
)
5413 aSign
= extractFloat128Sign( a
);
5414 bSign
= extractFloat128Sign( b
);
5415 if ( aSign
== bSign
) {
5416 return subFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5419 return addFloat128Sigs( a
, b
, aSign STATUS_VAR
);
5424 /*----------------------------------------------------------------------------
5425 | Returns the result of multiplying the quadruple-precision floating-point
5426 | values `a' and `b'. The operation is performed according to the IEC/IEEE
5427 | Standard for Binary Floating-Point Arithmetic.
5428 *----------------------------------------------------------------------------*/
5430 float128
float128_mul( float128 a
, float128 b STATUS_PARAM
)
5432 flag aSign
, bSign
, zSign
;
5433 int32 aExp
, bExp
, zExp
;
5434 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
, zSig3
;
5437 aSig1
= extractFloat128Frac1( a
);
5438 aSig0
= extractFloat128Frac0( a
);
5439 aExp
= extractFloat128Exp( a
);
5440 aSign
= extractFloat128Sign( a
);
5441 bSig1
= extractFloat128Frac1( b
);
5442 bSig0
= extractFloat128Frac0( b
);
5443 bExp
= extractFloat128Exp( b
);
5444 bSign
= extractFloat128Sign( b
);
5445 zSign
= aSign
^ bSign
;
5446 if ( aExp
== 0x7FFF ) {
5447 if ( ( aSig0
| aSig1
)
5448 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
5449 return propagateFloat128NaN( a
, b STATUS_VAR
);
5451 if ( ( bExp
| bSig0
| bSig1
) == 0 ) goto invalid
;
5452 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5454 if ( bExp
== 0x7FFF ) {
5455 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5456 if ( ( aExp
| aSig0
| aSig1
) == 0 ) {
5458 float_raise( float_flag_invalid STATUS_VAR
);
5459 z
.low
= float128_default_nan_low
;
5460 z
.high
= float128_default_nan_high
;
5463 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5466 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5467 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5470 if ( ( bSig0
| bSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5471 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5473 zExp
= aExp
+ bExp
- 0x4000;
5474 aSig0
|= LIT64( 0x0001000000000000 );
5475 shortShift128Left( bSig0
, bSig1
, 16, &bSig0
, &bSig1
);
5476 mul128To256( aSig0
, aSig1
, bSig0
, bSig1
, &zSig0
, &zSig1
, &zSig2
, &zSig3
);
5477 add128( zSig0
, zSig1
, aSig0
, aSig1
, &zSig0
, &zSig1
);
5478 zSig2
|= ( zSig3
!= 0 );
5479 if ( LIT64( 0x0002000000000000 ) <= zSig0
) {
5480 shift128ExtraRightJamming(
5481 zSig0
, zSig1
, zSig2
, 1, &zSig0
, &zSig1
, &zSig2
);
5484 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5488 /*----------------------------------------------------------------------------
5489 | Returns the result of dividing the quadruple-precision floating-point value
5490 | `a' by the corresponding value `b'. The operation is performed according to
5491 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5492 *----------------------------------------------------------------------------*/
5494 float128
float128_div( float128 a
, float128 b STATUS_PARAM
)
5496 flag aSign
, bSign
, zSign
;
5497 int32 aExp
, bExp
, zExp
;
5498 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, zSig0
, zSig1
, zSig2
;
5499 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
5502 aSig1
= extractFloat128Frac1( a
);
5503 aSig0
= extractFloat128Frac0( a
);
5504 aExp
= extractFloat128Exp( a
);
5505 aSign
= extractFloat128Sign( a
);
5506 bSig1
= extractFloat128Frac1( b
);
5507 bSig0
= extractFloat128Frac0( b
);
5508 bExp
= extractFloat128Exp( b
);
5509 bSign
= extractFloat128Sign( b
);
5510 zSign
= aSign
^ bSign
;
5511 if ( aExp
== 0x7FFF ) {
5512 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5513 if ( bExp
== 0x7FFF ) {
5514 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5517 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5519 if ( bExp
== 0x7FFF ) {
5520 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5521 return packFloat128( zSign
, 0, 0, 0 );
5524 if ( ( bSig0
| bSig1
) == 0 ) {
5525 if ( ( aExp
| aSig0
| aSig1
) == 0 ) {
5527 float_raise( float_flag_invalid STATUS_VAR
);
5528 z
.low
= float128_default_nan_low
;
5529 z
.high
= float128_default_nan_high
;
5532 float_raise( float_flag_divbyzero STATUS_VAR
);
5533 return packFloat128( zSign
, 0x7FFF, 0, 0 );
5535 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5538 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( zSign
, 0, 0, 0 );
5539 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5541 zExp
= aExp
- bExp
+ 0x3FFD;
5543 aSig0
| LIT64( 0x0001000000000000 ), aSig1
, 15, &aSig0
, &aSig1
);
5545 bSig0
| LIT64( 0x0001000000000000 ), bSig1
, 15, &bSig0
, &bSig1
);
5546 if ( le128( bSig0
, bSig1
, aSig0
, aSig1
) ) {
5547 shift128Right( aSig0
, aSig1
, 1, &aSig0
, &aSig1
);
5550 zSig0
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5551 mul128By64To192( bSig0
, bSig1
, zSig0
, &term0
, &term1
, &term2
);
5552 sub192( aSig0
, aSig1
, 0, term0
, term1
, term2
, &rem0
, &rem1
, &rem2
);
5553 while ( (int64_t) rem0
< 0 ) {
5555 add192( rem0
, rem1
, rem2
, 0, bSig0
, bSig1
, &rem0
, &rem1
, &rem2
);
5557 zSig1
= estimateDiv128To64( rem1
, rem2
, bSig0
);
5558 if ( ( zSig1
& 0x3FFF ) <= 4 ) {
5559 mul128By64To192( bSig0
, bSig1
, zSig1
, &term1
, &term2
, &term3
);
5560 sub192( rem1
, rem2
, 0, term1
, term2
, term3
, &rem1
, &rem2
, &rem3
);
5561 while ( (int64_t) rem1
< 0 ) {
5563 add192( rem1
, rem2
, rem3
, 0, bSig0
, bSig1
, &rem1
, &rem2
, &rem3
);
5565 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
5567 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 15, &zSig0
, &zSig1
, &zSig2
);
5568 return roundAndPackFloat128( zSign
, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5572 /*----------------------------------------------------------------------------
5573 | Returns the remainder of the quadruple-precision floating-point value `a'
5574 | with respect to the corresponding value `b'. The operation is performed
5575 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5576 *----------------------------------------------------------------------------*/
5578 float128
float128_rem( float128 a
, float128 b STATUS_PARAM
)
5581 int32 aExp
, bExp
, expDiff
;
5582 uint64_t aSig0
, aSig1
, bSig0
, bSig1
, q
, term0
, term1
, term2
;
5583 uint64_t allZero
, alternateASig0
, alternateASig1
, sigMean1
;
5587 aSig1
= extractFloat128Frac1( a
);
5588 aSig0
= extractFloat128Frac0( a
);
5589 aExp
= extractFloat128Exp( a
);
5590 aSign
= extractFloat128Sign( a
);
5591 bSig1
= extractFloat128Frac1( b
);
5592 bSig0
= extractFloat128Frac0( b
);
5593 bExp
= extractFloat128Exp( b
);
5594 if ( aExp
== 0x7FFF ) {
5595 if ( ( aSig0
| aSig1
)
5596 || ( ( bExp
== 0x7FFF ) && ( bSig0
| bSig1
) ) ) {
5597 return propagateFloat128NaN( a
, b STATUS_VAR
);
5601 if ( bExp
== 0x7FFF ) {
5602 if ( bSig0
| bSig1
) return propagateFloat128NaN( a
, b STATUS_VAR
);
5606 if ( ( bSig0
| bSig1
) == 0 ) {
5608 float_raise( float_flag_invalid STATUS_VAR
);
5609 z
.low
= float128_default_nan_low
;
5610 z
.high
= float128_default_nan_high
;
5613 normalizeFloat128Subnormal( bSig0
, bSig1
, &bExp
, &bSig0
, &bSig1
);
5616 if ( ( aSig0
| aSig1
) == 0 ) return a
;
5617 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5619 expDiff
= aExp
- bExp
;
5620 if ( expDiff
< -1 ) return a
;
5622 aSig0
| LIT64( 0x0001000000000000 ),
5624 15 - ( expDiff
< 0 ),
5629 bSig0
| LIT64( 0x0001000000000000 ), bSig1
, 15, &bSig0
, &bSig1
);
5630 q
= le128( bSig0
, bSig1
, aSig0
, aSig1
);
5631 if ( q
) sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
5633 while ( 0 < expDiff
) {
5634 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5635 q
= ( 4 < q
) ? q
- 4 : 0;
5636 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
5637 shortShift192Left( term0
, term1
, term2
, 61, &term1
, &term2
, &allZero
);
5638 shortShift128Left( aSig0
, aSig1
, 61, &aSig0
, &allZero
);
5639 sub128( aSig0
, 0, term1
, term2
, &aSig0
, &aSig1
);
5642 if ( -64 < expDiff
) {
5643 q
= estimateDiv128To64( aSig0
, aSig1
, bSig0
);
5644 q
= ( 4 < q
) ? q
- 4 : 0;
5646 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
5648 if ( expDiff
< 0 ) {
5649 shift128Right( aSig0
, aSig1
, - expDiff
, &aSig0
, &aSig1
);
5652 shortShift128Left( aSig0
, aSig1
, expDiff
, &aSig0
, &aSig1
);
5654 mul128By64To192( bSig0
, bSig1
, q
, &term0
, &term1
, &term2
);
5655 sub128( aSig0
, aSig1
, term1
, term2
, &aSig0
, &aSig1
);
5658 shift128Right( aSig0
, aSig1
, 12, &aSig0
, &aSig1
);
5659 shift128Right( bSig0
, bSig1
, 12, &bSig0
, &bSig1
);
5662 alternateASig0
= aSig0
;
5663 alternateASig1
= aSig1
;
5665 sub128( aSig0
, aSig1
, bSig0
, bSig1
, &aSig0
, &aSig1
);
5666 } while ( 0 <= (int64_t) aSig0
);
5668 aSig0
, aSig1
, alternateASig0
, alternateASig1
, (uint64_t *)&sigMean0
, &sigMean1
);
5669 if ( ( sigMean0
< 0 )
5670 || ( ( ( sigMean0
| sigMean1
) == 0 ) && ( q
& 1 ) ) ) {
5671 aSig0
= alternateASig0
;
5672 aSig1
= alternateASig1
;
5674 zSign
= ( (int64_t) aSig0
< 0 );
5675 if ( zSign
) sub128( 0, 0, aSig0
, aSig1
, &aSig0
, &aSig1
);
5677 normalizeRoundAndPackFloat128( aSign
^ zSign
, bExp
- 4, aSig0
, aSig1 STATUS_VAR
);
5681 /*----------------------------------------------------------------------------
5682 | Returns the square root of the quadruple-precision floating-point value `a'.
5683 | The operation is performed according to the IEC/IEEE Standard for Binary
5684 | Floating-Point Arithmetic.
5685 *----------------------------------------------------------------------------*/
5687 float128
float128_sqrt( float128 a STATUS_PARAM
)
5691 uint64_t aSig0
, aSig1
, zSig0
, zSig1
, zSig2
, doubleZSig0
;
5692 uint64_t rem0
, rem1
, rem2
, rem3
, term0
, term1
, term2
, term3
;
5695 aSig1
= extractFloat128Frac1( a
);
5696 aSig0
= extractFloat128Frac0( a
);
5697 aExp
= extractFloat128Exp( a
);
5698 aSign
= extractFloat128Sign( a
);
5699 if ( aExp
== 0x7FFF ) {
5700 if ( aSig0
| aSig1
) return propagateFloat128NaN( a
, a STATUS_VAR
);
5701 if ( ! aSign
) return a
;
5705 if ( ( aExp
| aSig0
| aSig1
) == 0 ) return a
;
5707 float_raise( float_flag_invalid STATUS_VAR
);
5708 z
.low
= float128_default_nan_low
;
5709 z
.high
= float128_default_nan_high
;
5713 if ( ( aSig0
| aSig1
) == 0 ) return packFloat128( 0, 0, 0, 0 );
5714 normalizeFloat128Subnormal( aSig0
, aSig1
, &aExp
, &aSig0
, &aSig1
);
5716 zExp
= ( ( aExp
- 0x3FFF )>>1 ) + 0x3FFE;
5717 aSig0
|= LIT64( 0x0001000000000000 );
5718 zSig0
= estimateSqrt32( aExp
, aSig0
>>17 );
5719 shortShift128Left( aSig0
, aSig1
, 13 - ( aExp
& 1 ), &aSig0
, &aSig1
);
5720 zSig0
= estimateDiv128To64( aSig0
, aSig1
, zSig0
<<32 ) + ( zSig0
<<30 );
5721 doubleZSig0
= zSig0
<<1;
5722 mul64To128( zSig0
, zSig0
, &term0
, &term1
);
5723 sub128( aSig0
, aSig1
, term0
, term1
, &rem0
, &rem1
);
5724 while ( (int64_t) rem0
< 0 ) {
5727 add128( rem0
, rem1
, zSig0
>>63, doubleZSig0
| 1, &rem0
, &rem1
);
5729 zSig1
= estimateDiv128To64( rem1
, 0, doubleZSig0
);
5730 if ( ( zSig1
& 0x1FFF ) <= 5 ) {
5731 if ( zSig1
== 0 ) zSig1
= 1;
5732 mul64To128( doubleZSig0
, zSig1
, &term1
, &term2
);
5733 sub128( rem1
, 0, term1
, term2
, &rem1
, &rem2
);
5734 mul64To128( zSig1
, zSig1
, &term2
, &term3
);
5735 sub192( rem1
, rem2
, 0, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
5736 while ( (int64_t) rem1
< 0 ) {
5738 shortShift128Left( 0, zSig1
, 1, &term2
, &term3
);
5740 term2
|= doubleZSig0
;
5741 add192( rem1
, rem2
, rem3
, 0, term2
, term3
, &rem1
, &rem2
, &rem3
);
5743 zSig1
|= ( ( rem1
| rem2
| rem3
) != 0 );
5745 shift128ExtraRightJamming( zSig0
, zSig1
, 0, 14, &zSig0
, &zSig1
, &zSig2
);
5746 return roundAndPackFloat128( 0, zExp
, zSig0
, zSig1
, zSig2 STATUS_VAR
);
5750 /*----------------------------------------------------------------------------
5751 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5752 | the corresponding value `b', and 0 otherwise. The invalid exception is
5753 | raised if either operand is a NaN. Otherwise, the comparison is performed
5754 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5755 *----------------------------------------------------------------------------*/
5757 int float128_eq( float128 a
, float128 b STATUS_PARAM
)
5760 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5761 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5762 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5763 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5765 float_raise( float_flag_invalid STATUS_VAR
);
5770 && ( ( a
.high
== b
.high
)
5772 && ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
5777 /*----------------------------------------------------------------------------
5778 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5779 | or equal to the corresponding value `b', and 0 otherwise. The comparison
5780 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5782 *----------------------------------------------------------------------------*/
5784 int float128_le( float128 a
, float128 b STATUS_PARAM
)
5788 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5789 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5790 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5791 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5793 float_raise( float_flag_invalid STATUS_VAR
);
5796 aSign
= extractFloat128Sign( a
);
5797 bSign
= extractFloat128Sign( b
);
5798 if ( aSign
!= bSign
) {
5801 || ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5805 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
5806 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
5810 /*----------------------------------------------------------------------------
5811 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5812 | the corresponding value `b', and 0 otherwise. The comparison is performed
5813 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5814 *----------------------------------------------------------------------------*/
5816 int float128_lt( float128 a
, float128 b STATUS_PARAM
)
5820 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5821 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5822 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5823 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5825 float_raise( float_flag_invalid STATUS_VAR
);
5828 aSign
= extractFloat128Sign( a
);
5829 bSign
= extractFloat128Sign( b
);
5830 if ( aSign
!= bSign
) {
5833 && ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5837 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
5838 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
5842 /*----------------------------------------------------------------------------
5843 | Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
5844 | be compared, and 0 otherwise. The comparison is performed according to the
5845 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5846 *----------------------------------------------------------------------------*/
5848 int float128_unordered( float128 a
, float128 b STATUS_PARAM
)
5850 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5851 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5852 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5853 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5855 float_raise( float_flag_invalid STATUS_VAR
);
5861 /*----------------------------------------------------------------------------
5862 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5863 | the corresponding value `b', and 0 otherwise. The comparison is performed
5864 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5865 *----------------------------------------------------------------------------*/
5867 int float128_eq_quiet( float128 a
, float128 b STATUS_PARAM
)
5870 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5871 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5872 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5873 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5875 if ( float128_is_signaling_nan( a
)
5876 || float128_is_signaling_nan( b
) ) {
5877 float_raise( float_flag_invalid STATUS_VAR
);
5883 && ( ( a
.high
== b
.high
)
5885 && ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) == 0 ) )
5890 /*----------------------------------------------------------------------------
5891 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5892 | or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
5893 | cause an exception. Otherwise, the comparison is performed according to the
5894 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5895 *----------------------------------------------------------------------------*/
5897 int float128_le_quiet( float128 a
, float128 b STATUS_PARAM
)
5901 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5902 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5903 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5904 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5906 if ( float128_is_signaling_nan( a
)
5907 || float128_is_signaling_nan( b
) ) {
5908 float_raise( float_flag_invalid STATUS_VAR
);
5912 aSign
= extractFloat128Sign( a
);
5913 bSign
= extractFloat128Sign( b
);
5914 if ( aSign
!= bSign
) {
5917 || ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5921 aSign
? le128( b
.high
, b
.low
, a
.high
, a
.low
)
5922 : le128( a
.high
, a
.low
, b
.high
, b
.low
);
5926 /*----------------------------------------------------------------------------
5927 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5928 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
5929 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
5930 | Standard for Binary Floating-Point Arithmetic.
5931 *----------------------------------------------------------------------------*/
5933 int float128_lt_quiet( float128 a
, float128 b STATUS_PARAM
)
5937 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5938 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5939 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5940 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5942 if ( float128_is_signaling_nan( a
)
5943 || float128_is_signaling_nan( b
) ) {
5944 float_raise( float_flag_invalid STATUS_VAR
);
5948 aSign
= extractFloat128Sign( a
);
5949 bSign
= extractFloat128Sign( b
);
5950 if ( aSign
!= bSign
) {
5953 && ( ( ( (uint64_t) ( ( a
.high
| b
.high
)<<1 ) ) | a
.low
| b
.low
)
5957 aSign
? lt128( b
.high
, b
.low
, a
.high
, a
.low
)
5958 : lt128( a
.high
, a
.low
, b
.high
, b
.low
);
5962 /*----------------------------------------------------------------------------
5963 | Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
5964 | be compared, and 0 otherwise. Quiet NaNs do not cause an exception. The
5965 | comparison is performed according to the IEC/IEEE Standard for Binary
5966 | Floating-Point Arithmetic.
5967 *----------------------------------------------------------------------------*/
5969 int float128_unordered_quiet( float128 a
, float128 b STATUS_PARAM
)
5971 if ( ( ( extractFloat128Exp( a
) == 0x7FFF )
5972 && ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) )
5973 || ( ( extractFloat128Exp( b
) == 0x7FFF )
5974 && ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )
5976 if ( float128_is_signaling_nan( a
)
5977 || float128_is_signaling_nan( b
) ) {
5978 float_raise( float_flag_invalid STATUS_VAR
);
5987 /* misc functions */
5988 float32
uint32_to_float32( unsigned int a STATUS_PARAM
)
5990 return int64_to_float32(a STATUS_VAR
);
5993 float64
uint32_to_float64( unsigned int a STATUS_PARAM
)
5995 return int64_to_float64(a STATUS_VAR
);
5998 unsigned int float32_to_uint32( float32 a STATUS_PARAM
)
6003 v
= float32_to_int64(a STATUS_VAR
);
6006 float_raise( float_flag_invalid STATUS_VAR
);
6007 } else if (v
> 0xffffffff) {
6009 float_raise( float_flag_invalid STATUS_VAR
);
6016 unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM
)
6021 v
= float32_to_int64_round_to_zero(a STATUS_VAR
);
6024 float_raise( float_flag_invalid STATUS_VAR
);
6025 } else if (v
> 0xffffffff) {
6027 float_raise( float_flag_invalid STATUS_VAR
);
6034 unsigned int float32_to_uint16_round_to_zero( float32 a STATUS_PARAM
)
6039 v
= float32_to_int64_round_to_zero(a STATUS_VAR
);
6042 float_raise( float_flag_invalid STATUS_VAR
);
6043 } else if (v
> 0xffff) {
6045 float_raise( float_flag_invalid STATUS_VAR
);
6052 unsigned int float64_to_uint32( float64 a STATUS_PARAM
)
6057 v
= float64_to_int64(a STATUS_VAR
);
6060 float_raise( float_flag_invalid STATUS_VAR
);
6061 } else if (v
> 0xffffffff) {
6063 float_raise( float_flag_invalid STATUS_VAR
);
6070 unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM
)
6075 v
= float64_to_int64_round_to_zero(a STATUS_VAR
);
6078 float_raise( float_flag_invalid STATUS_VAR
);
6079 } else if (v
> 0xffffffff) {
6081 float_raise( float_flag_invalid STATUS_VAR
);
6088 unsigned int float64_to_uint16_round_to_zero( float64 a STATUS_PARAM
)
6093 v
= float64_to_int64_round_to_zero(a STATUS_VAR
);
6096 float_raise( float_flag_invalid STATUS_VAR
);
6097 } else if (v
> 0xffff) {
6099 float_raise( float_flag_invalid STATUS_VAR
);
6106 /* FIXME: This looks broken. */
6107 uint64_t float64_to_uint64 (float64 a STATUS_PARAM
)
6111 v
= float64_val(int64_to_float64(INT64_MIN STATUS_VAR
));
6112 v
+= float64_val(a
);
6113 v
= float64_to_int64(make_float64(v
) STATUS_VAR
);
6115 return v
- INT64_MIN
;
6118 uint64_t float64_to_uint64_round_to_zero (float64 a STATUS_PARAM
)
6122 v
= float64_val(int64_to_float64(INT64_MIN STATUS_VAR
));
6123 v
+= float64_val(a
);
6124 v
= float64_to_int64_round_to_zero(make_float64(v
) STATUS_VAR
);
6126 return v
- INT64_MIN
;
6129 #define COMPARE(s, nan_exp) \
6130 INLINE int float ## s ## _compare_internal( float ## s a, float ## s b, \
6131 int is_quiet STATUS_PARAM ) \
6133 flag aSign, bSign; \
6134 uint ## s ## _t av, bv; \
6135 a = float ## s ## _squash_input_denormal(a STATUS_VAR); \
6136 b = float ## s ## _squash_input_denormal(b STATUS_VAR); \
6138 if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) && \
6139 extractFloat ## s ## Frac( a ) ) || \
6140 ( ( extractFloat ## s ## Exp( b ) == nan_exp ) && \
6141 extractFloat ## s ## Frac( b ) )) { \
6143 float ## s ## _is_signaling_nan( a ) || \
6144 float ## s ## _is_signaling_nan( b ) ) { \
6145 float_raise( float_flag_invalid STATUS_VAR); \
6147 return float_relation_unordered; \
6149 aSign = extractFloat ## s ## Sign( a ); \
6150 bSign = extractFloat ## s ## Sign( b ); \
6151 av = float ## s ## _val(a); \
6152 bv = float ## s ## _val(b); \
6153 if ( aSign != bSign ) { \
6154 if ( (uint ## s ## _t) ( ( av | bv )<<1 ) == 0 ) { \
6156 return float_relation_equal; \
6158 return 1 - (2 * aSign); \
6162 return float_relation_equal; \
6164 return 1 - 2 * (aSign ^ ( av < bv )); \
6169 int float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM ) \
6171 return float ## s ## _compare_internal(a, b, 0 STATUS_VAR); \
6174 int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
6176 return float ## s ## _compare_internal(a, b, 1 STATUS_VAR); \
6182 INLINE
int float128_compare_internal( float128 a
, float128 b
,
6183 int is_quiet STATUS_PARAM
)
6187 if (( ( extractFloat128Exp( a
) == 0x7fff ) &&
6188 ( extractFloat128Frac0( a
) | extractFloat128Frac1( a
) ) ) ||
6189 ( ( extractFloat128Exp( b
) == 0x7fff ) &&
6190 ( extractFloat128Frac0( b
) | extractFloat128Frac1( b
) ) )) {
6192 float128_is_signaling_nan( a
) ||
6193 float128_is_signaling_nan( b
) ) {
6194 float_raise( float_flag_invalid STATUS_VAR
);
6196 return float_relation_unordered
;
6198 aSign
= extractFloat128Sign( a
);
6199 bSign
= extractFloat128Sign( b
);
6200 if ( aSign
!= bSign
) {
6201 if ( ( ( ( a
.high
| b
.high
)<<1 ) | a
.low
| b
.low
) == 0 ) {
6203 return float_relation_equal
;
6205 return 1 - (2 * aSign
);
6208 if (a
.low
== b
.low
&& a
.high
== b
.high
) {
6209 return float_relation_equal
;
6211 return 1 - 2 * (aSign
^ ( lt128( a
.high
, a
.low
, b
.high
, b
.low
) ));
6216 int float128_compare( float128 a
, float128 b STATUS_PARAM
)
6218 return float128_compare_internal(a
, b
, 0 STATUS_VAR
);
6221 int float128_compare_quiet( float128 a
, float128 b STATUS_PARAM
)
6223 return float128_compare_internal(a
, b
, 1 STATUS_VAR
);
6226 /* min() and max() functions. These can't be implemented as
6227 * 'compare and pick one input' because that would mishandle
6228 * NaNs and +0 vs -0.
6230 #define MINMAX(s, nan_exp) \
6231 INLINE float ## s float ## s ## _minmax(float ## s a, float ## s b, \
6232 int ismin STATUS_PARAM ) \
6234 flag aSign, bSign; \
6235 uint ## s ## _t av, bv; \
6236 a = float ## s ## _squash_input_denormal(a STATUS_VAR); \
6237 b = float ## s ## _squash_input_denormal(b STATUS_VAR); \
6238 if (float ## s ## _is_any_nan(a) || \
6239 float ## s ## _is_any_nan(b)) { \
6240 return propagateFloat ## s ## NaN(a, b STATUS_VAR); \
6242 aSign = extractFloat ## s ## Sign(a); \
6243 bSign = extractFloat ## s ## Sign(b); \
6244 av = float ## s ## _val(a); \
6245 bv = float ## s ## _val(b); \
6246 if (aSign != bSign) { \
6248 return aSign ? a : b; \
6250 return aSign ? b : a; \
6254 return (aSign ^ (av < bv)) ? a : b; \
6256 return (aSign ^ (av < bv)) ? b : a; \
6261 float ## s float ## s ## _min(float ## s a, float ## s b STATUS_PARAM) \
6263 return float ## s ## _minmax(a, b, 1 STATUS_VAR); \
6266 float ## s float ## s ## _max(float ## s a, float ## s b STATUS_PARAM) \
6268 return float ## s ## _minmax(a, b, 0 STATUS_VAR); \
6275 /* Multiply A by 2 raised to the power N. */
6276 float32
float32_scalbn( float32 a
, int n STATUS_PARAM
)
6282 a
= float32_squash_input_denormal(a STATUS_VAR
);
6283 aSig
= extractFloat32Frac( a
);
6284 aExp
= extractFloat32Exp( a
);
6285 aSign
= extractFloat32Sign( a
);
6287 if ( aExp
== 0xFF ) {
6292 else if ( aSig
== 0 )
6297 return normalizeRoundAndPackFloat32( aSign
, aExp
, aSig STATUS_VAR
);
6300 float64
float64_scalbn( float64 a
, int n STATUS_PARAM
)
6306 a
= float64_squash_input_denormal(a STATUS_VAR
);
6307 aSig
= extractFloat64Frac( a
);
6308 aExp
= extractFloat64Exp( a
);
6309 aSign
= extractFloat64Sign( a
);
6311 if ( aExp
== 0x7FF ) {
6315 aSig
|= LIT64( 0x0010000000000000 );
6316 else if ( aSig
== 0 )
6321 return normalizeRoundAndPackFloat64( aSign
, aExp
, aSig STATUS_VAR
);
6325 floatx80
floatx80_scalbn( floatx80 a
, int n STATUS_PARAM
)
6331 aSig
= extractFloatx80Frac( a
);
6332 aExp
= extractFloatx80Exp( a
);
6333 aSign
= extractFloatx80Sign( a
);
6335 if ( aExp
== 0x7FF ) {
6338 if (aExp
== 0 && aSig
== 0)
6342 return normalizeRoundAndPackFloatx80( STATUS(floatx80_rounding_precision
),
6343 aSign
, aExp
, aSig
, 0 STATUS_VAR
);
6348 float128
float128_scalbn( float128 a
, int n STATUS_PARAM
)
6352 uint64_t aSig0
, aSig1
;
6354 aSig1
= extractFloat128Frac1( a
);
6355 aSig0
= extractFloat128Frac0( a
);
6356 aExp
= extractFloat128Exp( a
);
6357 aSign
= extractFloat128Sign( a
);
6358 if ( aExp
== 0x7FFF ) {
6362 aSig0
|= LIT64( 0x0001000000000000 );
6363 else if ( aSig0
== 0 && aSig1
== 0 )
6367 return normalizeRoundAndPackFloat128( aSign
, aExp
, aSig0
, aSig1