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