Fix most warnings (errors with -Werror) when debugging is enabled
[qemu-kvm/fedora.git] / fpu / softfloat.c
blob4d58744c87b6582363603c0162f19d18c60cd97e
2 /*============================================================================
4 This C source file is part of the SoftFloat IEC/IEEE Floating-point Arithmetic
5 Package, Release 2b.
7 Written by John R. Hauser. This work was made possible in part by the
8 International Computer Science Institute, located at Suite 600, 1947 Center
9 Street, Berkeley, California 94704. Funding was partially provided by the
10 National Science Foundation under grant MIP-9311980. The original version
11 of this code was written as part of a project to build a fixed-point vector
12 processor in collaboration with the University of California at Berkeley,
13 overseen by Profs. Nelson Morgan and John Wawrzynek. More information
14 is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
15 arithmetic/SoftFloat.html'.
17 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
18 been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
19 RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
20 AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
21 COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
22 EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
23 INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
24 OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
26 Derivative works are acceptable, even for commercial purposes, so long as
27 (1) the source code for the derivative work includes prominent notice that
28 the work is derivative, and (2) the source code includes prominent notice with
29 these four paragraphs for those parts of this code that are retained.
31 =============================================================================*/
33 /* FIXME: Flush-To-Zero only effects results. Denormal inputs should also
34 be flushed to zero. */
35 #include "softfloat.h"
37 /*----------------------------------------------------------------------------
38 | Primitive arithmetic functions, including multi-word arithmetic, and
39 | division and square root approximations. (Can be specialized to target if
40 | desired.)
41 *----------------------------------------------------------------------------*/
42 #include "softfloat-macros.h"
44 /*----------------------------------------------------------------------------
45 | Functions and definitions to determine: (1) whether tininess for underflow
46 | is detected before or after rounding by default, (2) what (if anything)
47 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
48 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
49 | are propagated from function inputs to output. These details are target-
50 | specific.
51 *----------------------------------------------------------------------------*/
52 #include "softfloat-specialize.h"
54 void set_float_rounding_mode(int val STATUS_PARAM)
56 STATUS(float_rounding_mode) = val;
59 void set_float_exception_flags(int val STATUS_PARAM)
61 STATUS(float_exception_flags) = val;
64 #ifdef FLOATX80
65 void set_floatx80_rounding_precision(int val STATUS_PARAM)
67 STATUS(floatx80_rounding_precision) = val;
69 #endif
71 /*----------------------------------------------------------------------------
72 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
73 | and 7, and returns the properly rounded 32-bit integer corresponding to the
74 | input. If `zSign' is 1, the input is negated before being converted to an
75 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
76 | is simply rounded to an integer, with the inexact exception raised if the
77 | input cannot be represented exactly as an integer. However, if the fixed-
78 | point input is too large, the invalid exception is raised and the largest
79 | positive or negative integer is returned.
80 *----------------------------------------------------------------------------*/
82 static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
84 int8 roundingMode;
85 flag roundNearestEven;
86 int8 roundIncrement, roundBits;
87 int32 z;
89 roundingMode = STATUS(float_rounding_mode);
90 roundNearestEven = ( roundingMode == float_round_nearest_even );
91 roundIncrement = 0x40;
92 if ( ! roundNearestEven ) {
93 if ( roundingMode == float_round_to_zero ) {
94 roundIncrement = 0;
96 else {
97 roundIncrement = 0x7F;
98 if ( zSign ) {
99 if ( roundingMode == float_round_up ) roundIncrement = 0;
101 else {
102 if ( roundingMode == float_round_down ) roundIncrement = 0;
106 roundBits = absZ & 0x7F;
107 absZ = ( absZ + roundIncrement )>>7;
108 absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
109 z = absZ;
110 if ( zSign ) z = - z;
111 if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
112 float_raise( float_flag_invalid STATUS_VAR);
113 return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
115 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
116 return z;
120 /*----------------------------------------------------------------------------
121 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
122 | `absZ1', with binary point between bits 63 and 64 (between the input words),
123 | and returns the properly rounded 64-bit integer corresponding to the input.
124 | If `zSign' is 1, the input is negated before being converted to an integer.
125 | Ordinarily, the fixed-point input is simply rounded to an integer, with
126 | the inexact exception raised if the input cannot be represented exactly as
127 | an integer. However, if the fixed-point input is too large, the invalid
128 | exception is raised and the largest positive or negative integer is
129 | returned.
130 *----------------------------------------------------------------------------*/
132 static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PARAM)
134 int8 roundingMode;
135 flag roundNearestEven, increment;
136 int64 z;
138 roundingMode = STATUS(float_rounding_mode);
139 roundNearestEven = ( roundingMode == float_round_nearest_even );
140 increment = ( (sbits64) absZ1 < 0 );
141 if ( ! roundNearestEven ) {
142 if ( roundingMode == float_round_to_zero ) {
143 increment = 0;
145 else {
146 if ( zSign ) {
147 increment = ( roundingMode == float_round_down ) && absZ1;
149 else {
150 increment = ( roundingMode == float_round_up ) && absZ1;
154 if ( increment ) {
155 ++absZ0;
156 if ( absZ0 == 0 ) goto overflow;
157 absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
159 z = absZ0;
160 if ( zSign ) z = - z;
161 if ( z && ( ( z < 0 ) ^ zSign ) ) {
162 overflow:
163 float_raise( float_flag_invalid STATUS_VAR);
164 return
165 zSign ? (sbits64) LIT64( 0x8000000000000000 )
166 : LIT64( 0x7FFFFFFFFFFFFFFF );
168 if ( absZ1 ) STATUS(float_exception_flags) |= float_flag_inexact;
169 return z;
173 /*----------------------------------------------------------------------------
174 | Returns the fraction bits of the single-precision floating-point value `a'.
175 *----------------------------------------------------------------------------*/
177 INLINE bits32 extractFloat32Frac( float32 a )
180 return float32_val(a) & 0x007FFFFF;
184 /*----------------------------------------------------------------------------
185 | Returns the exponent bits of the single-precision floating-point value `a'.
186 *----------------------------------------------------------------------------*/
188 INLINE int16 extractFloat32Exp( float32 a )
191 return ( float32_val(a)>>23 ) & 0xFF;
195 /*----------------------------------------------------------------------------
196 | Returns the sign bit of the single-precision floating-point value `a'.
197 *----------------------------------------------------------------------------*/
199 INLINE flag extractFloat32Sign( float32 a )
202 return float32_val(a)>>31;
206 /*----------------------------------------------------------------------------
207 | Normalizes the subnormal single-precision floating-point value represented
208 | by the denormalized significand `aSig'. The normalized exponent and
209 | significand are stored at the locations pointed to by `zExpPtr' and
210 | `zSigPtr', respectively.
211 *----------------------------------------------------------------------------*/
213 static void
214 normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
216 int8 shiftCount;
218 shiftCount = countLeadingZeros32( aSig ) - 8;
219 *zSigPtr = aSig<<shiftCount;
220 *zExpPtr = 1 - shiftCount;
224 /*----------------------------------------------------------------------------
225 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
226 | single-precision floating-point value, returning the result. After being
227 | shifted into the proper positions, the three fields are simply added
228 | together to form the result. This means that any integer portion of `zSig'
229 | will be added into the exponent. Since a properly normalized significand
230 | will have an integer portion equal to 1, the `zExp' input should be 1 less
231 | than the desired result exponent whenever `zSig' is a complete, normalized
232 | significand.
233 *----------------------------------------------------------------------------*/
235 INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
238 return make_float32(
239 ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig);
243 /*----------------------------------------------------------------------------
244 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
245 | and significand `zSig', and returns the proper single-precision floating-
246 | point value corresponding to the abstract input. Ordinarily, the abstract
247 | value is simply rounded and packed into the single-precision format, with
248 | the inexact exception raised if the abstract input cannot be represented
249 | exactly. However, if the abstract value is too large, the overflow and
250 | inexact exceptions are raised and an infinity or maximal finite value is
251 | returned. If the abstract value is too small, the input value is rounded to
252 | a subnormal number, and the underflow and inexact exceptions are raised if
253 | the abstract input cannot be represented exactly as a subnormal single-
254 | precision floating-point number.
255 | The input significand `zSig' has its binary point between bits 30
256 | and 29, which is 7 bits to the left of the usual location. This shifted
257 | significand must be normalized or smaller. If `zSig' is not normalized,
258 | `zExp' must be 0; in that case, the result returned is a subnormal number,
259 | and it must not require rounding. In the usual case that `zSig' is
260 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
261 | The handling of underflow and overflow follows the IEC/IEEE Standard for
262 | Binary Floating-Point Arithmetic.
263 *----------------------------------------------------------------------------*/
265 static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
267 int8 roundingMode;
268 flag roundNearestEven;
269 int8 roundIncrement, roundBits;
270 flag isTiny;
272 roundingMode = STATUS(float_rounding_mode);
273 roundNearestEven = ( roundingMode == float_round_nearest_even );
274 roundIncrement = 0x40;
275 if ( ! roundNearestEven ) {
276 if ( roundingMode == float_round_to_zero ) {
277 roundIncrement = 0;
279 else {
280 roundIncrement = 0x7F;
281 if ( zSign ) {
282 if ( roundingMode == float_round_up ) roundIncrement = 0;
284 else {
285 if ( roundingMode == float_round_down ) roundIncrement = 0;
289 roundBits = zSig & 0x7F;
290 if ( 0xFD <= (bits16) zExp ) {
291 if ( ( 0xFD < zExp )
292 || ( ( zExp == 0xFD )
293 && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
295 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
296 return packFloat32( zSign, 0xFF, - ( roundIncrement == 0 ));
298 if ( zExp < 0 ) {
299 if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
300 isTiny =
301 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
302 || ( zExp < -1 )
303 || ( zSig + roundIncrement < 0x80000000 );
304 shift32RightJamming( zSig, - zExp, &zSig );
305 zExp = 0;
306 roundBits = zSig & 0x7F;
307 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
310 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
311 zSig = ( zSig + roundIncrement )>>7;
312 zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
313 if ( zSig == 0 ) zExp = 0;
314 return packFloat32( zSign, zExp, zSig );
318 /*----------------------------------------------------------------------------
319 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
320 | and significand `zSig', and returns the proper single-precision floating-
321 | point value corresponding to the abstract input. This routine is just like
322 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
323 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
324 | floating-point exponent.
325 *----------------------------------------------------------------------------*/
327 static float32
328 normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
330 int8 shiftCount;
332 shiftCount = countLeadingZeros32( zSig ) - 1;
333 return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
337 /*----------------------------------------------------------------------------
338 | Returns the fraction bits of the double-precision floating-point value `a'.
339 *----------------------------------------------------------------------------*/
341 INLINE bits64 extractFloat64Frac( float64 a )
344 return float64_val(a) & LIT64( 0x000FFFFFFFFFFFFF );
348 /*----------------------------------------------------------------------------
349 | Returns the exponent bits of the double-precision floating-point value `a'.
350 *----------------------------------------------------------------------------*/
352 INLINE int16 extractFloat64Exp( float64 a )
355 return ( float64_val(a)>>52 ) & 0x7FF;
359 /*----------------------------------------------------------------------------
360 | Returns the sign bit of the double-precision floating-point value `a'.
361 *----------------------------------------------------------------------------*/
363 INLINE flag extractFloat64Sign( float64 a )
366 return float64_val(a)>>63;
370 /*----------------------------------------------------------------------------
371 | Normalizes the subnormal double-precision floating-point value represented
372 | by the denormalized significand `aSig'. The normalized exponent and
373 | significand are stored at the locations pointed to by `zExpPtr' and
374 | `zSigPtr', respectively.
375 *----------------------------------------------------------------------------*/
377 static void
378 normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
380 int8 shiftCount;
382 shiftCount = countLeadingZeros64( aSig ) - 11;
383 *zSigPtr = aSig<<shiftCount;
384 *zExpPtr = 1 - shiftCount;
388 /*----------------------------------------------------------------------------
389 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
390 | double-precision floating-point value, returning the result. After being
391 | shifted into the proper positions, the three fields are simply added
392 | together to form the result. This means that any integer portion of `zSig'
393 | will be added into the exponent. Since a properly normalized significand
394 | will have an integer portion equal to 1, the `zExp' input should be 1 less
395 | than the desired result exponent whenever `zSig' is a complete, normalized
396 | significand.
397 *----------------------------------------------------------------------------*/
399 INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
402 return make_float64(
403 ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig);
407 /*----------------------------------------------------------------------------
408 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
409 | and significand `zSig', and returns the proper double-precision floating-
410 | point value corresponding to the abstract input. Ordinarily, the abstract
411 | value is simply rounded and packed into the double-precision format, with
412 | the inexact exception raised if the abstract input cannot be represented
413 | exactly. However, if the abstract value is too large, the overflow and
414 | inexact exceptions are raised and an infinity or maximal finite value is
415 | returned. If the abstract value is too small, the input value is rounded
416 | to a subnormal number, and the underflow and inexact exceptions are raised
417 | if the abstract input cannot be represented exactly as a subnormal double-
418 | precision floating-point number.
419 | The input significand `zSig' has its binary point between bits 62
420 | and 61, which is 10 bits to the left of the usual location. This shifted
421 | significand must be normalized or smaller. If `zSig' is not normalized,
422 | `zExp' must be 0; in that case, the result returned is a subnormal number,
423 | and it must not require rounding. In the usual case that `zSig' is
424 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
425 | The handling of underflow and overflow follows the IEC/IEEE Standard for
426 | Binary Floating-Point Arithmetic.
427 *----------------------------------------------------------------------------*/
429 static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
431 int8 roundingMode;
432 flag roundNearestEven;
433 int16 roundIncrement, roundBits;
434 flag isTiny;
436 roundingMode = STATUS(float_rounding_mode);
437 roundNearestEven = ( roundingMode == float_round_nearest_even );
438 roundIncrement = 0x200;
439 if ( ! roundNearestEven ) {
440 if ( roundingMode == float_round_to_zero ) {
441 roundIncrement = 0;
443 else {
444 roundIncrement = 0x3FF;
445 if ( zSign ) {
446 if ( roundingMode == float_round_up ) roundIncrement = 0;
448 else {
449 if ( roundingMode == float_round_down ) roundIncrement = 0;
453 roundBits = zSig & 0x3FF;
454 if ( 0x7FD <= (bits16) zExp ) {
455 if ( ( 0x7FD < zExp )
456 || ( ( zExp == 0x7FD )
457 && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
459 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
460 return packFloat64( zSign, 0x7FF, - ( roundIncrement == 0 ));
462 if ( zExp < 0 ) {
463 if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
464 isTiny =
465 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
466 || ( zExp < -1 )
467 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
468 shift64RightJamming( zSig, - zExp, &zSig );
469 zExp = 0;
470 roundBits = zSig & 0x3FF;
471 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
474 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
475 zSig = ( zSig + roundIncrement )>>10;
476 zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
477 if ( zSig == 0 ) zExp = 0;
478 return packFloat64( zSign, zExp, zSig );
482 /*----------------------------------------------------------------------------
483 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
484 | and significand `zSig', and returns the proper double-precision floating-
485 | point value corresponding to the abstract input. This routine is just like
486 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
487 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
488 | floating-point exponent.
489 *----------------------------------------------------------------------------*/
491 static float64
492 normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
494 int8 shiftCount;
496 shiftCount = countLeadingZeros64( zSig ) - 1;
497 return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
501 #ifdef FLOATX80
503 /*----------------------------------------------------------------------------
504 | Returns the fraction bits of the extended double-precision floating-point
505 | value `a'.
506 *----------------------------------------------------------------------------*/
508 INLINE bits64 extractFloatx80Frac( floatx80 a )
511 return a.low;
515 /*----------------------------------------------------------------------------
516 | Returns the exponent bits of the extended double-precision floating-point
517 | value `a'.
518 *----------------------------------------------------------------------------*/
520 INLINE int32 extractFloatx80Exp( floatx80 a )
523 return a.high & 0x7FFF;
527 /*----------------------------------------------------------------------------
528 | Returns the sign bit of the extended double-precision floating-point value
529 | `a'.
530 *----------------------------------------------------------------------------*/
532 INLINE flag extractFloatx80Sign( floatx80 a )
535 return a.high>>15;
539 /*----------------------------------------------------------------------------
540 | Normalizes the subnormal extended double-precision floating-point value
541 | represented by the denormalized significand `aSig'. The normalized exponent
542 | and significand are stored at the locations pointed to by `zExpPtr' and
543 | `zSigPtr', respectively.
544 *----------------------------------------------------------------------------*/
546 static void
547 normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
549 int8 shiftCount;
551 shiftCount = countLeadingZeros64( aSig );
552 *zSigPtr = aSig<<shiftCount;
553 *zExpPtr = 1 - shiftCount;
557 /*----------------------------------------------------------------------------
558 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
559 | extended double-precision floating-point value, returning the result.
560 *----------------------------------------------------------------------------*/
562 INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
564 floatx80 z;
566 z.low = zSig;
567 z.high = ( ( (bits16) zSign )<<15 ) + zExp;
568 return z;
572 /*----------------------------------------------------------------------------
573 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
574 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
575 | and returns the proper extended double-precision floating-point value
576 | corresponding to the abstract input. Ordinarily, the abstract value is
577 | rounded and packed into the extended double-precision format, with the
578 | inexact exception raised if the abstract input cannot be represented
579 | exactly. However, if the abstract value is too large, the overflow and
580 | inexact exceptions are raised and an infinity or maximal finite value is
581 | returned. If the abstract value is too small, the input value is rounded to
582 | a subnormal number, and the underflow and inexact exceptions are raised if
583 | the abstract input cannot be represented exactly as a subnormal extended
584 | double-precision floating-point number.
585 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
586 | number of bits as single or double precision, respectively. Otherwise, the
587 | result is rounded to the full precision of the extended double-precision
588 | format.
589 | The input significand must be normalized or smaller. If the input
590 | significand is not normalized, `zExp' must be 0; in that case, the result
591 | returned is a subnormal number, and it must not require rounding. The
592 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
593 | Floating-Point Arithmetic.
594 *----------------------------------------------------------------------------*/
596 static floatx80
597 roundAndPackFloatx80(
598 int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
599 STATUS_PARAM)
601 int8 roundingMode;
602 flag roundNearestEven, increment, isTiny;
603 int64 roundIncrement, roundMask, roundBits;
605 roundingMode = STATUS(float_rounding_mode);
606 roundNearestEven = ( roundingMode == float_round_nearest_even );
607 if ( roundingPrecision == 80 ) goto precision80;
608 if ( roundingPrecision == 64 ) {
609 roundIncrement = LIT64( 0x0000000000000400 );
610 roundMask = LIT64( 0x00000000000007FF );
612 else if ( roundingPrecision == 32 ) {
613 roundIncrement = LIT64( 0x0000008000000000 );
614 roundMask = LIT64( 0x000000FFFFFFFFFF );
616 else {
617 goto precision80;
619 zSig0 |= ( zSig1 != 0 );
620 if ( ! roundNearestEven ) {
621 if ( roundingMode == float_round_to_zero ) {
622 roundIncrement = 0;
624 else {
625 roundIncrement = roundMask;
626 if ( zSign ) {
627 if ( roundingMode == float_round_up ) roundIncrement = 0;
629 else {
630 if ( roundingMode == float_round_down ) roundIncrement = 0;
634 roundBits = zSig0 & roundMask;
635 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
636 if ( ( 0x7FFE < zExp )
637 || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
639 goto overflow;
641 if ( zExp <= 0 ) {
642 if ( STATUS(flush_to_zero) ) return packFloatx80( zSign, 0, 0 );
643 isTiny =
644 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
645 || ( zExp < 0 )
646 || ( zSig0 <= zSig0 + roundIncrement );
647 shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
648 zExp = 0;
649 roundBits = zSig0 & roundMask;
650 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
651 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
652 zSig0 += roundIncrement;
653 if ( (sbits64) zSig0 < 0 ) zExp = 1;
654 roundIncrement = roundMask + 1;
655 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
656 roundMask |= roundIncrement;
658 zSig0 &= ~ roundMask;
659 return packFloatx80( zSign, zExp, zSig0 );
662 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
663 zSig0 += roundIncrement;
664 if ( zSig0 < roundIncrement ) {
665 ++zExp;
666 zSig0 = LIT64( 0x8000000000000000 );
668 roundIncrement = roundMask + 1;
669 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
670 roundMask |= roundIncrement;
672 zSig0 &= ~ roundMask;
673 if ( zSig0 == 0 ) zExp = 0;
674 return packFloatx80( zSign, zExp, zSig0 );
675 precision80:
676 increment = ( (sbits64) zSig1 < 0 );
677 if ( ! roundNearestEven ) {
678 if ( roundingMode == float_round_to_zero ) {
679 increment = 0;
681 else {
682 if ( zSign ) {
683 increment = ( roundingMode == float_round_down ) && zSig1;
685 else {
686 increment = ( roundingMode == float_round_up ) && zSig1;
690 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
691 if ( ( 0x7FFE < zExp )
692 || ( ( zExp == 0x7FFE )
693 && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
694 && increment
697 roundMask = 0;
698 overflow:
699 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
700 if ( ( roundingMode == float_round_to_zero )
701 || ( zSign && ( roundingMode == float_round_up ) )
702 || ( ! zSign && ( roundingMode == float_round_down ) )
704 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
706 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
708 if ( zExp <= 0 ) {
709 isTiny =
710 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
711 || ( zExp < 0 )
712 || ! increment
713 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
714 shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
715 zExp = 0;
716 if ( isTiny && zSig1 ) float_raise( float_flag_underflow STATUS_VAR);
717 if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
718 if ( roundNearestEven ) {
719 increment = ( (sbits64) zSig1 < 0 );
721 else {
722 if ( zSign ) {
723 increment = ( roundingMode == float_round_down ) && zSig1;
725 else {
726 increment = ( roundingMode == float_round_up ) && zSig1;
729 if ( increment ) {
730 ++zSig0;
731 zSig0 &=
732 ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
733 if ( (sbits64) zSig0 < 0 ) zExp = 1;
735 return packFloatx80( zSign, zExp, zSig0 );
738 if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
739 if ( increment ) {
740 ++zSig0;
741 if ( zSig0 == 0 ) {
742 ++zExp;
743 zSig0 = LIT64( 0x8000000000000000 );
745 else {
746 zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
749 else {
750 if ( zSig0 == 0 ) zExp = 0;
752 return packFloatx80( zSign, zExp, zSig0 );
756 /*----------------------------------------------------------------------------
757 | Takes an abstract floating-point value having sign `zSign', exponent
758 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
759 | and returns the proper extended double-precision floating-point value
760 | corresponding to the abstract input. This routine is just like
761 | `roundAndPackFloatx80' except that the input significand does not have to be
762 | normalized.
763 *----------------------------------------------------------------------------*/
765 static floatx80
766 normalizeRoundAndPackFloatx80(
767 int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
768 STATUS_PARAM)
770 int8 shiftCount;
772 if ( zSig0 == 0 ) {
773 zSig0 = zSig1;
774 zSig1 = 0;
775 zExp -= 64;
777 shiftCount = countLeadingZeros64( zSig0 );
778 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
779 zExp -= shiftCount;
780 return
781 roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 STATUS_VAR);
785 #endif
787 #ifdef FLOAT128
789 /*----------------------------------------------------------------------------
790 | Returns the least-significant 64 fraction bits of the quadruple-precision
791 | floating-point value `a'.
792 *----------------------------------------------------------------------------*/
794 INLINE bits64 extractFloat128Frac1( float128 a )
797 return a.low;
801 /*----------------------------------------------------------------------------
802 | Returns the most-significant 48 fraction bits of the quadruple-precision
803 | floating-point value `a'.
804 *----------------------------------------------------------------------------*/
806 INLINE bits64 extractFloat128Frac0( float128 a )
809 return a.high & LIT64( 0x0000FFFFFFFFFFFF );
813 /*----------------------------------------------------------------------------
814 | Returns the exponent bits of the quadruple-precision floating-point value
815 | `a'.
816 *----------------------------------------------------------------------------*/
818 INLINE int32 extractFloat128Exp( float128 a )
821 return ( a.high>>48 ) & 0x7FFF;
825 /*----------------------------------------------------------------------------
826 | Returns the sign bit of the quadruple-precision floating-point value `a'.
827 *----------------------------------------------------------------------------*/
829 INLINE flag extractFloat128Sign( float128 a )
832 return a.high>>63;
836 /*----------------------------------------------------------------------------
837 | Normalizes the subnormal quadruple-precision floating-point value
838 | represented by the denormalized significand formed by the concatenation of
839 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
840 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
841 | significand are stored at the location pointed to by `zSig0Ptr', and the
842 | least significant 64 bits of the normalized significand are stored at the
843 | location pointed to by `zSig1Ptr'.
844 *----------------------------------------------------------------------------*/
846 static void
847 normalizeFloat128Subnormal(
848 bits64 aSig0,
849 bits64 aSig1,
850 int32 *zExpPtr,
851 bits64 *zSig0Ptr,
852 bits64 *zSig1Ptr
855 int8 shiftCount;
857 if ( aSig0 == 0 ) {
858 shiftCount = countLeadingZeros64( aSig1 ) - 15;
859 if ( shiftCount < 0 ) {
860 *zSig0Ptr = aSig1>>( - shiftCount );
861 *zSig1Ptr = aSig1<<( shiftCount & 63 );
863 else {
864 *zSig0Ptr = aSig1<<shiftCount;
865 *zSig1Ptr = 0;
867 *zExpPtr = - shiftCount - 63;
869 else {
870 shiftCount = countLeadingZeros64( aSig0 ) - 15;
871 shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
872 *zExpPtr = 1 - shiftCount;
877 /*----------------------------------------------------------------------------
878 | Packs the sign `zSign', the exponent `zExp', and the significand formed
879 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
880 | floating-point value, returning the result. After being shifted into the
881 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
882 | added together to form the most significant 32 bits of the result. This
883 | means that any integer portion of `zSig0' will be added into the exponent.
884 | Since a properly normalized significand will have an integer portion equal
885 | to 1, the `zExp' input should be 1 less than the desired result exponent
886 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
887 | significand.
888 *----------------------------------------------------------------------------*/
890 INLINE float128
891 packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
893 float128 z;
895 z.low = zSig1;
896 z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
897 return z;
901 /*----------------------------------------------------------------------------
902 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
903 | and extended significand formed by the concatenation of `zSig0', `zSig1',
904 | and `zSig2', and returns the proper quadruple-precision floating-point value
905 | corresponding to the abstract input. Ordinarily, the abstract value is
906 | simply rounded and packed into the quadruple-precision format, with the
907 | inexact exception raised if the abstract input cannot be represented
908 | exactly. However, if the abstract value is too large, the overflow and
909 | inexact exceptions are raised and an infinity or maximal finite value is
910 | returned. If the abstract value is too small, the input value is rounded to
911 | a subnormal number, and the underflow and inexact exceptions are raised if
912 | the abstract input cannot be represented exactly as a subnormal quadruple-
913 | precision floating-point number.
914 | The input significand must be normalized or smaller. If the input
915 | significand is not normalized, `zExp' must be 0; in that case, the result
916 | returned is a subnormal number, and it must not require rounding. In the
917 | usual case that the input significand is normalized, `zExp' must be 1 less
918 | than the ``true'' floating-point exponent. The handling of underflow and
919 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
920 *----------------------------------------------------------------------------*/
922 static float128
923 roundAndPackFloat128(
924 flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 STATUS_PARAM)
926 int8 roundingMode;
927 flag roundNearestEven, increment, isTiny;
929 roundingMode = STATUS(float_rounding_mode);
930 roundNearestEven = ( roundingMode == float_round_nearest_even );
931 increment = ( (sbits64) zSig2 < 0 );
932 if ( ! roundNearestEven ) {
933 if ( roundingMode == float_round_to_zero ) {
934 increment = 0;
936 else {
937 if ( zSign ) {
938 increment = ( roundingMode == float_round_down ) && zSig2;
940 else {
941 increment = ( roundingMode == float_round_up ) && zSig2;
945 if ( 0x7FFD <= (bits32) zExp ) {
946 if ( ( 0x7FFD < zExp )
947 || ( ( zExp == 0x7FFD )
948 && eq128(
949 LIT64( 0x0001FFFFFFFFFFFF ),
950 LIT64( 0xFFFFFFFFFFFFFFFF ),
951 zSig0,
952 zSig1
954 && increment
957 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
958 if ( ( roundingMode == float_round_to_zero )
959 || ( zSign && ( roundingMode == float_round_up ) )
960 || ( ! zSign && ( roundingMode == float_round_down ) )
962 return
963 packFloat128(
964 zSign,
965 0x7FFE,
966 LIT64( 0x0000FFFFFFFFFFFF ),
967 LIT64( 0xFFFFFFFFFFFFFFFF )
970 return packFloat128( zSign, 0x7FFF, 0, 0 );
972 if ( zExp < 0 ) {
973 if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
974 isTiny =
975 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
976 || ( zExp < -1 )
977 || ! increment
978 || lt128(
979 zSig0,
980 zSig1,
981 LIT64( 0x0001FFFFFFFFFFFF ),
982 LIT64( 0xFFFFFFFFFFFFFFFF )
984 shift128ExtraRightJamming(
985 zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
986 zExp = 0;
987 if ( isTiny && zSig2 ) float_raise( float_flag_underflow STATUS_VAR);
988 if ( roundNearestEven ) {
989 increment = ( (sbits64) zSig2 < 0 );
991 else {
992 if ( zSign ) {
993 increment = ( roundingMode == float_round_down ) && zSig2;
995 else {
996 increment = ( roundingMode == float_round_up ) && zSig2;
1001 if ( zSig2 ) STATUS(float_exception_flags) |= float_flag_inexact;
1002 if ( increment ) {
1003 add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
1004 zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
1006 else {
1007 if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
1009 return packFloat128( zSign, zExp, zSig0, zSig1 );
1013 /*----------------------------------------------------------------------------
1014 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1015 | and significand formed by the concatenation of `zSig0' and `zSig1', and
1016 | returns the proper quadruple-precision floating-point value corresponding
1017 | to the abstract input. This routine is just like `roundAndPackFloat128'
1018 | except that the input significand has fewer bits and does not have to be
1019 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
1020 | point exponent.
1021 *----------------------------------------------------------------------------*/
1023 static float128
1024 normalizeRoundAndPackFloat128(
1025 flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 STATUS_PARAM)
1027 int8 shiftCount;
1028 bits64 zSig2;
1030 if ( zSig0 == 0 ) {
1031 zSig0 = zSig1;
1032 zSig1 = 0;
1033 zExp -= 64;
1035 shiftCount = countLeadingZeros64( zSig0 ) - 15;
1036 if ( 0 <= shiftCount ) {
1037 zSig2 = 0;
1038 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1040 else {
1041 shift128ExtraRightJamming(
1042 zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
1044 zExp -= shiftCount;
1045 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR);
1049 #endif
1051 /*----------------------------------------------------------------------------
1052 | Returns the result of converting the 32-bit two's complement integer `a'
1053 | to the single-precision floating-point format. The conversion is performed
1054 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1055 *----------------------------------------------------------------------------*/
1057 float32 int32_to_float32( int32 a STATUS_PARAM )
1059 flag zSign;
1061 if ( a == 0 ) return float32_zero;
1062 if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
1063 zSign = ( a < 0 );
1064 return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a STATUS_VAR );
1068 /*----------------------------------------------------------------------------
1069 | Returns the result of converting the 32-bit two's complement integer `a'
1070 | to the double-precision floating-point format. The conversion is performed
1071 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1072 *----------------------------------------------------------------------------*/
1074 float64 int32_to_float64( int32 a STATUS_PARAM )
1076 flag zSign;
1077 uint32 absA;
1078 int8 shiftCount;
1079 bits64 zSig;
1081 if ( a == 0 ) return float64_zero;
1082 zSign = ( a < 0 );
1083 absA = zSign ? - a : a;
1084 shiftCount = countLeadingZeros32( absA ) + 21;
1085 zSig = absA;
1086 return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );
1090 #ifdef FLOATX80
1092 /*----------------------------------------------------------------------------
1093 | Returns the result of converting the 32-bit two's complement integer `a'
1094 | to the extended double-precision floating-point format. The conversion
1095 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1096 | Arithmetic.
1097 *----------------------------------------------------------------------------*/
1099 floatx80 int32_to_floatx80( int32 a STATUS_PARAM )
1101 flag zSign;
1102 uint32 absA;
1103 int8 shiftCount;
1104 bits64 zSig;
1106 if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1107 zSign = ( a < 0 );
1108 absA = zSign ? - a : a;
1109 shiftCount = countLeadingZeros32( absA ) + 32;
1110 zSig = absA;
1111 return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
1115 #endif
1117 #ifdef FLOAT128
1119 /*----------------------------------------------------------------------------
1120 | Returns the result of converting the 32-bit two's complement integer `a' to
1121 | the quadruple-precision floating-point format. The conversion is performed
1122 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1123 *----------------------------------------------------------------------------*/
1125 float128 int32_to_float128( int32 a STATUS_PARAM )
1127 flag zSign;
1128 uint32 absA;
1129 int8 shiftCount;
1130 bits64 zSig0;
1132 if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1133 zSign = ( a < 0 );
1134 absA = zSign ? - a : a;
1135 shiftCount = countLeadingZeros32( absA ) + 17;
1136 zSig0 = absA;
1137 return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );
1141 #endif
1143 /*----------------------------------------------------------------------------
1144 | Returns the result of converting the 64-bit two's complement integer `a'
1145 | to the single-precision floating-point format. The conversion is performed
1146 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1147 *----------------------------------------------------------------------------*/
1149 float32 int64_to_float32( int64 a STATUS_PARAM )
1151 flag zSign;
1152 uint64 absA;
1153 int8 shiftCount;
1155 if ( a == 0 ) return float32_zero;
1156 zSign = ( a < 0 );
1157 absA = zSign ? - a : a;
1158 shiftCount = countLeadingZeros64( absA ) - 40;
1159 if ( 0 <= shiftCount ) {
1160 return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );
1162 else {
1163 shiftCount += 7;
1164 if ( shiftCount < 0 ) {
1165 shift64RightJamming( absA, - shiftCount, &absA );
1167 else {
1168 absA <<= shiftCount;
1170 return roundAndPackFloat32( zSign, 0x9C - shiftCount, absA STATUS_VAR );
1175 float32 uint64_to_float32( uint64 a STATUS_PARAM )
1177 int8 shiftCount;
1179 if ( a == 0 ) return float32_zero;
1180 shiftCount = countLeadingZeros64( a ) - 40;
1181 if ( 0 <= shiftCount ) {
1182 return packFloat32( 1 > 0, 0x95 - shiftCount, a<<shiftCount );
1184 else {
1185 shiftCount += 7;
1186 if ( shiftCount < 0 ) {
1187 shift64RightJamming( a, - shiftCount, &a );
1189 else {
1190 a <<= shiftCount;
1192 return roundAndPackFloat32( 1 > 0, 0x9C - shiftCount, a STATUS_VAR );
1196 /*----------------------------------------------------------------------------
1197 | Returns the result of converting the 64-bit two's complement integer `a'
1198 | to the double-precision floating-point format. The conversion is performed
1199 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1200 *----------------------------------------------------------------------------*/
1202 float64 int64_to_float64( int64 a STATUS_PARAM )
1204 flag zSign;
1206 if ( a == 0 ) return float64_zero;
1207 if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
1208 return packFloat64( 1, 0x43E, 0 );
1210 zSign = ( a < 0 );
1211 return normalizeRoundAndPackFloat64( zSign, 0x43C, zSign ? - a : a STATUS_VAR );
1215 float64 uint64_to_float64( uint64 a STATUS_PARAM )
1217 if ( a == 0 ) return float64_zero;
1218 return normalizeRoundAndPackFloat64( 0, 0x43C, a STATUS_VAR );
1222 #ifdef FLOATX80
1224 /*----------------------------------------------------------------------------
1225 | Returns the result of converting the 64-bit two's complement integer `a'
1226 | to the extended double-precision floating-point format. The conversion
1227 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1228 | Arithmetic.
1229 *----------------------------------------------------------------------------*/
1231 floatx80 int64_to_floatx80( int64 a STATUS_PARAM )
1233 flag zSign;
1234 uint64 absA;
1235 int8 shiftCount;
1237 if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1238 zSign = ( a < 0 );
1239 absA = zSign ? - a : a;
1240 shiftCount = countLeadingZeros64( absA );
1241 return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
1245 #endif
1247 #ifdef FLOAT128
1249 /*----------------------------------------------------------------------------
1250 | Returns the result of converting the 64-bit two's complement integer `a' to
1251 | the quadruple-precision floating-point format. The conversion is performed
1252 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1253 *----------------------------------------------------------------------------*/
1255 float128 int64_to_float128( int64 a STATUS_PARAM )
1257 flag zSign;
1258 uint64 absA;
1259 int8 shiftCount;
1260 int32 zExp;
1261 bits64 zSig0, zSig1;
1263 if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1264 zSign = ( a < 0 );
1265 absA = zSign ? - a : a;
1266 shiftCount = countLeadingZeros64( absA ) + 49;
1267 zExp = 0x406E - shiftCount;
1268 if ( 64 <= shiftCount ) {
1269 zSig1 = 0;
1270 zSig0 = absA;
1271 shiftCount -= 64;
1273 else {
1274 zSig1 = absA;
1275 zSig0 = 0;
1277 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1278 return packFloat128( zSign, zExp, zSig0, zSig1 );
1282 #endif
1284 /*----------------------------------------------------------------------------
1285 | Returns the result of converting the single-precision floating-point value
1286 | `a' to the 32-bit two's complement integer format. The conversion is
1287 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1288 | Arithmetic---which means in particular that the conversion is rounded
1289 | according to the current rounding mode. If `a' is a NaN, the largest
1290 | positive integer is returned. Otherwise, if the conversion overflows, the
1291 | largest integer with the same sign as `a' is returned.
1292 *----------------------------------------------------------------------------*/
1294 int32 float32_to_int32( float32 a STATUS_PARAM )
1296 flag aSign;
1297 int16 aExp, shiftCount;
1298 bits32 aSig;
1299 bits64 aSig64;
1301 aSig = extractFloat32Frac( a );
1302 aExp = extractFloat32Exp( a );
1303 aSign = extractFloat32Sign( a );
1304 if ( ( aExp == 0xFF ) && aSig ) aSign = 0;
1305 if ( aExp ) aSig |= 0x00800000;
1306 shiftCount = 0xAF - aExp;
1307 aSig64 = aSig;
1308 aSig64 <<= 32;
1309 if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );
1310 return roundAndPackInt32( aSign, aSig64 STATUS_VAR );
1314 /*----------------------------------------------------------------------------
1315 | Returns the result of converting the single-precision floating-point value
1316 | `a' to the 32-bit two's complement integer format. The conversion is
1317 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1318 | Arithmetic, except that the conversion is always rounded toward zero.
1319 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1320 | the conversion overflows, the largest integer with the same sign as `a' is
1321 | returned.
1322 *----------------------------------------------------------------------------*/
1324 int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
1326 flag aSign;
1327 int16 aExp, shiftCount;
1328 bits32 aSig;
1329 int32 z;
1331 aSig = extractFloat32Frac( a );
1332 aExp = extractFloat32Exp( a );
1333 aSign = extractFloat32Sign( a );
1334 shiftCount = aExp - 0x9E;
1335 if ( 0 <= shiftCount ) {
1336 if ( float32_val(a) != 0xCF000000 ) {
1337 float_raise( float_flag_invalid STATUS_VAR);
1338 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
1340 return (sbits32) 0x80000000;
1342 else if ( aExp <= 0x7E ) {
1343 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1344 return 0;
1346 aSig = ( aSig | 0x00800000 )<<8;
1347 z = aSig>>( - shiftCount );
1348 if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
1349 STATUS(float_exception_flags) |= float_flag_inexact;
1351 if ( aSign ) z = - z;
1352 return z;
1356 /*----------------------------------------------------------------------------
1357 | Returns the result of converting the single-precision floating-point value
1358 | `a' to the 64-bit two's complement integer format. The conversion is
1359 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1360 | Arithmetic---which means in particular that the conversion is rounded
1361 | according to the current rounding mode. If `a' is a NaN, the largest
1362 | positive integer is returned. Otherwise, if the conversion overflows, the
1363 | largest integer with the same sign as `a' is returned.
1364 *----------------------------------------------------------------------------*/
1366 int64 float32_to_int64( float32 a STATUS_PARAM )
1368 flag aSign;
1369 int16 aExp, shiftCount;
1370 bits32 aSig;
1371 bits64 aSig64, aSigExtra;
1373 aSig = extractFloat32Frac( a );
1374 aExp = extractFloat32Exp( a );
1375 aSign = extractFloat32Sign( a );
1376 shiftCount = 0xBE - aExp;
1377 if ( shiftCount < 0 ) {
1378 float_raise( float_flag_invalid STATUS_VAR);
1379 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1380 return LIT64( 0x7FFFFFFFFFFFFFFF );
1382 return (sbits64) LIT64( 0x8000000000000000 );
1384 if ( aExp ) aSig |= 0x00800000;
1385 aSig64 = aSig;
1386 aSig64 <<= 40;
1387 shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );
1388 return roundAndPackInt64( aSign, aSig64, aSigExtra STATUS_VAR );
1392 /*----------------------------------------------------------------------------
1393 | Returns the result of converting the single-precision floating-point value
1394 | `a' to the 64-bit two's complement integer format. The conversion is
1395 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1396 | Arithmetic, except that the conversion is always rounded toward zero. If
1397 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
1398 | conversion overflows, the largest integer with the same sign as `a' is
1399 | returned.
1400 *----------------------------------------------------------------------------*/
1402 int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
1404 flag aSign;
1405 int16 aExp, shiftCount;
1406 bits32 aSig;
1407 bits64 aSig64;
1408 int64 z;
1410 aSig = extractFloat32Frac( a );
1411 aExp = extractFloat32Exp( a );
1412 aSign = extractFloat32Sign( a );
1413 shiftCount = aExp - 0xBE;
1414 if ( 0 <= shiftCount ) {
1415 if ( float32_val(a) != 0xDF000000 ) {
1416 float_raise( float_flag_invalid STATUS_VAR);
1417 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1418 return LIT64( 0x7FFFFFFFFFFFFFFF );
1421 return (sbits64) LIT64( 0x8000000000000000 );
1423 else if ( aExp <= 0x7E ) {
1424 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1425 return 0;
1427 aSig64 = aSig | 0x00800000;
1428 aSig64 <<= 40;
1429 z = aSig64>>( - shiftCount );
1430 if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
1431 STATUS(float_exception_flags) |= float_flag_inexact;
1433 if ( aSign ) z = - z;
1434 return z;
1438 /*----------------------------------------------------------------------------
1439 | Returns the result of converting the single-precision floating-point value
1440 | `a' to the double-precision floating-point format. The conversion is
1441 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1442 | Arithmetic.
1443 *----------------------------------------------------------------------------*/
1445 float64 float32_to_float64( float32 a STATUS_PARAM )
1447 flag aSign;
1448 int16 aExp;
1449 bits32 aSig;
1451 aSig = extractFloat32Frac( a );
1452 aExp = extractFloat32Exp( a );
1453 aSign = extractFloat32Sign( a );
1454 if ( aExp == 0xFF ) {
1455 if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ));
1456 return packFloat64( aSign, 0x7FF, 0 );
1458 if ( aExp == 0 ) {
1459 if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
1460 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1461 --aExp;
1463 return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
1467 #ifdef FLOATX80
1469 /*----------------------------------------------------------------------------
1470 | Returns the result of converting the single-precision floating-point value
1471 | `a' to the extended double-precision floating-point format. The conversion
1472 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1473 | Arithmetic.
1474 *----------------------------------------------------------------------------*/
1476 floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
1478 flag aSign;
1479 int16 aExp;
1480 bits32 aSig;
1482 aSig = extractFloat32Frac( a );
1483 aExp = extractFloat32Exp( a );
1484 aSign = extractFloat32Sign( a );
1485 if ( aExp == 0xFF ) {
1486 if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) );
1487 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
1489 if ( aExp == 0 ) {
1490 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1491 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1493 aSig |= 0x00800000;
1494 return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
1498 #endif
1500 #ifdef FLOAT128
1502 /*----------------------------------------------------------------------------
1503 | Returns the result of converting the single-precision floating-point value
1504 | `a' to the double-precision floating-point format. The conversion is
1505 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1506 | Arithmetic.
1507 *----------------------------------------------------------------------------*/
1509 float128 float32_to_float128( float32 a STATUS_PARAM )
1511 flag aSign;
1512 int16 aExp;
1513 bits32 aSig;
1515 aSig = extractFloat32Frac( a );
1516 aExp = extractFloat32Exp( a );
1517 aSign = extractFloat32Sign( a );
1518 if ( aExp == 0xFF ) {
1519 if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) );
1520 return packFloat128( aSign, 0x7FFF, 0, 0 );
1522 if ( aExp == 0 ) {
1523 if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
1524 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1525 --aExp;
1527 return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
1531 #endif
1533 /*----------------------------------------------------------------------------
1534 | Rounds the single-precision floating-point value `a' to an integer, and
1535 | returns the result as a single-precision floating-point value. The
1536 | operation is performed according to the IEC/IEEE Standard for Binary
1537 | Floating-Point Arithmetic.
1538 *----------------------------------------------------------------------------*/
1540 float32 float32_round_to_int( float32 a STATUS_PARAM)
1542 flag aSign;
1543 int16 aExp;
1544 bits32 lastBitMask, roundBitsMask;
1545 int8 roundingMode;
1546 bits32 z;
1548 aExp = extractFloat32Exp( a );
1549 if ( 0x96 <= aExp ) {
1550 if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
1551 return propagateFloat32NaN( a, a STATUS_VAR );
1553 return a;
1555 if ( aExp <= 0x7E ) {
1556 if ( (bits32) ( float32_val(a)<<1 ) == 0 ) return a;
1557 STATUS(float_exception_flags) |= float_flag_inexact;
1558 aSign = extractFloat32Sign( a );
1559 switch ( STATUS(float_rounding_mode) ) {
1560 case float_round_nearest_even:
1561 if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
1562 return packFloat32( aSign, 0x7F, 0 );
1564 break;
1565 case float_round_down:
1566 return make_float32(aSign ? 0xBF800000 : 0);
1567 case float_round_up:
1568 return make_float32(aSign ? 0x80000000 : 0x3F800000);
1570 return packFloat32( aSign, 0, 0 );
1572 lastBitMask = 1;
1573 lastBitMask <<= 0x96 - aExp;
1574 roundBitsMask = lastBitMask - 1;
1575 z = float32_val(a);
1576 roundingMode = STATUS(float_rounding_mode);
1577 if ( roundingMode == float_round_nearest_even ) {
1578 z += lastBitMask>>1;
1579 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
1581 else if ( roundingMode != float_round_to_zero ) {
1582 if ( extractFloat32Sign( make_float32(z) ) ^ ( roundingMode == float_round_up ) ) {
1583 z += roundBitsMask;
1586 z &= ~ roundBitsMask;
1587 if ( z != float32_val(a) ) STATUS(float_exception_flags) |= float_flag_inexact;
1588 return make_float32(z);
1592 /*----------------------------------------------------------------------------
1593 | Returns the result of adding the absolute values of the single-precision
1594 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
1595 | before being returned. `zSign' is ignored if the result is a NaN.
1596 | The addition is performed according to the IEC/IEEE Standard for Binary
1597 | Floating-Point Arithmetic.
1598 *----------------------------------------------------------------------------*/
1600 static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1602 int16 aExp, bExp, zExp;
1603 bits32 aSig, bSig, zSig;
1604 int16 expDiff;
1606 aSig = extractFloat32Frac( a );
1607 aExp = extractFloat32Exp( a );
1608 bSig = extractFloat32Frac( b );
1609 bExp = extractFloat32Exp( b );
1610 expDiff = aExp - bExp;
1611 aSig <<= 6;
1612 bSig <<= 6;
1613 if ( 0 < expDiff ) {
1614 if ( aExp == 0xFF ) {
1615 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1616 return a;
1618 if ( bExp == 0 ) {
1619 --expDiff;
1621 else {
1622 bSig |= 0x20000000;
1624 shift32RightJamming( bSig, expDiff, &bSig );
1625 zExp = aExp;
1627 else if ( expDiff < 0 ) {
1628 if ( bExp == 0xFF ) {
1629 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1630 return packFloat32( zSign, 0xFF, 0 );
1632 if ( aExp == 0 ) {
1633 ++expDiff;
1635 else {
1636 aSig |= 0x20000000;
1638 shift32RightJamming( aSig, - expDiff, &aSig );
1639 zExp = bExp;
1641 else {
1642 if ( aExp == 0xFF ) {
1643 if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1644 return a;
1646 if ( aExp == 0 ) {
1647 if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
1648 return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
1650 zSig = 0x40000000 + aSig + bSig;
1651 zExp = aExp;
1652 goto roundAndPack;
1654 aSig |= 0x20000000;
1655 zSig = ( aSig + bSig )<<1;
1656 --zExp;
1657 if ( (sbits32) zSig < 0 ) {
1658 zSig = aSig + bSig;
1659 ++zExp;
1661 roundAndPack:
1662 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1666 /*----------------------------------------------------------------------------
1667 | Returns the result of subtracting the absolute values of the single-
1668 | precision floating-point values `a' and `b'. If `zSign' is 1, the
1669 | difference is negated before being returned. `zSign' is ignored if the
1670 | result is a NaN. The subtraction is performed according to the IEC/IEEE
1671 | Standard for Binary Floating-Point Arithmetic.
1672 *----------------------------------------------------------------------------*/
1674 static float32 subFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1676 int16 aExp, bExp, zExp;
1677 bits32 aSig, bSig, zSig;
1678 int16 expDiff;
1680 aSig = extractFloat32Frac( a );
1681 aExp = extractFloat32Exp( a );
1682 bSig = extractFloat32Frac( b );
1683 bExp = extractFloat32Exp( b );
1684 expDiff = aExp - bExp;
1685 aSig <<= 7;
1686 bSig <<= 7;
1687 if ( 0 < expDiff ) goto aExpBigger;
1688 if ( expDiff < 0 ) goto bExpBigger;
1689 if ( aExp == 0xFF ) {
1690 if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1691 float_raise( float_flag_invalid STATUS_VAR);
1692 return float32_default_nan;
1694 if ( aExp == 0 ) {
1695 aExp = 1;
1696 bExp = 1;
1698 if ( bSig < aSig ) goto aBigger;
1699 if ( aSig < bSig ) goto bBigger;
1700 return packFloat32( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
1701 bExpBigger:
1702 if ( bExp == 0xFF ) {
1703 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1704 return packFloat32( zSign ^ 1, 0xFF, 0 );
1706 if ( aExp == 0 ) {
1707 ++expDiff;
1709 else {
1710 aSig |= 0x40000000;
1712 shift32RightJamming( aSig, - expDiff, &aSig );
1713 bSig |= 0x40000000;
1714 bBigger:
1715 zSig = bSig - aSig;
1716 zExp = bExp;
1717 zSign ^= 1;
1718 goto normalizeRoundAndPack;
1719 aExpBigger:
1720 if ( aExp == 0xFF ) {
1721 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1722 return a;
1724 if ( bExp == 0 ) {
1725 --expDiff;
1727 else {
1728 bSig |= 0x40000000;
1730 shift32RightJamming( bSig, expDiff, &bSig );
1731 aSig |= 0x40000000;
1732 aBigger:
1733 zSig = aSig - bSig;
1734 zExp = aExp;
1735 normalizeRoundAndPack:
1736 --zExp;
1737 return normalizeRoundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1741 /*----------------------------------------------------------------------------
1742 | Returns the result of adding the single-precision floating-point values `a'
1743 | and `b'. The operation is performed according to the IEC/IEEE Standard for
1744 | Binary Floating-Point Arithmetic.
1745 *----------------------------------------------------------------------------*/
1747 float32 float32_add( float32 a, float32 b STATUS_PARAM )
1749 flag aSign, bSign;
1751 aSign = extractFloat32Sign( a );
1752 bSign = extractFloat32Sign( b );
1753 if ( aSign == bSign ) {
1754 return addFloat32Sigs( a, b, aSign STATUS_VAR);
1756 else {
1757 return subFloat32Sigs( a, b, aSign STATUS_VAR );
1762 /*----------------------------------------------------------------------------
1763 | Returns the result of subtracting the single-precision floating-point values
1764 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1765 | for Binary Floating-Point Arithmetic.
1766 *----------------------------------------------------------------------------*/
1768 float32 float32_sub( float32 a, float32 b STATUS_PARAM )
1770 flag aSign, bSign;
1772 aSign = extractFloat32Sign( a );
1773 bSign = extractFloat32Sign( b );
1774 if ( aSign == bSign ) {
1775 return subFloat32Sigs( a, b, aSign STATUS_VAR );
1777 else {
1778 return addFloat32Sigs( a, b, aSign STATUS_VAR );
1783 /*----------------------------------------------------------------------------
1784 | Returns the result of multiplying the single-precision floating-point values
1785 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1786 | for Binary Floating-Point Arithmetic.
1787 *----------------------------------------------------------------------------*/
1789 float32 float32_mul( float32 a, float32 b STATUS_PARAM )
1791 flag aSign, bSign, zSign;
1792 int16 aExp, bExp, zExp;
1793 bits32 aSig, bSig;
1794 bits64 zSig64;
1795 bits32 zSig;
1797 aSig = extractFloat32Frac( a );
1798 aExp = extractFloat32Exp( a );
1799 aSign = extractFloat32Sign( a );
1800 bSig = extractFloat32Frac( b );
1801 bExp = extractFloat32Exp( b );
1802 bSign = extractFloat32Sign( b );
1803 zSign = aSign ^ bSign;
1804 if ( aExp == 0xFF ) {
1805 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1806 return propagateFloat32NaN( a, b STATUS_VAR );
1808 if ( ( bExp | bSig ) == 0 ) {
1809 float_raise( float_flag_invalid STATUS_VAR);
1810 return float32_default_nan;
1812 return packFloat32( zSign, 0xFF, 0 );
1814 if ( bExp == 0xFF ) {
1815 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1816 if ( ( aExp | aSig ) == 0 ) {
1817 float_raise( float_flag_invalid STATUS_VAR);
1818 return float32_default_nan;
1820 return packFloat32( zSign, 0xFF, 0 );
1822 if ( aExp == 0 ) {
1823 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1824 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1826 if ( bExp == 0 ) {
1827 if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
1828 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1830 zExp = aExp + bExp - 0x7F;
1831 aSig = ( aSig | 0x00800000 )<<7;
1832 bSig = ( bSig | 0x00800000 )<<8;
1833 shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
1834 zSig = zSig64;
1835 if ( 0 <= (sbits32) ( zSig<<1 ) ) {
1836 zSig <<= 1;
1837 --zExp;
1839 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1843 /*----------------------------------------------------------------------------
1844 | Returns the result of dividing the single-precision floating-point value `a'
1845 | by the corresponding value `b'. The operation is performed according to the
1846 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1847 *----------------------------------------------------------------------------*/
1849 float32 float32_div( float32 a, float32 b STATUS_PARAM )
1851 flag aSign, bSign, zSign;
1852 int16 aExp, bExp, zExp;
1853 bits32 aSig, bSig, zSig;
1855 aSig = extractFloat32Frac( a );
1856 aExp = extractFloat32Exp( a );
1857 aSign = extractFloat32Sign( a );
1858 bSig = extractFloat32Frac( b );
1859 bExp = extractFloat32Exp( b );
1860 bSign = extractFloat32Sign( b );
1861 zSign = aSign ^ bSign;
1862 if ( aExp == 0xFF ) {
1863 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1864 if ( bExp == 0xFF ) {
1865 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1866 float_raise( float_flag_invalid STATUS_VAR);
1867 return float32_default_nan;
1869 return packFloat32( zSign, 0xFF, 0 );
1871 if ( bExp == 0xFF ) {
1872 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1873 return packFloat32( zSign, 0, 0 );
1875 if ( bExp == 0 ) {
1876 if ( bSig == 0 ) {
1877 if ( ( aExp | aSig ) == 0 ) {
1878 float_raise( float_flag_invalid STATUS_VAR);
1879 return float32_default_nan;
1881 float_raise( float_flag_divbyzero STATUS_VAR);
1882 return packFloat32( zSign, 0xFF, 0 );
1884 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1886 if ( aExp == 0 ) {
1887 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1888 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1890 zExp = aExp - bExp + 0x7D;
1891 aSig = ( aSig | 0x00800000 )<<7;
1892 bSig = ( bSig | 0x00800000 )<<8;
1893 if ( bSig <= ( aSig + aSig ) ) {
1894 aSig >>= 1;
1895 ++zExp;
1897 zSig = ( ( (bits64) aSig )<<32 ) / bSig;
1898 if ( ( zSig & 0x3F ) == 0 ) {
1899 zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
1901 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1905 /*----------------------------------------------------------------------------
1906 | Returns the remainder of the single-precision floating-point value `a'
1907 | with respect to the corresponding value `b'. The operation is performed
1908 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1909 *----------------------------------------------------------------------------*/
1911 float32 float32_rem( float32 a, float32 b STATUS_PARAM )
1913 flag aSign, bSign, zSign;
1914 int16 aExp, bExp, expDiff;
1915 bits32 aSig, bSig;
1916 bits32 q;
1917 bits64 aSig64, bSig64, q64;
1918 bits32 alternateASig;
1919 sbits32 sigMean;
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 if ( aExp == 0xFF ) {
1928 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1929 return propagateFloat32NaN( a, b STATUS_VAR );
1931 float_raise( float_flag_invalid STATUS_VAR);
1932 return float32_default_nan;
1934 if ( bExp == 0xFF ) {
1935 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1936 return a;
1938 if ( bExp == 0 ) {
1939 if ( bSig == 0 ) {
1940 float_raise( float_flag_invalid STATUS_VAR);
1941 return float32_default_nan;
1943 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1945 if ( aExp == 0 ) {
1946 if ( aSig == 0 ) return a;
1947 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1949 expDiff = aExp - bExp;
1950 aSig |= 0x00800000;
1951 bSig |= 0x00800000;
1952 if ( expDiff < 32 ) {
1953 aSig <<= 8;
1954 bSig <<= 8;
1955 if ( expDiff < 0 ) {
1956 if ( expDiff < -1 ) return a;
1957 aSig >>= 1;
1959 q = ( bSig <= aSig );
1960 if ( q ) aSig -= bSig;
1961 if ( 0 < expDiff ) {
1962 q = ( ( (bits64) aSig )<<32 ) / bSig;
1963 q >>= 32 - expDiff;
1964 bSig >>= 2;
1965 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
1967 else {
1968 aSig >>= 2;
1969 bSig >>= 2;
1972 else {
1973 if ( bSig <= aSig ) aSig -= bSig;
1974 aSig64 = ( (bits64) aSig )<<40;
1975 bSig64 = ( (bits64) bSig )<<40;
1976 expDiff -= 64;
1977 while ( 0 < expDiff ) {
1978 q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1979 q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1980 aSig64 = - ( ( bSig * q64 )<<38 );
1981 expDiff -= 62;
1983 expDiff += 64;
1984 q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1985 q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1986 q = q64>>( 64 - expDiff );
1987 bSig <<= 6;
1988 aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
1990 do {
1991 alternateASig = aSig;
1992 ++q;
1993 aSig -= bSig;
1994 } while ( 0 <= (sbits32) aSig );
1995 sigMean = aSig + alternateASig;
1996 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
1997 aSig = alternateASig;
1999 zSign = ( (sbits32) aSig < 0 );
2000 if ( zSign ) aSig = - aSig;
2001 return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig STATUS_VAR );
2005 /*----------------------------------------------------------------------------
2006 | Returns the square root of the single-precision floating-point value `a'.
2007 | The operation is performed according to the IEC/IEEE Standard for Binary
2008 | Floating-Point Arithmetic.
2009 *----------------------------------------------------------------------------*/
2011 float32 float32_sqrt( float32 a STATUS_PARAM )
2013 flag aSign;
2014 int16 aExp, zExp;
2015 bits32 aSig, zSig;
2016 bits64 rem, term;
2018 aSig = extractFloat32Frac( a );
2019 aExp = extractFloat32Exp( a );
2020 aSign = extractFloat32Sign( a );
2021 if ( aExp == 0xFF ) {
2022 if ( aSig ) return propagateFloat32NaN( a, float32_zero STATUS_VAR );
2023 if ( ! aSign ) return a;
2024 float_raise( float_flag_invalid STATUS_VAR);
2025 return float32_default_nan;
2027 if ( aSign ) {
2028 if ( ( aExp | aSig ) == 0 ) return a;
2029 float_raise( float_flag_invalid STATUS_VAR);
2030 return float32_default_nan;
2032 if ( aExp == 0 ) {
2033 if ( aSig == 0 ) return float32_zero;
2034 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
2036 zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
2037 aSig = ( aSig | 0x00800000 )<<8;
2038 zSig = estimateSqrt32( aExp, aSig ) + 2;
2039 if ( ( zSig & 0x7F ) <= 5 ) {
2040 if ( zSig < 2 ) {
2041 zSig = 0x7FFFFFFF;
2042 goto roundAndPack;
2044 aSig >>= aExp & 1;
2045 term = ( (bits64) zSig ) * zSig;
2046 rem = ( ( (bits64) aSig )<<32 ) - term;
2047 while ( (sbits64) rem < 0 ) {
2048 --zSig;
2049 rem += ( ( (bits64) zSig )<<1 ) | 1;
2051 zSig |= ( rem != 0 );
2053 shift32RightJamming( zSig, 1, &zSig );
2054 roundAndPack:
2055 return roundAndPackFloat32( 0, zExp, zSig STATUS_VAR );
2059 /*----------------------------------------------------------------------------
2060 | Returns the binary log of the single-precision floating-point value `a'.
2061 | The operation is performed according to the IEC/IEEE Standard for Binary
2062 | Floating-Point Arithmetic.
2063 *----------------------------------------------------------------------------*/
2064 float32 float32_log2( float32 a STATUS_PARAM )
2066 flag aSign, zSign;
2067 int16 aExp;
2068 bits32 aSig, zSig, i;
2070 aSig = extractFloat32Frac( a );
2071 aExp = extractFloat32Exp( a );
2072 aSign = extractFloat32Sign( a );
2074 if ( aExp == 0 ) {
2075 if ( aSig == 0 ) return packFloat32( 1, 0xFF, 0 );
2076 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
2078 if ( aSign ) {
2079 float_raise( float_flag_invalid STATUS_VAR);
2080 return float32_default_nan;
2082 if ( aExp == 0xFF ) {
2083 if ( aSig ) return propagateFloat32NaN( a, float32_zero STATUS_VAR );
2084 return a;
2087 aExp -= 0x7F;
2088 aSig |= 0x00800000;
2089 zSign = aExp < 0;
2090 zSig = aExp << 23;
2092 for (i = 1 << 22; i > 0; i >>= 1) {
2093 aSig = ( (bits64)aSig * aSig ) >> 23;
2094 if ( aSig & 0x01000000 ) {
2095 aSig >>= 1;
2096 zSig |= i;
2100 if ( zSign )
2101 zSig = -zSig;
2103 return normalizeRoundAndPackFloat32( zSign, 0x85, zSig STATUS_VAR );
2106 /*----------------------------------------------------------------------------
2107 | Returns 1 if the single-precision floating-point value `a' is equal to
2108 | the corresponding value `b', and 0 otherwise. The comparison is performed
2109 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2110 *----------------------------------------------------------------------------*/
2112 int float32_eq( float32 a, float32 b STATUS_PARAM )
2115 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2116 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2118 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2119 float_raise( float_flag_invalid STATUS_VAR);
2121 return 0;
2123 return ( float32_val(a) == float32_val(b) ) ||
2124 ( (bits32) ( ( float32_val(a) | float32_val(b) )<<1 ) == 0 );
2128 /*----------------------------------------------------------------------------
2129 | Returns 1 if the single-precision floating-point value `a' is less than
2130 | or equal to the corresponding value `b', and 0 otherwise. The comparison
2131 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2132 | Arithmetic.
2133 *----------------------------------------------------------------------------*/
2135 int float32_le( float32 a, float32 b STATUS_PARAM )
2137 flag aSign, bSign;
2138 bits32 av, bv;
2140 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2141 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2143 float_raise( float_flag_invalid STATUS_VAR);
2144 return 0;
2146 aSign = extractFloat32Sign( a );
2147 bSign = extractFloat32Sign( b );
2148 av = float32_val(a);
2149 bv = float32_val(b);
2150 if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
2151 return ( av == bv ) || ( aSign ^ ( av < bv ) );
2155 /*----------------------------------------------------------------------------
2156 | Returns 1 if the single-precision floating-point value `a' is less than
2157 | the corresponding value `b', and 0 otherwise. The comparison is performed
2158 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2159 *----------------------------------------------------------------------------*/
2161 int float32_lt( float32 a, float32 b STATUS_PARAM )
2163 flag aSign, bSign;
2164 bits32 av, bv;
2166 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2167 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2169 float_raise( float_flag_invalid STATUS_VAR);
2170 return 0;
2172 aSign = extractFloat32Sign( a );
2173 bSign = extractFloat32Sign( b );
2174 av = float32_val(a);
2175 bv = float32_val(b);
2176 if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
2177 return ( av != bv ) && ( aSign ^ ( av < bv ) );
2181 /*----------------------------------------------------------------------------
2182 | Returns 1 if the single-precision floating-point value `a' is equal to
2183 | the corresponding value `b', and 0 otherwise. The invalid exception is
2184 | raised if either operand is a NaN. Otherwise, the comparison is performed
2185 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2186 *----------------------------------------------------------------------------*/
2188 int float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
2190 bits32 av, bv;
2192 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2193 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2195 float_raise( float_flag_invalid STATUS_VAR);
2196 return 0;
2198 av = float32_val(a);
2199 bv = float32_val(b);
2200 return ( av == bv ) || ( (bits32) ( ( av | bv )<<1 ) == 0 );
2204 /*----------------------------------------------------------------------------
2205 | Returns 1 if the single-precision floating-point value `a' is less than or
2206 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
2207 | cause an exception. Otherwise, the comparison is performed according to the
2208 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2209 *----------------------------------------------------------------------------*/
2211 int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
2213 flag aSign, bSign;
2214 bits32 av, bv;
2216 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2217 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2219 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2220 float_raise( float_flag_invalid STATUS_VAR);
2222 return 0;
2224 aSign = extractFloat32Sign( a );
2225 bSign = extractFloat32Sign( b );
2226 av = float32_val(a);
2227 bv = float32_val(b);
2228 if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
2229 return ( av == bv ) || ( aSign ^ ( av < bv ) );
2233 /*----------------------------------------------------------------------------
2234 | Returns 1 if the single-precision floating-point value `a' is less than
2235 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2236 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
2237 | Standard for Binary Floating-Point Arithmetic.
2238 *----------------------------------------------------------------------------*/
2240 int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
2242 flag aSign, bSign;
2243 bits32 av, bv;
2245 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2246 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2248 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2249 float_raise( float_flag_invalid STATUS_VAR);
2251 return 0;
2253 aSign = extractFloat32Sign( a );
2254 bSign = extractFloat32Sign( b );
2255 av = float32_val(a);
2256 bv = float32_val(b);
2257 if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
2258 return ( av != bv ) && ( aSign ^ ( av < bv ) );
2262 /*----------------------------------------------------------------------------
2263 | Returns the result of converting the double-precision floating-point value
2264 | `a' to the 32-bit two's complement integer format. The conversion is
2265 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2266 | Arithmetic---which means in particular that the conversion is rounded
2267 | according to the current rounding mode. If `a' is a NaN, the largest
2268 | positive integer is returned. Otherwise, if the conversion overflows, the
2269 | largest integer with the same sign as `a' is returned.
2270 *----------------------------------------------------------------------------*/
2272 int32 float64_to_int32( float64 a STATUS_PARAM )
2274 flag aSign;
2275 int16 aExp, shiftCount;
2276 bits64 aSig;
2278 aSig = extractFloat64Frac( a );
2279 aExp = extractFloat64Exp( a );
2280 aSign = extractFloat64Sign( a );
2281 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2282 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2283 shiftCount = 0x42C - aExp;
2284 if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
2285 return roundAndPackInt32( aSign, aSig STATUS_VAR );
2289 /*----------------------------------------------------------------------------
2290 | Returns the result of converting the double-precision floating-point value
2291 | `a' to the 32-bit two's complement integer format. The conversion is
2292 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2293 | Arithmetic, except that the conversion is always rounded toward zero.
2294 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2295 | the conversion overflows, the largest integer with the same sign as `a' is
2296 | returned.
2297 *----------------------------------------------------------------------------*/
2299 int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
2301 flag aSign;
2302 int16 aExp, shiftCount;
2303 bits64 aSig, savedASig;
2304 int32 z;
2306 aSig = extractFloat64Frac( a );
2307 aExp = extractFloat64Exp( a );
2308 aSign = extractFloat64Sign( a );
2309 if ( 0x41E < aExp ) {
2310 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2311 goto invalid;
2313 else if ( aExp < 0x3FF ) {
2314 if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2315 return 0;
2317 aSig |= LIT64( 0x0010000000000000 );
2318 shiftCount = 0x433 - aExp;
2319 savedASig = aSig;
2320 aSig >>= shiftCount;
2321 z = aSig;
2322 if ( aSign ) z = - z;
2323 if ( ( z < 0 ) ^ aSign ) {
2324 invalid:
2325 float_raise( float_flag_invalid STATUS_VAR);
2326 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
2328 if ( ( aSig<<shiftCount ) != savedASig ) {
2329 STATUS(float_exception_flags) |= float_flag_inexact;
2331 return z;
2335 /*----------------------------------------------------------------------------
2336 | Returns the result of converting the double-precision floating-point value
2337 | `a' to the 64-bit two's complement integer format. The conversion is
2338 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2339 | Arithmetic---which means in particular that the conversion is rounded
2340 | according to the current rounding mode. If `a' is a NaN, the largest
2341 | positive integer is returned. Otherwise, if the conversion overflows, the
2342 | largest integer with the same sign as `a' is returned.
2343 *----------------------------------------------------------------------------*/
2345 int64 float64_to_int64( float64 a STATUS_PARAM )
2347 flag aSign;
2348 int16 aExp, shiftCount;
2349 bits64 aSig, aSigExtra;
2351 aSig = extractFloat64Frac( a );
2352 aExp = extractFloat64Exp( a );
2353 aSign = extractFloat64Sign( a );
2354 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2355 shiftCount = 0x433 - aExp;
2356 if ( shiftCount <= 0 ) {
2357 if ( 0x43E < aExp ) {
2358 float_raise( float_flag_invalid STATUS_VAR);
2359 if ( ! aSign
2360 || ( ( aExp == 0x7FF )
2361 && ( aSig != LIT64( 0x0010000000000000 ) ) )
2363 return LIT64( 0x7FFFFFFFFFFFFFFF );
2365 return (sbits64) LIT64( 0x8000000000000000 );
2367 aSigExtra = 0;
2368 aSig <<= - shiftCount;
2370 else {
2371 shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
2373 return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
2377 /*----------------------------------------------------------------------------
2378 | Returns the result of converting the double-precision floating-point value
2379 | `a' to the 64-bit two's complement integer format. The conversion is
2380 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2381 | Arithmetic, except that the conversion is always rounded toward zero.
2382 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2383 | the conversion overflows, the largest integer with the same sign as `a' is
2384 | returned.
2385 *----------------------------------------------------------------------------*/
2387 int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
2389 flag aSign;
2390 int16 aExp, shiftCount;
2391 bits64 aSig;
2392 int64 z;
2394 aSig = extractFloat64Frac( a );
2395 aExp = extractFloat64Exp( a );
2396 aSign = extractFloat64Sign( a );
2397 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2398 shiftCount = aExp - 0x433;
2399 if ( 0 <= shiftCount ) {
2400 if ( 0x43E <= aExp ) {
2401 if ( float64_val(a) != LIT64( 0xC3E0000000000000 ) ) {
2402 float_raise( float_flag_invalid STATUS_VAR);
2403 if ( ! aSign
2404 || ( ( aExp == 0x7FF )
2405 && ( aSig != LIT64( 0x0010000000000000 ) ) )
2407 return LIT64( 0x7FFFFFFFFFFFFFFF );
2410 return (sbits64) LIT64( 0x8000000000000000 );
2412 z = aSig<<shiftCount;
2414 else {
2415 if ( aExp < 0x3FE ) {
2416 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2417 return 0;
2419 z = aSig>>( - shiftCount );
2420 if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
2421 STATUS(float_exception_flags) |= float_flag_inexact;
2424 if ( aSign ) z = - z;
2425 return z;
2429 /*----------------------------------------------------------------------------
2430 | Returns the result of converting the double-precision floating-point value
2431 | `a' to the single-precision floating-point format. The conversion is
2432 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2433 | Arithmetic.
2434 *----------------------------------------------------------------------------*/
2436 float32 float64_to_float32( float64 a STATUS_PARAM )
2438 flag aSign;
2439 int16 aExp;
2440 bits64 aSig;
2441 bits32 zSig;
2443 aSig = extractFloat64Frac( a );
2444 aExp = extractFloat64Exp( a );
2445 aSign = extractFloat64Sign( a );
2446 if ( aExp == 0x7FF ) {
2447 if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR ) );
2448 return packFloat32( aSign, 0xFF, 0 );
2450 shift64RightJamming( aSig, 22, &aSig );
2451 zSig = aSig;
2452 if ( aExp || zSig ) {
2453 zSig |= 0x40000000;
2454 aExp -= 0x381;
2456 return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
2460 #ifdef FLOATX80
2462 /*----------------------------------------------------------------------------
2463 | Returns the result of converting the double-precision floating-point value
2464 | `a' to the extended double-precision floating-point format. The conversion
2465 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2466 | Arithmetic.
2467 *----------------------------------------------------------------------------*/
2469 floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
2471 flag aSign;
2472 int16 aExp;
2473 bits64 aSig;
2475 aSig = extractFloat64Frac( a );
2476 aExp = extractFloat64Exp( a );
2477 aSign = extractFloat64Sign( a );
2478 if ( aExp == 0x7FF ) {
2479 if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) );
2480 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
2482 if ( aExp == 0 ) {
2483 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
2484 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2486 return
2487 packFloatx80(
2488 aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
2492 #endif
2494 #ifdef FLOAT128
2496 /*----------------------------------------------------------------------------
2497 | Returns the result of converting the double-precision floating-point value
2498 | `a' to the quadruple-precision floating-point format. The conversion is
2499 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2500 | Arithmetic.
2501 *----------------------------------------------------------------------------*/
2503 float128 float64_to_float128( float64 a STATUS_PARAM )
2505 flag aSign;
2506 int16 aExp;
2507 bits64 aSig, zSig0, zSig1;
2509 aSig = extractFloat64Frac( a );
2510 aExp = extractFloat64Exp( a );
2511 aSign = extractFloat64Sign( a );
2512 if ( aExp == 0x7FF ) {
2513 if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR ) );
2514 return packFloat128( aSign, 0x7FFF, 0, 0 );
2516 if ( aExp == 0 ) {
2517 if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
2518 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2519 --aExp;
2521 shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
2522 return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
2526 #endif
2528 /*----------------------------------------------------------------------------
2529 | Rounds the double-precision floating-point value `a' to an integer, and
2530 | returns the result as a double-precision floating-point value. The
2531 | operation is performed according to the IEC/IEEE Standard for Binary
2532 | Floating-Point Arithmetic.
2533 *----------------------------------------------------------------------------*/
2535 float64 float64_round_to_int( float64 a STATUS_PARAM )
2537 flag aSign;
2538 int16 aExp;
2539 bits64 lastBitMask, roundBitsMask;
2540 int8 roundingMode;
2541 bits64 z;
2543 aExp = extractFloat64Exp( a );
2544 if ( 0x433 <= aExp ) {
2545 if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
2546 return propagateFloat64NaN( a, a STATUS_VAR );
2548 return a;
2550 if ( aExp < 0x3FF ) {
2551 if ( (bits64) ( float64_val(a)<<1 ) == 0 ) return a;
2552 STATUS(float_exception_flags) |= float_flag_inexact;
2553 aSign = extractFloat64Sign( a );
2554 switch ( STATUS(float_rounding_mode) ) {
2555 case float_round_nearest_even:
2556 if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
2557 return packFloat64( aSign, 0x3FF, 0 );
2559 break;
2560 case float_round_down:
2561 return make_float64(aSign ? LIT64( 0xBFF0000000000000 ) : 0);
2562 case float_round_up:
2563 return make_float64(
2564 aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 ));
2566 return packFloat64( aSign, 0, 0 );
2568 lastBitMask = 1;
2569 lastBitMask <<= 0x433 - aExp;
2570 roundBitsMask = lastBitMask - 1;
2571 z = float64_val(a);
2572 roundingMode = STATUS(float_rounding_mode);
2573 if ( roundingMode == float_round_nearest_even ) {
2574 z += lastBitMask>>1;
2575 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
2577 else if ( roundingMode != float_round_to_zero ) {
2578 if ( extractFloat64Sign( make_float64(z) ) ^ ( roundingMode == float_round_up ) ) {
2579 z += roundBitsMask;
2582 z &= ~ roundBitsMask;
2583 if ( z != float64_val(a) )
2584 STATUS(float_exception_flags) |= float_flag_inexact;
2585 return make_float64(z);
2589 float64 float64_trunc_to_int( float64 a STATUS_PARAM)
2591 int oldmode;
2592 float64 res;
2593 oldmode = STATUS(float_rounding_mode);
2594 STATUS(float_rounding_mode) = float_round_to_zero;
2595 res = float64_round_to_int(a STATUS_VAR);
2596 STATUS(float_rounding_mode) = oldmode;
2597 return res;
2600 /*----------------------------------------------------------------------------
2601 | Returns the result of adding the absolute values of the double-precision
2602 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
2603 | before being returned. `zSign' is ignored if the result is a NaN.
2604 | The addition is performed according to the IEC/IEEE Standard for Binary
2605 | Floating-Point Arithmetic.
2606 *----------------------------------------------------------------------------*/
2608 static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2610 int16 aExp, bExp, zExp;
2611 bits64 aSig, bSig, zSig;
2612 int16 expDiff;
2614 aSig = extractFloat64Frac( a );
2615 aExp = extractFloat64Exp( a );
2616 bSig = extractFloat64Frac( b );
2617 bExp = extractFloat64Exp( b );
2618 expDiff = aExp - bExp;
2619 aSig <<= 9;
2620 bSig <<= 9;
2621 if ( 0 < expDiff ) {
2622 if ( aExp == 0x7FF ) {
2623 if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2624 return a;
2626 if ( bExp == 0 ) {
2627 --expDiff;
2629 else {
2630 bSig |= LIT64( 0x2000000000000000 );
2632 shift64RightJamming( bSig, expDiff, &bSig );
2633 zExp = aExp;
2635 else if ( expDiff < 0 ) {
2636 if ( bExp == 0x7FF ) {
2637 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2638 return packFloat64( zSign, 0x7FF, 0 );
2640 if ( aExp == 0 ) {
2641 ++expDiff;
2643 else {
2644 aSig |= LIT64( 0x2000000000000000 );
2646 shift64RightJamming( aSig, - expDiff, &aSig );
2647 zExp = bExp;
2649 else {
2650 if ( aExp == 0x7FF ) {
2651 if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2652 return a;
2654 if ( aExp == 0 ) {
2655 if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
2656 return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
2658 zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
2659 zExp = aExp;
2660 goto roundAndPack;
2662 aSig |= LIT64( 0x2000000000000000 );
2663 zSig = ( aSig + bSig )<<1;
2664 --zExp;
2665 if ( (sbits64) zSig < 0 ) {
2666 zSig = aSig + bSig;
2667 ++zExp;
2669 roundAndPack:
2670 return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2674 /*----------------------------------------------------------------------------
2675 | Returns the result of subtracting the absolute values of the double-
2676 | precision floating-point values `a' and `b'. If `zSign' is 1, the
2677 | difference is negated before being returned. `zSign' is ignored if the
2678 | result is a NaN. The subtraction is performed according to the IEC/IEEE
2679 | Standard for Binary Floating-Point Arithmetic.
2680 *----------------------------------------------------------------------------*/
2682 static float64 subFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2684 int16 aExp, bExp, zExp;
2685 bits64 aSig, bSig, zSig;
2686 int16 expDiff;
2688 aSig = extractFloat64Frac( a );
2689 aExp = extractFloat64Exp( a );
2690 bSig = extractFloat64Frac( b );
2691 bExp = extractFloat64Exp( b );
2692 expDiff = aExp - bExp;
2693 aSig <<= 10;
2694 bSig <<= 10;
2695 if ( 0 < expDiff ) goto aExpBigger;
2696 if ( expDiff < 0 ) goto bExpBigger;
2697 if ( aExp == 0x7FF ) {
2698 if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2699 float_raise( float_flag_invalid STATUS_VAR);
2700 return float64_default_nan;
2702 if ( aExp == 0 ) {
2703 aExp = 1;
2704 bExp = 1;
2706 if ( bSig < aSig ) goto aBigger;
2707 if ( aSig < bSig ) goto bBigger;
2708 return packFloat64( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
2709 bExpBigger:
2710 if ( bExp == 0x7FF ) {
2711 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2712 return packFloat64( zSign ^ 1, 0x7FF, 0 );
2714 if ( aExp == 0 ) {
2715 ++expDiff;
2717 else {
2718 aSig |= LIT64( 0x4000000000000000 );
2720 shift64RightJamming( aSig, - expDiff, &aSig );
2721 bSig |= LIT64( 0x4000000000000000 );
2722 bBigger:
2723 zSig = bSig - aSig;
2724 zExp = bExp;
2725 zSign ^= 1;
2726 goto normalizeRoundAndPack;
2727 aExpBigger:
2728 if ( aExp == 0x7FF ) {
2729 if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2730 return a;
2732 if ( bExp == 0 ) {
2733 --expDiff;
2735 else {
2736 bSig |= LIT64( 0x4000000000000000 );
2738 shift64RightJamming( bSig, expDiff, &bSig );
2739 aSig |= LIT64( 0x4000000000000000 );
2740 aBigger:
2741 zSig = aSig - bSig;
2742 zExp = aExp;
2743 normalizeRoundAndPack:
2744 --zExp;
2745 return normalizeRoundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2749 /*----------------------------------------------------------------------------
2750 | Returns the result of adding the double-precision floating-point values `a'
2751 | and `b'. The operation is performed according to the IEC/IEEE Standard for
2752 | Binary Floating-Point Arithmetic.
2753 *----------------------------------------------------------------------------*/
2755 float64 float64_add( float64 a, float64 b STATUS_PARAM )
2757 flag aSign, bSign;
2759 aSign = extractFloat64Sign( a );
2760 bSign = extractFloat64Sign( b );
2761 if ( aSign == bSign ) {
2762 return addFloat64Sigs( a, b, aSign STATUS_VAR );
2764 else {
2765 return subFloat64Sigs( a, b, aSign STATUS_VAR );
2770 /*----------------------------------------------------------------------------
2771 | Returns the result of subtracting the double-precision floating-point values
2772 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2773 | for Binary Floating-Point Arithmetic.
2774 *----------------------------------------------------------------------------*/
2776 float64 float64_sub( float64 a, float64 b STATUS_PARAM )
2778 flag aSign, bSign;
2780 aSign = extractFloat64Sign( a );
2781 bSign = extractFloat64Sign( b );
2782 if ( aSign == bSign ) {
2783 return subFloat64Sigs( a, b, aSign STATUS_VAR );
2785 else {
2786 return addFloat64Sigs( a, b, aSign STATUS_VAR );
2791 /*----------------------------------------------------------------------------
2792 | Returns the result of multiplying the double-precision floating-point values
2793 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2794 | for Binary Floating-Point Arithmetic.
2795 *----------------------------------------------------------------------------*/
2797 float64 float64_mul( float64 a, float64 b STATUS_PARAM )
2799 flag aSign, bSign, zSign;
2800 int16 aExp, bExp, zExp;
2801 bits64 aSig, bSig, zSig0, zSig1;
2803 aSig = extractFloat64Frac( a );
2804 aExp = extractFloat64Exp( a );
2805 aSign = extractFloat64Sign( a );
2806 bSig = extractFloat64Frac( b );
2807 bExp = extractFloat64Exp( b );
2808 bSign = extractFloat64Sign( b );
2809 zSign = aSign ^ bSign;
2810 if ( aExp == 0x7FF ) {
2811 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2812 return propagateFloat64NaN( a, b STATUS_VAR );
2814 if ( ( bExp | bSig ) == 0 ) {
2815 float_raise( float_flag_invalid STATUS_VAR);
2816 return float64_default_nan;
2818 return packFloat64( zSign, 0x7FF, 0 );
2820 if ( bExp == 0x7FF ) {
2821 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2822 if ( ( aExp | aSig ) == 0 ) {
2823 float_raise( float_flag_invalid STATUS_VAR);
2824 return float64_default_nan;
2826 return packFloat64( zSign, 0x7FF, 0 );
2828 if ( aExp == 0 ) {
2829 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2830 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2832 if ( bExp == 0 ) {
2833 if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
2834 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2836 zExp = aExp + bExp - 0x3FF;
2837 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2838 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2839 mul64To128( aSig, bSig, &zSig0, &zSig1 );
2840 zSig0 |= ( zSig1 != 0 );
2841 if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
2842 zSig0 <<= 1;
2843 --zExp;
2845 return roundAndPackFloat64( zSign, zExp, zSig0 STATUS_VAR );
2849 /*----------------------------------------------------------------------------
2850 | Returns the result of dividing the double-precision floating-point value `a'
2851 | by the corresponding value `b'. The operation is performed according to
2852 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2853 *----------------------------------------------------------------------------*/
2855 float64 float64_div( float64 a, float64 b STATUS_PARAM )
2857 flag aSign, bSign, zSign;
2858 int16 aExp, bExp, zExp;
2859 bits64 aSig, bSig, zSig;
2860 bits64 rem0, rem1;
2861 bits64 term0, term1;
2863 aSig = extractFloat64Frac( a );
2864 aExp = extractFloat64Exp( a );
2865 aSign = extractFloat64Sign( a );
2866 bSig = extractFloat64Frac( b );
2867 bExp = extractFloat64Exp( b );
2868 bSign = extractFloat64Sign( b );
2869 zSign = aSign ^ bSign;
2870 if ( aExp == 0x7FF ) {
2871 if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2872 if ( bExp == 0x7FF ) {
2873 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2874 float_raise( float_flag_invalid STATUS_VAR);
2875 return float64_default_nan;
2877 return packFloat64( zSign, 0x7FF, 0 );
2879 if ( bExp == 0x7FF ) {
2880 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2881 return packFloat64( zSign, 0, 0 );
2883 if ( bExp == 0 ) {
2884 if ( bSig == 0 ) {
2885 if ( ( aExp | aSig ) == 0 ) {
2886 float_raise( float_flag_invalid STATUS_VAR);
2887 return float64_default_nan;
2889 float_raise( float_flag_divbyzero STATUS_VAR);
2890 return packFloat64( zSign, 0x7FF, 0 );
2892 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2894 if ( aExp == 0 ) {
2895 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2896 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2898 zExp = aExp - bExp + 0x3FD;
2899 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2900 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2901 if ( bSig <= ( aSig + aSig ) ) {
2902 aSig >>= 1;
2903 ++zExp;
2905 zSig = estimateDiv128To64( aSig, 0, bSig );
2906 if ( ( zSig & 0x1FF ) <= 2 ) {
2907 mul64To128( bSig, zSig, &term0, &term1 );
2908 sub128( aSig, 0, term0, term1, &rem0, &rem1 );
2909 while ( (sbits64) rem0 < 0 ) {
2910 --zSig;
2911 add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
2913 zSig |= ( rem1 != 0 );
2915 return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2919 /*----------------------------------------------------------------------------
2920 | Returns the remainder of the double-precision floating-point value `a'
2921 | with respect to the corresponding value `b'. The operation is performed
2922 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2923 *----------------------------------------------------------------------------*/
2925 float64 float64_rem( float64 a, float64 b STATUS_PARAM )
2927 flag aSign, bSign, zSign;
2928 int16 aExp, bExp, expDiff;
2929 bits64 aSig, bSig;
2930 bits64 q, alternateASig;
2931 sbits64 sigMean;
2933 aSig = extractFloat64Frac( a );
2934 aExp = extractFloat64Exp( a );
2935 aSign = extractFloat64Sign( a );
2936 bSig = extractFloat64Frac( b );
2937 bExp = extractFloat64Exp( b );
2938 bSign = extractFloat64Sign( b );
2939 if ( aExp == 0x7FF ) {
2940 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2941 return propagateFloat64NaN( a, b STATUS_VAR );
2943 float_raise( float_flag_invalid STATUS_VAR);
2944 return float64_default_nan;
2946 if ( bExp == 0x7FF ) {
2947 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2948 return a;
2950 if ( bExp == 0 ) {
2951 if ( bSig == 0 ) {
2952 float_raise( float_flag_invalid STATUS_VAR);
2953 return float64_default_nan;
2955 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2957 if ( aExp == 0 ) {
2958 if ( aSig == 0 ) return a;
2959 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2961 expDiff = aExp - bExp;
2962 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
2963 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2964 if ( expDiff < 0 ) {
2965 if ( expDiff < -1 ) return a;
2966 aSig >>= 1;
2968 q = ( bSig <= aSig );
2969 if ( q ) aSig -= bSig;
2970 expDiff -= 64;
2971 while ( 0 < expDiff ) {
2972 q = estimateDiv128To64( aSig, 0, bSig );
2973 q = ( 2 < q ) ? q - 2 : 0;
2974 aSig = - ( ( bSig>>2 ) * q );
2975 expDiff -= 62;
2977 expDiff += 64;
2978 if ( 0 < expDiff ) {
2979 q = estimateDiv128To64( aSig, 0, bSig );
2980 q = ( 2 < q ) ? q - 2 : 0;
2981 q >>= 64 - expDiff;
2982 bSig >>= 2;
2983 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
2985 else {
2986 aSig >>= 2;
2987 bSig >>= 2;
2989 do {
2990 alternateASig = aSig;
2991 ++q;
2992 aSig -= bSig;
2993 } while ( 0 <= (sbits64) aSig );
2994 sigMean = aSig + alternateASig;
2995 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
2996 aSig = alternateASig;
2998 zSign = ( (sbits64) aSig < 0 );
2999 if ( zSign ) aSig = - aSig;
3000 return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig STATUS_VAR );
3004 /*----------------------------------------------------------------------------
3005 | Returns the square root of the double-precision floating-point value `a'.
3006 | The operation is performed according to the IEC/IEEE Standard for Binary
3007 | Floating-Point Arithmetic.
3008 *----------------------------------------------------------------------------*/
3010 float64 float64_sqrt( float64 a STATUS_PARAM )
3012 flag aSign;
3013 int16 aExp, zExp;
3014 bits64 aSig, zSig, doubleZSig;
3015 bits64 rem0, rem1, term0, term1;
3017 aSig = extractFloat64Frac( a );
3018 aExp = extractFloat64Exp( a );
3019 aSign = extractFloat64Sign( a );
3020 if ( aExp == 0x7FF ) {
3021 if ( aSig ) return propagateFloat64NaN( a, a STATUS_VAR );
3022 if ( ! aSign ) return a;
3023 float_raise( float_flag_invalid STATUS_VAR);
3024 return float64_default_nan;
3026 if ( aSign ) {
3027 if ( ( aExp | aSig ) == 0 ) return a;
3028 float_raise( float_flag_invalid STATUS_VAR);
3029 return float64_default_nan;
3031 if ( aExp == 0 ) {
3032 if ( aSig == 0 ) return float64_zero;
3033 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
3035 zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
3036 aSig |= LIT64( 0x0010000000000000 );
3037 zSig = estimateSqrt32( aExp, aSig>>21 );
3038 aSig <<= 9 - ( aExp & 1 );
3039 zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
3040 if ( ( zSig & 0x1FF ) <= 5 ) {
3041 doubleZSig = zSig<<1;
3042 mul64To128( zSig, zSig, &term0, &term1 );
3043 sub128( aSig, 0, term0, term1, &rem0, &rem1 );
3044 while ( (sbits64) rem0 < 0 ) {
3045 --zSig;
3046 doubleZSig -= 2;
3047 add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
3049 zSig |= ( ( rem0 | rem1 ) != 0 );
3051 return roundAndPackFloat64( 0, zExp, zSig STATUS_VAR );
3055 /*----------------------------------------------------------------------------
3056 | Returns the binary log of the double-precision floating-point value `a'.
3057 | The operation is performed according to the IEC/IEEE Standard for Binary
3058 | Floating-Point Arithmetic.
3059 *----------------------------------------------------------------------------*/
3060 float64 float64_log2( float64 a STATUS_PARAM )
3062 flag aSign, zSign;
3063 int16 aExp;
3064 bits64 aSig, aSig0, aSig1, zSig, i;
3066 aSig = extractFloat64Frac( a );
3067 aExp = extractFloat64Exp( a );
3068 aSign = extractFloat64Sign( a );
3070 if ( aExp == 0 ) {
3071 if ( aSig == 0 ) return packFloat64( 1, 0x7FF, 0 );
3072 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
3074 if ( aSign ) {
3075 float_raise( float_flag_invalid STATUS_VAR);
3076 return float64_default_nan;
3078 if ( aExp == 0x7FF ) {
3079 if ( aSig ) return propagateFloat64NaN( a, float64_zero STATUS_VAR );
3080 return a;
3083 aExp -= 0x3FF;
3084 aSig |= LIT64( 0x0010000000000000 );
3085 zSign = aExp < 0;
3086 zSig = (bits64)aExp << 52;
3087 for (i = 1LL << 51; i > 0; i >>= 1) {
3088 mul64To128( aSig, aSig, &aSig0, &aSig1 );
3089 aSig = ( aSig0 << 12 ) | ( aSig1 >> 52 );
3090 if ( aSig & LIT64( 0x0020000000000000 ) ) {
3091 aSig >>= 1;
3092 zSig |= i;
3096 if ( zSign )
3097 zSig = -zSig;
3098 return normalizeRoundAndPackFloat64( zSign, 0x408, zSig STATUS_VAR );
3101 /*----------------------------------------------------------------------------
3102 | Returns 1 if the double-precision floating-point value `a' is equal to the
3103 | corresponding value `b', and 0 otherwise. The comparison is performed
3104 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3105 *----------------------------------------------------------------------------*/
3107 int float64_eq( float64 a, float64 b STATUS_PARAM )
3109 bits64 av, bv;
3111 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3112 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3114 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3115 float_raise( float_flag_invalid STATUS_VAR);
3117 return 0;
3119 av = float64_val(a);
3120 bv = float64_val(b);
3121 return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
3125 /*----------------------------------------------------------------------------
3126 | Returns 1 if the double-precision floating-point value `a' is less than or
3127 | equal to the corresponding value `b', and 0 otherwise. The comparison is
3128 | performed according to the IEC/IEEE Standard for Binary Floating-Point
3129 | Arithmetic.
3130 *----------------------------------------------------------------------------*/
3132 int float64_le( float64 a, float64 b STATUS_PARAM )
3134 flag aSign, bSign;
3135 bits64 av, bv;
3137 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3138 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3140 float_raise( float_flag_invalid STATUS_VAR);
3141 return 0;
3143 aSign = extractFloat64Sign( a );
3144 bSign = extractFloat64Sign( b );
3145 av = float64_val(a);
3146 bv = float64_val(b);
3147 if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
3148 return ( av == bv ) || ( aSign ^ ( av < bv ) );
3152 /*----------------------------------------------------------------------------
3153 | Returns 1 if the double-precision floating-point value `a' is less than
3154 | the corresponding value `b', and 0 otherwise. The comparison is performed
3155 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3156 *----------------------------------------------------------------------------*/
3158 int float64_lt( float64 a, float64 b STATUS_PARAM )
3160 flag aSign, bSign;
3161 bits64 av, bv;
3163 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3164 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3166 float_raise( float_flag_invalid STATUS_VAR);
3167 return 0;
3169 aSign = extractFloat64Sign( a );
3170 bSign = extractFloat64Sign( b );
3171 av = float64_val(a);
3172 bv = float64_val(b);
3173 if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
3174 return ( av != bv ) && ( aSign ^ ( av < bv ) );
3178 /*----------------------------------------------------------------------------
3179 | Returns 1 if the double-precision floating-point value `a' is equal to the
3180 | corresponding value `b', and 0 otherwise. The invalid exception is raised
3181 | if either operand is a NaN. Otherwise, the comparison is performed
3182 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3183 *----------------------------------------------------------------------------*/
3185 int float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
3187 bits64 av, bv;
3189 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3190 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3192 float_raise( float_flag_invalid STATUS_VAR);
3193 return 0;
3195 av = float64_val(a);
3196 bv = float64_val(b);
3197 return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
3201 /*----------------------------------------------------------------------------
3202 | Returns 1 if the double-precision floating-point value `a' is less than or
3203 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
3204 | cause an exception. Otherwise, the comparison is performed according to the
3205 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3206 *----------------------------------------------------------------------------*/
3208 int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
3210 flag aSign, bSign;
3211 bits64 av, bv;
3213 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3214 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3216 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3217 float_raise( float_flag_invalid STATUS_VAR);
3219 return 0;
3221 aSign = extractFloat64Sign( a );
3222 bSign = extractFloat64Sign( b );
3223 av = float64_val(a);
3224 bv = float64_val(b);
3225 if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
3226 return ( av == bv ) || ( aSign ^ ( av < bv ) );
3230 /*----------------------------------------------------------------------------
3231 | Returns 1 if the double-precision floating-point value `a' is less than
3232 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
3233 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
3234 | Standard for Binary Floating-Point Arithmetic.
3235 *----------------------------------------------------------------------------*/
3237 int float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
3239 flag aSign, bSign;
3240 bits64 av, bv;
3242 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3243 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3245 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3246 float_raise( float_flag_invalid STATUS_VAR);
3248 return 0;
3250 aSign = extractFloat64Sign( a );
3251 bSign = extractFloat64Sign( b );
3252 av = float64_val(a);
3253 bv = float64_val(b);
3254 if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
3255 return ( av != bv ) && ( aSign ^ ( av < bv ) );
3259 #ifdef FLOATX80
3261 /*----------------------------------------------------------------------------
3262 | Returns the result of converting the extended double-precision floating-
3263 | point value `a' to the 32-bit two's complement integer format. The
3264 | conversion is performed according to the IEC/IEEE Standard for Binary
3265 | Floating-Point Arithmetic---which means in particular that the conversion
3266 | is rounded according to the current rounding mode. If `a' is a NaN, the
3267 | largest positive integer is returned. Otherwise, if the conversion
3268 | overflows, the largest integer with the same sign as `a' is returned.
3269 *----------------------------------------------------------------------------*/
3271 int32 floatx80_to_int32( floatx80 a STATUS_PARAM )
3273 flag aSign;
3274 int32 aExp, shiftCount;
3275 bits64 aSig;
3277 aSig = extractFloatx80Frac( a );
3278 aExp = extractFloatx80Exp( a );
3279 aSign = extractFloatx80Sign( a );
3280 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3281 shiftCount = 0x4037 - aExp;
3282 if ( shiftCount <= 0 ) shiftCount = 1;
3283 shift64RightJamming( aSig, shiftCount, &aSig );
3284 return roundAndPackInt32( aSign, aSig STATUS_VAR );
3288 /*----------------------------------------------------------------------------
3289 | Returns the result of converting the extended double-precision floating-
3290 | point value `a' to the 32-bit two's complement integer format. The
3291 | conversion is performed according to the IEC/IEEE Standard for Binary
3292 | Floating-Point Arithmetic, except that the conversion is always rounded
3293 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3294 | Otherwise, if the conversion overflows, the largest integer with the same
3295 | sign as `a' is returned.
3296 *----------------------------------------------------------------------------*/
3298 int32 floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM )
3300 flag aSign;
3301 int32 aExp, shiftCount;
3302 bits64 aSig, savedASig;
3303 int32 z;
3305 aSig = extractFloatx80Frac( a );
3306 aExp = extractFloatx80Exp( a );
3307 aSign = extractFloatx80Sign( a );
3308 if ( 0x401E < aExp ) {
3309 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3310 goto invalid;
3312 else if ( aExp < 0x3FFF ) {
3313 if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3314 return 0;
3316 shiftCount = 0x403E - aExp;
3317 savedASig = aSig;
3318 aSig >>= shiftCount;
3319 z = aSig;
3320 if ( aSign ) z = - z;
3321 if ( ( z < 0 ) ^ aSign ) {
3322 invalid:
3323 float_raise( float_flag_invalid STATUS_VAR);
3324 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
3326 if ( ( aSig<<shiftCount ) != savedASig ) {
3327 STATUS(float_exception_flags) |= float_flag_inexact;
3329 return z;
3333 /*----------------------------------------------------------------------------
3334 | Returns the result of converting the extended double-precision floating-
3335 | point value `a' to the 64-bit two's complement integer format. The
3336 | conversion is performed according to the IEC/IEEE Standard for Binary
3337 | Floating-Point Arithmetic---which means in particular that the conversion
3338 | is rounded according to the current rounding mode. If `a' is a NaN,
3339 | the largest positive integer is returned. Otherwise, if the conversion
3340 | overflows, the largest integer with the same sign as `a' is returned.
3341 *----------------------------------------------------------------------------*/
3343 int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
3345 flag aSign;
3346 int32 aExp, shiftCount;
3347 bits64 aSig, aSigExtra;
3349 aSig = extractFloatx80Frac( a );
3350 aExp = extractFloatx80Exp( a );
3351 aSign = extractFloatx80Sign( a );
3352 shiftCount = 0x403E - aExp;
3353 if ( shiftCount <= 0 ) {
3354 if ( shiftCount ) {
3355 float_raise( float_flag_invalid STATUS_VAR);
3356 if ( ! aSign
3357 || ( ( aExp == 0x7FFF )
3358 && ( aSig != LIT64( 0x8000000000000000 ) ) )
3360 return LIT64( 0x7FFFFFFFFFFFFFFF );
3362 return (sbits64) LIT64( 0x8000000000000000 );
3364 aSigExtra = 0;
3366 else {
3367 shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
3369 return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
3373 /*----------------------------------------------------------------------------
3374 | Returns the result of converting the extended double-precision floating-
3375 | point value `a' to the 64-bit two's complement integer format. The
3376 | conversion is performed according to the IEC/IEEE Standard for Binary
3377 | Floating-Point Arithmetic, except that the conversion is always rounded
3378 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3379 | Otherwise, if the conversion overflows, the largest integer with the same
3380 | sign as `a' is returned.
3381 *----------------------------------------------------------------------------*/
3383 int64 floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM )
3385 flag aSign;
3386 int32 aExp, shiftCount;
3387 bits64 aSig;
3388 int64 z;
3390 aSig = extractFloatx80Frac( a );
3391 aExp = extractFloatx80Exp( a );
3392 aSign = extractFloatx80Sign( a );
3393 shiftCount = aExp - 0x403E;
3394 if ( 0 <= shiftCount ) {
3395 aSig &= LIT64( 0x7FFFFFFFFFFFFFFF );
3396 if ( ( a.high != 0xC03E ) || aSig ) {
3397 float_raise( float_flag_invalid STATUS_VAR);
3398 if ( ! aSign || ( ( aExp == 0x7FFF ) && aSig ) ) {
3399 return LIT64( 0x7FFFFFFFFFFFFFFF );
3402 return (sbits64) LIT64( 0x8000000000000000 );
3404 else if ( aExp < 0x3FFF ) {
3405 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3406 return 0;
3408 z = aSig>>( - shiftCount );
3409 if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
3410 STATUS(float_exception_flags) |= float_flag_inexact;
3412 if ( aSign ) z = - z;
3413 return z;
3417 /*----------------------------------------------------------------------------
3418 | Returns the result of converting the extended double-precision floating-
3419 | point value `a' to the single-precision floating-point format. The
3420 | conversion is performed according to the IEC/IEEE Standard for Binary
3421 | Floating-Point Arithmetic.
3422 *----------------------------------------------------------------------------*/
3424 float32 floatx80_to_float32( floatx80 a STATUS_PARAM )
3426 flag aSign;
3427 int32 aExp;
3428 bits64 aSig;
3430 aSig = extractFloatx80Frac( a );
3431 aExp = extractFloatx80Exp( a );
3432 aSign = extractFloatx80Sign( a );
3433 if ( aExp == 0x7FFF ) {
3434 if ( (bits64) ( aSig<<1 ) ) {
3435 return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) );
3437 return packFloat32( aSign, 0xFF, 0 );
3439 shift64RightJamming( aSig, 33, &aSig );
3440 if ( aExp || aSig ) aExp -= 0x3F81;
3441 return roundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
3445 /*----------------------------------------------------------------------------
3446 | Returns the result of converting the extended double-precision floating-
3447 | point value `a' to the double-precision floating-point format. The
3448 | conversion is performed according to the IEC/IEEE Standard for Binary
3449 | Floating-Point Arithmetic.
3450 *----------------------------------------------------------------------------*/
3452 float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
3454 flag aSign;
3455 int32 aExp;
3456 bits64 aSig, zSig;
3458 aSig = extractFloatx80Frac( a );
3459 aExp = extractFloatx80Exp( a );
3460 aSign = extractFloatx80Sign( a );
3461 if ( aExp == 0x7FFF ) {
3462 if ( (bits64) ( aSig<<1 ) ) {
3463 return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) );
3465 return packFloat64( aSign, 0x7FF, 0 );
3467 shift64RightJamming( aSig, 1, &zSig );
3468 if ( aExp || aSig ) aExp -= 0x3C01;
3469 return roundAndPackFloat64( aSign, aExp, zSig STATUS_VAR );
3473 #ifdef FLOAT128
3475 /*----------------------------------------------------------------------------
3476 | Returns the result of converting the extended double-precision floating-
3477 | point value `a' to the quadruple-precision floating-point format. The
3478 | conversion is performed according to the IEC/IEEE Standard for Binary
3479 | Floating-Point Arithmetic.
3480 *----------------------------------------------------------------------------*/
3482 float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
3484 flag aSign;
3485 int16 aExp;
3486 bits64 aSig, zSig0, zSig1;
3488 aSig = extractFloatx80Frac( a );
3489 aExp = extractFloatx80Exp( a );
3490 aSign = extractFloatx80Sign( a );
3491 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
3492 return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) );
3494 shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
3495 return packFloat128( aSign, aExp, zSig0, zSig1 );
3499 #endif
3501 /*----------------------------------------------------------------------------
3502 | Rounds the extended double-precision floating-point value `a' to an integer,
3503 | and returns the result as an extended quadruple-precision floating-point
3504 | value. The operation is performed according to the IEC/IEEE Standard for
3505 | Binary Floating-Point Arithmetic.
3506 *----------------------------------------------------------------------------*/
3508 floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM )
3510 flag aSign;
3511 int32 aExp;
3512 bits64 lastBitMask, roundBitsMask;
3513 int8 roundingMode;
3514 floatx80 z;
3516 aExp = extractFloatx80Exp( a );
3517 if ( 0x403E <= aExp ) {
3518 if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
3519 return propagateFloatx80NaN( a, a STATUS_VAR );
3521 return a;
3523 if ( aExp < 0x3FFF ) {
3524 if ( ( aExp == 0 )
3525 && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
3526 return a;
3528 STATUS(float_exception_flags) |= float_flag_inexact;
3529 aSign = extractFloatx80Sign( a );
3530 switch ( STATUS(float_rounding_mode) ) {
3531 case float_round_nearest_even:
3532 if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
3534 return
3535 packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3537 break;
3538 case float_round_down:
3539 return
3540 aSign ?
3541 packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
3542 : packFloatx80( 0, 0, 0 );
3543 case float_round_up:
3544 return
3545 aSign ? packFloatx80( 1, 0, 0 )
3546 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
3548 return packFloatx80( aSign, 0, 0 );
3550 lastBitMask = 1;
3551 lastBitMask <<= 0x403E - aExp;
3552 roundBitsMask = lastBitMask - 1;
3553 z = a;
3554 roundingMode = STATUS(float_rounding_mode);
3555 if ( roundingMode == float_round_nearest_even ) {
3556 z.low += lastBitMask>>1;
3557 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
3559 else if ( roundingMode != float_round_to_zero ) {
3560 if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
3561 z.low += roundBitsMask;
3564 z.low &= ~ roundBitsMask;
3565 if ( z.low == 0 ) {
3566 ++z.high;
3567 z.low = LIT64( 0x8000000000000000 );
3569 if ( z.low != a.low ) STATUS(float_exception_flags) |= float_flag_inexact;
3570 return z;
3574 /*----------------------------------------------------------------------------
3575 | Returns the result of adding the absolute values of the extended double-
3576 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
3577 | negated before being returned. `zSign' is ignored if the result is a NaN.
3578 | The addition is performed according to the IEC/IEEE Standard for Binary
3579 | Floating-Point Arithmetic.
3580 *----------------------------------------------------------------------------*/
3582 static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM)
3584 int32 aExp, bExp, zExp;
3585 bits64 aSig, bSig, zSig0, zSig1;
3586 int32 expDiff;
3588 aSig = extractFloatx80Frac( a );
3589 aExp = extractFloatx80Exp( a );
3590 bSig = extractFloatx80Frac( b );
3591 bExp = extractFloatx80Exp( b );
3592 expDiff = aExp - bExp;
3593 if ( 0 < expDiff ) {
3594 if ( aExp == 0x7FFF ) {
3595 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3596 return a;
3598 if ( bExp == 0 ) --expDiff;
3599 shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3600 zExp = aExp;
3602 else if ( expDiff < 0 ) {
3603 if ( bExp == 0x7FFF ) {
3604 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3605 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3607 if ( aExp == 0 ) ++expDiff;
3608 shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3609 zExp = bExp;
3611 else {
3612 if ( aExp == 0x7FFF ) {
3613 if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3614 return propagateFloatx80NaN( a, b STATUS_VAR );
3616 return a;
3618 zSig1 = 0;
3619 zSig0 = aSig + bSig;
3620 if ( aExp == 0 ) {
3621 normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
3622 goto roundAndPack;
3624 zExp = aExp;
3625 goto shiftRight1;
3627 zSig0 = aSig + bSig;
3628 if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
3629 shiftRight1:
3630 shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
3631 zSig0 |= LIT64( 0x8000000000000000 );
3632 ++zExp;
3633 roundAndPack:
3634 return
3635 roundAndPackFloatx80(
3636 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3640 /*----------------------------------------------------------------------------
3641 | Returns the result of subtracting the absolute values of the extended
3642 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
3643 | difference is negated before being returned. `zSign' is ignored if the
3644 | result is a NaN. The subtraction is performed according to the IEC/IEEE
3645 | Standard for Binary Floating-Point Arithmetic.
3646 *----------------------------------------------------------------------------*/
3648 static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM )
3650 int32 aExp, bExp, zExp;
3651 bits64 aSig, bSig, zSig0, zSig1;
3652 int32 expDiff;
3653 floatx80 z;
3655 aSig = extractFloatx80Frac( a );
3656 aExp = extractFloatx80Exp( a );
3657 bSig = extractFloatx80Frac( b );
3658 bExp = extractFloatx80Exp( b );
3659 expDiff = aExp - bExp;
3660 if ( 0 < expDiff ) goto aExpBigger;
3661 if ( expDiff < 0 ) goto bExpBigger;
3662 if ( aExp == 0x7FFF ) {
3663 if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3664 return propagateFloatx80NaN( a, b STATUS_VAR );
3666 float_raise( float_flag_invalid STATUS_VAR);
3667 z.low = floatx80_default_nan_low;
3668 z.high = floatx80_default_nan_high;
3669 return z;
3671 if ( aExp == 0 ) {
3672 aExp = 1;
3673 bExp = 1;
3675 zSig1 = 0;
3676 if ( bSig < aSig ) goto aBigger;
3677 if ( aSig < bSig ) goto bBigger;
3678 return packFloatx80( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
3679 bExpBigger:
3680 if ( bExp == 0x7FFF ) {
3681 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3682 return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
3684 if ( aExp == 0 ) ++expDiff;
3685 shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3686 bBigger:
3687 sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
3688 zExp = bExp;
3689 zSign ^= 1;
3690 goto normalizeRoundAndPack;
3691 aExpBigger:
3692 if ( aExp == 0x7FFF ) {
3693 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3694 return a;
3696 if ( bExp == 0 ) --expDiff;
3697 shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3698 aBigger:
3699 sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
3700 zExp = aExp;
3701 normalizeRoundAndPack:
3702 return
3703 normalizeRoundAndPackFloatx80(
3704 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3708 /*----------------------------------------------------------------------------
3709 | Returns the result of adding the extended double-precision floating-point
3710 | values `a' and `b'. The operation is performed according to the IEC/IEEE
3711 | Standard for Binary Floating-Point Arithmetic.
3712 *----------------------------------------------------------------------------*/
3714 floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM )
3716 flag aSign, bSign;
3718 aSign = extractFloatx80Sign( a );
3719 bSign = extractFloatx80Sign( b );
3720 if ( aSign == bSign ) {
3721 return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3723 else {
3724 return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3729 /*----------------------------------------------------------------------------
3730 | Returns the result of subtracting the extended double-precision floating-
3731 | point values `a' and `b'. The operation is performed according to the
3732 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3733 *----------------------------------------------------------------------------*/
3735 floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM )
3737 flag aSign, bSign;
3739 aSign = extractFloatx80Sign( a );
3740 bSign = extractFloatx80Sign( b );
3741 if ( aSign == bSign ) {
3742 return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3744 else {
3745 return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3750 /*----------------------------------------------------------------------------
3751 | Returns the result of multiplying the extended double-precision floating-
3752 | point values `a' and `b'. The operation is performed according to the
3753 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3754 *----------------------------------------------------------------------------*/
3756 floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
3758 flag aSign, bSign, zSign;
3759 int32 aExp, bExp, zExp;
3760 bits64 aSig, bSig, zSig0, zSig1;
3761 floatx80 z;
3763 aSig = extractFloatx80Frac( a );
3764 aExp = extractFloatx80Exp( a );
3765 aSign = extractFloatx80Sign( a );
3766 bSig = extractFloatx80Frac( b );
3767 bExp = extractFloatx80Exp( b );
3768 bSign = extractFloatx80Sign( b );
3769 zSign = aSign ^ bSign;
3770 if ( aExp == 0x7FFF ) {
3771 if ( (bits64) ( aSig<<1 )
3772 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3773 return propagateFloatx80NaN( a, b STATUS_VAR );
3775 if ( ( bExp | bSig ) == 0 ) goto invalid;
3776 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3778 if ( bExp == 0x7FFF ) {
3779 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3780 if ( ( aExp | aSig ) == 0 ) {
3781 invalid:
3782 float_raise( float_flag_invalid STATUS_VAR);
3783 z.low = floatx80_default_nan_low;
3784 z.high = floatx80_default_nan_high;
3785 return z;
3787 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3789 if ( aExp == 0 ) {
3790 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3791 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3793 if ( bExp == 0 ) {
3794 if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
3795 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3797 zExp = aExp + bExp - 0x3FFE;
3798 mul64To128( aSig, bSig, &zSig0, &zSig1 );
3799 if ( 0 < (sbits64) zSig0 ) {
3800 shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
3801 --zExp;
3803 return
3804 roundAndPackFloatx80(
3805 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3809 /*----------------------------------------------------------------------------
3810 | Returns the result of dividing the extended double-precision floating-point
3811 | value `a' by the corresponding value `b'. The operation is performed
3812 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3813 *----------------------------------------------------------------------------*/
3815 floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
3817 flag aSign, bSign, zSign;
3818 int32 aExp, bExp, zExp;
3819 bits64 aSig, bSig, zSig0, zSig1;
3820 bits64 rem0, rem1, rem2, term0, term1, term2;
3821 floatx80 z;
3823 aSig = extractFloatx80Frac( a );
3824 aExp = extractFloatx80Exp( a );
3825 aSign = extractFloatx80Sign( a );
3826 bSig = extractFloatx80Frac( b );
3827 bExp = extractFloatx80Exp( b );
3828 bSign = extractFloatx80Sign( b );
3829 zSign = aSign ^ bSign;
3830 if ( aExp == 0x7FFF ) {
3831 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3832 if ( bExp == 0x7FFF ) {
3833 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3834 goto invalid;
3836 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3838 if ( bExp == 0x7FFF ) {
3839 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3840 return packFloatx80( zSign, 0, 0 );
3842 if ( bExp == 0 ) {
3843 if ( bSig == 0 ) {
3844 if ( ( aExp | aSig ) == 0 ) {
3845 invalid:
3846 float_raise( float_flag_invalid STATUS_VAR);
3847 z.low = floatx80_default_nan_low;
3848 z.high = floatx80_default_nan_high;
3849 return z;
3851 float_raise( float_flag_divbyzero STATUS_VAR);
3852 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3854 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3856 if ( aExp == 0 ) {
3857 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3858 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3860 zExp = aExp - bExp + 0x3FFE;
3861 rem1 = 0;
3862 if ( bSig <= aSig ) {
3863 shift128Right( aSig, 0, 1, &aSig, &rem1 );
3864 ++zExp;
3866 zSig0 = estimateDiv128To64( aSig, rem1, bSig );
3867 mul64To128( bSig, zSig0, &term0, &term1 );
3868 sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
3869 while ( (sbits64) rem0 < 0 ) {
3870 --zSig0;
3871 add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
3873 zSig1 = estimateDiv128To64( rem1, 0, bSig );
3874 if ( (bits64) ( zSig1<<1 ) <= 8 ) {
3875 mul64To128( bSig, zSig1, &term1, &term2 );
3876 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3877 while ( (sbits64) rem1 < 0 ) {
3878 --zSig1;
3879 add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
3881 zSig1 |= ( ( rem1 | rem2 ) != 0 );
3883 return
3884 roundAndPackFloatx80(
3885 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3889 /*----------------------------------------------------------------------------
3890 | Returns the remainder of the extended double-precision floating-point value
3891 | `a' with respect to the corresponding value `b'. The operation is performed
3892 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3893 *----------------------------------------------------------------------------*/
3895 floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
3897 flag aSign, bSign, zSign;
3898 int32 aExp, bExp, expDiff;
3899 bits64 aSig0, aSig1, bSig;
3900 bits64 q, term0, term1, alternateASig0, alternateASig1;
3901 floatx80 z;
3903 aSig0 = extractFloatx80Frac( a );
3904 aExp = extractFloatx80Exp( a );
3905 aSign = extractFloatx80Sign( a );
3906 bSig = extractFloatx80Frac( b );
3907 bExp = extractFloatx80Exp( b );
3908 bSign = extractFloatx80Sign( b );
3909 if ( aExp == 0x7FFF ) {
3910 if ( (bits64) ( aSig0<<1 )
3911 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3912 return propagateFloatx80NaN( a, b STATUS_VAR );
3914 goto invalid;
3916 if ( bExp == 0x7FFF ) {
3917 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3918 return a;
3920 if ( bExp == 0 ) {
3921 if ( bSig == 0 ) {
3922 invalid:
3923 float_raise( float_flag_invalid STATUS_VAR);
3924 z.low = floatx80_default_nan_low;
3925 z.high = floatx80_default_nan_high;
3926 return z;
3928 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3930 if ( aExp == 0 ) {
3931 if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
3932 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3934 bSig |= LIT64( 0x8000000000000000 );
3935 zSign = aSign;
3936 expDiff = aExp - bExp;
3937 aSig1 = 0;
3938 if ( expDiff < 0 ) {
3939 if ( expDiff < -1 ) return a;
3940 shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
3941 expDiff = 0;
3943 q = ( bSig <= aSig0 );
3944 if ( q ) aSig0 -= bSig;
3945 expDiff -= 64;
3946 while ( 0 < expDiff ) {
3947 q = estimateDiv128To64( aSig0, aSig1, bSig );
3948 q = ( 2 < q ) ? q - 2 : 0;
3949 mul64To128( bSig, q, &term0, &term1 );
3950 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3951 shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
3952 expDiff -= 62;
3954 expDiff += 64;
3955 if ( 0 < expDiff ) {
3956 q = estimateDiv128To64( aSig0, aSig1, bSig );
3957 q = ( 2 < q ) ? q - 2 : 0;
3958 q >>= 64 - expDiff;
3959 mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
3960 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3961 shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
3962 while ( le128( term0, term1, aSig0, aSig1 ) ) {
3963 ++q;
3964 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3967 else {
3968 term1 = 0;
3969 term0 = bSig;
3971 sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
3972 if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
3973 || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
3974 && ( q & 1 ) )
3976 aSig0 = alternateASig0;
3977 aSig1 = alternateASig1;
3978 zSign = ! zSign;
3980 return
3981 normalizeRoundAndPackFloatx80(
3982 80, zSign, bExp + expDiff, aSig0, aSig1 STATUS_VAR );
3986 /*----------------------------------------------------------------------------
3987 | Returns the square root of the extended double-precision floating-point
3988 | value `a'. The operation is performed according to the IEC/IEEE Standard
3989 | for Binary Floating-Point Arithmetic.
3990 *----------------------------------------------------------------------------*/
3992 floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
3994 flag aSign;
3995 int32 aExp, zExp;
3996 bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
3997 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
3998 floatx80 z;
4000 aSig0 = extractFloatx80Frac( a );
4001 aExp = extractFloatx80Exp( a );
4002 aSign = extractFloatx80Sign( a );
4003 if ( aExp == 0x7FFF ) {
4004 if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a STATUS_VAR );
4005 if ( ! aSign ) return a;
4006 goto invalid;
4008 if ( aSign ) {
4009 if ( ( aExp | aSig0 ) == 0 ) return a;
4010 invalid:
4011 float_raise( float_flag_invalid STATUS_VAR);
4012 z.low = floatx80_default_nan_low;
4013 z.high = floatx80_default_nan_high;
4014 return z;
4016 if ( aExp == 0 ) {
4017 if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
4018 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
4020 zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
4021 zSig0 = estimateSqrt32( aExp, aSig0>>32 );
4022 shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
4023 zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
4024 doubleZSig0 = zSig0<<1;
4025 mul64To128( zSig0, zSig0, &term0, &term1 );
4026 sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
4027 while ( (sbits64) rem0 < 0 ) {
4028 --zSig0;
4029 doubleZSig0 -= 2;
4030 add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
4032 zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
4033 if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
4034 if ( zSig1 == 0 ) zSig1 = 1;
4035 mul64To128( doubleZSig0, zSig1, &term1, &term2 );
4036 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
4037 mul64To128( zSig1, zSig1, &term2, &term3 );
4038 sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
4039 while ( (sbits64) rem1 < 0 ) {
4040 --zSig1;
4041 shortShift128Left( 0, zSig1, 1, &term2, &term3 );
4042 term3 |= 1;
4043 term2 |= doubleZSig0;
4044 add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
4046 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
4048 shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
4049 zSig0 |= doubleZSig0;
4050 return
4051 roundAndPackFloatx80(
4052 STATUS(floatx80_rounding_precision), 0, zExp, zSig0, zSig1 STATUS_VAR );
4056 /*----------------------------------------------------------------------------
4057 | Returns 1 if the extended double-precision floating-point value `a' is
4058 | equal to the corresponding value `b', and 0 otherwise. The comparison is
4059 | performed according to the IEC/IEEE Standard for Binary Floating-Point
4060 | Arithmetic.
4061 *----------------------------------------------------------------------------*/
4063 int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
4066 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4067 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4068 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4069 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4071 if ( floatx80_is_signaling_nan( a )
4072 || floatx80_is_signaling_nan( b ) ) {
4073 float_raise( float_flag_invalid STATUS_VAR);
4075 return 0;
4077 return
4078 ( a.low == b.low )
4079 && ( ( a.high == b.high )
4080 || ( ( a.low == 0 )
4081 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
4086 /*----------------------------------------------------------------------------
4087 | Returns 1 if the extended double-precision floating-point value `a' is
4088 | less than or equal to the corresponding value `b', and 0 otherwise. The
4089 | comparison is performed according to the IEC/IEEE Standard for Binary
4090 | Floating-Point Arithmetic.
4091 *----------------------------------------------------------------------------*/
4093 int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
4095 flag aSign, bSign;
4097 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4098 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4099 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4100 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4102 float_raise( float_flag_invalid STATUS_VAR);
4103 return 0;
4105 aSign = extractFloatx80Sign( a );
4106 bSign = extractFloatx80Sign( b );
4107 if ( aSign != bSign ) {
4108 return
4109 aSign
4110 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4111 == 0 );
4113 return
4114 aSign ? le128( b.high, b.low, a.high, a.low )
4115 : le128( a.high, a.low, b.high, b.low );
4119 /*----------------------------------------------------------------------------
4120 | Returns 1 if the extended double-precision floating-point value `a' is
4121 | less than the corresponding value `b', and 0 otherwise. The comparison
4122 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4123 | Arithmetic.
4124 *----------------------------------------------------------------------------*/
4126 int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
4128 flag aSign, bSign;
4130 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4131 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4132 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4133 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4135 float_raise( float_flag_invalid STATUS_VAR);
4136 return 0;
4138 aSign = extractFloatx80Sign( a );
4139 bSign = extractFloatx80Sign( b );
4140 if ( aSign != bSign ) {
4141 return
4142 aSign
4143 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4144 != 0 );
4146 return
4147 aSign ? lt128( b.high, b.low, a.high, a.low )
4148 : lt128( a.high, a.low, b.high, b.low );
4152 /*----------------------------------------------------------------------------
4153 | Returns 1 if the extended double-precision floating-point value `a' is equal
4154 | to the corresponding value `b', and 0 otherwise. The invalid exception is
4155 | raised if either operand is a NaN. Otherwise, the comparison is performed
4156 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4157 *----------------------------------------------------------------------------*/
4159 int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
4162 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4163 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4164 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4165 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4167 float_raise( float_flag_invalid STATUS_VAR);
4168 return 0;
4170 return
4171 ( a.low == b.low )
4172 && ( ( a.high == b.high )
4173 || ( ( a.low == 0 )
4174 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
4179 /*----------------------------------------------------------------------------
4180 | Returns 1 if the extended double-precision floating-point value `a' is less
4181 | than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
4182 | do not cause an exception. Otherwise, the comparison is performed according
4183 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4184 *----------------------------------------------------------------------------*/
4186 int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4188 flag aSign, bSign;
4190 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4191 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4192 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4193 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4195 if ( floatx80_is_signaling_nan( a )
4196 || floatx80_is_signaling_nan( b ) ) {
4197 float_raise( float_flag_invalid STATUS_VAR);
4199 return 0;
4201 aSign = extractFloatx80Sign( a );
4202 bSign = extractFloatx80Sign( b );
4203 if ( aSign != bSign ) {
4204 return
4205 aSign
4206 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4207 == 0 );
4209 return
4210 aSign ? le128( b.high, b.low, a.high, a.low )
4211 : le128( a.high, a.low, b.high, b.low );
4215 /*----------------------------------------------------------------------------
4216 | Returns 1 if the extended double-precision floating-point value `a' is less
4217 | than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
4218 | an exception. Otherwise, the comparison is performed according to the
4219 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4220 *----------------------------------------------------------------------------*/
4222 int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4224 flag aSign, bSign;
4226 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4227 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4228 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4229 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4231 if ( floatx80_is_signaling_nan( a )
4232 || floatx80_is_signaling_nan( b ) ) {
4233 float_raise( float_flag_invalid STATUS_VAR);
4235 return 0;
4237 aSign = extractFloatx80Sign( a );
4238 bSign = extractFloatx80Sign( b );
4239 if ( aSign != bSign ) {
4240 return
4241 aSign
4242 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4243 != 0 );
4245 return
4246 aSign ? lt128( b.high, b.low, a.high, a.low )
4247 : lt128( a.high, a.low, b.high, b.low );
4251 #endif
4253 #ifdef FLOAT128
4255 /*----------------------------------------------------------------------------
4256 | Returns the result of converting the quadruple-precision floating-point
4257 | value `a' to the 32-bit two's complement integer format. The conversion
4258 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4259 | Arithmetic---which means in particular that the conversion is rounded
4260 | according to the current rounding mode. If `a' is a NaN, the largest
4261 | positive integer is returned. Otherwise, if the conversion overflows, the
4262 | largest integer with the same sign as `a' is returned.
4263 *----------------------------------------------------------------------------*/
4265 int32 float128_to_int32( float128 a STATUS_PARAM )
4267 flag aSign;
4268 int32 aExp, shiftCount;
4269 bits64 aSig0, aSig1;
4271 aSig1 = extractFloat128Frac1( a );
4272 aSig0 = extractFloat128Frac0( a );
4273 aExp = extractFloat128Exp( a );
4274 aSign = extractFloat128Sign( a );
4275 if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
4276 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4277 aSig0 |= ( aSig1 != 0 );
4278 shiftCount = 0x4028 - aExp;
4279 if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
4280 return roundAndPackInt32( aSign, aSig0 STATUS_VAR );
4284 /*----------------------------------------------------------------------------
4285 | Returns the result of converting the quadruple-precision floating-point
4286 | value `a' to the 32-bit two's complement integer format. The conversion
4287 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4288 | Arithmetic, except that the conversion is always rounded toward zero. If
4289 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
4290 | conversion overflows, the largest integer with the same sign as `a' is
4291 | returned.
4292 *----------------------------------------------------------------------------*/
4294 int32 float128_to_int32_round_to_zero( float128 a STATUS_PARAM )
4296 flag aSign;
4297 int32 aExp, shiftCount;
4298 bits64 aSig0, aSig1, savedASig;
4299 int32 z;
4301 aSig1 = extractFloat128Frac1( a );
4302 aSig0 = extractFloat128Frac0( a );
4303 aExp = extractFloat128Exp( a );
4304 aSign = extractFloat128Sign( a );
4305 aSig0 |= ( aSig1 != 0 );
4306 if ( 0x401E < aExp ) {
4307 if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
4308 goto invalid;
4310 else if ( aExp < 0x3FFF ) {
4311 if ( aExp || aSig0 ) STATUS(float_exception_flags) |= float_flag_inexact;
4312 return 0;
4314 aSig0 |= LIT64( 0x0001000000000000 );
4315 shiftCount = 0x402F - aExp;
4316 savedASig = aSig0;
4317 aSig0 >>= shiftCount;
4318 z = aSig0;
4319 if ( aSign ) z = - z;
4320 if ( ( z < 0 ) ^ aSign ) {
4321 invalid:
4322 float_raise( float_flag_invalid STATUS_VAR);
4323 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
4325 if ( ( aSig0<<shiftCount ) != savedASig ) {
4326 STATUS(float_exception_flags) |= float_flag_inexact;
4328 return z;
4332 /*----------------------------------------------------------------------------
4333 | Returns the result of converting the quadruple-precision floating-point
4334 | value `a' to the 64-bit two's complement integer format. The conversion
4335 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4336 | Arithmetic---which means in particular that the conversion is rounded
4337 | according to the current rounding mode. If `a' is a NaN, the largest
4338 | positive integer is returned. Otherwise, if the conversion overflows, the
4339 | largest integer with the same sign as `a' is returned.
4340 *----------------------------------------------------------------------------*/
4342 int64 float128_to_int64( float128 a STATUS_PARAM )
4344 flag aSign;
4345 int32 aExp, shiftCount;
4346 bits64 aSig0, aSig1;
4348 aSig1 = extractFloat128Frac1( a );
4349 aSig0 = extractFloat128Frac0( a );
4350 aExp = extractFloat128Exp( a );
4351 aSign = extractFloat128Sign( a );
4352 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4353 shiftCount = 0x402F - aExp;
4354 if ( shiftCount <= 0 ) {
4355 if ( 0x403E < aExp ) {
4356 float_raise( float_flag_invalid STATUS_VAR);
4357 if ( ! aSign
4358 || ( ( aExp == 0x7FFF )
4359 && ( aSig1 || ( aSig0 != LIT64( 0x0001000000000000 ) ) )
4362 return LIT64( 0x7FFFFFFFFFFFFFFF );
4364 return (sbits64) LIT64( 0x8000000000000000 );
4366 shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
4368 else {
4369 shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
4371 return roundAndPackInt64( aSign, aSig0, aSig1 STATUS_VAR );
4375 /*----------------------------------------------------------------------------
4376 | Returns the result of converting the quadruple-precision floating-point
4377 | value `a' to the 64-bit two's complement integer format. The conversion
4378 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4379 | Arithmetic, except that the conversion is always rounded toward zero.
4380 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
4381 | the conversion overflows, the largest integer with the same sign as `a' is
4382 | returned.
4383 *----------------------------------------------------------------------------*/
4385 int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
4387 flag aSign;
4388 int32 aExp, shiftCount;
4389 bits64 aSig0, aSig1;
4390 int64 z;
4392 aSig1 = extractFloat128Frac1( a );
4393 aSig0 = extractFloat128Frac0( a );
4394 aExp = extractFloat128Exp( a );
4395 aSign = extractFloat128Sign( a );
4396 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4397 shiftCount = aExp - 0x402F;
4398 if ( 0 < shiftCount ) {
4399 if ( 0x403E <= aExp ) {
4400 aSig0 &= LIT64( 0x0000FFFFFFFFFFFF );
4401 if ( ( a.high == LIT64( 0xC03E000000000000 ) )
4402 && ( aSig1 < LIT64( 0x0002000000000000 ) ) ) {
4403 if ( aSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
4405 else {
4406 float_raise( float_flag_invalid STATUS_VAR);
4407 if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
4408 return LIT64( 0x7FFFFFFFFFFFFFFF );
4411 return (sbits64) LIT64( 0x8000000000000000 );
4413 z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
4414 if ( (bits64) ( aSig1<<shiftCount ) ) {
4415 STATUS(float_exception_flags) |= float_flag_inexact;
4418 else {
4419 if ( aExp < 0x3FFF ) {
4420 if ( aExp | aSig0 | aSig1 ) {
4421 STATUS(float_exception_flags) |= float_flag_inexact;
4423 return 0;
4425 z = aSig0>>( - shiftCount );
4426 if ( aSig1
4427 || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
4428 STATUS(float_exception_flags) |= float_flag_inexact;
4431 if ( aSign ) z = - z;
4432 return z;
4436 /*----------------------------------------------------------------------------
4437 | Returns the result of converting the quadruple-precision floating-point
4438 | value `a' to the single-precision floating-point format. The conversion
4439 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4440 | Arithmetic.
4441 *----------------------------------------------------------------------------*/
4443 float32 float128_to_float32( float128 a STATUS_PARAM )
4445 flag aSign;
4446 int32 aExp;
4447 bits64 aSig0, aSig1;
4448 bits32 zSig;
4450 aSig1 = extractFloat128Frac1( a );
4451 aSig0 = extractFloat128Frac0( a );
4452 aExp = extractFloat128Exp( a );
4453 aSign = extractFloat128Sign( a );
4454 if ( aExp == 0x7FFF ) {
4455 if ( aSig0 | aSig1 ) {
4456 return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR ) );
4458 return packFloat32( aSign, 0xFF, 0 );
4460 aSig0 |= ( aSig1 != 0 );
4461 shift64RightJamming( aSig0, 18, &aSig0 );
4462 zSig = aSig0;
4463 if ( aExp || zSig ) {
4464 zSig |= 0x40000000;
4465 aExp -= 0x3F81;
4467 return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
4471 /*----------------------------------------------------------------------------
4472 | Returns the result of converting the quadruple-precision floating-point
4473 | value `a' to the double-precision floating-point format. The conversion
4474 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4475 | Arithmetic.
4476 *----------------------------------------------------------------------------*/
4478 float64 float128_to_float64( float128 a STATUS_PARAM )
4480 flag aSign;
4481 int32 aExp;
4482 bits64 aSig0, aSig1;
4484 aSig1 = extractFloat128Frac1( a );
4485 aSig0 = extractFloat128Frac0( a );
4486 aExp = extractFloat128Exp( a );
4487 aSign = extractFloat128Sign( a );
4488 if ( aExp == 0x7FFF ) {
4489 if ( aSig0 | aSig1 ) {
4490 return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR ) );
4492 return packFloat64( aSign, 0x7FF, 0 );
4494 shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4495 aSig0 |= ( aSig1 != 0 );
4496 if ( aExp || aSig0 ) {
4497 aSig0 |= LIT64( 0x4000000000000000 );
4498 aExp -= 0x3C01;
4500 return roundAndPackFloat64( aSign, aExp, aSig0 STATUS_VAR );
4504 #ifdef FLOATX80
4506 /*----------------------------------------------------------------------------
4507 | Returns the result of converting the quadruple-precision floating-point
4508 | value `a' to the extended double-precision floating-point format. The
4509 | conversion is performed according to the IEC/IEEE Standard for Binary
4510 | Floating-Point Arithmetic.
4511 *----------------------------------------------------------------------------*/
4513 floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
4515 flag aSign;
4516 int32 aExp;
4517 bits64 aSig0, aSig1;
4519 aSig1 = extractFloat128Frac1( a );
4520 aSig0 = extractFloat128Frac0( a );
4521 aExp = extractFloat128Exp( a );
4522 aSign = extractFloat128Sign( a );
4523 if ( aExp == 0x7FFF ) {
4524 if ( aSig0 | aSig1 ) {
4525 return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR ) );
4527 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
4529 if ( aExp == 0 ) {
4530 if ( ( aSig0 | aSig1 ) == 0 ) return packFloatx80( aSign, 0, 0 );
4531 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4533 else {
4534 aSig0 |= LIT64( 0x0001000000000000 );
4536 shortShift128Left( aSig0, aSig1, 15, &aSig0, &aSig1 );
4537 return roundAndPackFloatx80( 80, aSign, aExp, aSig0, aSig1 STATUS_VAR );
4541 #endif
4543 /*----------------------------------------------------------------------------
4544 | Rounds the quadruple-precision floating-point value `a' to an integer, and
4545 | returns the result as a quadruple-precision floating-point value. The
4546 | operation is performed according to the IEC/IEEE Standard for Binary
4547 | Floating-Point Arithmetic.
4548 *----------------------------------------------------------------------------*/
4550 float128 float128_round_to_int( float128 a STATUS_PARAM )
4552 flag aSign;
4553 int32 aExp;
4554 bits64 lastBitMask, roundBitsMask;
4555 int8 roundingMode;
4556 float128 z;
4558 aExp = extractFloat128Exp( a );
4559 if ( 0x402F <= aExp ) {
4560 if ( 0x406F <= aExp ) {
4561 if ( ( aExp == 0x7FFF )
4562 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
4564 return propagateFloat128NaN( a, a STATUS_VAR );
4566 return a;
4568 lastBitMask = 1;
4569 lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
4570 roundBitsMask = lastBitMask - 1;
4571 z = a;
4572 roundingMode = STATUS(float_rounding_mode);
4573 if ( roundingMode == float_round_nearest_even ) {
4574 if ( lastBitMask ) {
4575 add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
4576 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
4578 else {
4579 if ( (sbits64) z.low < 0 ) {
4580 ++z.high;
4581 if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
4585 else if ( roundingMode != float_round_to_zero ) {
4586 if ( extractFloat128Sign( z )
4587 ^ ( roundingMode == float_round_up ) ) {
4588 add128( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
4591 z.low &= ~ roundBitsMask;
4593 else {
4594 if ( aExp < 0x3FFF ) {
4595 if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
4596 STATUS(float_exception_flags) |= float_flag_inexact;
4597 aSign = extractFloat128Sign( a );
4598 switch ( STATUS(float_rounding_mode) ) {
4599 case float_round_nearest_even:
4600 if ( ( aExp == 0x3FFE )
4601 && ( extractFloat128Frac0( a )
4602 | extractFloat128Frac1( a ) )
4604 return packFloat128( aSign, 0x3FFF, 0, 0 );
4606 break;
4607 case float_round_down:
4608 return
4609 aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
4610 : packFloat128( 0, 0, 0, 0 );
4611 case float_round_up:
4612 return
4613 aSign ? packFloat128( 1, 0, 0, 0 )
4614 : packFloat128( 0, 0x3FFF, 0, 0 );
4616 return packFloat128( aSign, 0, 0, 0 );
4618 lastBitMask = 1;
4619 lastBitMask <<= 0x402F - aExp;
4620 roundBitsMask = lastBitMask - 1;
4621 z.low = 0;
4622 z.high = a.high;
4623 roundingMode = STATUS(float_rounding_mode);
4624 if ( roundingMode == float_round_nearest_even ) {
4625 z.high += lastBitMask>>1;
4626 if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
4627 z.high &= ~ lastBitMask;
4630 else if ( roundingMode != float_round_to_zero ) {
4631 if ( extractFloat128Sign( z )
4632 ^ ( roundingMode == float_round_up ) ) {
4633 z.high |= ( a.low != 0 );
4634 z.high += roundBitsMask;
4637 z.high &= ~ roundBitsMask;
4639 if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
4640 STATUS(float_exception_flags) |= float_flag_inexact;
4642 return z;
4646 /*----------------------------------------------------------------------------
4647 | Returns the result of adding the absolute values of the quadruple-precision
4648 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
4649 | before being returned. `zSign' is ignored if the result is a NaN.
4650 | The addition is performed according to the IEC/IEEE Standard for Binary
4651 | Floating-Point Arithmetic.
4652 *----------------------------------------------------------------------------*/
4654 static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4656 int32 aExp, bExp, zExp;
4657 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4658 int32 expDiff;
4660 aSig1 = extractFloat128Frac1( a );
4661 aSig0 = extractFloat128Frac0( a );
4662 aExp = extractFloat128Exp( a );
4663 bSig1 = extractFloat128Frac1( b );
4664 bSig0 = extractFloat128Frac0( b );
4665 bExp = extractFloat128Exp( b );
4666 expDiff = aExp - bExp;
4667 if ( 0 < expDiff ) {
4668 if ( aExp == 0x7FFF ) {
4669 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4670 return a;
4672 if ( bExp == 0 ) {
4673 --expDiff;
4675 else {
4676 bSig0 |= LIT64( 0x0001000000000000 );
4678 shift128ExtraRightJamming(
4679 bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
4680 zExp = aExp;
4682 else if ( expDiff < 0 ) {
4683 if ( bExp == 0x7FFF ) {
4684 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4685 return packFloat128( zSign, 0x7FFF, 0, 0 );
4687 if ( aExp == 0 ) {
4688 ++expDiff;
4690 else {
4691 aSig0 |= LIT64( 0x0001000000000000 );
4693 shift128ExtraRightJamming(
4694 aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
4695 zExp = bExp;
4697 else {
4698 if ( aExp == 0x7FFF ) {
4699 if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4700 return propagateFloat128NaN( a, b STATUS_VAR );
4702 return a;
4704 add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4705 if ( aExp == 0 ) {
4706 if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
4707 return packFloat128( zSign, 0, zSig0, zSig1 );
4709 zSig2 = 0;
4710 zSig0 |= LIT64( 0x0002000000000000 );
4711 zExp = aExp;
4712 goto shiftRight1;
4714 aSig0 |= LIT64( 0x0001000000000000 );
4715 add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4716 --zExp;
4717 if ( zSig0 < LIT64( 0x0002000000000000 ) ) goto roundAndPack;
4718 ++zExp;
4719 shiftRight1:
4720 shift128ExtraRightJamming(
4721 zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4722 roundAndPack:
4723 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4727 /*----------------------------------------------------------------------------
4728 | Returns the result of subtracting the absolute values of the quadruple-
4729 | precision floating-point values `a' and `b'. If `zSign' is 1, the
4730 | difference is negated before being returned. `zSign' is ignored if the
4731 | result is a NaN. The subtraction is performed according to the IEC/IEEE
4732 | Standard for Binary Floating-Point Arithmetic.
4733 *----------------------------------------------------------------------------*/
4735 static float128 subFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4737 int32 aExp, bExp, zExp;
4738 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
4739 int32 expDiff;
4740 float128 z;
4742 aSig1 = extractFloat128Frac1( a );
4743 aSig0 = extractFloat128Frac0( a );
4744 aExp = extractFloat128Exp( a );
4745 bSig1 = extractFloat128Frac1( b );
4746 bSig0 = extractFloat128Frac0( b );
4747 bExp = extractFloat128Exp( b );
4748 expDiff = aExp - bExp;
4749 shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4750 shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
4751 if ( 0 < expDiff ) goto aExpBigger;
4752 if ( expDiff < 0 ) goto bExpBigger;
4753 if ( aExp == 0x7FFF ) {
4754 if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4755 return propagateFloat128NaN( a, b STATUS_VAR );
4757 float_raise( float_flag_invalid STATUS_VAR);
4758 z.low = float128_default_nan_low;
4759 z.high = float128_default_nan_high;
4760 return z;
4762 if ( aExp == 0 ) {
4763 aExp = 1;
4764 bExp = 1;
4766 if ( bSig0 < aSig0 ) goto aBigger;
4767 if ( aSig0 < bSig0 ) goto bBigger;
4768 if ( bSig1 < aSig1 ) goto aBigger;
4769 if ( aSig1 < bSig1 ) goto bBigger;
4770 return packFloat128( STATUS(float_rounding_mode) == float_round_down, 0, 0, 0 );
4771 bExpBigger:
4772 if ( bExp == 0x7FFF ) {
4773 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4774 return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
4776 if ( aExp == 0 ) {
4777 ++expDiff;
4779 else {
4780 aSig0 |= LIT64( 0x4000000000000000 );
4782 shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
4783 bSig0 |= LIT64( 0x4000000000000000 );
4784 bBigger:
4785 sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
4786 zExp = bExp;
4787 zSign ^= 1;
4788 goto normalizeRoundAndPack;
4789 aExpBigger:
4790 if ( aExp == 0x7FFF ) {
4791 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4792 return a;
4794 if ( bExp == 0 ) {
4795 --expDiff;
4797 else {
4798 bSig0 |= LIT64( 0x4000000000000000 );
4800 shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
4801 aSig0 |= LIT64( 0x4000000000000000 );
4802 aBigger:
4803 sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4804 zExp = aExp;
4805 normalizeRoundAndPack:
4806 --zExp;
4807 return normalizeRoundAndPackFloat128( zSign, zExp - 14, zSig0, zSig1 STATUS_VAR );
4811 /*----------------------------------------------------------------------------
4812 | Returns the result of adding the quadruple-precision floating-point values
4813 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
4814 | for Binary Floating-Point Arithmetic.
4815 *----------------------------------------------------------------------------*/
4817 float128 float128_add( float128 a, float128 b STATUS_PARAM )
4819 flag aSign, bSign;
4821 aSign = extractFloat128Sign( a );
4822 bSign = extractFloat128Sign( b );
4823 if ( aSign == bSign ) {
4824 return addFloat128Sigs( a, b, aSign STATUS_VAR );
4826 else {
4827 return subFloat128Sigs( a, b, aSign STATUS_VAR );
4832 /*----------------------------------------------------------------------------
4833 | Returns the result of subtracting the quadruple-precision floating-point
4834 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4835 | Standard for Binary Floating-Point Arithmetic.
4836 *----------------------------------------------------------------------------*/
4838 float128 float128_sub( float128 a, float128 b STATUS_PARAM )
4840 flag aSign, bSign;
4842 aSign = extractFloat128Sign( a );
4843 bSign = extractFloat128Sign( b );
4844 if ( aSign == bSign ) {
4845 return subFloat128Sigs( a, b, aSign STATUS_VAR );
4847 else {
4848 return addFloat128Sigs( a, b, aSign STATUS_VAR );
4853 /*----------------------------------------------------------------------------
4854 | Returns the result of multiplying the quadruple-precision floating-point
4855 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4856 | Standard for Binary Floating-Point Arithmetic.
4857 *----------------------------------------------------------------------------*/
4859 float128 float128_mul( float128 a, float128 b STATUS_PARAM )
4861 flag aSign, bSign, zSign;
4862 int32 aExp, bExp, zExp;
4863 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
4864 float128 z;
4866 aSig1 = extractFloat128Frac1( a );
4867 aSig0 = extractFloat128Frac0( a );
4868 aExp = extractFloat128Exp( a );
4869 aSign = extractFloat128Sign( a );
4870 bSig1 = extractFloat128Frac1( b );
4871 bSig0 = extractFloat128Frac0( b );
4872 bExp = extractFloat128Exp( b );
4873 bSign = extractFloat128Sign( b );
4874 zSign = aSign ^ bSign;
4875 if ( aExp == 0x7FFF ) {
4876 if ( ( aSig0 | aSig1 )
4877 || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
4878 return propagateFloat128NaN( a, b STATUS_VAR );
4880 if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
4881 return packFloat128( zSign, 0x7FFF, 0, 0 );
4883 if ( bExp == 0x7FFF ) {
4884 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4885 if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4886 invalid:
4887 float_raise( float_flag_invalid STATUS_VAR);
4888 z.low = float128_default_nan_low;
4889 z.high = float128_default_nan_high;
4890 return z;
4892 return packFloat128( zSign, 0x7FFF, 0, 0 );
4894 if ( aExp == 0 ) {
4895 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4896 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4898 if ( bExp == 0 ) {
4899 if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4900 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4902 zExp = aExp + bExp - 0x4000;
4903 aSig0 |= LIT64( 0x0001000000000000 );
4904 shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
4905 mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
4906 add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
4907 zSig2 |= ( zSig3 != 0 );
4908 if ( LIT64( 0x0002000000000000 ) <= zSig0 ) {
4909 shift128ExtraRightJamming(
4910 zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4911 ++zExp;
4913 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4917 /*----------------------------------------------------------------------------
4918 | Returns the result of dividing the quadruple-precision floating-point value
4919 | `a' by the corresponding value `b'. The operation is performed according to
4920 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4921 *----------------------------------------------------------------------------*/
4923 float128 float128_div( float128 a, float128 b STATUS_PARAM )
4925 flag aSign, bSign, zSign;
4926 int32 aExp, bExp, zExp;
4927 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4928 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
4929 float128 z;
4931 aSig1 = extractFloat128Frac1( a );
4932 aSig0 = extractFloat128Frac0( a );
4933 aExp = extractFloat128Exp( a );
4934 aSign = extractFloat128Sign( a );
4935 bSig1 = extractFloat128Frac1( b );
4936 bSig0 = extractFloat128Frac0( b );
4937 bExp = extractFloat128Exp( b );
4938 bSign = extractFloat128Sign( b );
4939 zSign = aSign ^ bSign;
4940 if ( aExp == 0x7FFF ) {
4941 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4942 if ( bExp == 0x7FFF ) {
4943 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4944 goto invalid;
4946 return packFloat128( zSign, 0x7FFF, 0, 0 );
4948 if ( bExp == 0x7FFF ) {
4949 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4950 return packFloat128( zSign, 0, 0, 0 );
4952 if ( bExp == 0 ) {
4953 if ( ( bSig0 | bSig1 ) == 0 ) {
4954 if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4955 invalid:
4956 float_raise( float_flag_invalid STATUS_VAR);
4957 z.low = float128_default_nan_low;
4958 z.high = float128_default_nan_high;
4959 return z;
4961 float_raise( float_flag_divbyzero STATUS_VAR);
4962 return packFloat128( zSign, 0x7FFF, 0, 0 );
4964 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4966 if ( aExp == 0 ) {
4967 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4968 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4970 zExp = aExp - bExp + 0x3FFD;
4971 shortShift128Left(
4972 aSig0 | LIT64( 0x0001000000000000 ), aSig1, 15, &aSig0, &aSig1 );
4973 shortShift128Left(
4974 bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
4975 if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
4976 shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
4977 ++zExp;
4979 zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
4980 mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
4981 sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
4982 while ( (sbits64) rem0 < 0 ) {
4983 --zSig0;
4984 add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
4986 zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
4987 if ( ( zSig1 & 0x3FFF ) <= 4 ) {
4988 mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
4989 sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
4990 while ( (sbits64) rem1 < 0 ) {
4991 --zSig1;
4992 add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
4994 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
4996 shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
4997 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
5001 /*----------------------------------------------------------------------------
5002 | Returns the remainder of the quadruple-precision floating-point value `a'
5003 | with respect to the corresponding value `b'. The operation is performed
5004 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5005 *----------------------------------------------------------------------------*/
5007 float128 float128_rem( float128 a, float128 b STATUS_PARAM )
5009 flag aSign, bSign, zSign;
5010 int32 aExp, bExp, expDiff;
5011 bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
5012 bits64 allZero, alternateASig0, alternateASig1, sigMean1;
5013 sbits64 sigMean0;
5014 float128 z;
5016 aSig1 = extractFloat128Frac1( a );
5017 aSig0 = extractFloat128Frac0( a );
5018 aExp = extractFloat128Exp( a );
5019 aSign = extractFloat128Sign( a );
5020 bSig1 = extractFloat128Frac1( b );
5021 bSig0 = extractFloat128Frac0( b );
5022 bExp = extractFloat128Exp( b );
5023 bSign = extractFloat128Sign( b );
5024 if ( aExp == 0x7FFF ) {
5025 if ( ( aSig0 | aSig1 )
5026 || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
5027 return propagateFloat128NaN( a, b STATUS_VAR );
5029 goto invalid;
5031 if ( bExp == 0x7FFF ) {
5032 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
5033 return a;
5035 if ( bExp == 0 ) {
5036 if ( ( bSig0 | bSig1 ) == 0 ) {
5037 invalid:
5038 float_raise( float_flag_invalid STATUS_VAR);
5039 z.low = float128_default_nan_low;
5040 z.high = float128_default_nan_high;
5041 return z;
5043 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
5045 if ( aExp == 0 ) {
5046 if ( ( aSig0 | aSig1 ) == 0 ) return a;
5047 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
5049 expDiff = aExp - bExp;
5050 if ( expDiff < -1 ) return a;
5051 shortShift128Left(
5052 aSig0 | LIT64( 0x0001000000000000 ),
5053 aSig1,
5054 15 - ( expDiff < 0 ),
5055 &aSig0,
5056 &aSig1
5058 shortShift128Left(
5059 bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
5060 q = le128( bSig0, bSig1, aSig0, aSig1 );
5061 if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
5062 expDiff -= 64;
5063 while ( 0 < expDiff ) {
5064 q = estimateDiv128To64( aSig0, aSig1, bSig0 );
5065 q = ( 4 < q ) ? q - 4 : 0;
5066 mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
5067 shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZero );
5068 shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero );
5069 sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 );
5070 expDiff -= 61;
5072 if ( -64 < expDiff ) {
5073 q = estimateDiv128To64( aSig0, aSig1, bSig0 );
5074 q = ( 4 < q ) ? q - 4 : 0;
5075 q >>= - expDiff;
5076 shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
5077 expDiff += 52;
5078 if ( expDiff < 0 ) {
5079 shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
5081 else {
5082 shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
5084 mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
5085 sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
5087 else {
5088 shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 );
5089 shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
5091 do {
5092 alternateASig0 = aSig0;
5093 alternateASig1 = aSig1;
5094 ++q;
5095 sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
5096 } while ( 0 <= (sbits64) aSig0 );
5097 add128(
5098 aSig0, aSig1, alternateASig0, alternateASig1, (bits64 *)&sigMean0, &sigMean1 );
5099 if ( ( sigMean0 < 0 )
5100 || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
5101 aSig0 = alternateASig0;
5102 aSig1 = alternateASig1;
5104 zSign = ( (sbits64) aSig0 < 0 );
5105 if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
5106 return
5107 normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 STATUS_VAR );
5111 /*----------------------------------------------------------------------------
5112 | Returns the square root of the quadruple-precision floating-point value `a'.
5113 | The operation is performed according to the IEC/IEEE Standard for Binary
5114 | Floating-Point Arithmetic.
5115 *----------------------------------------------------------------------------*/
5117 float128 float128_sqrt( float128 a STATUS_PARAM )
5119 flag aSign;
5120 int32 aExp, zExp;
5121 bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
5122 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
5123 float128 z;
5125 aSig1 = extractFloat128Frac1( a );
5126 aSig0 = extractFloat128Frac0( a );
5127 aExp = extractFloat128Exp( a );
5128 aSign = extractFloat128Sign( a );
5129 if ( aExp == 0x7FFF ) {
5130 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, a STATUS_VAR );
5131 if ( ! aSign ) return a;
5132 goto invalid;
5134 if ( aSign ) {
5135 if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
5136 invalid:
5137 float_raise( float_flag_invalid STATUS_VAR);
5138 z.low = float128_default_nan_low;
5139 z.high = float128_default_nan_high;
5140 return z;
5142 if ( aExp == 0 ) {
5143 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( 0, 0, 0, 0 );
5144 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
5146 zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFE;
5147 aSig0 |= LIT64( 0x0001000000000000 );
5148 zSig0 = estimateSqrt32( aExp, aSig0>>17 );
5149 shortShift128Left( aSig0, aSig1, 13 - ( aExp & 1 ), &aSig0, &aSig1 );
5150 zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
5151 doubleZSig0 = zSig0<<1;
5152 mul64To128( zSig0, zSig0, &term0, &term1 );
5153 sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
5154 while ( (sbits64) rem0 < 0 ) {
5155 --zSig0;
5156 doubleZSig0 -= 2;
5157 add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
5159 zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
5160 if ( ( zSig1 & 0x1FFF ) <= 5 ) {
5161 if ( zSig1 == 0 ) zSig1 = 1;
5162 mul64To128( doubleZSig0, zSig1, &term1, &term2 );
5163 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
5164 mul64To128( zSig1, zSig1, &term2, &term3 );
5165 sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
5166 while ( (sbits64) rem1 < 0 ) {
5167 --zSig1;
5168 shortShift128Left( 0, zSig1, 1, &term2, &term3 );
5169 term3 |= 1;
5170 term2 |= doubleZSig0;
5171 add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
5173 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
5175 shift128ExtraRightJamming( zSig0, zSig1, 0, 14, &zSig0, &zSig1, &zSig2 );
5176 return roundAndPackFloat128( 0, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
5180 /*----------------------------------------------------------------------------
5181 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5182 | the corresponding value `b', and 0 otherwise. The comparison is performed
5183 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5184 *----------------------------------------------------------------------------*/
5186 int float128_eq( float128 a, float128 b STATUS_PARAM )
5189 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5190 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5191 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5192 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5194 if ( float128_is_signaling_nan( a )
5195 || float128_is_signaling_nan( b ) ) {
5196 float_raise( float_flag_invalid STATUS_VAR);
5198 return 0;
5200 return
5201 ( a.low == b.low )
5202 && ( ( a.high == b.high )
5203 || ( ( a.low == 0 )
5204 && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5209 /*----------------------------------------------------------------------------
5210 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5211 | or equal to the corresponding value `b', and 0 otherwise. The comparison
5212 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5213 | Arithmetic.
5214 *----------------------------------------------------------------------------*/
5216 int float128_le( float128 a, float128 b STATUS_PARAM )
5218 flag aSign, bSign;
5220 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5221 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5222 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5223 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5225 float_raise( float_flag_invalid STATUS_VAR);
5226 return 0;
5228 aSign = extractFloat128Sign( a );
5229 bSign = extractFloat128Sign( b );
5230 if ( aSign != bSign ) {
5231 return
5232 aSign
5233 || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5234 == 0 );
5236 return
5237 aSign ? le128( b.high, b.low, a.high, a.low )
5238 : le128( a.high, a.low, b.high, b.low );
5242 /*----------------------------------------------------------------------------
5243 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5244 | the corresponding value `b', and 0 otherwise. The comparison is performed
5245 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5246 *----------------------------------------------------------------------------*/
5248 int float128_lt( float128 a, float128 b STATUS_PARAM )
5250 flag aSign, bSign;
5252 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5253 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5254 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5255 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5257 float_raise( float_flag_invalid STATUS_VAR);
5258 return 0;
5260 aSign = extractFloat128Sign( a );
5261 bSign = extractFloat128Sign( b );
5262 if ( aSign != bSign ) {
5263 return
5264 aSign
5265 && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5266 != 0 );
5268 return
5269 aSign ? lt128( b.high, b.low, a.high, a.low )
5270 : lt128( a.high, a.low, b.high, b.low );
5274 /*----------------------------------------------------------------------------
5275 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5276 | the corresponding value `b', and 0 otherwise. The invalid exception is
5277 | raised if either operand is a NaN. Otherwise, the comparison is performed
5278 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5279 *----------------------------------------------------------------------------*/
5281 int float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
5284 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5285 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5286 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5287 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5289 float_raise( float_flag_invalid STATUS_VAR);
5290 return 0;
5292 return
5293 ( a.low == b.low )
5294 && ( ( a.high == b.high )
5295 || ( ( a.low == 0 )
5296 && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5301 /*----------------------------------------------------------------------------
5302 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5303 | or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
5304 | cause an exception. Otherwise, the comparison is performed according to the
5305 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5306 *----------------------------------------------------------------------------*/
5308 int float128_le_quiet( float128 a, float128 b STATUS_PARAM )
5310 flag aSign, bSign;
5312 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5313 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5314 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5315 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5317 if ( float128_is_signaling_nan( a )
5318 || float128_is_signaling_nan( b ) ) {
5319 float_raise( float_flag_invalid STATUS_VAR);
5321 return 0;
5323 aSign = extractFloat128Sign( a );
5324 bSign = extractFloat128Sign( b );
5325 if ( aSign != bSign ) {
5326 return
5327 aSign
5328 || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5329 == 0 );
5331 return
5332 aSign ? le128( b.high, b.low, a.high, a.low )
5333 : le128( a.high, a.low, b.high, b.low );
5337 /*----------------------------------------------------------------------------
5338 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5339 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
5340 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
5341 | Standard for Binary Floating-Point Arithmetic.
5342 *----------------------------------------------------------------------------*/
5344 int float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
5346 flag aSign, bSign;
5348 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5349 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5350 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5351 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5353 if ( float128_is_signaling_nan( a )
5354 || float128_is_signaling_nan( b ) ) {
5355 float_raise( float_flag_invalid STATUS_VAR);
5357 return 0;
5359 aSign = extractFloat128Sign( a );
5360 bSign = extractFloat128Sign( b );
5361 if ( aSign != bSign ) {
5362 return
5363 aSign
5364 && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5365 != 0 );
5367 return
5368 aSign ? lt128( b.high, b.low, a.high, a.low )
5369 : lt128( a.high, a.low, b.high, b.low );
5373 #endif
5375 /* misc functions */
5376 float32 uint32_to_float32( unsigned int a STATUS_PARAM )
5378 return int64_to_float32(a STATUS_VAR);
5381 float64 uint32_to_float64( unsigned int a STATUS_PARAM )
5383 return int64_to_float64(a STATUS_VAR);
5386 unsigned int float32_to_uint32( float32 a STATUS_PARAM )
5388 int64_t v;
5389 unsigned int res;
5391 v = float32_to_int64(a STATUS_VAR);
5392 if (v < 0) {
5393 res = 0;
5394 float_raise( float_flag_invalid STATUS_VAR);
5395 } else if (v > 0xffffffff) {
5396 res = 0xffffffff;
5397 float_raise( float_flag_invalid STATUS_VAR);
5398 } else {
5399 res = v;
5401 return res;
5404 unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
5406 int64_t v;
5407 unsigned int res;
5409 v = float32_to_int64_round_to_zero(a STATUS_VAR);
5410 if (v < 0) {
5411 res = 0;
5412 float_raise( float_flag_invalid STATUS_VAR);
5413 } else if (v > 0xffffffff) {
5414 res = 0xffffffff;
5415 float_raise( float_flag_invalid STATUS_VAR);
5416 } else {
5417 res = v;
5419 return res;
5422 unsigned int float64_to_uint32( float64 a STATUS_PARAM )
5424 int64_t v;
5425 unsigned int res;
5427 v = float64_to_int64(a STATUS_VAR);
5428 if (v < 0) {
5429 res = 0;
5430 float_raise( float_flag_invalid STATUS_VAR);
5431 } else if (v > 0xffffffff) {
5432 res = 0xffffffff;
5433 float_raise( float_flag_invalid STATUS_VAR);
5434 } else {
5435 res = v;
5437 return res;
5440 unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
5442 int64_t v;
5443 unsigned int res;
5445 v = float64_to_int64_round_to_zero(a STATUS_VAR);
5446 if (v < 0) {
5447 res = 0;
5448 float_raise( float_flag_invalid STATUS_VAR);
5449 } else if (v > 0xffffffff) {
5450 res = 0xffffffff;
5451 float_raise( float_flag_invalid STATUS_VAR);
5452 } else {
5453 res = v;
5455 return res;
5458 /* FIXME: This looks broken. */
5459 uint64_t float64_to_uint64 (float64 a STATUS_PARAM)
5461 int64_t v;
5463 v = float64_val(int64_to_float64(INT64_MIN STATUS_VAR));
5464 v += float64_val(a);
5465 v = float64_to_int64(make_float64(v) STATUS_VAR);
5467 return v - INT64_MIN;
5470 uint64_t float64_to_uint64_round_to_zero (float64 a STATUS_PARAM)
5472 int64_t v;
5474 v = float64_val(int64_to_float64(INT64_MIN STATUS_VAR));
5475 v += float64_val(a);
5476 v = float64_to_int64_round_to_zero(make_float64(v) STATUS_VAR);
5478 return v - INT64_MIN;
5481 #define COMPARE(s, nan_exp) \
5482 INLINE int float ## s ## _compare_internal( float ## s a, float ## s b, \
5483 int is_quiet STATUS_PARAM ) \
5485 flag aSign, bSign; \
5486 bits ## s av, bv; \
5488 if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) && \
5489 extractFloat ## s ## Frac( a ) ) || \
5490 ( ( extractFloat ## s ## Exp( b ) == nan_exp ) && \
5491 extractFloat ## s ## Frac( b ) )) { \
5492 if (!is_quiet || \
5493 float ## s ## _is_signaling_nan( a ) || \
5494 float ## s ## _is_signaling_nan( b ) ) { \
5495 float_raise( float_flag_invalid STATUS_VAR); \
5497 return float_relation_unordered; \
5499 aSign = extractFloat ## s ## Sign( a ); \
5500 bSign = extractFloat ## s ## Sign( b ); \
5501 av = float ## s ## _val(a); \
5502 bv = float ## s ## _val(b); \
5503 if ( aSign != bSign ) { \
5504 if ( (bits ## s) ( ( av | bv )<<1 ) == 0 ) { \
5505 /* zero case */ \
5506 return float_relation_equal; \
5507 } else { \
5508 return 1 - (2 * aSign); \
5510 } else { \
5511 if (av == bv) { \
5512 return float_relation_equal; \
5513 } else { \
5514 return 1 - 2 * (aSign ^ ( av < bv )); \
5519 int float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM ) \
5521 return float ## s ## _compare_internal(a, b, 0 STATUS_VAR); \
5524 int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
5526 return float ## s ## _compare_internal(a, b, 1 STATUS_VAR); \
5529 COMPARE(32, 0xff)
5530 COMPARE(64, 0x7ff)
5532 INLINE int float128_compare_internal( float128 a, float128 b,
5533 int is_quiet STATUS_PARAM )
5535 flag aSign, bSign;
5537 if (( ( extractFloat128Exp( a ) == 0x7fff ) &&
5538 ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) ) ||
5539 ( ( extractFloat128Exp( b ) == 0x7fff ) &&
5540 ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )) {
5541 if (!is_quiet ||
5542 float128_is_signaling_nan( a ) ||
5543 float128_is_signaling_nan( b ) ) {
5544 float_raise( float_flag_invalid STATUS_VAR);
5546 return float_relation_unordered;
5548 aSign = extractFloat128Sign( a );
5549 bSign = extractFloat128Sign( b );
5550 if ( aSign != bSign ) {
5551 if ( ( ( ( a.high | b.high )<<1 ) | a.low | b.low ) == 0 ) {
5552 /* zero case */
5553 return float_relation_equal;
5554 } else {
5555 return 1 - (2 * aSign);
5557 } else {
5558 if (a.low == b.low && a.high == b.high) {
5559 return float_relation_equal;
5560 } else {
5561 return 1 - 2 * (aSign ^ ( lt128( a.high, a.low, b.high, b.low ) ));
5566 int float128_compare( float128 a, float128 b STATUS_PARAM )
5568 return float128_compare_internal(a, b, 0 STATUS_VAR);
5571 int float128_compare_quiet( float128 a, float128 b STATUS_PARAM )
5573 return float128_compare_internal(a, b, 1 STATUS_VAR);
5576 /* Multiply A by 2 raised to the power N. */
5577 float32 float32_scalbn( float32 a, int n STATUS_PARAM )
5579 flag aSign;
5580 int16 aExp;
5581 bits32 aSig;
5583 aSig = extractFloat32Frac( a );
5584 aExp = extractFloat32Exp( a );
5585 aSign = extractFloat32Sign( a );
5587 if ( aExp == 0xFF ) {
5588 return a;
5590 if ( aExp != 0 )
5591 aSig |= 0x00800000;
5592 else if ( aSig == 0 )
5593 return a;
5595 aExp += n - 1;
5596 aSig <<= 7;
5597 return normalizeRoundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
5600 float64 float64_scalbn( float64 a, int n STATUS_PARAM )
5602 flag aSign;
5603 int16 aExp;
5604 bits64 aSig;
5606 aSig = extractFloat64Frac( a );
5607 aExp = extractFloat64Exp( a );
5608 aSign = extractFloat64Sign( a );
5610 if ( aExp == 0x7FF ) {
5611 return a;
5613 if ( aExp != 0 )
5614 aSig |= LIT64( 0x0010000000000000 );
5615 else if ( aSig == 0 )
5616 return a;
5618 aExp += n - 1;
5619 aSig <<= 10;
5620 return normalizeRoundAndPackFloat64( aSign, aExp, aSig STATUS_VAR );
5623 #ifdef FLOATX80
5624 floatx80 floatx80_scalbn( floatx80 a, int n STATUS_PARAM )
5626 flag aSign;
5627 int16 aExp;
5628 bits64 aSig;
5630 aSig = extractFloatx80Frac( a );
5631 aExp = extractFloatx80Exp( a );
5632 aSign = extractFloatx80Sign( a );
5634 if ( aExp == 0x7FF ) {
5635 return a;
5637 if (aExp == 0 && aSig == 0)
5638 return a;
5640 aExp += n;
5641 return normalizeRoundAndPackFloatx80( STATUS(floatx80_rounding_precision),
5642 aSign, aExp, aSig, 0 STATUS_VAR );
5644 #endif
5646 #ifdef FLOAT128
5647 float128 float128_scalbn( float128 a, int n STATUS_PARAM )
5649 flag aSign;
5650 int32 aExp;
5651 bits64 aSig0, aSig1;
5653 aSig1 = extractFloat128Frac1( a );
5654 aSig0 = extractFloat128Frac0( a );
5655 aExp = extractFloat128Exp( a );
5656 aSign = extractFloat128Sign( a );
5657 if ( aExp == 0x7FFF ) {
5658 return a;
5660 if ( aExp != 0 )
5661 aSig0 |= LIT64( 0x0001000000000000 );
5662 else if ( aSig0 == 0 && aSig1 == 0 )
5663 return a;
5665 aExp += n - 1;
5666 return normalizeRoundAndPackFloat128( aSign, aExp, aSig0, aSig1
5667 STATUS_VAR );
5670 #endif