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