debuggers.hg

view tools/ioemu/fpu/softfloat.c @ 0:7d21f7218375

Exact replica of unstable on 051908 + README-this
author Mukesh Rathor
date Mon May 19 15:34:57 2008 -0700 (2008-05-19)
parents
children
line source
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 #include "softfloat.h"
35 /*----------------------------------------------------------------------------
36 | Primitive arithmetic functions, including multi-word arithmetic, and
37 | division and square root approximations. (Can be specialized to target if
38 | desired.)
39 *----------------------------------------------------------------------------*/
40 #include "softfloat-macros.h"
42 /*----------------------------------------------------------------------------
43 | Functions and definitions to determine: (1) whether tininess for underflow
44 | is detected before or after rounding by default, (2) what (if anything)
45 | happens when exceptions are raised, (3) how signaling NaNs are distinguished
46 | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
47 | are propagated from function inputs to output. These details are target-
48 | specific.
49 *----------------------------------------------------------------------------*/
50 #include "softfloat-specialize.h"
52 void set_float_rounding_mode(int val STATUS_PARAM)
53 {
54 STATUS(float_rounding_mode) = val;
55 }
57 void set_float_exception_flags(int val STATUS_PARAM)
58 {
59 STATUS(float_exception_flags) = val;
60 }
62 #ifdef FLOATX80
63 void set_floatx80_rounding_precision(int val STATUS_PARAM)
64 {
65 STATUS(floatx80_rounding_precision) = val;
66 }
67 #endif
69 /*----------------------------------------------------------------------------
70 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
71 | and 7, and returns the properly rounded 32-bit integer corresponding to the
72 | input. If `zSign' is 1, the input is negated before being converted to an
73 | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
74 | is simply rounded to an integer, with the inexact exception raised if the
75 | input cannot be represented exactly as an integer. However, if the fixed-
76 | point input is too large, the invalid exception is raised and the largest
77 | positive or negative integer is returned.
78 *----------------------------------------------------------------------------*/
80 static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
81 {
82 int8 roundingMode;
83 flag roundNearestEven;
84 int8 roundIncrement, roundBits;
85 int32 z;
87 roundingMode = STATUS(float_rounding_mode);
88 roundNearestEven = ( roundingMode == float_round_nearest_even );
89 roundIncrement = 0x40;
90 if ( ! roundNearestEven ) {
91 if ( roundingMode == float_round_to_zero ) {
92 roundIncrement = 0;
93 }
94 else {
95 roundIncrement = 0x7F;
96 if ( zSign ) {
97 if ( roundingMode == float_round_up ) roundIncrement = 0;
98 }
99 else {
100 if ( roundingMode == float_round_down ) roundIncrement = 0;
101 }
102 }
103 }
104 roundBits = absZ & 0x7F;
105 absZ = ( absZ + roundIncrement )>>7;
106 absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
107 z = absZ;
108 if ( zSign ) z = - z;
109 if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
110 float_raise( float_flag_invalid STATUS_VAR);
111 return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
112 }
113 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
114 return z;
116 }
118 /*----------------------------------------------------------------------------
119 | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
120 | `absZ1', with binary point between bits 63 and 64 (between the input words),
121 | and returns the properly rounded 64-bit integer corresponding to the input.
122 | If `zSign' is 1, the input is negated before being converted to an integer.
123 | Ordinarily, the fixed-point input is simply rounded to an integer, with
124 | the inexact exception raised if the input cannot be represented exactly as
125 | an integer. However, if the fixed-point input is too large, the invalid
126 | exception is raised and the largest positive or negative integer is
127 | returned.
128 *----------------------------------------------------------------------------*/
130 static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PARAM)
131 {
132 int8 roundingMode;
133 flag roundNearestEven, increment;
134 int64 z;
136 roundingMode = STATUS(float_rounding_mode);
137 roundNearestEven = ( roundingMode == float_round_nearest_even );
138 increment = ( (sbits64) absZ1 < 0 );
139 if ( ! roundNearestEven ) {
140 if ( roundingMode == float_round_to_zero ) {
141 increment = 0;
142 }
143 else {
144 if ( zSign ) {
145 increment = ( roundingMode == float_round_down ) && absZ1;
146 }
147 else {
148 increment = ( roundingMode == float_round_up ) && absZ1;
149 }
150 }
151 }
152 if ( increment ) {
153 ++absZ0;
154 if ( absZ0 == 0 ) goto overflow;
155 absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
156 }
157 z = absZ0;
158 if ( zSign ) z = - z;
159 if ( z && ( ( z < 0 ) ^ zSign ) ) {
160 overflow:
161 float_raise( float_flag_invalid STATUS_VAR);
162 return
163 zSign ? (sbits64) LIT64( 0x8000000000000000 )
164 : LIT64( 0x7FFFFFFFFFFFFFFF );
165 }
166 if ( absZ1 ) STATUS(float_exception_flags) |= float_flag_inexact;
167 return z;
169 }
171 /*----------------------------------------------------------------------------
172 | Returns the fraction bits of the single-precision floating-point value `a'.
173 *----------------------------------------------------------------------------*/
175 INLINE bits32 extractFloat32Frac( float32 a )
176 {
178 return a & 0x007FFFFF;
180 }
182 /*----------------------------------------------------------------------------
183 | Returns the exponent bits of the single-precision floating-point value `a'.
184 *----------------------------------------------------------------------------*/
186 INLINE int16 extractFloat32Exp( float32 a )
187 {
189 return ( a>>23 ) & 0xFF;
191 }
193 /*----------------------------------------------------------------------------
194 | Returns the sign bit of the single-precision floating-point value `a'.
195 *----------------------------------------------------------------------------*/
197 INLINE flag extractFloat32Sign( float32 a )
198 {
200 return a>>31;
202 }
204 /*----------------------------------------------------------------------------
205 | Normalizes the subnormal single-precision floating-point value represented
206 | by the denormalized significand `aSig'. The normalized exponent and
207 | significand are stored at the locations pointed to by `zExpPtr' and
208 | `zSigPtr', respectively.
209 *----------------------------------------------------------------------------*/
211 static void
212 normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
213 {
214 int8 shiftCount;
216 shiftCount = countLeadingZeros32( aSig ) - 8;
217 *zSigPtr = aSig<<shiftCount;
218 *zExpPtr = 1 - shiftCount;
220 }
222 /*----------------------------------------------------------------------------
223 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
224 | single-precision floating-point value, returning the result. After being
225 | shifted into the proper positions, the three fields are simply added
226 | together to form the result. This means that any integer portion of `zSig'
227 | will be added into the exponent. Since a properly normalized significand
228 | will have an integer portion equal to 1, the `zExp' input should be 1 less
229 | than the desired result exponent whenever `zSig' is a complete, normalized
230 | significand.
231 *----------------------------------------------------------------------------*/
233 INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
234 {
236 return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
238 }
240 /*----------------------------------------------------------------------------
241 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
242 | and significand `zSig', and returns the proper single-precision floating-
243 | point value corresponding to the abstract input. Ordinarily, the abstract
244 | value is simply rounded and packed into the single-precision format, with
245 | the inexact exception raised if the abstract input cannot be represented
246 | exactly. However, if the abstract value is too large, the overflow and
247 | inexact exceptions are raised and an infinity or maximal finite value is
248 | returned. If the abstract value is too small, the input value is rounded to
249 | a subnormal number, and the underflow and inexact exceptions are raised if
250 | the abstract input cannot be represented exactly as a subnormal single-
251 | precision floating-point number.
252 | The input significand `zSig' has its binary point between bits 30
253 | and 29, which is 7 bits to the left of the usual location. This shifted
254 | significand must be normalized or smaller. If `zSig' is not normalized,
255 | `zExp' must be 0; in that case, the result returned is a subnormal number,
256 | and it must not require rounding. In the usual case that `zSig' is
257 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
258 | The handling of underflow and overflow follows the IEC/IEEE Standard for
259 | Binary Floating-Point Arithmetic.
260 *----------------------------------------------------------------------------*/
262 static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
263 {
264 int8 roundingMode;
265 flag roundNearestEven;
266 int8 roundIncrement, roundBits;
267 flag isTiny;
269 roundingMode = STATUS(float_rounding_mode);
270 roundNearestEven = ( roundingMode == float_round_nearest_even );
271 roundIncrement = 0x40;
272 if ( ! roundNearestEven ) {
273 if ( roundingMode == float_round_to_zero ) {
274 roundIncrement = 0;
275 }
276 else {
277 roundIncrement = 0x7F;
278 if ( zSign ) {
279 if ( roundingMode == float_round_up ) roundIncrement = 0;
280 }
281 else {
282 if ( roundingMode == float_round_down ) roundIncrement = 0;
283 }
284 }
285 }
286 roundBits = zSig & 0x7F;
287 if ( 0xFD <= (bits16) zExp ) {
288 if ( ( 0xFD < zExp )
289 || ( ( zExp == 0xFD )
290 && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
291 ) {
292 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
293 return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
294 }
295 if ( zExp < 0 ) {
296 isTiny =
297 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
298 || ( zExp < -1 )
299 || ( zSig + roundIncrement < 0x80000000 );
300 shift32RightJamming( zSig, - zExp, &zSig );
301 zExp = 0;
302 roundBits = zSig & 0x7F;
303 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
304 }
305 }
306 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
307 zSig = ( zSig + roundIncrement )>>7;
308 zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
309 if ( zSig == 0 ) zExp = 0;
310 return packFloat32( zSign, zExp, zSig );
312 }
314 /*----------------------------------------------------------------------------
315 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
316 | and significand `zSig', and returns the proper single-precision floating-
317 | point value corresponding to the abstract input. This routine is just like
318 | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
319 | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
320 | floating-point exponent.
321 *----------------------------------------------------------------------------*/
323 static float32
324 normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
325 {
326 int8 shiftCount;
328 shiftCount = countLeadingZeros32( zSig ) - 1;
329 return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
331 }
333 /*----------------------------------------------------------------------------
334 | Returns the fraction bits of the double-precision floating-point value `a'.
335 *----------------------------------------------------------------------------*/
337 INLINE bits64 extractFloat64Frac( float64 a )
338 {
340 return a & LIT64( 0x000FFFFFFFFFFFFF );
342 }
344 /*----------------------------------------------------------------------------
345 | Returns the exponent bits of the double-precision floating-point value `a'.
346 *----------------------------------------------------------------------------*/
348 INLINE int16 extractFloat64Exp( float64 a )
349 {
351 return ( a>>52 ) & 0x7FF;
353 }
355 /*----------------------------------------------------------------------------
356 | Returns the sign bit of the double-precision floating-point value `a'.
357 *----------------------------------------------------------------------------*/
359 INLINE flag extractFloat64Sign( float64 a )
360 {
362 return a>>63;
364 }
366 /*----------------------------------------------------------------------------
367 | Normalizes the subnormal double-precision floating-point value represented
368 | by the denormalized significand `aSig'. The normalized exponent and
369 | significand are stored at the locations pointed to by `zExpPtr' and
370 | `zSigPtr', respectively.
371 *----------------------------------------------------------------------------*/
373 static void
374 normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
375 {
376 int8 shiftCount;
378 shiftCount = countLeadingZeros64( aSig ) - 11;
379 *zSigPtr = aSig<<shiftCount;
380 *zExpPtr = 1 - shiftCount;
382 }
384 /*----------------------------------------------------------------------------
385 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
386 | double-precision floating-point value, returning the result. After being
387 | shifted into the proper positions, the three fields are simply added
388 | together to form the result. This means that any integer portion of `zSig'
389 | will be added into the exponent. Since a properly normalized significand
390 | will have an integer portion equal to 1, the `zExp' input should be 1 less
391 | than the desired result exponent whenever `zSig' is a complete, normalized
392 | significand.
393 *----------------------------------------------------------------------------*/
395 INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
396 {
398 return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;
400 }
402 /*----------------------------------------------------------------------------
403 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
404 | and significand `zSig', and returns the proper double-precision floating-
405 | point value corresponding to the abstract input. Ordinarily, the abstract
406 | value is simply rounded and packed into the double-precision format, with
407 | the inexact exception raised if the abstract input cannot be represented
408 | exactly. However, if the abstract value is too large, the overflow and
409 | inexact exceptions are raised and an infinity or maximal finite value is
410 | returned. If the abstract value is too small, the input value is rounded
411 | to a subnormal number, and the underflow and inexact exceptions are raised
412 | if the abstract input cannot be represented exactly as a subnormal double-
413 | precision floating-point number.
414 | The input significand `zSig' has its binary point between bits 62
415 | and 61, which is 10 bits to the left of the usual location. This shifted
416 | significand must be normalized or smaller. If `zSig' is not normalized,
417 | `zExp' must be 0; in that case, the result returned is a subnormal number,
418 | and it must not require rounding. In the usual case that `zSig' is
419 | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
420 | The handling of underflow and overflow follows the IEC/IEEE Standard for
421 | Binary Floating-Point Arithmetic.
422 *----------------------------------------------------------------------------*/
424 static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
425 {
426 int8 roundingMode;
427 flag roundNearestEven;
428 int16 roundIncrement, roundBits;
429 flag isTiny;
431 roundingMode = STATUS(float_rounding_mode);
432 roundNearestEven = ( roundingMode == float_round_nearest_even );
433 roundIncrement = 0x200;
434 if ( ! roundNearestEven ) {
435 if ( roundingMode == float_round_to_zero ) {
436 roundIncrement = 0;
437 }
438 else {
439 roundIncrement = 0x3FF;
440 if ( zSign ) {
441 if ( roundingMode == float_round_up ) roundIncrement = 0;
442 }
443 else {
444 if ( roundingMode == float_round_down ) roundIncrement = 0;
445 }
446 }
447 }
448 roundBits = zSig & 0x3FF;
449 if ( 0x7FD <= (bits16) zExp ) {
450 if ( ( 0x7FD < zExp )
451 || ( ( zExp == 0x7FD )
452 && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
453 ) {
454 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
455 return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );
456 }
457 if ( zExp < 0 ) {
458 isTiny =
459 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
460 || ( zExp < -1 )
461 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
462 shift64RightJamming( zSig, - zExp, &zSig );
463 zExp = 0;
464 roundBits = zSig & 0x3FF;
465 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
466 }
467 }
468 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
469 zSig = ( zSig + roundIncrement )>>10;
470 zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
471 if ( zSig == 0 ) zExp = 0;
472 return packFloat64( zSign, zExp, zSig );
474 }
476 /*----------------------------------------------------------------------------
477 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
478 | and significand `zSig', and returns the proper double-precision floating-
479 | point value corresponding to the abstract input. This routine is just like
480 | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
481 | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
482 | floating-point exponent.
483 *----------------------------------------------------------------------------*/
485 static float64
486 normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
487 {
488 int8 shiftCount;
490 shiftCount = countLeadingZeros64( zSig ) - 1;
491 return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
493 }
495 #ifdef FLOATX80
497 /*----------------------------------------------------------------------------
498 | Returns the fraction bits of the extended double-precision floating-point
499 | value `a'.
500 *----------------------------------------------------------------------------*/
502 INLINE bits64 extractFloatx80Frac( floatx80 a )
503 {
505 return a.low;
507 }
509 /*----------------------------------------------------------------------------
510 | Returns the exponent bits of the extended double-precision floating-point
511 | value `a'.
512 *----------------------------------------------------------------------------*/
514 INLINE int32 extractFloatx80Exp( floatx80 a )
515 {
517 return a.high & 0x7FFF;
519 }
521 /*----------------------------------------------------------------------------
522 | Returns the sign bit of the extended double-precision floating-point value
523 | `a'.
524 *----------------------------------------------------------------------------*/
526 INLINE flag extractFloatx80Sign( floatx80 a )
527 {
529 return a.high>>15;
531 }
533 /*----------------------------------------------------------------------------
534 | Normalizes the subnormal extended double-precision floating-point value
535 | represented by the denormalized significand `aSig'. The normalized exponent
536 | and significand are stored at the locations pointed to by `zExpPtr' and
537 | `zSigPtr', respectively.
538 *----------------------------------------------------------------------------*/
540 static void
541 normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
542 {
543 int8 shiftCount;
545 shiftCount = countLeadingZeros64( aSig );
546 *zSigPtr = aSig<<shiftCount;
547 *zExpPtr = 1 - shiftCount;
549 }
551 /*----------------------------------------------------------------------------
552 | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
553 | extended double-precision floating-point value, returning the result.
554 *----------------------------------------------------------------------------*/
556 INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
557 {
558 floatx80 z;
560 z.low = zSig;
561 z.high = ( ( (bits16) zSign )<<15 ) + zExp;
562 return z;
564 }
566 /*----------------------------------------------------------------------------
567 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
568 | and extended significand formed by the concatenation of `zSig0' and `zSig1',
569 | and returns the proper extended double-precision floating-point value
570 | corresponding to the abstract input. Ordinarily, the abstract value is
571 | rounded and packed into the extended double-precision format, with the
572 | inexact exception raised if the abstract input cannot be represented
573 | exactly. However, if the abstract value is too large, the overflow and
574 | inexact exceptions are raised and an infinity or maximal finite value is
575 | returned. If the abstract value is too small, the input value is rounded to
576 | a subnormal number, and the underflow and inexact exceptions are raised if
577 | the abstract input cannot be represented exactly as a subnormal extended
578 | double-precision floating-point number.
579 | If `roundingPrecision' is 32 or 64, the result is rounded to the same
580 | number of bits as single or double precision, respectively. Otherwise, the
581 | result is rounded to the full precision of the extended double-precision
582 | format.
583 | The input significand must be normalized or smaller. If the input
584 | significand is not normalized, `zExp' must be 0; in that case, the result
585 | returned is a subnormal number, and it must not require rounding. The
586 | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
587 | Floating-Point Arithmetic.
588 *----------------------------------------------------------------------------*/
590 static floatx80
591 roundAndPackFloatx80(
592 int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
593 STATUS_PARAM)
594 {
595 int8 roundingMode;
596 flag roundNearestEven, increment, isTiny;
597 int64 roundIncrement, roundMask, roundBits;
599 roundingMode = STATUS(float_rounding_mode);
600 roundNearestEven = ( roundingMode == float_round_nearest_even );
601 if ( roundingPrecision == 80 ) goto precision80;
602 if ( roundingPrecision == 64 ) {
603 roundIncrement = LIT64( 0x0000000000000400 );
604 roundMask = LIT64( 0x00000000000007FF );
605 }
606 else if ( roundingPrecision == 32 ) {
607 roundIncrement = LIT64( 0x0000008000000000 );
608 roundMask = LIT64( 0x000000FFFFFFFFFF );
609 }
610 else {
611 goto precision80;
612 }
613 zSig0 |= ( zSig1 != 0 );
614 if ( ! roundNearestEven ) {
615 if ( roundingMode == float_round_to_zero ) {
616 roundIncrement = 0;
617 }
618 else {
619 roundIncrement = roundMask;
620 if ( zSign ) {
621 if ( roundingMode == float_round_up ) roundIncrement = 0;
622 }
623 else {
624 if ( roundingMode == float_round_down ) roundIncrement = 0;
625 }
626 }
627 }
628 roundBits = zSig0 & roundMask;
629 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
630 if ( ( 0x7FFE < zExp )
631 || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
632 ) {
633 goto overflow;
634 }
635 if ( zExp <= 0 ) {
636 isTiny =
637 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
638 || ( zExp < 0 )
639 || ( zSig0 <= zSig0 + roundIncrement );
640 shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
641 zExp = 0;
642 roundBits = zSig0 & roundMask;
643 if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
644 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
645 zSig0 += roundIncrement;
646 if ( (sbits64) zSig0 < 0 ) zExp = 1;
647 roundIncrement = roundMask + 1;
648 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
649 roundMask |= roundIncrement;
650 }
651 zSig0 &= ~ roundMask;
652 return packFloatx80( zSign, zExp, zSig0 );
653 }
654 }
655 if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
656 zSig0 += roundIncrement;
657 if ( zSig0 < roundIncrement ) {
658 ++zExp;
659 zSig0 = LIT64( 0x8000000000000000 );
660 }
661 roundIncrement = roundMask + 1;
662 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
663 roundMask |= roundIncrement;
664 }
665 zSig0 &= ~ roundMask;
666 if ( zSig0 == 0 ) zExp = 0;
667 return packFloatx80( zSign, zExp, zSig0 );
668 precision80:
669 increment = ( (sbits64) zSig1 < 0 );
670 if ( ! roundNearestEven ) {
671 if ( roundingMode == float_round_to_zero ) {
672 increment = 0;
673 }
674 else {
675 if ( zSign ) {
676 increment = ( roundingMode == float_round_down ) && zSig1;
677 }
678 else {
679 increment = ( roundingMode == float_round_up ) && zSig1;
680 }
681 }
682 }
683 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
684 if ( ( 0x7FFE < zExp )
685 || ( ( zExp == 0x7FFE )
686 && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
687 && increment
688 )
689 ) {
690 roundMask = 0;
691 overflow:
692 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
693 if ( ( roundingMode == float_round_to_zero )
694 || ( zSign && ( roundingMode == float_round_up ) )
695 || ( ! zSign && ( roundingMode == float_round_down ) )
696 ) {
697 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
698 }
699 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
700 }
701 if ( zExp <= 0 ) {
702 isTiny =
703 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
704 || ( zExp < 0 )
705 || ! increment
706 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
707 shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
708 zExp = 0;
709 if ( isTiny && zSig1 ) float_raise( float_flag_underflow STATUS_VAR);
710 if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
711 if ( roundNearestEven ) {
712 increment = ( (sbits64) zSig1 < 0 );
713 }
714 else {
715 if ( zSign ) {
716 increment = ( roundingMode == float_round_down ) && zSig1;
717 }
718 else {
719 increment = ( roundingMode == float_round_up ) && zSig1;
720 }
721 }
722 if ( increment ) {
723 ++zSig0;
724 zSig0 &=
725 ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
726 if ( (sbits64) zSig0 < 0 ) zExp = 1;
727 }
728 return packFloatx80( zSign, zExp, zSig0 );
729 }
730 }
731 if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
732 if ( increment ) {
733 ++zSig0;
734 if ( zSig0 == 0 ) {
735 ++zExp;
736 zSig0 = LIT64( 0x8000000000000000 );
737 }
738 else {
739 zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
740 }
741 }
742 else {
743 if ( zSig0 == 0 ) zExp = 0;
744 }
745 return packFloatx80( zSign, zExp, zSig0 );
747 }
749 /*----------------------------------------------------------------------------
750 | Takes an abstract floating-point value having sign `zSign', exponent
751 | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
752 | and returns the proper extended double-precision floating-point value
753 | corresponding to the abstract input. This routine is just like
754 | `roundAndPackFloatx80' except that the input significand does not have to be
755 | normalized.
756 *----------------------------------------------------------------------------*/
758 static floatx80
759 normalizeRoundAndPackFloatx80(
760 int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
761 STATUS_PARAM)
762 {
763 int8 shiftCount;
765 if ( zSig0 == 0 ) {
766 zSig0 = zSig1;
767 zSig1 = 0;
768 zExp -= 64;
769 }
770 shiftCount = countLeadingZeros64( zSig0 );
771 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
772 zExp -= shiftCount;
773 return
774 roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 STATUS_VAR);
776 }
778 #endif
780 #ifdef FLOAT128
782 /*----------------------------------------------------------------------------
783 | Returns the least-significant 64 fraction bits of the quadruple-precision
784 | floating-point value `a'.
785 *----------------------------------------------------------------------------*/
787 INLINE bits64 extractFloat128Frac1( float128 a )
788 {
790 return a.low;
792 }
794 /*----------------------------------------------------------------------------
795 | Returns the most-significant 48 fraction bits of the quadruple-precision
796 | floating-point value `a'.
797 *----------------------------------------------------------------------------*/
799 INLINE bits64 extractFloat128Frac0( float128 a )
800 {
802 return a.high & LIT64( 0x0000FFFFFFFFFFFF );
804 }
806 /*----------------------------------------------------------------------------
807 | Returns the exponent bits of the quadruple-precision floating-point value
808 | `a'.
809 *----------------------------------------------------------------------------*/
811 INLINE int32 extractFloat128Exp( float128 a )
812 {
814 return ( a.high>>48 ) & 0x7FFF;
816 }
818 /*----------------------------------------------------------------------------
819 | Returns the sign bit of the quadruple-precision floating-point value `a'.
820 *----------------------------------------------------------------------------*/
822 INLINE flag extractFloat128Sign( float128 a )
823 {
825 return a.high>>63;
827 }
829 /*----------------------------------------------------------------------------
830 | Normalizes the subnormal quadruple-precision floating-point value
831 | represented by the denormalized significand formed by the concatenation of
832 | `aSig0' and `aSig1'. The normalized exponent is stored at the location
833 | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
834 | significand are stored at the location pointed to by `zSig0Ptr', and the
835 | least significant 64 bits of the normalized significand are stored at the
836 | location pointed to by `zSig1Ptr'.
837 *----------------------------------------------------------------------------*/
839 static void
840 normalizeFloat128Subnormal(
841 bits64 aSig0,
842 bits64 aSig1,
843 int32 *zExpPtr,
844 bits64 *zSig0Ptr,
845 bits64 *zSig1Ptr
846 )
847 {
848 int8 shiftCount;
850 if ( aSig0 == 0 ) {
851 shiftCount = countLeadingZeros64( aSig1 ) - 15;
852 if ( shiftCount < 0 ) {
853 *zSig0Ptr = aSig1>>( - shiftCount );
854 *zSig1Ptr = aSig1<<( shiftCount & 63 );
855 }
856 else {
857 *zSig0Ptr = aSig1<<shiftCount;
858 *zSig1Ptr = 0;
859 }
860 *zExpPtr = - shiftCount - 63;
861 }
862 else {
863 shiftCount = countLeadingZeros64( aSig0 ) - 15;
864 shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
865 *zExpPtr = 1 - shiftCount;
866 }
868 }
870 /*----------------------------------------------------------------------------
871 | Packs the sign `zSign', the exponent `zExp', and the significand formed
872 | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
873 | floating-point value, returning the result. After being shifted into the
874 | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
875 | added together to form the most significant 32 bits of the result. This
876 | means that any integer portion of `zSig0' will be added into the exponent.
877 | Since a properly normalized significand will have an integer portion equal
878 | to 1, the `zExp' input should be 1 less than the desired result exponent
879 | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
880 | significand.
881 *----------------------------------------------------------------------------*/
883 INLINE float128
884 packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
885 {
886 float128 z;
888 z.low = zSig1;
889 z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
890 return z;
892 }
894 /*----------------------------------------------------------------------------
895 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
896 | and extended significand formed by the concatenation of `zSig0', `zSig1',
897 | and `zSig2', and returns the proper quadruple-precision floating-point value
898 | corresponding to the abstract input. Ordinarily, the abstract value is
899 | simply rounded and packed into the quadruple-precision format, with the
900 | inexact exception raised if the abstract input cannot be represented
901 | exactly. However, if the abstract value is too large, the overflow and
902 | inexact exceptions are raised and an infinity or maximal finite value is
903 | returned. If the abstract value is too small, the input value is rounded to
904 | a subnormal number, and the underflow and inexact exceptions are raised if
905 | the abstract input cannot be represented exactly as a subnormal quadruple-
906 | precision floating-point number.
907 | The input significand must be normalized or smaller. If the input
908 | significand is not normalized, `zExp' must be 0; in that case, the result
909 | returned is a subnormal number, and it must not require rounding. In the
910 | usual case that the input significand is normalized, `zExp' must be 1 less
911 | than the ``true'' floating-point exponent. The handling of underflow and
912 | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
913 *----------------------------------------------------------------------------*/
915 static float128
916 roundAndPackFloat128(
917 flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 STATUS_PARAM)
918 {
919 int8 roundingMode;
920 flag roundNearestEven, increment, isTiny;
922 roundingMode = STATUS(float_rounding_mode);
923 roundNearestEven = ( roundingMode == float_round_nearest_even );
924 increment = ( (sbits64) zSig2 < 0 );
925 if ( ! roundNearestEven ) {
926 if ( roundingMode == float_round_to_zero ) {
927 increment = 0;
928 }
929 else {
930 if ( zSign ) {
931 increment = ( roundingMode == float_round_down ) && zSig2;
932 }
933 else {
934 increment = ( roundingMode == float_round_up ) && zSig2;
935 }
936 }
937 }
938 if ( 0x7FFD <= (bits32) zExp ) {
939 if ( ( 0x7FFD < zExp )
940 || ( ( zExp == 0x7FFD )
941 && eq128(
942 LIT64( 0x0001FFFFFFFFFFFF ),
943 LIT64( 0xFFFFFFFFFFFFFFFF ),
944 zSig0,
945 zSig1
946 )
947 && increment
948 )
949 ) {
950 float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
951 if ( ( roundingMode == float_round_to_zero )
952 || ( zSign && ( roundingMode == float_round_up ) )
953 || ( ! zSign && ( roundingMode == float_round_down ) )
954 ) {
955 return
956 packFloat128(
957 zSign,
958 0x7FFE,
959 LIT64( 0x0000FFFFFFFFFFFF ),
960 LIT64( 0xFFFFFFFFFFFFFFFF )
961 );
962 }
963 return packFloat128( zSign, 0x7FFF, 0, 0 );
964 }
965 if ( zExp < 0 ) {
966 isTiny =
967 ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
968 || ( zExp < -1 )
969 || ! increment
970 || lt128(
971 zSig0,
972 zSig1,
973 LIT64( 0x0001FFFFFFFFFFFF ),
974 LIT64( 0xFFFFFFFFFFFFFFFF )
975 );
976 shift128ExtraRightJamming(
977 zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
978 zExp = 0;
979 if ( isTiny && zSig2 ) float_raise( float_flag_underflow STATUS_VAR);
980 if ( roundNearestEven ) {
981 increment = ( (sbits64) zSig2 < 0 );
982 }
983 else {
984 if ( zSign ) {
985 increment = ( roundingMode == float_round_down ) && zSig2;
986 }
987 else {
988 increment = ( roundingMode == float_round_up ) && zSig2;
989 }
990 }
991 }
992 }
993 if ( zSig2 ) STATUS(float_exception_flags) |= float_flag_inexact;
994 if ( increment ) {
995 add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
996 zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
997 }
998 else {
999 if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
1001 return packFloat128( zSign, zExp, zSig0, zSig1 );
1005 /*----------------------------------------------------------------------------
1006 | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1007 | and significand formed by the concatenation of `zSig0' and `zSig1', and
1008 | returns the proper quadruple-precision floating-point value corresponding
1009 | to the abstract input. This routine is just like `roundAndPackFloat128'
1010 | except that the input significand has fewer bits and does not have to be
1011 | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
1012 | point exponent.
1013 *----------------------------------------------------------------------------*/
1015 static float128
1016 normalizeRoundAndPackFloat128(
1017 flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 STATUS_PARAM)
1019 int8 shiftCount;
1020 bits64 zSig2;
1022 if ( zSig0 == 0 ) {
1023 zSig0 = zSig1;
1024 zSig1 = 0;
1025 zExp -= 64;
1027 shiftCount = countLeadingZeros64( zSig0 ) - 15;
1028 if ( 0 <= shiftCount ) {
1029 zSig2 = 0;
1030 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1032 else {
1033 shift128ExtraRightJamming(
1034 zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
1036 zExp -= shiftCount;
1037 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR);
1041 #endif
1043 /*----------------------------------------------------------------------------
1044 | Returns the result of converting the 32-bit two's complement integer `a'
1045 | to the single-precision floating-point format. The conversion is performed
1046 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1047 *----------------------------------------------------------------------------*/
1049 float32 int32_to_float32( int32 a STATUS_PARAM )
1051 flag zSign;
1053 if ( a == 0 ) return 0;
1054 if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
1055 zSign = ( a < 0 );
1056 return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a STATUS_VAR );
1060 /*----------------------------------------------------------------------------
1061 | Returns the result of converting the 32-bit two's complement integer `a'
1062 | to the double-precision floating-point format. The conversion is performed
1063 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1064 *----------------------------------------------------------------------------*/
1066 float64 int32_to_float64( int32 a STATUS_PARAM )
1068 flag zSign;
1069 uint32 absA;
1070 int8 shiftCount;
1071 bits64 zSig;
1073 if ( a == 0 ) return 0;
1074 zSign = ( a < 0 );
1075 absA = zSign ? - a : a;
1076 shiftCount = countLeadingZeros32( absA ) + 21;
1077 zSig = absA;
1078 return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );
1082 #ifdef FLOATX80
1084 /*----------------------------------------------------------------------------
1085 | Returns the result of converting the 32-bit two's complement integer `a'
1086 | to the extended double-precision floating-point format. The conversion
1087 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1088 | Arithmetic.
1089 *----------------------------------------------------------------------------*/
1091 floatx80 int32_to_floatx80( int32 a STATUS_PARAM )
1093 flag zSign;
1094 uint32 absA;
1095 int8 shiftCount;
1096 bits64 zSig;
1098 if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1099 zSign = ( a < 0 );
1100 absA = zSign ? - a : a;
1101 shiftCount = countLeadingZeros32( absA ) + 32;
1102 zSig = absA;
1103 return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
1107 #endif
1109 #ifdef FLOAT128
1111 /*----------------------------------------------------------------------------
1112 | Returns the result of converting the 32-bit two's complement integer `a' to
1113 | the quadruple-precision floating-point format. The conversion is performed
1114 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1115 *----------------------------------------------------------------------------*/
1117 float128 int32_to_float128( int32 a STATUS_PARAM )
1119 flag zSign;
1120 uint32 absA;
1121 int8 shiftCount;
1122 bits64 zSig0;
1124 if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1125 zSign = ( a < 0 );
1126 absA = zSign ? - a : a;
1127 shiftCount = countLeadingZeros32( absA ) + 17;
1128 zSig0 = absA;
1129 return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );
1133 #endif
1135 /*----------------------------------------------------------------------------
1136 | Returns the result of converting the 64-bit two's complement integer `a'
1137 | to the single-precision floating-point format. The conversion is performed
1138 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1139 *----------------------------------------------------------------------------*/
1141 float32 int64_to_float32( int64 a STATUS_PARAM )
1143 flag zSign;
1144 uint64 absA;
1145 int8 shiftCount;
1147 if ( a == 0 ) return 0;
1148 zSign = ( a < 0 );
1149 absA = zSign ? - a : a;
1150 shiftCount = countLeadingZeros64( absA ) - 40;
1151 if ( 0 <= shiftCount ) {
1152 return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );
1154 else {
1155 shiftCount += 7;
1156 if ( shiftCount < 0 ) {
1157 shift64RightJamming( absA, - shiftCount, &absA );
1159 else {
1160 absA <<= shiftCount;
1162 return roundAndPackFloat32( zSign, 0x9C - shiftCount, absA STATUS_VAR );
1167 /*----------------------------------------------------------------------------
1168 | Returns the result of converting the 64-bit two's complement integer `a'
1169 | to the double-precision floating-point format. The conversion is performed
1170 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1171 *----------------------------------------------------------------------------*/
1173 float64 int64_to_float64( int64 a STATUS_PARAM )
1175 flag zSign;
1177 if ( a == 0 ) return 0;
1178 if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
1179 return packFloat64( 1, 0x43E, 0 );
1181 zSign = ( a < 0 );
1182 return normalizeRoundAndPackFloat64( zSign, 0x43C, zSign ? - a : a STATUS_VAR );
1186 #ifdef FLOATX80
1188 /*----------------------------------------------------------------------------
1189 | Returns the result of converting the 64-bit two's complement integer `a'
1190 | to the extended double-precision floating-point format. The conversion
1191 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1192 | Arithmetic.
1193 *----------------------------------------------------------------------------*/
1195 floatx80 int64_to_floatx80( int64 a STATUS_PARAM )
1197 flag zSign;
1198 uint64 absA;
1199 int8 shiftCount;
1201 if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1202 zSign = ( a < 0 );
1203 absA = zSign ? - a : a;
1204 shiftCount = countLeadingZeros64( absA );
1205 return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
1209 #endif
1211 #ifdef FLOAT128
1213 /*----------------------------------------------------------------------------
1214 | Returns the result of converting the 64-bit two's complement integer `a' to
1215 | the quadruple-precision floating-point format. The conversion is performed
1216 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1217 *----------------------------------------------------------------------------*/
1219 float128 int64_to_float128( int64 a STATUS_PARAM )
1221 flag zSign;
1222 uint64 absA;
1223 int8 shiftCount;
1224 int32 zExp;
1225 bits64 zSig0, zSig1;
1227 if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1228 zSign = ( a < 0 );
1229 absA = zSign ? - a : a;
1230 shiftCount = countLeadingZeros64( absA ) + 49;
1231 zExp = 0x406E - shiftCount;
1232 if ( 64 <= shiftCount ) {
1233 zSig1 = 0;
1234 zSig0 = absA;
1235 shiftCount -= 64;
1237 else {
1238 zSig1 = absA;
1239 zSig0 = 0;
1241 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1242 return packFloat128( zSign, zExp, zSig0, zSig1 );
1246 #endif
1248 /*----------------------------------------------------------------------------
1249 | Returns the result of converting the single-precision floating-point value
1250 | `a' to the 32-bit two's complement integer format. The conversion is
1251 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1252 | Arithmetic---which means in particular that the conversion is rounded
1253 | according to the current rounding mode. If `a' is a NaN, the largest
1254 | positive integer is returned. Otherwise, if the conversion overflows, the
1255 | largest integer with the same sign as `a' is returned.
1256 *----------------------------------------------------------------------------*/
1258 int32 float32_to_int32( float32 a STATUS_PARAM )
1260 flag aSign;
1261 int16 aExp, shiftCount;
1262 bits32 aSig;
1263 bits64 aSig64;
1265 aSig = extractFloat32Frac( a );
1266 aExp = extractFloat32Exp( a );
1267 aSign = extractFloat32Sign( a );
1268 if ( ( aExp == 0xFF ) && aSig ) aSign = 0;
1269 if ( aExp ) aSig |= 0x00800000;
1270 shiftCount = 0xAF - aExp;
1271 aSig64 = aSig;
1272 aSig64 <<= 32;
1273 if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );
1274 return roundAndPackInt32( aSign, aSig64 STATUS_VAR );
1278 /*----------------------------------------------------------------------------
1279 | Returns the result of converting the single-precision floating-point value
1280 | `a' to the 32-bit two's complement integer format. The conversion is
1281 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1282 | Arithmetic, except that the conversion is always rounded toward zero.
1283 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1284 | the conversion overflows, the largest integer with the same sign as `a' is
1285 | returned.
1286 *----------------------------------------------------------------------------*/
1288 int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
1290 flag aSign;
1291 int16 aExp, shiftCount;
1292 bits32 aSig;
1293 int32 z;
1295 aSig = extractFloat32Frac( a );
1296 aExp = extractFloat32Exp( a );
1297 aSign = extractFloat32Sign( a );
1298 shiftCount = aExp - 0x9E;
1299 if ( 0 <= shiftCount ) {
1300 if ( a != 0xCF000000 ) {
1301 float_raise( float_flag_invalid STATUS_VAR);
1302 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
1304 return (sbits32) 0x80000000;
1306 else if ( aExp <= 0x7E ) {
1307 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1308 return 0;
1310 aSig = ( aSig | 0x00800000 )<<8;
1311 z = aSig>>( - shiftCount );
1312 if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
1313 STATUS(float_exception_flags) |= float_flag_inexact;
1315 if ( aSign ) z = - z;
1316 return z;
1320 /*----------------------------------------------------------------------------
1321 | Returns the result of converting the single-precision floating-point value
1322 | `a' to the 64-bit two's complement integer format. The conversion is
1323 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1324 | Arithmetic---which means in particular that the conversion is rounded
1325 | according to the current rounding mode. If `a' is a NaN, the largest
1326 | positive integer is returned. Otherwise, if the conversion overflows, the
1327 | largest integer with the same sign as `a' is returned.
1328 *----------------------------------------------------------------------------*/
1330 int64 float32_to_int64( float32 a STATUS_PARAM )
1332 flag aSign;
1333 int16 aExp, shiftCount;
1334 bits32 aSig;
1335 bits64 aSig64, aSigExtra;
1337 aSig = extractFloat32Frac( a );
1338 aExp = extractFloat32Exp( a );
1339 aSign = extractFloat32Sign( a );
1340 shiftCount = 0xBE - aExp;
1341 if ( shiftCount < 0 ) {
1342 float_raise( float_flag_invalid STATUS_VAR);
1343 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1344 return LIT64( 0x7FFFFFFFFFFFFFFF );
1346 return (sbits64) LIT64( 0x8000000000000000 );
1348 if ( aExp ) aSig |= 0x00800000;
1349 aSig64 = aSig;
1350 aSig64 <<= 40;
1351 shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );
1352 return roundAndPackInt64( aSign, aSig64, aSigExtra STATUS_VAR );
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, except that the conversion is always rounded toward zero. If
1361 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
1362 | conversion overflows, the largest integer with the same sign as `a' is
1363 | returned.
1364 *----------------------------------------------------------------------------*/
1366 int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
1368 flag aSign;
1369 int16 aExp, shiftCount;
1370 bits32 aSig;
1371 bits64 aSig64;
1372 int64 z;
1374 aSig = extractFloat32Frac( a );
1375 aExp = extractFloat32Exp( a );
1376 aSign = extractFloat32Sign( a );
1377 shiftCount = aExp - 0xBE;
1378 if ( 0 <= shiftCount ) {
1379 if ( a != 0xDF000000 ) {
1380 float_raise( float_flag_invalid STATUS_VAR);
1381 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1382 return LIT64( 0x7FFFFFFFFFFFFFFF );
1385 return (sbits64) LIT64( 0x8000000000000000 );
1387 else if ( aExp <= 0x7E ) {
1388 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1389 return 0;
1391 aSig64 = aSig | 0x00800000;
1392 aSig64 <<= 40;
1393 z = aSig64>>( - shiftCount );
1394 if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
1395 STATUS(float_exception_flags) |= float_flag_inexact;
1397 if ( aSign ) z = - z;
1398 return z;
1402 /*----------------------------------------------------------------------------
1403 | Returns the result of converting the single-precision floating-point value
1404 | `a' to the double-precision floating-point format. The conversion is
1405 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1406 | Arithmetic.
1407 *----------------------------------------------------------------------------*/
1409 float64 float32_to_float64( float32 a STATUS_PARAM )
1411 flag aSign;
1412 int16 aExp;
1413 bits32 aSig;
1415 aSig = extractFloat32Frac( a );
1416 aExp = extractFloat32Exp( a );
1417 aSign = extractFloat32Sign( a );
1418 if ( aExp == 0xFF ) {
1419 if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ));
1420 return packFloat64( aSign, 0x7FF, 0 );
1422 if ( aExp == 0 ) {
1423 if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
1424 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1425 --aExp;
1427 return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
1431 #ifdef FLOATX80
1433 /*----------------------------------------------------------------------------
1434 | Returns the result of converting the single-precision floating-point value
1435 | `a' to the extended double-precision floating-point format. The conversion
1436 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1437 | Arithmetic.
1438 *----------------------------------------------------------------------------*/
1440 floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
1442 flag aSign;
1443 int16 aExp;
1444 bits32 aSig;
1446 aSig = extractFloat32Frac( a );
1447 aExp = extractFloat32Exp( a );
1448 aSign = extractFloat32Sign( a );
1449 if ( aExp == 0xFF ) {
1450 if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) );
1451 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
1453 if ( aExp == 0 ) {
1454 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1455 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1457 aSig |= 0x00800000;
1458 return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
1462 #endif
1464 #ifdef FLOAT128
1466 /*----------------------------------------------------------------------------
1467 | Returns the result of converting the single-precision floating-point value
1468 | `a' to the double-precision floating-point format. The conversion is
1469 | performed according to the IEC/IEEE Standard for Binary Floating-Point
1470 | Arithmetic.
1471 *----------------------------------------------------------------------------*/
1473 float128 float32_to_float128( float32 a STATUS_PARAM )
1475 flag aSign;
1476 int16 aExp;
1477 bits32 aSig;
1479 aSig = extractFloat32Frac( a );
1480 aExp = extractFloat32Exp( a );
1481 aSign = extractFloat32Sign( a );
1482 if ( aExp == 0xFF ) {
1483 if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) );
1484 return packFloat128( aSign, 0x7FFF, 0, 0 );
1486 if ( aExp == 0 ) {
1487 if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
1488 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1489 --aExp;
1491 return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
1495 #endif
1497 /*----------------------------------------------------------------------------
1498 | Rounds the single-precision floating-point value `a' to an integer, and
1499 | returns the result as a single-precision floating-point value. The
1500 | operation is performed according to the IEC/IEEE Standard for Binary
1501 | Floating-Point Arithmetic.
1502 *----------------------------------------------------------------------------*/
1504 float32 float32_round_to_int( float32 a STATUS_PARAM)
1506 flag aSign;
1507 int16 aExp;
1508 bits32 lastBitMask, roundBitsMask;
1509 int8 roundingMode;
1510 float32 z;
1512 aExp = extractFloat32Exp( a );
1513 if ( 0x96 <= aExp ) {
1514 if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
1515 return propagateFloat32NaN( a, a STATUS_VAR );
1517 return a;
1519 if ( aExp <= 0x7E ) {
1520 if ( (bits32) ( a<<1 ) == 0 ) return a;
1521 STATUS(float_exception_flags) |= float_flag_inexact;
1522 aSign = extractFloat32Sign( a );
1523 switch ( STATUS(float_rounding_mode) ) {
1524 case float_round_nearest_even:
1525 if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
1526 return packFloat32( aSign, 0x7F, 0 );
1528 break;
1529 case float_round_down:
1530 return aSign ? 0xBF800000 : 0;
1531 case float_round_up:
1532 return aSign ? 0x80000000 : 0x3F800000;
1534 return packFloat32( aSign, 0, 0 );
1536 lastBitMask = 1;
1537 lastBitMask <<= 0x96 - aExp;
1538 roundBitsMask = lastBitMask - 1;
1539 z = a;
1540 roundingMode = STATUS(float_rounding_mode);
1541 if ( roundingMode == float_round_nearest_even ) {
1542 z += lastBitMask>>1;
1543 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
1545 else if ( roundingMode != float_round_to_zero ) {
1546 if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
1547 z += roundBitsMask;
1550 z &= ~ roundBitsMask;
1551 if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
1552 return z;
1556 /*----------------------------------------------------------------------------
1557 | Returns the result of adding the absolute values of the single-precision
1558 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
1559 | before being returned. `zSign' is ignored if the result is a NaN.
1560 | The addition is performed according to the IEC/IEEE Standard for Binary
1561 | Floating-Point Arithmetic.
1562 *----------------------------------------------------------------------------*/
1564 static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1566 int16 aExp, bExp, zExp;
1567 bits32 aSig, bSig, zSig;
1568 int16 expDiff;
1570 aSig = extractFloat32Frac( a );
1571 aExp = extractFloat32Exp( a );
1572 bSig = extractFloat32Frac( b );
1573 bExp = extractFloat32Exp( b );
1574 expDiff = aExp - bExp;
1575 aSig <<= 6;
1576 bSig <<= 6;
1577 if ( 0 < expDiff ) {
1578 if ( aExp == 0xFF ) {
1579 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1580 return a;
1582 if ( bExp == 0 ) {
1583 --expDiff;
1585 else {
1586 bSig |= 0x20000000;
1588 shift32RightJamming( bSig, expDiff, &bSig );
1589 zExp = aExp;
1591 else if ( expDiff < 0 ) {
1592 if ( bExp == 0xFF ) {
1593 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1594 return packFloat32( zSign, 0xFF, 0 );
1596 if ( aExp == 0 ) {
1597 ++expDiff;
1599 else {
1600 aSig |= 0x20000000;
1602 shift32RightJamming( aSig, - expDiff, &aSig );
1603 zExp = bExp;
1605 else {
1606 if ( aExp == 0xFF ) {
1607 if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1608 return a;
1610 if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
1611 zSig = 0x40000000 + aSig + bSig;
1612 zExp = aExp;
1613 goto roundAndPack;
1615 aSig |= 0x20000000;
1616 zSig = ( aSig + bSig )<<1;
1617 --zExp;
1618 if ( (sbits32) zSig < 0 ) {
1619 zSig = aSig + bSig;
1620 ++zExp;
1622 roundAndPack:
1623 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1627 /*----------------------------------------------------------------------------
1628 | Returns the result of subtracting the absolute values of the single-
1629 | precision floating-point values `a' and `b'. If `zSign' is 1, the
1630 | difference is negated before being returned. `zSign' is ignored if the
1631 | result is a NaN. The subtraction is performed according to the IEC/IEEE
1632 | Standard for Binary Floating-Point Arithmetic.
1633 *----------------------------------------------------------------------------*/
1635 static float32 subFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1637 int16 aExp, bExp, zExp;
1638 bits32 aSig, bSig, zSig;
1639 int16 expDiff;
1641 aSig = extractFloat32Frac( a );
1642 aExp = extractFloat32Exp( a );
1643 bSig = extractFloat32Frac( b );
1644 bExp = extractFloat32Exp( b );
1645 expDiff = aExp - bExp;
1646 aSig <<= 7;
1647 bSig <<= 7;
1648 if ( 0 < expDiff ) goto aExpBigger;
1649 if ( expDiff < 0 ) goto bExpBigger;
1650 if ( aExp == 0xFF ) {
1651 if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1652 float_raise( float_flag_invalid STATUS_VAR);
1653 return float32_default_nan;
1655 if ( aExp == 0 ) {
1656 aExp = 1;
1657 bExp = 1;
1659 if ( bSig < aSig ) goto aBigger;
1660 if ( aSig < bSig ) goto bBigger;
1661 return packFloat32( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
1662 bExpBigger:
1663 if ( bExp == 0xFF ) {
1664 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1665 return packFloat32( zSign ^ 1, 0xFF, 0 );
1667 if ( aExp == 0 ) {
1668 ++expDiff;
1670 else {
1671 aSig |= 0x40000000;
1673 shift32RightJamming( aSig, - expDiff, &aSig );
1674 bSig |= 0x40000000;
1675 bBigger:
1676 zSig = bSig - aSig;
1677 zExp = bExp;
1678 zSign ^= 1;
1679 goto normalizeRoundAndPack;
1680 aExpBigger:
1681 if ( aExp == 0xFF ) {
1682 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1683 return a;
1685 if ( bExp == 0 ) {
1686 --expDiff;
1688 else {
1689 bSig |= 0x40000000;
1691 shift32RightJamming( bSig, expDiff, &bSig );
1692 aSig |= 0x40000000;
1693 aBigger:
1694 zSig = aSig - bSig;
1695 zExp = aExp;
1696 normalizeRoundAndPack:
1697 --zExp;
1698 return normalizeRoundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1702 /*----------------------------------------------------------------------------
1703 | Returns the result of adding the single-precision floating-point values `a'
1704 | and `b'. The operation is performed according to the IEC/IEEE Standard for
1705 | Binary Floating-Point Arithmetic.
1706 *----------------------------------------------------------------------------*/
1708 float32 float32_add( float32 a, float32 b STATUS_PARAM )
1710 flag aSign, bSign;
1712 aSign = extractFloat32Sign( a );
1713 bSign = extractFloat32Sign( b );
1714 if ( aSign == bSign ) {
1715 return addFloat32Sigs( a, b, aSign STATUS_VAR);
1717 else {
1718 return subFloat32Sigs( a, b, aSign STATUS_VAR );
1723 /*----------------------------------------------------------------------------
1724 | Returns the result of subtracting the single-precision floating-point values
1725 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1726 | for Binary Floating-Point Arithmetic.
1727 *----------------------------------------------------------------------------*/
1729 float32 float32_sub( float32 a, float32 b STATUS_PARAM )
1731 flag aSign, bSign;
1733 aSign = extractFloat32Sign( a );
1734 bSign = extractFloat32Sign( b );
1735 if ( aSign == bSign ) {
1736 return subFloat32Sigs( a, b, aSign STATUS_VAR );
1738 else {
1739 return addFloat32Sigs( a, b, aSign STATUS_VAR );
1744 /*----------------------------------------------------------------------------
1745 | Returns the result of multiplying the single-precision floating-point values
1746 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1747 | for Binary Floating-Point Arithmetic.
1748 *----------------------------------------------------------------------------*/
1750 float32 float32_mul( float32 a, float32 b STATUS_PARAM )
1752 flag aSign, bSign, zSign;
1753 int16 aExp, bExp, zExp;
1754 bits32 aSig, bSig;
1755 bits64 zSig64;
1756 bits32 zSig;
1758 aSig = extractFloat32Frac( a );
1759 aExp = extractFloat32Exp( a );
1760 aSign = extractFloat32Sign( a );
1761 bSig = extractFloat32Frac( b );
1762 bExp = extractFloat32Exp( b );
1763 bSign = extractFloat32Sign( b );
1764 zSign = aSign ^ bSign;
1765 if ( aExp == 0xFF ) {
1766 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1767 return propagateFloat32NaN( a, b STATUS_VAR );
1769 if ( ( bExp | bSig ) == 0 ) {
1770 float_raise( float_flag_invalid STATUS_VAR);
1771 return float32_default_nan;
1773 return packFloat32( zSign, 0xFF, 0 );
1775 if ( bExp == 0xFF ) {
1776 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1777 if ( ( aExp | aSig ) == 0 ) {
1778 float_raise( float_flag_invalid STATUS_VAR);
1779 return float32_default_nan;
1781 return packFloat32( zSign, 0xFF, 0 );
1783 if ( aExp == 0 ) {
1784 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1785 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1787 if ( bExp == 0 ) {
1788 if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
1789 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1791 zExp = aExp + bExp - 0x7F;
1792 aSig = ( aSig | 0x00800000 )<<7;
1793 bSig = ( bSig | 0x00800000 )<<8;
1794 shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
1795 zSig = zSig64;
1796 if ( 0 <= (sbits32) ( zSig<<1 ) ) {
1797 zSig <<= 1;
1798 --zExp;
1800 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1804 /*----------------------------------------------------------------------------
1805 | Returns the result of dividing the single-precision floating-point value `a'
1806 | by the corresponding value `b'. The operation is performed according to the
1807 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1808 *----------------------------------------------------------------------------*/
1810 float32 float32_div( float32 a, float32 b STATUS_PARAM )
1812 flag aSign, bSign, zSign;
1813 int16 aExp, bExp, zExp;
1814 bits32 aSig, bSig, zSig;
1816 aSig = extractFloat32Frac( a );
1817 aExp = extractFloat32Exp( a );
1818 aSign = extractFloat32Sign( a );
1819 bSig = extractFloat32Frac( b );
1820 bExp = extractFloat32Exp( b );
1821 bSign = extractFloat32Sign( b );
1822 zSign = aSign ^ bSign;
1823 if ( aExp == 0xFF ) {
1824 if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1825 if ( bExp == 0xFF ) {
1826 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1827 float_raise( float_flag_invalid STATUS_VAR);
1828 return float32_default_nan;
1830 return packFloat32( zSign, 0xFF, 0 );
1832 if ( bExp == 0xFF ) {
1833 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1834 return packFloat32( zSign, 0, 0 );
1836 if ( bExp == 0 ) {
1837 if ( bSig == 0 ) {
1838 if ( ( aExp | aSig ) == 0 ) {
1839 float_raise( float_flag_invalid STATUS_VAR);
1840 return float32_default_nan;
1842 float_raise( float_flag_divbyzero STATUS_VAR);
1843 return packFloat32( zSign, 0xFF, 0 );
1845 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1847 if ( aExp == 0 ) {
1848 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1849 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1851 zExp = aExp - bExp + 0x7D;
1852 aSig = ( aSig | 0x00800000 )<<7;
1853 bSig = ( bSig | 0x00800000 )<<8;
1854 if ( bSig <= ( aSig + aSig ) ) {
1855 aSig >>= 1;
1856 ++zExp;
1858 zSig = ( ( (bits64) aSig )<<32 ) / bSig;
1859 if ( ( zSig & 0x3F ) == 0 ) {
1860 zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
1862 return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1866 /*----------------------------------------------------------------------------
1867 | Returns the remainder of the single-precision floating-point value `a'
1868 | with respect to the corresponding value `b'. The operation is performed
1869 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1870 *----------------------------------------------------------------------------*/
1872 float32 float32_rem( float32 a, float32 b STATUS_PARAM )
1874 flag aSign, bSign, zSign;
1875 int16 aExp, bExp, expDiff;
1876 bits32 aSig, bSig;
1877 bits32 q;
1878 bits64 aSig64, bSig64, q64;
1879 bits32 alternateASig;
1880 sbits32 sigMean;
1882 aSig = extractFloat32Frac( a );
1883 aExp = extractFloat32Exp( a );
1884 aSign = extractFloat32Sign( a );
1885 bSig = extractFloat32Frac( b );
1886 bExp = extractFloat32Exp( b );
1887 bSign = extractFloat32Sign( b );
1888 if ( aExp == 0xFF ) {
1889 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1890 return propagateFloat32NaN( a, b STATUS_VAR );
1892 float_raise( float_flag_invalid STATUS_VAR);
1893 return float32_default_nan;
1895 if ( bExp == 0xFF ) {
1896 if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1897 return a;
1899 if ( bExp == 0 ) {
1900 if ( bSig == 0 ) {
1901 float_raise( float_flag_invalid STATUS_VAR);
1902 return float32_default_nan;
1904 normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1906 if ( aExp == 0 ) {
1907 if ( aSig == 0 ) return a;
1908 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1910 expDiff = aExp - bExp;
1911 aSig |= 0x00800000;
1912 bSig |= 0x00800000;
1913 if ( expDiff < 32 ) {
1914 aSig <<= 8;
1915 bSig <<= 8;
1916 if ( expDiff < 0 ) {
1917 if ( expDiff < -1 ) return a;
1918 aSig >>= 1;
1920 q = ( bSig <= aSig );
1921 if ( q ) aSig -= bSig;
1922 if ( 0 < expDiff ) {
1923 q = ( ( (bits64) aSig )<<32 ) / bSig;
1924 q >>= 32 - expDiff;
1925 bSig >>= 2;
1926 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
1928 else {
1929 aSig >>= 2;
1930 bSig >>= 2;
1933 else {
1934 if ( bSig <= aSig ) aSig -= bSig;
1935 aSig64 = ( (bits64) aSig )<<40;
1936 bSig64 = ( (bits64) bSig )<<40;
1937 expDiff -= 64;
1938 while ( 0 < expDiff ) {
1939 q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1940 q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1941 aSig64 = - ( ( bSig * q64 )<<38 );
1942 expDiff -= 62;
1944 expDiff += 64;
1945 q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1946 q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1947 q = q64>>( 64 - expDiff );
1948 bSig <<= 6;
1949 aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
1951 do {
1952 alternateASig = aSig;
1953 ++q;
1954 aSig -= bSig;
1955 } while ( 0 <= (sbits32) aSig );
1956 sigMean = aSig + alternateASig;
1957 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
1958 aSig = alternateASig;
1960 zSign = ( (sbits32) aSig < 0 );
1961 if ( zSign ) aSig = - aSig;
1962 return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig STATUS_VAR );
1966 /*----------------------------------------------------------------------------
1967 | Returns the square root of the single-precision floating-point value `a'.
1968 | The operation is performed according to the IEC/IEEE Standard for Binary
1969 | Floating-Point Arithmetic.
1970 *----------------------------------------------------------------------------*/
1972 float32 float32_sqrt( float32 a STATUS_PARAM )
1974 flag aSign;
1975 int16 aExp, zExp;
1976 bits32 aSig, zSig;
1977 bits64 rem, term;
1979 aSig = extractFloat32Frac( a );
1980 aExp = extractFloat32Exp( a );
1981 aSign = extractFloat32Sign( a );
1982 if ( aExp == 0xFF ) {
1983 if ( aSig ) return propagateFloat32NaN( a, 0 STATUS_VAR );
1984 if ( ! aSign ) return a;
1985 float_raise( float_flag_invalid STATUS_VAR);
1986 return float32_default_nan;
1988 if ( aSign ) {
1989 if ( ( aExp | aSig ) == 0 ) return a;
1990 float_raise( float_flag_invalid STATUS_VAR);
1991 return float32_default_nan;
1993 if ( aExp == 0 ) {
1994 if ( aSig == 0 ) return 0;
1995 normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1997 zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
1998 aSig = ( aSig | 0x00800000 )<<8;
1999 zSig = estimateSqrt32( aExp, aSig ) + 2;
2000 if ( ( zSig & 0x7F ) <= 5 ) {
2001 if ( zSig < 2 ) {
2002 zSig = 0x7FFFFFFF;
2003 goto roundAndPack;
2005 aSig >>= aExp & 1;
2006 term = ( (bits64) zSig ) * zSig;
2007 rem = ( ( (bits64) aSig )<<32 ) - term;
2008 while ( (sbits64) rem < 0 ) {
2009 --zSig;
2010 rem += ( ( (bits64) zSig )<<1 ) | 1;
2012 zSig |= ( rem != 0 );
2014 shift32RightJamming( zSig, 1, &zSig );
2015 roundAndPack:
2016 return roundAndPackFloat32( 0, zExp, zSig STATUS_VAR );
2020 /*----------------------------------------------------------------------------
2021 | Returns 1 if the single-precision floating-point value `a' is equal to
2022 | the corresponding value `b', and 0 otherwise. The comparison is performed
2023 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2024 *----------------------------------------------------------------------------*/
2026 int float32_eq( float32 a, float32 b STATUS_PARAM )
2029 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2030 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2031 ) {
2032 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2033 float_raise( float_flag_invalid STATUS_VAR);
2035 return 0;
2037 return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
2041 /*----------------------------------------------------------------------------
2042 | Returns 1 if the single-precision floating-point value `a' is less than
2043 | or equal to the corresponding value `b', and 0 otherwise. The comparison
2044 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2045 | Arithmetic.
2046 *----------------------------------------------------------------------------*/
2048 int float32_le( float32 a, float32 b STATUS_PARAM )
2050 flag aSign, bSign;
2052 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2053 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2054 ) {
2055 float_raise( float_flag_invalid STATUS_VAR);
2056 return 0;
2058 aSign = extractFloat32Sign( a );
2059 bSign = extractFloat32Sign( b );
2060 if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
2061 return ( a == b ) || ( aSign ^ ( a < b ) );
2065 /*----------------------------------------------------------------------------
2066 | Returns 1 if the single-precision floating-point value `a' is less than
2067 | the corresponding value `b', and 0 otherwise. The comparison is performed
2068 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2069 *----------------------------------------------------------------------------*/
2071 int float32_lt( float32 a, float32 b STATUS_PARAM )
2073 flag aSign, bSign;
2075 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2076 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2077 ) {
2078 float_raise( float_flag_invalid STATUS_VAR);
2079 return 0;
2081 aSign = extractFloat32Sign( a );
2082 bSign = extractFloat32Sign( b );
2083 if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
2084 return ( a != b ) && ( aSign ^ ( a < b ) );
2088 /*----------------------------------------------------------------------------
2089 | Returns 1 if the single-precision floating-point value `a' is equal to
2090 | the corresponding value `b', and 0 otherwise. The invalid exception is
2091 | raised if either operand is a NaN. Otherwise, the comparison is performed
2092 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2093 *----------------------------------------------------------------------------*/
2095 int float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
2098 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2099 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2100 ) {
2101 float_raise( float_flag_invalid STATUS_VAR);
2102 return 0;
2104 return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
2108 /*----------------------------------------------------------------------------
2109 | Returns 1 if the single-precision floating-point value `a' is less than or
2110 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
2111 | cause an exception. Otherwise, the comparison is performed according to the
2112 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2113 *----------------------------------------------------------------------------*/
2115 int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
2117 flag aSign, bSign;
2119 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2120 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2121 ) {
2122 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2123 float_raise( float_flag_invalid STATUS_VAR);
2125 return 0;
2127 aSign = extractFloat32Sign( a );
2128 bSign = extractFloat32Sign( b );
2129 if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
2130 return ( a == b ) || ( aSign ^ ( a < b ) );
2134 /*----------------------------------------------------------------------------
2135 | Returns 1 if the single-precision floating-point value `a' is less than
2136 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2137 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
2138 | Standard for Binary Floating-Point Arithmetic.
2139 *----------------------------------------------------------------------------*/
2141 int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
2143 flag aSign, bSign;
2145 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2146 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2147 ) {
2148 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2149 float_raise( float_flag_invalid STATUS_VAR);
2151 return 0;
2153 aSign = extractFloat32Sign( a );
2154 bSign = extractFloat32Sign( b );
2155 if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
2156 return ( a != b ) && ( aSign ^ ( a < b ) );
2160 /*----------------------------------------------------------------------------
2161 | Returns the result of converting the double-precision floating-point value
2162 | `a' to the 32-bit two's complement integer format. The conversion is
2163 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2164 | Arithmetic---which means in particular that the conversion is rounded
2165 | according to the current rounding mode. If `a' is a NaN, the largest
2166 | positive integer is returned. Otherwise, if the conversion overflows, the
2167 | largest integer with the same sign as `a' is returned.
2168 *----------------------------------------------------------------------------*/
2170 int32 float64_to_int32( float64 a STATUS_PARAM )
2172 flag aSign;
2173 int16 aExp, shiftCount;
2174 bits64 aSig;
2176 aSig = extractFloat64Frac( a );
2177 aExp = extractFloat64Exp( a );
2178 aSign = extractFloat64Sign( a );
2179 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2180 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2181 shiftCount = 0x42C - aExp;
2182 if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
2183 return roundAndPackInt32( aSign, aSig STATUS_VAR );
2187 /*----------------------------------------------------------------------------
2188 | Returns the result of converting the double-precision floating-point value
2189 | `a' to the 32-bit two's complement integer format. The conversion is
2190 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2191 | Arithmetic, except that the conversion is always rounded toward zero.
2192 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2193 | the conversion overflows, the largest integer with the same sign as `a' is
2194 | returned.
2195 *----------------------------------------------------------------------------*/
2197 int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
2199 flag aSign;
2200 int16 aExp, shiftCount;
2201 bits64 aSig, savedASig;
2202 int32 z;
2204 aSig = extractFloat64Frac( a );
2205 aExp = extractFloat64Exp( a );
2206 aSign = extractFloat64Sign( a );
2207 if ( 0x41E < aExp ) {
2208 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2209 goto invalid;
2211 else if ( aExp < 0x3FF ) {
2212 if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2213 return 0;
2215 aSig |= LIT64( 0x0010000000000000 );
2216 shiftCount = 0x433 - aExp;
2217 savedASig = aSig;
2218 aSig >>= shiftCount;
2219 z = aSig;
2220 if ( aSign ) z = - z;
2221 if ( ( z < 0 ) ^ aSign ) {
2222 invalid:
2223 float_raise( float_flag_invalid STATUS_VAR);
2224 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
2226 if ( ( aSig<<shiftCount ) != savedASig ) {
2227 STATUS(float_exception_flags) |= float_flag_inexact;
2229 return z;
2233 /*----------------------------------------------------------------------------
2234 | Returns the result of converting the double-precision floating-point value
2235 | `a' to the 64-bit two's complement integer format. The conversion is
2236 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2237 | Arithmetic---which means in particular that the conversion is rounded
2238 | according to the current rounding mode. If `a' is a NaN, the largest
2239 | positive integer is returned. Otherwise, if the conversion overflows, the
2240 | largest integer with the same sign as `a' is returned.
2241 *----------------------------------------------------------------------------*/
2243 int64 float64_to_int64( float64 a STATUS_PARAM )
2245 flag aSign;
2246 int16 aExp, shiftCount;
2247 bits64 aSig, aSigExtra;
2249 aSig = extractFloat64Frac( a );
2250 aExp = extractFloat64Exp( a );
2251 aSign = extractFloat64Sign( a );
2252 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2253 shiftCount = 0x433 - aExp;
2254 if ( shiftCount <= 0 ) {
2255 if ( 0x43E < aExp ) {
2256 float_raise( float_flag_invalid STATUS_VAR);
2257 if ( ! aSign
2258 || ( ( aExp == 0x7FF )
2259 && ( aSig != LIT64( 0x0010000000000000 ) ) )
2260 ) {
2261 return LIT64( 0x7FFFFFFFFFFFFFFF );
2263 return (sbits64) LIT64( 0x8000000000000000 );
2265 aSigExtra = 0;
2266 aSig <<= - shiftCount;
2268 else {
2269 shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
2271 return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
2275 /*----------------------------------------------------------------------------
2276 | Returns the result of converting the double-precision floating-point value
2277 | `a' to the 64-bit two's complement integer format. The conversion is
2278 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2279 | Arithmetic, except that the conversion is always rounded toward zero.
2280 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2281 | the conversion overflows, the largest integer with the same sign as `a' is
2282 | returned.
2283 *----------------------------------------------------------------------------*/
2285 int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
2287 flag aSign;
2288 int16 aExp, shiftCount;
2289 bits64 aSig;
2290 int64 z;
2292 aSig = extractFloat64Frac( a );
2293 aExp = extractFloat64Exp( a );
2294 aSign = extractFloat64Sign( a );
2295 if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2296 shiftCount = aExp - 0x433;
2297 if ( 0 <= shiftCount ) {
2298 if ( 0x43E <= aExp ) {
2299 if ( a != LIT64( 0xC3E0000000000000 ) ) {
2300 float_raise( float_flag_invalid STATUS_VAR);
2301 if ( ! aSign
2302 || ( ( aExp == 0x7FF )
2303 && ( aSig != LIT64( 0x0010000000000000 ) ) )
2304 ) {
2305 return LIT64( 0x7FFFFFFFFFFFFFFF );
2308 return (sbits64) LIT64( 0x8000000000000000 );
2310 z = aSig<<shiftCount;
2312 else {
2313 if ( aExp < 0x3FE ) {
2314 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2315 return 0;
2317 z = aSig>>( - shiftCount );
2318 if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
2319 STATUS(float_exception_flags) |= float_flag_inexact;
2322 if ( aSign ) z = - z;
2323 return z;
2327 /*----------------------------------------------------------------------------
2328 | Returns the result of converting the double-precision floating-point value
2329 | `a' to the single-precision floating-point format. The conversion is
2330 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2331 | Arithmetic.
2332 *----------------------------------------------------------------------------*/
2334 float32 float64_to_float32( float64 a STATUS_PARAM )
2336 flag aSign;
2337 int16 aExp;
2338 bits64 aSig;
2339 bits32 zSig;
2341 aSig = extractFloat64Frac( a );
2342 aExp = extractFloat64Exp( a );
2343 aSign = extractFloat64Sign( a );
2344 if ( aExp == 0x7FF ) {
2345 if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR ) );
2346 return packFloat32( aSign, 0xFF, 0 );
2348 shift64RightJamming( aSig, 22, &aSig );
2349 zSig = aSig;
2350 if ( aExp || zSig ) {
2351 zSig |= 0x40000000;
2352 aExp -= 0x381;
2354 return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
2358 #ifdef FLOATX80
2360 /*----------------------------------------------------------------------------
2361 | Returns the result of converting the double-precision floating-point value
2362 | `a' to the extended double-precision floating-point format. The conversion
2363 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2364 | Arithmetic.
2365 *----------------------------------------------------------------------------*/
2367 floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
2369 flag aSign;
2370 int16 aExp;
2371 bits64 aSig;
2373 aSig = extractFloat64Frac( a );
2374 aExp = extractFloat64Exp( a );
2375 aSign = extractFloat64Sign( a );
2376 if ( aExp == 0x7FF ) {
2377 if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) );
2378 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
2380 if ( aExp == 0 ) {
2381 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
2382 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2384 return
2385 packFloatx80(
2386 aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
2390 #endif
2392 #ifdef FLOAT128
2394 /*----------------------------------------------------------------------------
2395 | Returns the result of converting the double-precision floating-point value
2396 | `a' to the quadruple-precision floating-point format. The conversion is
2397 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2398 | Arithmetic.
2399 *----------------------------------------------------------------------------*/
2401 float128 float64_to_float128( float64 a STATUS_PARAM )
2403 flag aSign;
2404 int16 aExp;
2405 bits64 aSig, zSig0, zSig1;
2407 aSig = extractFloat64Frac( a );
2408 aExp = extractFloat64Exp( a );
2409 aSign = extractFloat64Sign( a );
2410 if ( aExp == 0x7FF ) {
2411 if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR ) );
2412 return packFloat128( aSign, 0x7FFF, 0, 0 );
2414 if ( aExp == 0 ) {
2415 if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
2416 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2417 --aExp;
2419 shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
2420 return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
2424 #endif
2426 /*----------------------------------------------------------------------------
2427 | Rounds the double-precision floating-point value `a' to an integer, and
2428 | returns the result as a double-precision floating-point value. The
2429 | operation is performed according to the IEC/IEEE Standard for Binary
2430 | Floating-Point Arithmetic.
2431 *----------------------------------------------------------------------------*/
2433 float64 float64_round_to_int( float64 a STATUS_PARAM )
2435 flag aSign;
2436 int16 aExp;
2437 bits64 lastBitMask, roundBitsMask;
2438 int8 roundingMode;
2439 float64 z;
2441 aExp = extractFloat64Exp( a );
2442 if ( 0x433 <= aExp ) {
2443 if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
2444 return propagateFloat64NaN( a, a STATUS_VAR );
2446 return a;
2448 if ( aExp < 0x3FF ) {
2449 if ( (bits64) ( a<<1 ) == 0 ) return a;
2450 STATUS(float_exception_flags) |= float_flag_inexact;
2451 aSign = extractFloat64Sign( a );
2452 switch ( STATUS(float_rounding_mode) ) {
2453 case float_round_nearest_even:
2454 if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
2455 return packFloat64( aSign, 0x3FF, 0 );
2457 break;
2458 case float_round_down:
2459 return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
2460 case float_round_up:
2461 return
2462 aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
2464 return packFloat64( aSign, 0, 0 );
2466 lastBitMask = 1;
2467 lastBitMask <<= 0x433 - aExp;
2468 roundBitsMask = lastBitMask - 1;
2469 z = a;
2470 roundingMode = STATUS(float_rounding_mode);
2471 if ( roundingMode == float_round_nearest_even ) {
2472 z += lastBitMask>>1;
2473 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
2475 else if ( roundingMode != float_round_to_zero ) {
2476 if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
2477 z += roundBitsMask;
2480 z &= ~ roundBitsMask;
2481 if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
2482 return z;
2486 float64 float64_trunc_to_int( float64 a STATUS_PARAM)
2488 int oldmode;
2489 float64 res;
2490 oldmode = STATUS(float_rounding_mode);
2491 STATUS(float_rounding_mode) = float_round_to_zero;
2492 res = float64_round_to_int(a STATUS_VAR);
2493 STATUS(float_rounding_mode) = oldmode;
2494 return res;
2497 /*----------------------------------------------------------------------------
2498 | Returns the result of adding the absolute values of the double-precision
2499 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
2500 | before being returned. `zSign' is ignored if the result is a NaN.
2501 | The addition is performed according to the IEC/IEEE Standard for Binary
2502 | Floating-Point Arithmetic.
2503 *----------------------------------------------------------------------------*/
2505 static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2507 int16 aExp, bExp, zExp;
2508 bits64 aSig, bSig, zSig;
2509 int16 expDiff;
2511 aSig = extractFloat64Frac( a );
2512 aExp = extractFloat64Exp( a );
2513 bSig = extractFloat64Frac( b );
2514 bExp = extractFloat64Exp( b );
2515 expDiff = aExp - bExp;
2516 aSig <<= 9;
2517 bSig <<= 9;
2518 if ( 0 < expDiff ) {
2519 if ( aExp == 0x7FF ) {
2520 if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2521 return a;
2523 if ( bExp == 0 ) {
2524 --expDiff;
2526 else {
2527 bSig |= LIT64( 0x2000000000000000 );
2529 shift64RightJamming( bSig, expDiff, &bSig );
2530 zExp = aExp;
2532 else if ( expDiff < 0 ) {
2533 if ( bExp == 0x7FF ) {
2534 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2535 return packFloat64( zSign, 0x7FF, 0 );
2537 if ( aExp == 0 ) {
2538 ++expDiff;
2540 else {
2541 aSig |= LIT64( 0x2000000000000000 );
2543 shift64RightJamming( aSig, - expDiff, &aSig );
2544 zExp = bExp;
2546 else {
2547 if ( aExp == 0x7FF ) {
2548 if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2549 return a;
2551 if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
2552 zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
2553 zExp = aExp;
2554 goto roundAndPack;
2556 aSig |= LIT64( 0x2000000000000000 );
2557 zSig = ( aSig + bSig )<<1;
2558 --zExp;
2559 if ( (sbits64) zSig < 0 ) {
2560 zSig = aSig + bSig;
2561 ++zExp;
2563 roundAndPack:
2564 return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2568 /*----------------------------------------------------------------------------
2569 | Returns the result of subtracting the absolute values of the double-
2570 | precision floating-point values `a' and `b'. If `zSign' is 1, the
2571 | difference is negated before being returned. `zSign' is ignored if the
2572 | result is a NaN. The subtraction is performed according to the IEC/IEEE
2573 | Standard for Binary Floating-Point Arithmetic.
2574 *----------------------------------------------------------------------------*/
2576 static float64 subFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2578 int16 aExp, bExp, zExp;
2579 bits64 aSig, bSig, zSig;
2580 int16 expDiff;
2582 aSig = extractFloat64Frac( a );
2583 aExp = extractFloat64Exp( a );
2584 bSig = extractFloat64Frac( b );
2585 bExp = extractFloat64Exp( b );
2586 expDiff = aExp - bExp;
2587 aSig <<= 10;
2588 bSig <<= 10;
2589 if ( 0 < expDiff ) goto aExpBigger;
2590 if ( expDiff < 0 ) goto bExpBigger;
2591 if ( aExp == 0x7FF ) {
2592 if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2593 float_raise( float_flag_invalid STATUS_VAR);
2594 return float64_default_nan;
2596 if ( aExp == 0 ) {
2597 aExp = 1;
2598 bExp = 1;
2600 if ( bSig < aSig ) goto aBigger;
2601 if ( aSig < bSig ) goto bBigger;
2602 return packFloat64( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
2603 bExpBigger:
2604 if ( bExp == 0x7FF ) {
2605 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2606 return packFloat64( zSign ^ 1, 0x7FF, 0 );
2608 if ( aExp == 0 ) {
2609 ++expDiff;
2611 else {
2612 aSig |= LIT64( 0x4000000000000000 );
2614 shift64RightJamming( aSig, - expDiff, &aSig );
2615 bSig |= LIT64( 0x4000000000000000 );
2616 bBigger:
2617 zSig = bSig - aSig;
2618 zExp = bExp;
2619 zSign ^= 1;
2620 goto normalizeRoundAndPack;
2621 aExpBigger:
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( 0x4000000000000000 );
2632 shift64RightJamming( bSig, expDiff, &bSig );
2633 aSig |= LIT64( 0x4000000000000000 );
2634 aBigger:
2635 zSig = aSig - bSig;
2636 zExp = aExp;
2637 normalizeRoundAndPack:
2638 --zExp;
2639 return normalizeRoundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2643 /*----------------------------------------------------------------------------
2644 | Returns the result of adding the double-precision floating-point values `a'
2645 | and `b'. The operation is performed according to the IEC/IEEE Standard for
2646 | Binary Floating-Point Arithmetic.
2647 *----------------------------------------------------------------------------*/
2649 float64 float64_add( float64 a, float64 b STATUS_PARAM )
2651 flag aSign, bSign;
2653 aSign = extractFloat64Sign( a );
2654 bSign = extractFloat64Sign( b );
2655 if ( aSign == bSign ) {
2656 return addFloat64Sigs( a, b, aSign STATUS_VAR );
2658 else {
2659 return subFloat64Sigs( a, b, aSign STATUS_VAR );
2664 /*----------------------------------------------------------------------------
2665 | Returns the result of subtracting the double-precision floating-point values
2666 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2667 | for Binary Floating-Point Arithmetic.
2668 *----------------------------------------------------------------------------*/
2670 float64 float64_sub( float64 a, float64 b STATUS_PARAM )
2672 flag aSign, bSign;
2674 aSign = extractFloat64Sign( a );
2675 bSign = extractFloat64Sign( b );
2676 if ( aSign == bSign ) {
2677 return subFloat64Sigs( a, b, aSign STATUS_VAR );
2679 else {
2680 return addFloat64Sigs( a, b, aSign STATUS_VAR );
2685 /*----------------------------------------------------------------------------
2686 | Returns the result of multiplying the double-precision floating-point values
2687 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2688 | for Binary Floating-Point Arithmetic.
2689 *----------------------------------------------------------------------------*/
2691 float64 float64_mul( float64 a, float64 b STATUS_PARAM )
2693 flag aSign, bSign, zSign;
2694 int16 aExp, bExp, zExp;
2695 bits64 aSig, bSig, zSig0, zSig1;
2697 aSig = extractFloat64Frac( a );
2698 aExp = extractFloat64Exp( a );
2699 aSign = extractFloat64Sign( a );
2700 bSig = extractFloat64Frac( b );
2701 bExp = extractFloat64Exp( b );
2702 bSign = extractFloat64Sign( b );
2703 zSign = aSign ^ bSign;
2704 if ( aExp == 0x7FF ) {
2705 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2706 return propagateFloat64NaN( a, b STATUS_VAR );
2708 if ( ( bExp | bSig ) == 0 ) {
2709 float_raise( float_flag_invalid STATUS_VAR);
2710 return float64_default_nan;
2712 return packFloat64( zSign, 0x7FF, 0 );
2714 if ( bExp == 0x7FF ) {
2715 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2716 if ( ( aExp | aSig ) == 0 ) {
2717 float_raise( float_flag_invalid STATUS_VAR);
2718 return float64_default_nan;
2720 return packFloat64( zSign, 0x7FF, 0 );
2722 if ( aExp == 0 ) {
2723 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2724 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2726 if ( bExp == 0 ) {
2727 if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
2728 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2730 zExp = aExp + bExp - 0x3FF;
2731 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2732 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2733 mul64To128( aSig, bSig, &zSig0, &zSig1 );
2734 zSig0 |= ( zSig1 != 0 );
2735 if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
2736 zSig0 <<= 1;
2737 --zExp;
2739 return roundAndPackFloat64( zSign, zExp, zSig0 STATUS_VAR );
2743 /*----------------------------------------------------------------------------
2744 | Returns the result of dividing the double-precision floating-point value `a'
2745 | by the corresponding value `b'. The operation is performed according to
2746 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2747 *----------------------------------------------------------------------------*/
2749 float64 float64_div( float64 a, float64 b STATUS_PARAM )
2751 flag aSign, bSign, zSign;
2752 int16 aExp, bExp, zExp;
2753 bits64 aSig, bSig, zSig;
2754 bits64 rem0, rem1;
2755 bits64 term0, term1;
2757 aSig = extractFloat64Frac( a );
2758 aExp = extractFloat64Exp( a );
2759 aSign = extractFloat64Sign( a );
2760 bSig = extractFloat64Frac( b );
2761 bExp = extractFloat64Exp( b );
2762 bSign = extractFloat64Sign( b );
2763 zSign = aSign ^ bSign;
2764 if ( aExp == 0x7FF ) {
2765 if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2766 if ( bExp == 0x7FF ) {
2767 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2768 float_raise( float_flag_invalid STATUS_VAR);
2769 return float64_default_nan;
2771 return packFloat64( zSign, 0x7FF, 0 );
2773 if ( bExp == 0x7FF ) {
2774 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2775 return packFloat64( zSign, 0, 0 );
2777 if ( bExp == 0 ) {
2778 if ( bSig == 0 ) {
2779 if ( ( aExp | aSig ) == 0 ) {
2780 float_raise( float_flag_invalid STATUS_VAR);
2781 return float64_default_nan;
2783 float_raise( float_flag_divbyzero STATUS_VAR);
2784 return packFloat64( zSign, 0x7FF, 0 );
2786 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2788 if ( aExp == 0 ) {
2789 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2790 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2792 zExp = aExp - bExp + 0x3FD;
2793 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2794 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2795 if ( bSig <= ( aSig + aSig ) ) {
2796 aSig >>= 1;
2797 ++zExp;
2799 zSig = estimateDiv128To64( aSig, 0, bSig );
2800 if ( ( zSig & 0x1FF ) <= 2 ) {
2801 mul64To128( bSig, zSig, &term0, &term1 );
2802 sub128( aSig, 0, term0, term1, &rem0, &rem1 );
2803 while ( (sbits64) rem0 < 0 ) {
2804 --zSig;
2805 add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
2807 zSig |= ( rem1 != 0 );
2809 return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2813 /*----------------------------------------------------------------------------
2814 | Returns the remainder of the double-precision floating-point value `a'
2815 | with respect to the corresponding value `b'. The operation is performed
2816 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2817 *----------------------------------------------------------------------------*/
2819 float64 float64_rem( float64 a, float64 b STATUS_PARAM )
2821 flag aSign, bSign, zSign;
2822 int16 aExp, bExp, expDiff;
2823 bits64 aSig, bSig;
2824 bits64 q, alternateASig;
2825 sbits64 sigMean;
2827 aSig = extractFloat64Frac( a );
2828 aExp = extractFloat64Exp( a );
2829 aSign = extractFloat64Sign( a );
2830 bSig = extractFloat64Frac( b );
2831 bExp = extractFloat64Exp( b );
2832 bSign = extractFloat64Sign( b );
2833 if ( aExp == 0x7FF ) {
2834 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2835 return propagateFloat64NaN( a, b STATUS_VAR );
2837 float_raise( float_flag_invalid STATUS_VAR);
2838 return float64_default_nan;
2840 if ( bExp == 0x7FF ) {
2841 if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2842 return a;
2844 if ( bExp == 0 ) {
2845 if ( bSig == 0 ) {
2846 float_raise( float_flag_invalid STATUS_VAR);
2847 return float64_default_nan;
2849 normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2851 if ( aExp == 0 ) {
2852 if ( aSig == 0 ) return a;
2853 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2855 expDiff = aExp - bExp;
2856 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
2857 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2858 if ( expDiff < 0 ) {
2859 if ( expDiff < -1 ) return a;
2860 aSig >>= 1;
2862 q = ( bSig <= aSig );
2863 if ( q ) aSig -= bSig;
2864 expDiff -= 64;
2865 while ( 0 < expDiff ) {
2866 q = estimateDiv128To64( aSig, 0, bSig );
2867 q = ( 2 < q ) ? q - 2 : 0;
2868 aSig = - ( ( bSig>>2 ) * q );
2869 expDiff -= 62;
2871 expDiff += 64;
2872 if ( 0 < expDiff ) {
2873 q = estimateDiv128To64( aSig, 0, bSig );
2874 q = ( 2 < q ) ? q - 2 : 0;
2875 q >>= 64 - expDiff;
2876 bSig >>= 2;
2877 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
2879 else {
2880 aSig >>= 2;
2881 bSig >>= 2;
2883 do {
2884 alternateASig = aSig;
2885 ++q;
2886 aSig -= bSig;
2887 } while ( 0 <= (sbits64) aSig );
2888 sigMean = aSig + alternateASig;
2889 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
2890 aSig = alternateASig;
2892 zSign = ( (sbits64) aSig < 0 );
2893 if ( zSign ) aSig = - aSig;
2894 return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig STATUS_VAR );
2898 /*----------------------------------------------------------------------------
2899 | Returns the square root of the double-precision floating-point value `a'.
2900 | The operation is performed according to the IEC/IEEE Standard for Binary
2901 | Floating-Point Arithmetic.
2902 *----------------------------------------------------------------------------*/
2904 float64 float64_sqrt( float64 a STATUS_PARAM )
2906 flag aSign;
2907 int16 aExp, zExp;
2908 bits64 aSig, zSig, doubleZSig;
2909 bits64 rem0, rem1, term0, term1;
2911 aSig = extractFloat64Frac( a );
2912 aExp = extractFloat64Exp( a );
2913 aSign = extractFloat64Sign( a );
2914 if ( aExp == 0x7FF ) {
2915 if ( aSig ) return propagateFloat64NaN( a, a STATUS_VAR );
2916 if ( ! aSign ) return a;
2917 float_raise( float_flag_invalid STATUS_VAR);
2918 return float64_default_nan;
2920 if ( aSign ) {
2921 if ( ( aExp | aSig ) == 0 ) return a;
2922 float_raise( float_flag_invalid STATUS_VAR);
2923 return float64_default_nan;
2925 if ( aExp == 0 ) {
2926 if ( aSig == 0 ) return 0;
2927 normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2929 zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
2930 aSig |= LIT64( 0x0010000000000000 );
2931 zSig = estimateSqrt32( aExp, aSig>>21 );
2932 aSig <<= 9 - ( aExp & 1 );
2933 zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
2934 if ( ( zSig & 0x1FF ) <= 5 ) {
2935 doubleZSig = zSig<<1;
2936 mul64To128( zSig, zSig, &term0, &term1 );
2937 sub128( aSig, 0, term0, term1, &rem0, &rem1 );
2938 while ( (sbits64) rem0 < 0 ) {
2939 --zSig;
2940 doubleZSig -= 2;
2941 add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
2943 zSig |= ( ( rem0 | rem1 ) != 0 );
2945 return roundAndPackFloat64( 0, zExp, zSig STATUS_VAR );
2949 /*----------------------------------------------------------------------------
2950 | Returns 1 if the double-precision floating-point value `a' is equal to the
2951 | corresponding value `b', and 0 otherwise. The comparison is performed
2952 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2953 *----------------------------------------------------------------------------*/
2955 int float64_eq( float64 a, float64 b STATUS_PARAM )
2958 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
2959 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
2960 ) {
2961 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
2962 float_raise( float_flag_invalid STATUS_VAR);
2964 return 0;
2966 return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
2970 /*----------------------------------------------------------------------------
2971 | Returns 1 if the double-precision floating-point value `a' is less than or
2972 | equal to the corresponding value `b', and 0 otherwise. The comparison is
2973 | performed according to the IEC/IEEE Standard for Binary Floating-Point
2974 | Arithmetic.
2975 *----------------------------------------------------------------------------*/
2977 int float64_le( float64 a, float64 b STATUS_PARAM )
2979 flag aSign, bSign;
2981 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
2982 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
2983 ) {
2984 float_raise( float_flag_invalid STATUS_VAR);
2985 return 0;
2987 aSign = extractFloat64Sign( a );
2988 bSign = extractFloat64Sign( b );
2989 if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
2990 return ( a == b ) || ( aSign ^ ( a < b ) );
2994 /*----------------------------------------------------------------------------
2995 | Returns 1 if the double-precision floating-point value `a' is less than
2996 | the corresponding value `b', and 0 otherwise. The comparison is performed
2997 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2998 *----------------------------------------------------------------------------*/
3000 int float64_lt( float64 a, float64 b STATUS_PARAM )
3002 flag aSign, bSign;
3004 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3005 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3006 ) {
3007 float_raise( float_flag_invalid STATUS_VAR);
3008 return 0;
3010 aSign = extractFloat64Sign( a );
3011 bSign = extractFloat64Sign( b );
3012 if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
3013 return ( a != b ) && ( aSign ^ ( a < b ) );
3017 /*----------------------------------------------------------------------------
3018 | Returns 1 if the double-precision floating-point value `a' is equal to the
3019 | corresponding value `b', and 0 otherwise. The invalid exception is raised
3020 | if either operand is a NaN. Otherwise, the comparison is performed
3021 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3022 *----------------------------------------------------------------------------*/
3024 int float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
3027 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3028 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3029 ) {
3030 float_raise( float_flag_invalid STATUS_VAR);
3031 return 0;
3033 return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
3037 /*----------------------------------------------------------------------------
3038 | Returns 1 if the double-precision floating-point value `a' is less than or
3039 | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
3040 | cause an exception. Otherwise, the comparison is performed according to the
3041 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3042 *----------------------------------------------------------------------------*/
3044 int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
3046 flag aSign, bSign;
3048 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3049 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3050 ) {
3051 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3052 float_raise( float_flag_invalid STATUS_VAR);
3054 return 0;
3056 aSign = extractFloat64Sign( a );
3057 bSign = extractFloat64Sign( b );
3058 if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
3059 return ( a == b ) || ( aSign ^ ( a < b ) );
3063 /*----------------------------------------------------------------------------
3064 | Returns 1 if the double-precision floating-point value `a' is less than
3065 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
3066 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
3067 | Standard for Binary Floating-Point Arithmetic.
3068 *----------------------------------------------------------------------------*/
3070 int float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
3072 flag aSign, bSign;
3074 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3075 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3076 ) {
3077 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3078 float_raise( float_flag_invalid STATUS_VAR);
3080 return 0;
3082 aSign = extractFloat64Sign( a );
3083 bSign = extractFloat64Sign( b );
3084 if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
3085 return ( a != b ) && ( aSign ^ ( a < b ) );
3089 #ifdef FLOATX80
3091 /*----------------------------------------------------------------------------
3092 | Returns the result of converting the extended double-precision floating-
3093 | point value `a' to the 32-bit two's complement integer format. The
3094 | conversion is performed according to the IEC/IEEE Standard for Binary
3095 | Floating-Point Arithmetic---which means in particular that the conversion
3096 | is rounded according to the current rounding mode. If `a' is a NaN, the
3097 | largest positive integer is returned. Otherwise, if the conversion
3098 | overflows, the largest integer with the same sign as `a' is returned.
3099 *----------------------------------------------------------------------------*/
3101 int32 floatx80_to_int32( floatx80 a STATUS_PARAM )
3103 flag aSign;
3104 int32 aExp, shiftCount;
3105 bits64 aSig;
3107 aSig = extractFloatx80Frac( a );
3108 aExp = extractFloatx80Exp( a );
3109 aSign = extractFloatx80Sign( a );
3110 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3111 shiftCount = 0x4037 - aExp;
3112 if ( shiftCount <= 0 ) shiftCount = 1;
3113 shift64RightJamming( aSig, shiftCount, &aSig );
3114 return roundAndPackInt32( aSign, aSig STATUS_VAR );
3118 /*----------------------------------------------------------------------------
3119 | Returns the result of converting the extended double-precision floating-
3120 | point value `a' to the 32-bit two's complement integer format. The
3121 | conversion is performed according to the IEC/IEEE Standard for Binary
3122 | Floating-Point Arithmetic, except that the conversion is always rounded
3123 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3124 | Otherwise, if the conversion overflows, the largest integer with the same
3125 | sign as `a' is returned.
3126 *----------------------------------------------------------------------------*/
3128 int32 floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM )
3130 flag aSign;
3131 int32 aExp, shiftCount;
3132 bits64 aSig, savedASig;
3133 int32 z;
3135 aSig = extractFloatx80Frac( a );
3136 aExp = extractFloatx80Exp( a );
3137 aSign = extractFloatx80Sign( a );
3138 if ( 0x401E < aExp ) {
3139 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3140 goto invalid;
3142 else if ( aExp < 0x3FFF ) {
3143 if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3144 return 0;
3146 shiftCount = 0x403E - aExp;
3147 savedASig = aSig;
3148 aSig >>= shiftCount;
3149 z = aSig;
3150 if ( aSign ) z = - z;
3151 if ( ( z < 0 ) ^ aSign ) {
3152 invalid:
3153 float_raise( float_flag_invalid STATUS_VAR);
3154 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
3156 if ( ( aSig<<shiftCount ) != savedASig ) {
3157 STATUS(float_exception_flags) |= float_flag_inexact;
3159 return z;
3163 /*----------------------------------------------------------------------------
3164 | Returns the result of converting the extended double-precision floating-
3165 | point value `a' to the 64-bit two's complement integer format. The
3166 | conversion is performed according to the IEC/IEEE Standard for Binary
3167 | Floating-Point Arithmetic---which means in particular that the conversion
3168 | is rounded according to the current rounding mode. If `a' is a NaN,
3169 | the largest positive integer is returned. Otherwise, if the conversion
3170 | overflows, the largest integer with the same sign as `a' is returned.
3171 *----------------------------------------------------------------------------*/
3173 int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
3175 flag aSign;
3176 int32 aExp, shiftCount;
3177 bits64 aSig, aSigExtra;
3179 aSig = extractFloatx80Frac( a );
3180 aExp = extractFloatx80Exp( a );
3181 aSign = extractFloatx80Sign( a );
3182 shiftCount = 0x403E - aExp;
3183 if ( shiftCount <= 0 ) {
3184 if ( shiftCount ) {
3185 float_raise( float_flag_invalid STATUS_VAR);
3186 if ( ! aSign
3187 || ( ( aExp == 0x7FFF )
3188 && ( aSig != LIT64( 0x8000000000000000 ) ) )
3189 ) {
3190 return LIT64( 0x7FFFFFFFFFFFFFFF );
3192 return (sbits64) LIT64( 0x8000000000000000 );
3194 aSigExtra = 0;
3196 else {
3197 shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
3199 return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
3203 /*----------------------------------------------------------------------------
3204 | Returns the result of converting the extended double-precision floating-
3205 | point value `a' to the 64-bit two's complement integer format. The
3206 | conversion is performed according to the IEC/IEEE Standard for Binary
3207 | Floating-Point Arithmetic, except that the conversion is always rounded
3208 | toward zero. If `a' is a NaN, the largest positive integer is returned.
3209 | Otherwise, if the conversion overflows, the largest integer with the same
3210 | sign as `a' is returned.
3211 *----------------------------------------------------------------------------*/
3213 int64 floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM )
3215 flag aSign;
3216 int32 aExp, shiftCount;
3217 bits64 aSig;
3218 int64 z;
3220 aSig = extractFloatx80Frac( a );
3221 aExp = extractFloatx80Exp( a );
3222 aSign = extractFloatx80Sign( a );
3223 shiftCount = aExp - 0x403E;
3224 if ( 0 <= shiftCount ) {
3225 aSig &= LIT64( 0x7FFFFFFFFFFFFFFF );
3226 if ( ( a.high != 0xC03E ) || aSig ) {
3227 float_raise( float_flag_invalid STATUS_VAR);
3228 if ( ! aSign || ( ( aExp == 0x7FFF ) && aSig ) ) {
3229 return LIT64( 0x7FFFFFFFFFFFFFFF );
3232 return (sbits64) LIT64( 0x8000000000000000 );
3234 else if ( aExp < 0x3FFF ) {
3235 if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3236 return 0;
3238 z = aSig>>( - shiftCount );
3239 if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
3240 STATUS(float_exception_flags) |= float_flag_inexact;
3242 if ( aSign ) z = - z;
3243 return z;
3247 /*----------------------------------------------------------------------------
3248 | Returns the result of converting the extended double-precision floating-
3249 | point value `a' to the single-precision floating-point format. The
3250 | conversion is performed according to the IEC/IEEE Standard for Binary
3251 | Floating-Point Arithmetic.
3252 *----------------------------------------------------------------------------*/
3254 float32 floatx80_to_float32( floatx80 a STATUS_PARAM )
3256 flag aSign;
3257 int32 aExp;
3258 bits64 aSig;
3260 aSig = extractFloatx80Frac( a );
3261 aExp = extractFloatx80Exp( a );
3262 aSign = extractFloatx80Sign( a );
3263 if ( aExp == 0x7FFF ) {
3264 if ( (bits64) ( aSig<<1 ) ) {
3265 return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) );
3267 return packFloat32( aSign, 0xFF, 0 );
3269 shift64RightJamming( aSig, 33, &aSig );
3270 if ( aExp || aSig ) aExp -= 0x3F81;
3271 return roundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
3275 /*----------------------------------------------------------------------------
3276 | Returns the result of converting the extended double-precision floating-
3277 | point value `a' to the double-precision floating-point format. The
3278 | conversion is performed according to the IEC/IEEE Standard for Binary
3279 | Floating-Point Arithmetic.
3280 *----------------------------------------------------------------------------*/
3282 float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
3284 flag aSign;
3285 int32 aExp;
3286 bits64 aSig, zSig;
3288 aSig = extractFloatx80Frac( a );
3289 aExp = extractFloatx80Exp( a );
3290 aSign = extractFloatx80Sign( a );
3291 if ( aExp == 0x7FFF ) {
3292 if ( (bits64) ( aSig<<1 ) ) {
3293 return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) );
3295 return packFloat64( aSign, 0x7FF, 0 );
3297 shift64RightJamming( aSig, 1, &zSig );
3298 if ( aExp || aSig ) aExp -= 0x3C01;
3299 return roundAndPackFloat64( aSign, aExp, zSig STATUS_VAR );
3303 #ifdef FLOAT128
3305 /*----------------------------------------------------------------------------
3306 | Returns the result of converting the extended double-precision floating-
3307 | point value `a' to the quadruple-precision floating-point format. The
3308 | conversion is performed according to the IEC/IEEE Standard for Binary
3309 | Floating-Point Arithmetic.
3310 *----------------------------------------------------------------------------*/
3312 float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
3314 flag aSign;
3315 int16 aExp;
3316 bits64 aSig, zSig0, zSig1;
3318 aSig = extractFloatx80Frac( a );
3319 aExp = extractFloatx80Exp( a );
3320 aSign = extractFloatx80Sign( a );
3321 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
3322 return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) );
3324 shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
3325 return packFloat128( aSign, aExp, zSig0, zSig1 );
3329 #endif
3331 /*----------------------------------------------------------------------------
3332 | Rounds the extended double-precision floating-point value `a' to an integer,
3333 | and returns the result as an extended quadruple-precision floating-point
3334 | value. The operation is performed according to the IEC/IEEE Standard for
3335 | Binary Floating-Point Arithmetic.
3336 *----------------------------------------------------------------------------*/
3338 floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM )
3340 flag aSign;
3341 int32 aExp;
3342 bits64 lastBitMask, roundBitsMask;
3343 int8 roundingMode;
3344 floatx80 z;
3346 aExp = extractFloatx80Exp( a );
3347 if ( 0x403E <= aExp ) {
3348 if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
3349 return propagateFloatx80NaN( a, a STATUS_VAR );
3351 return a;
3353 if ( aExp < 0x3FFF ) {
3354 if ( ( aExp == 0 )
3355 && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
3356 return a;
3358 STATUS(float_exception_flags) |= float_flag_inexact;
3359 aSign = extractFloatx80Sign( a );
3360 switch ( STATUS(float_rounding_mode) ) {
3361 case float_round_nearest_even:
3362 if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
3363 ) {
3364 return
3365 packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3367 break;
3368 case float_round_down:
3369 return
3370 aSign ?
3371 packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
3372 : packFloatx80( 0, 0, 0 );
3373 case float_round_up:
3374 return
3375 aSign ? packFloatx80( 1, 0, 0 )
3376 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
3378 return packFloatx80( aSign, 0, 0 );
3380 lastBitMask = 1;
3381 lastBitMask <<= 0x403E - aExp;
3382 roundBitsMask = lastBitMask - 1;
3383 z = a;
3384 roundingMode = STATUS(float_rounding_mode);
3385 if ( roundingMode == float_round_nearest_even ) {
3386 z.low += lastBitMask>>1;
3387 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
3389 else if ( roundingMode != float_round_to_zero ) {
3390 if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
3391 z.low += roundBitsMask;
3394 z.low &= ~ roundBitsMask;
3395 if ( z.low == 0 ) {
3396 ++z.high;
3397 z.low = LIT64( 0x8000000000000000 );
3399 if ( z.low != a.low ) STATUS(float_exception_flags) |= float_flag_inexact;
3400 return z;
3404 /*----------------------------------------------------------------------------
3405 | Returns the result of adding the absolute values of the extended double-
3406 | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
3407 | negated before being returned. `zSign' is ignored if the result is a NaN.
3408 | The addition is performed according to the IEC/IEEE Standard for Binary
3409 | Floating-Point Arithmetic.
3410 *----------------------------------------------------------------------------*/
3412 static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM)
3414 int32 aExp, bExp, zExp;
3415 bits64 aSig, bSig, zSig0, zSig1;
3416 int32 expDiff;
3418 aSig = extractFloatx80Frac( a );
3419 aExp = extractFloatx80Exp( a );
3420 bSig = extractFloatx80Frac( b );
3421 bExp = extractFloatx80Exp( b );
3422 expDiff = aExp - bExp;
3423 if ( 0 < expDiff ) {
3424 if ( aExp == 0x7FFF ) {
3425 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3426 return a;
3428 if ( bExp == 0 ) --expDiff;
3429 shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3430 zExp = aExp;
3432 else if ( expDiff < 0 ) {
3433 if ( bExp == 0x7FFF ) {
3434 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3435 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3437 if ( aExp == 0 ) ++expDiff;
3438 shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3439 zExp = bExp;
3441 else {
3442 if ( aExp == 0x7FFF ) {
3443 if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3444 return propagateFloatx80NaN( a, b STATUS_VAR );
3446 return a;
3448 zSig1 = 0;
3449 zSig0 = aSig + bSig;
3450 if ( aExp == 0 ) {
3451 normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
3452 goto roundAndPack;
3454 zExp = aExp;
3455 goto shiftRight1;
3457 zSig0 = aSig + bSig;
3458 if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
3459 shiftRight1:
3460 shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
3461 zSig0 |= LIT64( 0x8000000000000000 );
3462 ++zExp;
3463 roundAndPack:
3464 return
3465 roundAndPackFloatx80(
3466 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3470 /*----------------------------------------------------------------------------
3471 | Returns the result of subtracting the absolute values of the extended
3472 | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
3473 | difference is negated before being returned. `zSign' is ignored if the
3474 | result is a NaN. The subtraction is performed according to the IEC/IEEE
3475 | Standard for Binary Floating-Point Arithmetic.
3476 *----------------------------------------------------------------------------*/
3478 static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM )
3480 int32 aExp, bExp, zExp;
3481 bits64 aSig, bSig, zSig0, zSig1;
3482 int32 expDiff;
3483 floatx80 z;
3485 aSig = extractFloatx80Frac( a );
3486 aExp = extractFloatx80Exp( a );
3487 bSig = extractFloatx80Frac( b );
3488 bExp = extractFloatx80Exp( b );
3489 expDiff = aExp - bExp;
3490 if ( 0 < expDiff ) goto aExpBigger;
3491 if ( expDiff < 0 ) goto bExpBigger;
3492 if ( aExp == 0x7FFF ) {
3493 if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3494 return propagateFloatx80NaN( a, b STATUS_VAR );
3496 float_raise( float_flag_invalid STATUS_VAR);
3497 z.low = floatx80_default_nan_low;
3498 z.high = floatx80_default_nan_high;
3499 return z;
3501 if ( aExp == 0 ) {
3502 aExp = 1;
3503 bExp = 1;
3505 zSig1 = 0;
3506 if ( bSig < aSig ) goto aBigger;
3507 if ( aSig < bSig ) goto bBigger;
3508 return packFloatx80( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
3509 bExpBigger:
3510 if ( bExp == 0x7FFF ) {
3511 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3512 return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
3514 if ( aExp == 0 ) ++expDiff;
3515 shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3516 bBigger:
3517 sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
3518 zExp = bExp;
3519 zSign ^= 1;
3520 goto normalizeRoundAndPack;
3521 aExpBigger:
3522 if ( aExp == 0x7FFF ) {
3523 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3524 return a;
3526 if ( bExp == 0 ) --expDiff;
3527 shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3528 aBigger:
3529 sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
3530 zExp = aExp;
3531 normalizeRoundAndPack:
3532 return
3533 normalizeRoundAndPackFloatx80(
3534 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3538 /*----------------------------------------------------------------------------
3539 | Returns the result of adding the extended double-precision floating-point
3540 | values `a' and `b'. The operation is performed according to the IEC/IEEE
3541 | Standard for Binary Floating-Point Arithmetic.
3542 *----------------------------------------------------------------------------*/
3544 floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM )
3546 flag aSign, bSign;
3548 aSign = extractFloatx80Sign( a );
3549 bSign = extractFloatx80Sign( b );
3550 if ( aSign == bSign ) {
3551 return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3553 else {
3554 return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3559 /*----------------------------------------------------------------------------
3560 | Returns the result of subtracting the extended double-precision floating-
3561 | point values `a' and `b'. The operation is performed according to the
3562 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3563 *----------------------------------------------------------------------------*/
3565 floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM )
3567 flag aSign, bSign;
3569 aSign = extractFloatx80Sign( a );
3570 bSign = extractFloatx80Sign( b );
3571 if ( aSign == bSign ) {
3572 return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3574 else {
3575 return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3580 /*----------------------------------------------------------------------------
3581 | Returns the result of multiplying the extended double-precision floating-
3582 | point values `a' and `b'. The operation is performed according to the
3583 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3584 *----------------------------------------------------------------------------*/
3586 floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
3588 flag aSign, bSign, zSign;
3589 int32 aExp, bExp, zExp;
3590 bits64 aSig, bSig, zSig0, zSig1;
3591 floatx80 z;
3593 aSig = extractFloatx80Frac( a );
3594 aExp = extractFloatx80Exp( a );
3595 aSign = extractFloatx80Sign( a );
3596 bSig = extractFloatx80Frac( b );
3597 bExp = extractFloatx80Exp( b );
3598 bSign = extractFloatx80Sign( b );
3599 zSign = aSign ^ bSign;
3600 if ( aExp == 0x7FFF ) {
3601 if ( (bits64) ( aSig<<1 )
3602 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3603 return propagateFloatx80NaN( a, b STATUS_VAR );
3605 if ( ( bExp | bSig ) == 0 ) goto invalid;
3606 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3608 if ( bExp == 0x7FFF ) {
3609 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3610 if ( ( aExp | aSig ) == 0 ) {
3611 invalid:
3612 float_raise( float_flag_invalid STATUS_VAR);
3613 z.low = floatx80_default_nan_low;
3614 z.high = floatx80_default_nan_high;
3615 return z;
3617 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3619 if ( aExp == 0 ) {
3620 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3621 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3623 if ( bExp == 0 ) {
3624 if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
3625 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3627 zExp = aExp + bExp - 0x3FFE;
3628 mul64To128( aSig, bSig, &zSig0, &zSig1 );
3629 if ( 0 < (sbits64) zSig0 ) {
3630 shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
3631 --zExp;
3633 return
3634 roundAndPackFloatx80(
3635 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3639 /*----------------------------------------------------------------------------
3640 | Returns the result of dividing the extended double-precision floating-point
3641 | value `a' by the corresponding value `b'. The operation is performed
3642 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3643 *----------------------------------------------------------------------------*/
3645 floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
3647 flag aSign, bSign, zSign;
3648 int32 aExp, bExp, zExp;
3649 bits64 aSig, bSig, zSig0, zSig1;
3650 bits64 rem0, rem1, rem2, term0, term1, term2;
3651 floatx80 z;
3653 aSig = extractFloatx80Frac( a );
3654 aExp = extractFloatx80Exp( a );
3655 aSign = extractFloatx80Sign( a );
3656 bSig = extractFloatx80Frac( b );
3657 bExp = extractFloatx80Exp( b );
3658 bSign = extractFloatx80Sign( b );
3659 zSign = aSign ^ bSign;
3660 if ( aExp == 0x7FFF ) {
3661 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3662 if ( bExp == 0x7FFF ) {
3663 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3664 goto invalid;
3666 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3668 if ( bExp == 0x7FFF ) {
3669 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3670 return packFloatx80( zSign, 0, 0 );
3672 if ( bExp == 0 ) {
3673 if ( bSig == 0 ) {
3674 if ( ( aExp | aSig ) == 0 ) {
3675 invalid:
3676 float_raise( float_flag_invalid STATUS_VAR);
3677 z.low = floatx80_default_nan_low;
3678 z.high = floatx80_default_nan_high;
3679 return z;
3681 float_raise( float_flag_divbyzero STATUS_VAR);
3682 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3684 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3686 if ( aExp == 0 ) {
3687 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3688 normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3690 zExp = aExp - bExp + 0x3FFE;
3691 rem1 = 0;
3692 if ( bSig <= aSig ) {
3693 shift128Right( aSig, 0, 1, &aSig, &rem1 );
3694 ++zExp;
3696 zSig0 = estimateDiv128To64( aSig, rem1, bSig );
3697 mul64To128( bSig, zSig0, &term0, &term1 );
3698 sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
3699 while ( (sbits64) rem0 < 0 ) {
3700 --zSig0;
3701 add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
3703 zSig1 = estimateDiv128To64( rem1, 0, bSig );
3704 if ( (bits64) ( zSig1<<1 ) <= 8 ) {
3705 mul64To128( bSig, zSig1, &term1, &term2 );
3706 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3707 while ( (sbits64) rem1 < 0 ) {
3708 --zSig1;
3709 add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
3711 zSig1 |= ( ( rem1 | rem2 ) != 0 );
3713 return
3714 roundAndPackFloatx80(
3715 STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3719 /*----------------------------------------------------------------------------
3720 | Returns the remainder of the extended double-precision floating-point value
3721 | `a' with respect to the corresponding value `b'. The operation is performed
3722 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3723 *----------------------------------------------------------------------------*/
3725 floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
3727 flag aSign, bSign, zSign;
3728 int32 aExp, bExp, expDiff;
3729 bits64 aSig0, aSig1, bSig;
3730 bits64 q, term0, term1, alternateASig0, alternateASig1;
3731 floatx80 z;
3733 aSig0 = extractFloatx80Frac( a );
3734 aExp = extractFloatx80Exp( a );
3735 aSign = extractFloatx80Sign( a );
3736 bSig = extractFloatx80Frac( b );
3737 bExp = extractFloatx80Exp( b );
3738 bSign = extractFloatx80Sign( b );
3739 if ( aExp == 0x7FFF ) {
3740 if ( (bits64) ( aSig0<<1 )
3741 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3742 return propagateFloatx80NaN( a, b STATUS_VAR );
3744 goto invalid;
3746 if ( bExp == 0x7FFF ) {
3747 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3748 return a;
3750 if ( bExp == 0 ) {
3751 if ( bSig == 0 ) {
3752 invalid:
3753 float_raise( float_flag_invalid STATUS_VAR);
3754 z.low = floatx80_default_nan_low;
3755 z.high = floatx80_default_nan_high;
3756 return z;
3758 normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3760 if ( aExp == 0 ) {
3761 if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
3762 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3764 bSig |= LIT64( 0x8000000000000000 );
3765 zSign = aSign;
3766 expDiff = aExp - bExp;
3767 aSig1 = 0;
3768 if ( expDiff < 0 ) {
3769 if ( expDiff < -1 ) return a;
3770 shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
3771 expDiff = 0;
3773 q = ( bSig <= aSig0 );
3774 if ( q ) aSig0 -= bSig;
3775 expDiff -= 64;
3776 while ( 0 < expDiff ) {
3777 q = estimateDiv128To64( aSig0, aSig1, bSig );
3778 q = ( 2 < q ) ? q - 2 : 0;
3779 mul64To128( bSig, q, &term0, &term1 );
3780 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3781 shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
3782 expDiff -= 62;
3784 expDiff += 64;
3785 if ( 0 < expDiff ) {
3786 q = estimateDiv128To64( aSig0, aSig1, bSig );
3787 q = ( 2 < q ) ? q - 2 : 0;
3788 q >>= 64 - expDiff;
3789 mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
3790 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3791 shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
3792 while ( le128( term0, term1, aSig0, aSig1 ) ) {
3793 ++q;
3794 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3797 else {
3798 term1 = 0;
3799 term0 = bSig;
3801 sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
3802 if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
3803 || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
3804 && ( q & 1 ) )
3805 ) {
3806 aSig0 = alternateASig0;
3807 aSig1 = alternateASig1;
3808 zSign = ! zSign;
3810 return
3811 normalizeRoundAndPackFloatx80(
3812 80, zSign, bExp + expDiff, aSig0, aSig1 STATUS_VAR );
3816 /*----------------------------------------------------------------------------
3817 | Returns the square root of the extended double-precision floating-point
3818 | value `a'. The operation is performed according to the IEC/IEEE Standard
3819 | for Binary Floating-Point Arithmetic.
3820 *----------------------------------------------------------------------------*/
3822 floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
3824 flag aSign;
3825 int32 aExp, zExp;
3826 bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
3827 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
3828 floatx80 z;
3830 aSig0 = extractFloatx80Frac( a );
3831 aExp = extractFloatx80Exp( a );
3832 aSign = extractFloatx80Sign( a );
3833 if ( aExp == 0x7FFF ) {
3834 if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a STATUS_VAR );
3835 if ( ! aSign ) return a;
3836 goto invalid;
3838 if ( aSign ) {
3839 if ( ( aExp | aSig0 ) == 0 ) return a;
3840 invalid:
3841 float_raise( float_flag_invalid STATUS_VAR);
3842 z.low = floatx80_default_nan_low;
3843 z.high = floatx80_default_nan_high;
3844 return z;
3846 if ( aExp == 0 ) {
3847 if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
3848 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3850 zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
3851 zSig0 = estimateSqrt32( aExp, aSig0>>32 );
3852 shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
3853 zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
3854 doubleZSig0 = zSig0<<1;
3855 mul64To128( zSig0, zSig0, &term0, &term1 );
3856 sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
3857 while ( (sbits64) rem0 < 0 ) {
3858 --zSig0;
3859 doubleZSig0 -= 2;
3860 add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
3862 zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
3863 if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
3864 if ( zSig1 == 0 ) zSig1 = 1;
3865 mul64To128( doubleZSig0, zSig1, &term1, &term2 );
3866 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3867 mul64To128( zSig1, zSig1, &term2, &term3 );
3868 sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
3869 while ( (sbits64) rem1 < 0 ) {
3870 --zSig1;
3871 shortShift128Left( 0, zSig1, 1, &term2, &term3 );
3872 term3 |= 1;
3873 term2 |= doubleZSig0;
3874 add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
3876 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
3878 shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
3879 zSig0 |= doubleZSig0;
3880 return
3881 roundAndPackFloatx80(
3882 STATUS(floatx80_rounding_precision), 0, zExp, zSig0, zSig1 STATUS_VAR );
3886 /*----------------------------------------------------------------------------
3887 | Returns 1 if the extended double-precision floating-point value `a' is
3888 | equal to the corresponding value `b', and 0 otherwise. The comparison is
3889 | performed according to the IEC/IEEE Standard for Binary Floating-Point
3890 | Arithmetic.
3891 *----------------------------------------------------------------------------*/
3893 int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
3896 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3897 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3898 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3899 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3900 ) {
3901 if ( floatx80_is_signaling_nan( a )
3902 || floatx80_is_signaling_nan( b ) ) {
3903 float_raise( float_flag_invalid STATUS_VAR);
3905 return 0;
3907 return
3908 ( a.low == b.low )
3909 && ( ( a.high == b.high )
3910 || ( ( a.low == 0 )
3911 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
3912 );
3916 /*----------------------------------------------------------------------------
3917 | Returns 1 if the extended double-precision floating-point value `a' is
3918 | less than or equal to the corresponding value `b', and 0 otherwise. The
3919 | comparison is performed according to the IEC/IEEE Standard for Binary
3920 | Floating-Point Arithmetic.
3921 *----------------------------------------------------------------------------*/
3923 int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
3925 flag aSign, bSign;
3927 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3928 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3929 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3930 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3931 ) {
3932 float_raise( float_flag_invalid STATUS_VAR);
3933 return 0;
3935 aSign = extractFloatx80Sign( a );
3936 bSign = extractFloatx80Sign( b );
3937 if ( aSign != bSign ) {
3938 return
3939 aSign
3940 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3941 == 0 );
3943 return
3944 aSign ? le128( b.high, b.low, a.high, a.low )
3945 : le128( a.high, a.low, b.high, b.low );
3949 /*----------------------------------------------------------------------------
3950 | Returns 1 if the extended double-precision floating-point value `a' is
3951 | less than the corresponding value `b', and 0 otherwise. The comparison
3952 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
3953 | Arithmetic.
3954 *----------------------------------------------------------------------------*/
3956 int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
3958 flag aSign, bSign;
3960 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3961 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3962 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3963 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3964 ) {
3965 float_raise( float_flag_invalid STATUS_VAR);
3966 return 0;
3968 aSign = extractFloatx80Sign( a );
3969 bSign = extractFloatx80Sign( b );
3970 if ( aSign != bSign ) {
3971 return
3972 aSign
3973 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3974 != 0 );
3976 return
3977 aSign ? lt128( b.high, b.low, a.high, a.low )
3978 : lt128( a.high, a.low, b.high, b.low );
3982 /*----------------------------------------------------------------------------
3983 | Returns 1 if the extended double-precision floating-point value `a' is equal
3984 | to the corresponding value `b', and 0 otherwise. The invalid exception is
3985 | raised if either operand is a NaN. Otherwise, the comparison is performed
3986 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3987 *----------------------------------------------------------------------------*/
3989 int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
3992 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3993 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3994 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3995 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3996 ) {
3997 float_raise( float_flag_invalid STATUS_VAR);
3998 return 0;
4000 return
4001 ( a.low == b.low )
4002 && ( ( a.high == b.high )
4003 || ( ( a.low == 0 )
4004 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
4005 );
4009 /*----------------------------------------------------------------------------
4010 | Returns 1 if the extended double-precision floating-point value `a' is less
4011 | than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
4012 | do not cause an exception. Otherwise, the comparison is performed according
4013 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4014 *----------------------------------------------------------------------------*/
4016 int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4018 flag aSign, bSign;
4020 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4021 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4022 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4023 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4024 ) {
4025 if ( floatx80_is_signaling_nan( a )
4026 || floatx80_is_signaling_nan( b ) ) {
4027 float_raise( float_flag_invalid STATUS_VAR);
4029 return 0;
4031 aSign = extractFloatx80Sign( a );
4032 bSign = extractFloatx80Sign( b );
4033 if ( aSign != bSign ) {
4034 return
4035 aSign
4036 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4037 == 0 );
4039 return
4040 aSign ? le128( b.high, b.low, a.high, a.low )
4041 : le128( a.high, a.low, b.high, b.low );
4045 /*----------------------------------------------------------------------------
4046 | Returns 1 if the extended double-precision floating-point value `a' is less
4047 | than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
4048 | an exception. Otherwise, the comparison is performed according to the
4049 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4050 *----------------------------------------------------------------------------*/
4052 int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4054 flag aSign, bSign;
4056 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4057 && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4058 || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4059 && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4060 ) {
4061 if ( floatx80_is_signaling_nan( a )
4062 || floatx80_is_signaling_nan( b ) ) {
4063 float_raise( float_flag_invalid STATUS_VAR);
4065 return 0;
4067 aSign = extractFloatx80Sign( a );
4068 bSign = extractFloatx80Sign( b );
4069 if ( aSign != bSign ) {
4070 return
4071 aSign
4072 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4073 != 0 );
4075 return
4076 aSign ? lt128( b.high, b.low, a.high, a.low )
4077 : lt128( a.high, a.low, b.high, b.low );
4081 #endif
4083 #ifdef FLOAT128
4085 /*----------------------------------------------------------------------------
4086 | Returns the result of converting the quadruple-precision floating-point
4087 | value `a' to the 32-bit two's complement integer format. The conversion
4088 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4089 | Arithmetic---which means in particular that the conversion is rounded
4090 | according to the current rounding mode. If `a' is a NaN, the largest
4091 | positive integer is returned. Otherwise, if the conversion overflows, the
4092 | largest integer with the same sign as `a' is returned.
4093 *----------------------------------------------------------------------------*/
4095 int32 float128_to_int32( float128 a STATUS_PARAM )
4097 flag aSign;
4098 int32 aExp, shiftCount;
4099 bits64 aSig0, aSig1;
4101 aSig1 = extractFloat128Frac1( a );
4102 aSig0 = extractFloat128Frac0( a );
4103 aExp = extractFloat128Exp( a );
4104 aSign = extractFloat128Sign( a );
4105 if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
4106 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4107 aSig0 |= ( aSig1 != 0 );
4108 shiftCount = 0x4028 - aExp;
4109 if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
4110 return roundAndPackInt32( aSign, aSig0 STATUS_VAR );
4114 /*----------------------------------------------------------------------------
4115 | Returns the result of converting the quadruple-precision floating-point
4116 | value `a' to the 32-bit two's complement integer format. The conversion
4117 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4118 | Arithmetic, except that the conversion is always rounded toward zero. If
4119 | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
4120 | conversion overflows, the largest integer with the same sign as `a' is
4121 | returned.
4122 *----------------------------------------------------------------------------*/
4124 int32 float128_to_int32_round_to_zero( float128 a STATUS_PARAM )
4126 flag aSign;
4127 int32 aExp, shiftCount;
4128 bits64 aSig0, aSig1, savedASig;
4129 int32 z;
4131 aSig1 = extractFloat128Frac1( a );
4132 aSig0 = extractFloat128Frac0( a );
4133 aExp = extractFloat128Exp( a );
4134 aSign = extractFloat128Sign( a );
4135 aSig0 |= ( aSig1 != 0 );
4136 if ( 0x401E < aExp ) {
4137 if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
4138 goto invalid;
4140 else if ( aExp < 0x3FFF ) {
4141 if ( aExp || aSig0 ) STATUS(float_exception_flags) |= float_flag_inexact;
4142 return 0;
4144 aSig0 |= LIT64( 0x0001000000000000 );
4145 shiftCount = 0x402F - aExp;
4146 savedASig = aSig0;
4147 aSig0 >>= shiftCount;
4148 z = aSig0;
4149 if ( aSign ) z = - z;
4150 if ( ( z < 0 ) ^ aSign ) {
4151 invalid:
4152 float_raise( float_flag_invalid STATUS_VAR);
4153 return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
4155 if ( ( aSig0<<shiftCount ) != savedASig ) {
4156 STATUS(float_exception_flags) |= float_flag_inexact;
4158 return z;
4162 /*----------------------------------------------------------------------------
4163 | Returns the result of converting the quadruple-precision floating-point
4164 | value `a' to the 64-bit two's complement integer format. The conversion
4165 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4166 | Arithmetic---which means in particular that the conversion is rounded
4167 | according to the current rounding mode. If `a' is a NaN, the largest
4168 | positive integer is returned. Otherwise, if the conversion overflows, the
4169 | largest integer with the same sign as `a' is returned.
4170 *----------------------------------------------------------------------------*/
4172 int64 float128_to_int64( float128 a STATUS_PARAM )
4174 flag aSign;
4175 int32 aExp, shiftCount;
4176 bits64 aSig0, aSig1;
4178 aSig1 = extractFloat128Frac1( a );
4179 aSig0 = extractFloat128Frac0( a );
4180 aExp = extractFloat128Exp( a );
4181 aSign = extractFloat128Sign( a );
4182 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4183 shiftCount = 0x402F - aExp;
4184 if ( shiftCount <= 0 ) {
4185 if ( 0x403E < aExp ) {
4186 float_raise( float_flag_invalid STATUS_VAR);
4187 if ( ! aSign
4188 || ( ( aExp == 0x7FFF )
4189 && ( aSig1 || ( aSig0 != LIT64( 0x0001000000000000 ) ) )
4191 ) {
4192 return LIT64( 0x7FFFFFFFFFFFFFFF );
4194 return (sbits64) LIT64( 0x8000000000000000 );
4196 shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
4198 else {
4199 shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
4201 return roundAndPackInt64( aSign, aSig0, aSig1 STATUS_VAR );
4205 /*----------------------------------------------------------------------------
4206 | Returns the result of converting the quadruple-precision floating-point
4207 | value `a' to the 64-bit two's complement integer format. The conversion
4208 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4209 | Arithmetic, except that the conversion is always rounded toward zero.
4210 | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
4211 | the conversion overflows, the largest integer with the same sign as `a' is
4212 | returned.
4213 *----------------------------------------------------------------------------*/
4215 int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
4217 flag aSign;
4218 int32 aExp, shiftCount;
4219 bits64 aSig0, aSig1;
4220 int64 z;
4222 aSig1 = extractFloat128Frac1( a );
4223 aSig0 = extractFloat128Frac0( a );
4224 aExp = extractFloat128Exp( a );
4225 aSign = extractFloat128Sign( a );
4226 if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4227 shiftCount = aExp - 0x402F;
4228 if ( 0 < shiftCount ) {
4229 if ( 0x403E <= aExp ) {
4230 aSig0 &= LIT64( 0x0000FFFFFFFFFFFF );
4231 if ( ( a.high == LIT64( 0xC03E000000000000 ) )
4232 && ( aSig1 < LIT64( 0x0002000000000000 ) ) ) {
4233 if ( aSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
4235 else {
4236 float_raise( float_flag_invalid STATUS_VAR);
4237 if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
4238 return LIT64( 0x7FFFFFFFFFFFFFFF );
4241 return (sbits64) LIT64( 0x8000000000000000 );
4243 z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
4244 if ( (bits64) ( aSig1<<shiftCount ) ) {
4245 STATUS(float_exception_flags) |= float_flag_inexact;
4248 else {
4249 if ( aExp < 0x3FFF ) {
4250 if ( aExp | aSig0 | aSig1 ) {
4251 STATUS(float_exception_flags) |= float_flag_inexact;
4253 return 0;
4255 z = aSig0>>( - shiftCount );
4256 if ( aSig1
4257 || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
4258 STATUS(float_exception_flags) |= float_flag_inexact;
4261 if ( aSign ) z = - z;
4262 return z;
4266 /*----------------------------------------------------------------------------
4267 | Returns the result of converting the quadruple-precision floating-point
4268 | value `a' to the single-precision floating-point format. The conversion
4269 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4270 | Arithmetic.
4271 *----------------------------------------------------------------------------*/
4273 float32 float128_to_float32( float128 a STATUS_PARAM )
4275 flag aSign;
4276 int32 aExp;
4277 bits64 aSig0, aSig1;
4278 bits32 zSig;
4280 aSig1 = extractFloat128Frac1( a );
4281 aSig0 = extractFloat128Frac0( a );
4282 aExp = extractFloat128Exp( a );
4283 aSign = extractFloat128Sign( a );
4284 if ( aExp == 0x7FFF ) {
4285 if ( aSig0 | aSig1 ) {
4286 return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR ) );
4288 return packFloat32( aSign, 0xFF, 0 );
4290 aSig0 |= ( aSig1 != 0 );
4291 shift64RightJamming( aSig0, 18, &aSig0 );
4292 zSig = aSig0;
4293 if ( aExp || zSig ) {
4294 zSig |= 0x40000000;
4295 aExp -= 0x3F81;
4297 return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
4301 /*----------------------------------------------------------------------------
4302 | Returns the result of converting the quadruple-precision floating-point
4303 | value `a' to the double-precision floating-point format. The conversion
4304 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4305 | Arithmetic.
4306 *----------------------------------------------------------------------------*/
4308 float64 float128_to_float64( float128 a STATUS_PARAM )
4310 flag aSign;
4311 int32 aExp;
4312 bits64 aSig0, aSig1;
4314 aSig1 = extractFloat128Frac1( a );
4315 aSig0 = extractFloat128Frac0( a );
4316 aExp = extractFloat128Exp( a );
4317 aSign = extractFloat128Sign( a );
4318 if ( aExp == 0x7FFF ) {
4319 if ( aSig0 | aSig1 ) {
4320 return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR ) );
4322 return packFloat64( aSign, 0x7FF, 0 );
4324 shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4325 aSig0 |= ( aSig1 != 0 );
4326 if ( aExp || aSig0 ) {
4327 aSig0 |= LIT64( 0x4000000000000000 );
4328 aExp -= 0x3C01;
4330 return roundAndPackFloat64( aSign, aExp, aSig0 STATUS_VAR );
4334 #ifdef FLOATX80
4336 /*----------------------------------------------------------------------------
4337 | Returns the result of converting the quadruple-precision floating-point
4338 | value `a' to the extended double-precision floating-point format. The
4339 | conversion is performed according to the IEC/IEEE Standard for Binary
4340 | Floating-Point Arithmetic.
4341 *----------------------------------------------------------------------------*/
4343 floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
4345 flag aSign;
4346 int32 aExp;
4347 bits64 aSig0, aSig1;
4349 aSig1 = extractFloat128Frac1( a );
4350 aSig0 = extractFloat128Frac0( a );
4351 aExp = extractFloat128Exp( a );
4352 aSign = extractFloat128Sign( a );
4353 if ( aExp == 0x7FFF ) {
4354 if ( aSig0 | aSig1 ) {
4355 return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR ) );
4357 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
4359 if ( aExp == 0 ) {
4360 if ( ( aSig0 | aSig1 ) == 0 ) return packFloatx80( aSign, 0, 0 );
4361 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4363 else {
4364 aSig0 |= LIT64( 0x0001000000000000 );
4366 shortShift128Left( aSig0, aSig1, 15, &aSig0, &aSig1 );
4367 return roundAndPackFloatx80( 80, aSign, aExp, aSig0, aSig1 STATUS_VAR );
4371 #endif
4373 /*----------------------------------------------------------------------------
4374 | Rounds the quadruple-precision floating-point value `a' to an integer, and
4375 | returns the result as a quadruple-precision floating-point value. The
4376 | operation is performed according to the IEC/IEEE Standard for Binary
4377 | Floating-Point Arithmetic.
4378 *----------------------------------------------------------------------------*/
4380 float128 float128_round_to_int( float128 a STATUS_PARAM )
4382 flag aSign;
4383 int32 aExp;
4384 bits64 lastBitMask, roundBitsMask;
4385 int8 roundingMode;
4386 float128 z;
4388 aExp = extractFloat128Exp( a );
4389 if ( 0x402F <= aExp ) {
4390 if ( 0x406F <= aExp ) {
4391 if ( ( aExp == 0x7FFF )
4392 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
4393 ) {
4394 return propagateFloat128NaN( a, a STATUS_VAR );
4396 return a;
4398 lastBitMask = 1;
4399 lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
4400 roundBitsMask = lastBitMask - 1;
4401 z = a;
4402 roundingMode = STATUS(float_rounding_mode);
4403 if ( roundingMode == float_round_nearest_even ) {
4404 if ( lastBitMask ) {
4405 add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
4406 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
4408 else {
4409 if ( (sbits64) z.low < 0 ) {
4410 ++z.high;
4411 if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
4415 else if ( roundingMode != float_round_to_zero ) {
4416 if ( extractFloat128Sign( z )
4417 ^ ( roundingMode == float_round_up ) ) {
4418 add128( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
4421 z.low &= ~ roundBitsMask;
4423 else {
4424 if ( aExp < 0x3FFF ) {
4425 if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
4426 STATUS(float_exception_flags) |= float_flag_inexact;
4427 aSign = extractFloat128Sign( a );
4428 switch ( STATUS(float_rounding_mode) ) {
4429 case float_round_nearest_even:
4430 if ( ( aExp == 0x3FFE )
4431 && ( extractFloat128Frac0( a )
4432 | extractFloat128Frac1( a ) )
4433 ) {
4434 return packFloat128( aSign, 0x3FFF, 0, 0 );
4436 break;
4437 case float_round_down:
4438 return
4439 aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
4440 : packFloat128( 0, 0, 0, 0 );
4441 case float_round_up:
4442 return
4443 aSign ? packFloat128( 1, 0, 0, 0 )
4444 : packFloat128( 0, 0x3FFF, 0, 0 );
4446 return packFloat128( aSign, 0, 0, 0 );
4448 lastBitMask = 1;
4449 lastBitMask <<= 0x402F - aExp;
4450 roundBitsMask = lastBitMask - 1;
4451 z.low = 0;
4452 z.high = a.high;
4453 roundingMode = STATUS(float_rounding_mode);
4454 if ( roundingMode == float_round_nearest_even ) {
4455 z.high += lastBitMask>>1;
4456 if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
4457 z.high &= ~ lastBitMask;
4460 else if ( roundingMode != float_round_to_zero ) {
4461 if ( extractFloat128Sign( z )
4462 ^ ( roundingMode == float_round_up ) ) {
4463 z.high |= ( a.low != 0 );
4464 z.high += roundBitsMask;
4467 z.high &= ~ roundBitsMask;
4469 if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
4470 STATUS(float_exception_flags) |= float_flag_inexact;
4472 return z;
4476 /*----------------------------------------------------------------------------
4477 | Returns the result of adding the absolute values of the quadruple-precision
4478 | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
4479 | before being returned. `zSign' is ignored if the result is a NaN.
4480 | The addition is performed according to the IEC/IEEE Standard for Binary
4481 | Floating-Point Arithmetic.
4482 *----------------------------------------------------------------------------*/
4484 static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4486 int32 aExp, bExp, zExp;
4487 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4488 int32 expDiff;
4490 aSig1 = extractFloat128Frac1( a );
4491 aSig0 = extractFloat128Frac0( a );
4492 aExp = extractFloat128Exp( a );
4493 bSig1 = extractFloat128Frac1( b );
4494 bSig0 = extractFloat128Frac0( b );
4495 bExp = extractFloat128Exp( b );
4496 expDiff = aExp - bExp;
4497 if ( 0 < expDiff ) {
4498 if ( aExp == 0x7FFF ) {
4499 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4500 return a;
4502 if ( bExp == 0 ) {
4503 --expDiff;
4505 else {
4506 bSig0 |= LIT64( 0x0001000000000000 );
4508 shift128ExtraRightJamming(
4509 bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
4510 zExp = aExp;
4512 else if ( expDiff < 0 ) {
4513 if ( bExp == 0x7FFF ) {
4514 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4515 return packFloat128( zSign, 0x7FFF, 0, 0 );
4517 if ( aExp == 0 ) {
4518 ++expDiff;
4520 else {
4521 aSig0 |= LIT64( 0x0001000000000000 );
4523 shift128ExtraRightJamming(
4524 aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
4525 zExp = bExp;
4527 else {
4528 if ( aExp == 0x7FFF ) {
4529 if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4530 return propagateFloat128NaN( a, b STATUS_VAR );
4532 return a;
4534 add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4535 if ( aExp == 0 ) return packFloat128( zSign, 0, zSig0, zSig1 );
4536 zSig2 = 0;
4537 zSig0 |= LIT64( 0x0002000000000000 );
4538 zExp = aExp;
4539 goto shiftRight1;
4541 aSig0 |= LIT64( 0x0001000000000000 );
4542 add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4543 --zExp;
4544 if ( zSig0 < LIT64( 0x0002000000000000 ) ) goto roundAndPack;
4545 ++zExp;
4546 shiftRight1:
4547 shift128ExtraRightJamming(
4548 zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4549 roundAndPack:
4550 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4554 /*----------------------------------------------------------------------------
4555 | Returns the result of subtracting the absolute values of the quadruple-
4556 | precision floating-point values `a' and `b'. If `zSign' is 1, the
4557 | difference is negated before being returned. `zSign' is ignored if the
4558 | result is a NaN. The subtraction is performed according to the IEC/IEEE
4559 | Standard for Binary Floating-Point Arithmetic.
4560 *----------------------------------------------------------------------------*/
4562 static float128 subFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4564 int32 aExp, bExp, zExp;
4565 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
4566 int32 expDiff;
4567 float128 z;
4569 aSig1 = extractFloat128Frac1( a );
4570 aSig0 = extractFloat128Frac0( a );
4571 aExp = extractFloat128Exp( a );
4572 bSig1 = extractFloat128Frac1( b );
4573 bSig0 = extractFloat128Frac0( b );
4574 bExp = extractFloat128Exp( b );
4575 expDiff = aExp - bExp;
4576 shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4577 shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
4578 if ( 0 < expDiff ) goto aExpBigger;
4579 if ( expDiff < 0 ) goto bExpBigger;
4580 if ( aExp == 0x7FFF ) {
4581 if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4582 return propagateFloat128NaN( a, b STATUS_VAR );
4584 float_raise( float_flag_invalid STATUS_VAR);
4585 z.low = float128_default_nan_low;
4586 z.high = float128_default_nan_high;
4587 return z;
4589 if ( aExp == 0 ) {
4590 aExp = 1;
4591 bExp = 1;
4593 if ( bSig0 < aSig0 ) goto aBigger;
4594 if ( aSig0 < bSig0 ) goto bBigger;
4595 if ( bSig1 < aSig1 ) goto aBigger;
4596 if ( aSig1 < bSig1 ) goto bBigger;
4597 return packFloat128( STATUS(float_rounding_mode) == float_round_down, 0, 0, 0 );
4598 bExpBigger:
4599 if ( bExp == 0x7FFF ) {
4600 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4601 return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
4603 if ( aExp == 0 ) {
4604 ++expDiff;
4606 else {
4607 aSig0 |= LIT64( 0x4000000000000000 );
4609 shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
4610 bSig0 |= LIT64( 0x4000000000000000 );
4611 bBigger:
4612 sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
4613 zExp = bExp;
4614 zSign ^= 1;
4615 goto normalizeRoundAndPack;
4616 aExpBigger:
4617 if ( aExp == 0x7FFF ) {
4618 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4619 return a;
4621 if ( bExp == 0 ) {
4622 --expDiff;
4624 else {
4625 bSig0 |= LIT64( 0x4000000000000000 );
4627 shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
4628 aSig0 |= LIT64( 0x4000000000000000 );
4629 aBigger:
4630 sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4631 zExp = aExp;
4632 normalizeRoundAndPack:
4633 --zExp;
4634 return normalizeRoundAndPackFloat128( zSign, zExp - 14, zSig0, zSig1 STATUS_VAR );
4638 /*----------------------------------------------------------------------------
4639 | Returns the result of adding the quadruple-precision floating-point values
4640 | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
4641 | for Binary Floating-Point Arithmetic.
4642 *----------------------------------------------------------------------------*/
4644 float128 float128_add( float128 a, float128 b STATUS_PARAM )
4646 flag aSign, bSign;
4648 aSign = extractFloat128Sign( a );
4649 bSign = extractFloat128Sign( b );
4650 if ( aSign == bSign ) {
4651 return addFloat128Sigs( a, b, aSign STATUS_VAR );
4653 else {
4654 return subFloat128Sigs( a, b, aSign STATUS_VAR );
4659 /*----------------------------------------------------------------------------
4660 | Returns the result of subtracting the quadruple-precision floating-point
4661 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4662 | Standard for Binary Floating-Point Arithmetic.
4663 *----------------------------------------------------------------------------*/
4665 float128 float128_sub( float128 a, float128 b STATUS_PARAM )
4667 flag aSign, bSign;
4669 aSign = extractFloat128Sign( a );
4670 bSign = extractFloat128Sign( b );
4671 if ( aSign == bSign ) {
4672 return subFloat128Sigs( a, b, aSign STATUS_VAR );
4674 else {
4675 return addFloat128Sigs( a, b, aSign STATUS_VAR );
4680 /*----------------------------------------------------------------------------
4681 | Returns the result of multiplying the quadruple-precision floating-point
4682 | values `a' and `b'. The operation is performed according to the IEC/IEEE
4683 | Standard for Binary Floating-Point Arithmetic.
4684 *----------------------------------------------------------------------------*/
4686 float128 float128_mul( float128 a, float128 b STATUS_PARAM )
4688 flag aSign, bSign, zSign;
4689 int32 aExp, bExp, zExp;
4690 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
4691 float128 z;
4693 aSig1 = extractFloat128Frac1( a );
4694 aSig0 = extractFloat128Frac0( a );
4695 aExp = extractFloat128Exp( a );
4696 aSign = extractFloat128Sign( a );
4697 bSig1 = extractFloat128Frac1( b );
4698 bSig0 = extractFloat128Frac0( b );
4699 bExp = extractFloat128Exp( b );
4700 bSign = extractFloat128Sign( b );
4701 zSign = aSign ^ bSign;
4702 if ( aExp == 0x7FFF ) {
4703 if ( ( aSig0 | aSig1 )
4704 || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
4705 return propagateFloat128NaN( a, b STATUS_VAR );
4707 if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
4708 return packFloat128( zSign, 0x7FFF, 0, 0 );
4710 if ( bExp == 0x7FFF ) {
4711 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4712 if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4713 invalid:
4714 float_raise( float_flag_invalid STATUS_VAR);
4715 z.low = float128_default_nan_low;
4716 z.high = float128_default_nan_high;
4717 return z;
4719 return packFloat128( zSign, 0x7FFF, 0, 0 );
4721 if ( aExp == 0 ) {
4722 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4723 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4725 if ( bExp == 0 ) {
4726 if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4727 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4729 zExp = aExp + bExp - 0x4000;
4730 aSig0 |= LIT64( 0x0001000000000000 );
4731 shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
4732 mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
4733 add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
4734 zSig2 |= ( zSig3 != 0 );
4735 if ( LIT64( 0x0002000000000000 ) <= zSig0 ) {
4736 shift128ExtraRightJamming(
4737 zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4738 ++zExp;
4740 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4744 /*----------------------------------------------------------------------------
4745 | Returns the result of dividing the quadruple-precision floating-point value
4746 | `a' by the corresponding value `b'. The operation is performed according to
4747 | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4748 *----------------------------------------------------------------------------*/
4750 float128 float128_div( float128 a, float128 b STATUS_PARAM )
4752 flag aSign, bSign, zSign;
4753 int32 aExp, bExp, zExp;
4754 bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4755 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
4756 float128 z;
4758 aSig1 = extractFloat128Frac1( a );
4759 aSig0 = extractFloat128Frac0( a );
4760 aExp = extractFloat128Exp( a );
4761 aSign = extractFloat128Sign( a );
4762 bSig1 = extractFloat128Frac1( b );
4763 bSig0 = extractFloat128Frac0( b );
4764 bExp = extractFloat128Exp( b );
4765 bSign = extractFloat128Sign( b );
4766 zSign = aSign ^ bSign;
4767 if ( aExp == 0x7FFF ) {
4768 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4769 if ( bExp == 0x7FFF ) {
4770 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4771 goto invalid;
4773 return packFloat128( zSign, 0x7FFF, 0, 0 );
4775 if ( bExp == 0x7FFF ) {
4776 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4777 return packFloat128( zSign, 0, 0, 0 );
4779 if ( bExp == 0 ) {
4780 if ( ( bSig0 | bSig1 ) == 0 ) {
4781 if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4782 invalid:
4783 float_raise( float_flag_invalid STATUS_VAR);
4784 z.low = float128_default_nan_low;
4785 z.high = float128_default_nan_high;
4786 return z;
4788 float_raise( float_flag_divbyzero STATUS_VAR);
4789 return packFloat128( zSign, 0x7FFF, 0, 0 );
4791 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4793 if ( aExp == 0 ) {
4794 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4795 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4797 zExp = aExp - bExp + 0x3FFD;
4798 shortShift128Left(
4799 aSig0 | LIT64( 0x0001000000000000 ), aSig1, 15, &aSig0, &aSig1 );
4800 shortShift128Left(
4801 bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
4802 if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
4803 shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
4804 ++zExp;
4806 zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
4807 mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
4808 sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
4809 while ( (sbits64) rem0 < 0 ) {
4810 --zSig0;
4811 add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
4813 zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
4814 if ( ( zSig1 & 0x3FFF ) <= 4 ) {
4815 mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
4816 sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
4817 while ( (sbits64) rem1 < 0 ) {
4818 --zSig1;
4819 add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
4821 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
4823 shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
4824 return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4828 /*----------------------------------------------------------------------------
4829 | Returns the remainder of the quadruple-precision floating-point value `a'
4830 | with respect to the corresponding value `b'. The operation is performed
4831 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4832 *----------------------------------------------------------------------------*/
4834 float128 float128_rem( float128 a, float128 b STATUS_PARAM )
4836 flag aSign, bSign, zSign;
4837 int32 aExp, bExp, expDiff;
4838 bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
4839 bits64 allZero, alternateASig0, alternateASig1, sigMean1;
4840 sbits64 sigMean0;
4841 float128 z;
4843 aSig1 = extractFloat128Frac1( a );
4844 aSig0 = extractFloat128Frac0( a );
4845 aExp = extractFloat128Exp( a );
4846 aSign = extractFloat128Sign( a );
4847 bSig1 = extractFloat128Frac1( b );
4848 bSig0 = extractFloat128Frac0( b );
4849 bExp = extractFloat128Exp( b );
4850 bSign = extractFloat128Sign( b );
4851 if ( aExp == 0x7FFF ) {
4852 if ( ( aSig0 | aSig1 )
4853 || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
4854 return propagateFloat128NaN( a, b STATUS_VAR );
4856 goto invalid;
4858 if ( bExp == 0x7FFF ) {
4859 if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4860 return a;
4862 if ( bExp == 0 ) {
4863 if ( ( bSig0 | bSig1 ) == 0 ) {
4864 invalid:
4865 float_raise( float_flag_invalid STATUS_VAR);
4866 z.low = float128_default_nan_low;
4867 z.high = float128_default_nan_high;
4868 return z;
4870 normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4872 if ( aExp == 0 ) {
4873 if ( ( aSig0 | aSig1 ) == 0 ) return a;
4874 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4876 expDiff = aExp - bExp;
4877 if ( expDiff < -1 ) return a;
4878 shortShift128Left(
4879 aSig0 | LIT64( 0x0001000000000000 ),
4880 aSig1,
4881 15 - ( expDiff < 0 ),
4882 &aSig0,
4883 &aSig1
4884 );
4885 shortShift128Left(
4886 bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
4887 q = le128( bSig0, bSig1, aSig0, aSig1 );
4888 if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
4889 expDiff -= 64;
4890 while ( 0 < expDiff ) {
4891 q = estimateDiv128To64( aSig0, aSig1, bSig0 );
4892 q = ( 4 < q ) ? q - 4 : 0;
4893 mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
4894 shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZero );
4895 shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero );
4896 sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 );
4897 expDiff -= 61;
4899 if ( -64 < expDiff ) {
4900 q = estimateDiv128To64( aSig0, aSig1, bSig0 );
4901 q = ( 4 < q ) ? q - 4 : 0;
4902 q >>= - expDiff;
4903 shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
4904 expDiff += 52;
4905 if ( expDiff < 0 ) {
4906 shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
4908 else {
4909 shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
4911 mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
4912 sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
4914 else {
4915 shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 );
4916 shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
4918 do {
4919 alternateASig0 = aSig0;
4920 alternateASig1 = aSig1;
4921 ++q;
4922 sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
4923 } while ( 0 <= (sbits64) aSig0 );
4924 add128(
4925 aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
4926 if ( ( sigMean0 < 0 )
4927 || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
4928 aSig0 = alternateASig0;
4929 aSig1 = alternateASig1;
4931 zSign = ( (sbits64) aSig0 < 0 );
4932 if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
4933 return
4934 normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 STATUS_VAR );
4938 /*----------------------------------------------------------------------------
4939 | Returns the square root of the quadruple-precision floating-point value `a'.
4940 | The operation is performed according to the IEC/IEEE Standard for Binary
4941 | Floating-Point Arithmetic.
4942 *----------------------------------------------------------------------------*/
4944 float128 float128_sqrt( float128 a STATUS_PARAM )
4946 flag aSign;
4947 int32 aExp, zExp;
4948 bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
4949 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
4950 float128 z;
4952 aSig1 = extractFloat128Frac1( a );
4953 aSig0 = extractFloat128Frac0( a );
4954 aExp = extractFloat128Exp( a );
4955 aSign = extractFloat128Sign( a );
4956 if ( aExp == 0x7FFF ) {
4957 if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, a STATUS_VAR );
4958 if ( ! aSign ) return a;
4959 goto invalid;
4961 if ( aSign ) {
4962 if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
4963 invalid:
4964 float_raise( float_flag_invalid STATUS_VAR);
4965 z.low = float128_default_nan_low;
4966 z.high = float128_default_nan_high;
4967 return z;
4969 if ( aExp == 0 ) {
4970 if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( 0, 0, 0, 0 );
4971 normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4973 zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFE;
4974 aSig0 |= LIT64( 0x0001000000000000 );
4975 zSig0 = estimateSqrt32( aExp, aSig0>>17 );
4976 shortShift128Left( aSig0, aSig1, 13 - ( aExp & 1 ), &aSig0, &aSig1 );
4977 zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
4978 doubleZSig0 = zSig0<<1;
4979 mul64To128( zSig0, zSig0, &term0, &term1 );
4980 sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
4981 while ( (sbits64) rem0 < 0 ) {
4982 --zSig0;
4983 doubleZSig0 -= 2;
4984 add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
4986 zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
4987 if ( ( zSig1 & 0x1FFF ) <= 5 ) {
4988 if ( zSig1 == 0 ) zSig1 = 1;
4989 mul64To128( doubleZSig0, zSig1, &term1, &term2 );
4990 sub128( rem1, 0, term1, term2, &rem1, &rem2 );
4991 mul64To128( zSig1, zSig1, &term2, &term3 );
4992 sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
4993 while ( (sbits64) rem1 < 0 ) {
4994 --zSig1;
4995 shortShift128Left( 0, zSig1, 1, &term2, &term3 );
4996 term3 |= 1;
4997 term2 |= doubleZSig0;
4998 add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
5000 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
5002 shift128ExtraRightJamming( zSig0, zSig1, 0, 14, &zSig0, &zSig1, &zSig2 );
5003 return roundAndPackFloat128( 0, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
5007 /*----------------------------------------------------------------------------
5008 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5009 | the corresponding value `b', and 0 otherwise. The comparison is performed
5010 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5011 *----------------------------------------------------------------------------*/
5013 int float128_eq( float128 a, float128 b STATUS_PARAM )
5016 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5017 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5018 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5019 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5020 ) {
5021 if ( float128_is_signaling_nan( a )
5022 || float128_is_signaling_nan( b ) ) {
5023 float_raise( float_flag_invalid STATUS_VAR);
5025 return 0;
5027 return
5028 ( a.low == b.low )
5029 && ( ( a.high == b.high )
5030 || ( ( a.low == 0 )
5031 && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5032 );
5036 /*----------------------------------------------------------------------------
5037 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5038 | or equal to the corresponding value `b', and 0 otherwise. The comparison
5039 | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5040 | Arithmetic.
5041 *----------------------------------------------------------------------------*/
5043 int float128_le( float128 a, float128 b STATUS_PARAM )
5045 flag aSign, bSign;
5047 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5048 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5049 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5050 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5051 ) {
5052 float_raise( float_flag_invalid STATUS_VAR);
5053 return 0;
5055 aSign = extractFloat128Sign( a );
5056 bSign = extractFloat128Sign( b );
5057 if ( aSign != bSign ) {
5058 return
5059 aSign
5060 || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5061 == 0 );
5063 return
5064 aSign ? le128( b.high, b.low, a.high, a.low )
5065 : le128( a.high, a.low, b.high, b.low );
5069 /*----------------------------------------------------------------------------
5070 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5071 | the corresponding value `b', and 0 otherwise. The comparison is performed
5072 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5073 *----------------------------------------------------------------------------*/
5075 int float128_lt( float128 a, float128 b STATUS_PARAM )
5077 flag aSign, bSign;
5079 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5080 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5081 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5082 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5083 ) {
5084 float_raise( float_flag_invalid STATUS_VAR);
5085 return 0;
5087 aSign = extractFloat128Sign( a );
5088 bSign = extractFloat128Sign( b );
5089 if ( aSign != bSign ) {
5090 return
5091 aSign
5092 && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5093 != 0 );
5095 return
5096 aSign ? lt128( b.high, b.low, a.high, a.low )
5097 : lt128( a.high, a.low, b.high, b.low );
5101 /*----------------------------------------------------------------------------
5102 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5103 | the corresponding value `b', and 0 otherwise. The invalid exception is
5104 | raised if either operand is a NaN. Otherwise, the comparison is performed
5105 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5106 *----------------------------------------------------------------------------*/
5108 int float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
5111 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5112 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5113 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5114 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5115 ) {
5116 float_raise( float_flag_invalid STATUS_VAR);
5117 return 0;
5119 return
5120 ( a.low == b.low )
5121 && ( ( a.high == b.high )
5122 || ( ( a.low == 0 )
5123 && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5124 );
5128 /*----------------------------------------------------------------------------
5129 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5130 | or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
5131 | cause an exception. Otherwise, the comparison is performed according to the
5132 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5133 *----------------------------------------------------------------------------*/
5135 int float128_le_quiet( float128 a, float128 b STATUS_PARAM )
5137 flag aSign, bSign;
5139 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5140 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5141 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5142 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5143 ) {
5144 if ( float128_is_signaling_nan( a )
5145 || float128_is_signaling_nan( b ) ) {
5146 float_raise( float_flag_invalid STATUS_VAR);
5148 return 0;
5150 aSign = extractFloat128Sign( a );
5151 bSign = extractFloat128Sign( b );
5152 if ( aSign != bSign ) {
5153 return
5154 aSign
5155 || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5156 == 0 );
5158 return
5159 aSign ? le128( b.high, b.low, a.high, a.low )
5160 : le128( a.high, a.low, b.high, b.low );
5164 /*----------------------------------------------------------------------------
5165 | Returns 1 if the quadruple-precision floating-point value `a' is less than
5166 | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
5167 | exception. Otherwise, the comparison is performed according to the IEC/IEEE
5168 | Standard for Binary Floating-Point Arithmetic.
5169 *----------------------------------------------------------------------------*/
5171 int float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
5173 flag aSign, bSign;
5175 if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5176 && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5177 || ( ( extractFloat128Exp( b ) == 0x7FFF )
5178 && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5179 ) {
5180 if ( float128_is_signaling_nan( a )
5181 || float128_is_signaling_nan( b ) ) {
5182 float_raise( float_flag_invalid STATUS_VAR);
5184 return 0;
5186 aSign = extractFloat128Sign( a );
5187 bSign = extractFloat128Sign( b );
5188 if ( aSign != bSign ) {
5189 return
5190 aSign
5191 && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5192 != 0 );
5194 return
5195 aSign ? lt128( b.high, b.low, a.high, a.low )
5196 : lt128( a.high, a.low, b.high, b.low );
5200 #endif
5202 /* misc functions */
5203 float32 uint32_to_float32( unsigned int a STATUS_PARAM )
5205 return int64_to_float32(a STATUS_VAR);
5208 float64 uint32_to_float64( unsigned int a STATUS_PARAM )
5210 return int64_to_float64(a STATUS_VAR);
5213 unsigned int float32_to_uint32( float32 a STATUS_PARAM )
5215 int64_t v;
5216 unsigned int res;
5218 v = float32_to_int64(a STATUS_VAR);
5219 if (v < 0) {
5220 res = 0;
5221 float_raise( float_flag_invalid STATUS_VAR);
5222 } else if (v > 0xffffffff) {
5223 res = 0xffffffff;
5224 float_raise( float_flag_invalid STATUS_VAR);
5225 } else {
5226 res = v;
5228 return res;
5231 unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
5233 int64_t v;
5234 unsigned int res;
5236 v = float32_to_int64_round_to_zero(a STATUS_VAR);
5237 if (v < 0) {
5238 res = 0;
5239 float_raise( float_flag_invalid STATUS_VAR);
5240 } else if (v > 0xffffffff) {
5241 res = 0xffffffff;
5242 float_raise( float_flag_invalid STATUS_VAR);
5243 } else {
5244 res = v;
5246 return res;
5249 unsigned int float64_to_uint32( float64 a STATUS_PARAM )
5251 int64_t v;
5252 unsigned int res;
5254 v = float64_to_int64(a STATUS_VAR);
5255 if (v < 0) {
5256 res = 0;
5257 float_raise( float_flag_invalid STATUS_VAR);
5258 } else if (v > 0xffffffff) {
5259 res = 0xffffffff;
5260 float_raise( float_flag_invalid STATUS_VAR);
5261 } else {
5262 res = v;
5264 return res;
5267 unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
5269 int64_t v;
5270 unsigned int res;
5272 v = float64_to_int64_round_to_zero(a STATUS_VAR);
5273 if (v < 0) {
5274 res = 0;
5275 float_raise( float_flag_invalid STATUS_VAR);
5276 } else if (v > 0xffffffff) {
5277 res = 0xffffffff;
5278 float_raise( float_flag_invalid STATUS_VAR);
5279 } else {
5280 res = v;
5282 return res;
5285 #define COMPARE(s, nan_exp) \
5286 INLINE int float ## s ## _compare_internal( float ## s a, float ## s b, \
5287 int is_quiet STATUS_PARAM ) \
5288 { \
5289 flag aSign, bSign; \
5291 if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) && \
5292 extractFloat ## s ## Frac( a ) ) || \
5293 ( ( extractFloat ## s ## Exp( b ) == nan_exp ) && \
5294 extractFloat ## s ## Frac( b ) )) { \
5295 if (!is_quiet || \
5296 float ## s ## _is_signaling_nan( a ) || \
5297 float ## s ## _is_signaling_nan( b ) ) { \
5298 float_raise( float_flag_invalid STATUS_VAR); \
5299 } \
5300 return float_relation_unordered; \
5301 } \
5302 aSign = extractFloat ## s ## Sign( a ); \
5303 bSign = extractFloat ## s ## Sign( b ); \
5304 if ( aSign != bSign ) { \
5305 if ( (bits ## s) ( ( a | b )<<1 ) == 0 ) { \
5306 /* zero case */ \
5307 return float_relation_equal; \
5308 } else { \
5309 return 1 - (2 * aSign); \
5310 } \
5311 } else { \
5312 if (a == b) { \
5313 return float_relation_equal; \
5314 } else { \
5315 return 1 - 2 * (aSign ^ ( a < b )); \
5316 } \
5317 } \
5318 } \
5320 int float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM ) \
5321 { \
5322 return float ## s ## _compare_internal(a, b, 0 STATUS_VAR); \
5323 } \
5325 int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
5326 { \
5327 return float ## s ## _compare_internal(a, b, 1 STATUS_VAR); \
5330 COMPARE(32, 0xff)
5331 COMPARE(64, 0x7ff)