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