![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to softfloat.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Qemu by Fabrice Bellard] / qemu / fpu|
Return to softfloat.c CVS log | Up to [Qemu by Fabrice Bellard] / qemu / fpu |
1.1 root 1:
2: /*============================================================================
3:
4: This C source file is part of the SoftFloat IEC/IEEE Floating-point Arithmetic
5: Package, Release 2b.
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 `http://www.cs.berkeley.edu/~jhauser/
15: arithmetic/SoftFloat.html'.
16:
17: THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
18: been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
19: RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
20: AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
21: COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
22: EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
23: INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
24: OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
25:
26: Derivative works are acceptable, even for commercial purposes, so long as
27: (1) the source code for the derivative work includes prominent notice that
28: the work is derivative, and (2) the source code includes prominent notice with
29: these four paragraphs for those parts of this code that are retained.
30:
31: =============================================================================*/
32:
33: #include "softfloat.h"
34:
35: /*----------------------------------------------------------------------------
36: | Primitive arithmetic functions, including multi-word arithmetic, and
37: | division and square root approximations. (Can be specialized to target if
38: | desired.)
39: *----------------------------------------------------------------------------*/
40: #include "softfloat-macros.h"
41:
42: /*----------------------------------------------------------------------------
43: | Functions and definitions to determine: (1) whether tininess for underflow
44: | is detected before or after rounding by default, (2) what (if anything)
45: | happens when exceptions are raised, (3) how signaling NaNs are distinguished
46: | from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
47: | are propagated from function inputs to output. These details are target-
48: | specific.
49: *----------------------------------------------------------------------------*/
50: #include "softfloat-specialize.h"
51:
52: void set_float_rounding_mode(int val STATUS_PARAM)
53: {
54: STATUS(float_rounding_mode) = val;
55: }
56:
57: void set_float_exception_flags(int val STATUS_PARAM)
58: {
59: STATUS(float_exception_flags) = val;
60: }
61:
62: #ifdef FLOATX80
63: void set_floatx80_rounding_precision(int val STATUS_PARAM)
64: {
65: STATUS(floatx80_rounding_precision) = val;
66: }
67: #endif
68:
69: /*----------------------------------------------------------------------------
70: | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
71: | and 7, and returns the properly rounded 32-bit integer corresponding to the
72: | input. If `zSign' is 1, the input is negated before being converted to an
73: | integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point input
74: | is simply rounded to an integer, with the inexact exception raised if the
75: | input cannot be represented exactly as an integer. However, if the fixed-
76: | point input is too large, the invalid exception is raised and the largest
77: | positive or negative integer is returned.
78: *----------------------------------------------------------------------------*/
79:
80: static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
81: {
82: int8 roundingMode;
83: flag roundNearestEven;
84: int8 roundIncrement, roundBits;
85: int32 z;
86:
87: roundingMode = STATUS(float_rounding_mode);
88: roundNearestEven = ( roundingMode == float_round_nearest_even );
89: roundIncrement = 0x40;
90: if ( ! roundNearestEven ) {
91: if ( roundingMode == float_round_to_zero ) {
92: roundIncrement = 0;
93: }
94: else {
95: roundIncrement = 0x7F;
96: if ( zSign ) {
97: if ( roundingMode == float_round_up ) roundIncrement = 0;
98: }
99: else {
100: if ( roundingMode == float_round_down ) roundIncrement = 0;
101: }
102: }
103: }
104: roundBits = absZ & 0x7F;
105: absZ = ( absZ + roundIncrement )>>7;
106: absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
107: z = absZ;
108: if ( zSign ) z = - z;
109: if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
110: float_raise( float_flag_invalid STATUS_VAR);
111: return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
112: }
113: if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
114: return z;
115:
116: }
117:
118: /*----------------------------------------------------------------------------
119: | Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
120: | `absZ1', with binary point between bits 63 and 64 (between the input words),
121: | and returns the properly rounded 64-bit integer corresponding to the input.
122: | If `zSign' is 1, the input is negated before being converted to an integer.
123: | Ordinarily, the fixed-point input is simply rounded to an integer, with
124: | the inexact exception raised if the input cannot be represented exactly as
125: | an integer. However, if the fixed-point input is too large, the invalid
126: | exception is raised and the largest positive or negative integer is
127: | returned.
128: *----------------------------------------------------------------------------*/
129:
130: static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PARAM)
131: {
132: int8 roundingMode;
133: flag roundNearestEven, increment;
134: int64 z;
135:
136: roundingMode = STATUS(float_rounding_mode);
137: roundNearestEven = ( roundingMode == float_round_nearest_even );
138: increment = ( (sbits64) absZ1 < 0 );
139: if ( ! roundNearestEven ) {
140: if ( roundingMode == float_round_to_zero ) {
141: increment = 0;
142: }
143: else {
144: if ( zSign ) {
145: increment = ( roundingMode == float_round_down ) && absZ1;
146: }
147: else {
148: increment = ( roundingMode == float_round_up ) && absZ1;
149: }
150: }
151: }
152: if ( increment ) {
153: ++absZ0;
154: if ( absZ0 == 0 ) goto overflow;
155: absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
156: }
157: z = absZ0;
158: if ( zSign ) z = - z;
159: if ( z && ( ( z < 0 ) ^ zSign ) ) {
160: overflow:
161: float_raise( float_flag_invalid STATUS_VAR);
162: return
163: zSign ? (sbits64) LIT64( 0x8000000000000000 )
164: : LIT64( 0x7FFFFFFFFFFFFFFF );
165: }
166: if ( absZ1 ) STATUS(float_exception_flags) |= float_flag_inexact;
167: return z;
168:
169: }
170:
171: /*----------------------------------------------------------------------------
172: | Returns the fraction bits of the single-precision floating-point value `a'.
173: *----------------------------------------------------------------------------*/
174:
175: INLINE bits32 extractFloat32Frac( float32 a )
176: {
177:
178: return a & 0x007FFFFF;
179:
180: }
181:
182: /*----------------------------------------------------------------------------
183: | Returns the exponent bits of the single-precision floating-point value `a'.
184: *----------------------------------------------------------------------------*/
185:
186: INLINE int16 extractFloat32Exp( float32 a )
187: {
188:
189: return ( a>>23 ) & 0xFF;
190:
191: }
192:
193: /*----------------------------------------------------------------------------
194: | Returns the sign bit of the single-precision floating-point value `a'.
195: *----------------------------------------------------------------------------*/
196:
197: INLINE flag extractFloat32Sign( float32 a )
198: {
199:
200: return a>>31;
201:
202: }
203:
204: /*----------------------------------------------------------------------------
205: | Normalizes the subnormal single-precision floating-point value represented
206: | by the denormalized significand `aSig'. The normalized exponent and
207: | significand are stored at the locations pointed to by `zExpPtr' and
208: | `zSigPtr', respectively.
209: *----------------------------------------------------------------------------*/
210:
211: static void
212: normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
213: {
214: int8 shiftCount;
215:
216: shiftCount = countLeadingZeros32( aSig ) - 8;
217: *zSigPtr = aSig<<shiftCount;
218: *zExpPtr = 1 - shiftCount;
219:
220: }
221:
222: /*----------------------------------------------------------------------------
223: | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
224: | single-precision floating-point value, returning the result. After being
225: | shifted into the proper positions, the three fields are simply added
226: | together to form the result. This means that any integer portion of `zSig'
227: | will be added into the exponent. Since a properly normalized significand
228: | will have an integer portion equal to 1, the `zExp' input should be 1 less
229: | than the desired result exponent whenever `zSig' is a complete, normalized
230: | significand.
231: *----------------------------------------------------------------------------*/
232:
233: INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
234: {
235:
236: return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
237:
238: }
239:
240: /*----------------------------------------------------------------------------
241: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
242: | and significand `zSig', and returns the proper single-precision floating-
243: | point value corresponding to the abstract input. Ordinarily, the abstract
244: | value is simply rounded and packed into the single-precision format, with
245: | the inexact exception raised if the abstract input cannot be represented
246: | exactly. However, if the abstract value is too large, the overflow and
247: | inexact exceptions are raised and an infinity or maximal finite value is
248: | returned. If the abstract value is too small, the input value is rounded to
249: | a subnormal number, and the underflow and inexact exceptions are raised if
250: | the abstract input cannot be represented exactly as a subnormal single-
251: | precision floating-point number.
252: | The input significand `zSig' has its binary point between bits 30
253: | and 29, which is 7 bits to the left of the usual location. This shifted
254: | significand must be normalized or smaller. If `zSig' is not normalized,
255: | `zExp' must be 0; in that case, the result returned is a subnormal number,
256: | and it must not require rounding. In the usual case that `zSig' is
257: | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
258: | The handling of underflow and overflow follows the IEC/IEEE Standard for
259: | Binary Floating-Point Arithmetic.
260: *----------------------------------------------------------------------------*/
261:
262: static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
263: {
264: int8 roundingMode;
265: flag roundNearestEven;
266: int8 roundIncrement, roundBits;
267: flag isTiny;
268:
269: roundingMode = STATUS(float_rounding_mode);
270: roundNearestEven = ( roundingMode == float_round_nearest_even );
271: roundIncrement = 0x40;
272: if ( ! roundNearestEven ) {
273: if ( roundingMode == float_round_to_zero ) {
274: roundIncrement = 0;
275: }
276: else {
277: roundIncrement = 0x7F;
278: if ( zSign ) {
279: if ( roundingMode == float_round_up ) roundIncrement = 0;
280: }
281: else {
282: if ( roundingMode == float_round_down ) roundIncrement = 0;
283: }
284: }
285: }
286: roundBits = zSig & 0x7F;
287: if ( 0xFD <= (bits16) zExp ) {
288: if ( ( 0xFD < zExp )
289: || ( ( zExp == 0xFD )
290: && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
291: ) {
292: float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
293: return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
294: }
295: if ( zExp < 0 ) {
296: isTiny =
297: ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
298: || ( zExp < -1 )
299: || ( zSig + roundIncrement < 0x80000000 );
300: shift32RightJamming( zSig, - zExp, &zSig );
301: zExp = 0;
302: roundBits = zSig & 0x7F;
303: if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
304: }
305: }
306: if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
307: zSig = ( zSig + roundIncrement )>>7;
308: zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
309: if ( zSig == 0 ) zExp = 0;
310: return packFloat32( zSign, zExp, zSig );
311:
312: }
313:
314: /*----------------------------------------------------------------------------
315: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
316: | and significand `zSig', and returns the proper single-precision floating-
317: | point value corresponding to the abstract input. This routine is just like
318: | `roundAndPackFloat32' except that `zSig' does not have to be normalized.
319: | Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
320: | floating-point exponent.
321: *----------------------------------------------------------------------------*/
322:
323: static float32
324: normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
325: {
326: int8 shiftCount;
327:
328: shiftCount = countLeadingZeros32( zSig ) - 1;
329: return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
330:
331: }
332:
333: /*----------------------------------------------------------------------------
334: | Returns the fraction bits of the double-precision floating-point value `a'.
335: *----------------------------------------------------------------------------*/
336:
337: INLINE bits64 extractFloat64Frac( float64 a )
338: {
339:
340: return a & LIT64( 0x000FFFFFFFFFFFFF );
341:
342: }
343:
344: /*----------------------------------------------------------------------------
345: | Returns the exponent bits of the double-precision floating-point value `a'.
346: *----------------------------------------------------------------------------*/
347:
348: INLINE int16 extractFloat64Exp( float64 a )
349: {
350:
351: return ( a>>52 ) & 0x7FF;
352:
353: }
354:
355: /*----------------------------------------------------------------------------
356: | Returns the sign bit of the double-precision floating-point value `a'.
357: *----------------------------------------------------------------------------*/
358:
359: INLINE flag extractFloat64Sign( float64 a )
360: {
361:
362: return a>>63;
363:
364: }
365:
366: /*----------------------------------------------------------------------------
367: | Normalizes the subnormal double-precision floating-point value represented
368: | by the denormalized significand `aSig'. The normalized exponent and
369: | significand are stored at the locations pointed to by `zExpPtr' and
370: | `zSigPtr', respectively.
371: *----------------------------------------------------------------------------*/
372:
373: static void
374: normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
375: {
376: int8 shiftCount;
377:
378: shiftCount = countLeadingZeros64( aSig ) - 11;
379: *zSigPtr = aSig<<shiftCount;
380: *zExpPtr = 1 - shiftCount;
381:
382: }
383:
384: /*----------------------------------------------------------------------------
385: | Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
386: | double-precision floating-point value, returning the result. After being
387: | shifted into the proper positions, the three fields are simply added
388: | together to form the result. This means that any integer portion of `zSig'
389: | will be added into the exponent. Since a properly normalized significand
390: | will have an integer portion equal to 1, the `zExp' input should be 1 less
391: | than the desired result exponent whenever `zSig' is a complete, normalized
392: | significand.
393: *----------------------------------------------------------------------------*/
394:
395: INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
396: {
397:
398: return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;
399:
400: }
401:
402: /*----------------------------------------------------------------------------
403: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
404: | and significand `zSig', and returns the proper double-precision floating-
405: | point value corresponding to the abstract input. Ordinarily, the abstract
406: | value is simply rounded and packed into the double-precision format, with
407: | the inexact exception raised if the abstract input cannot be represented
408: | exactly. However, if the abstract value is too large, the overflow and
409: | inexact exceptions are raised and an infinity or maximal finite value is
410: | returned. If the abstract value is too small, the input value is rounded
411: | to a subnormal number, and the underflow and inexact exceptions are raised
412: | if the abstract input cannot be represented exactly as a subnormal double-
413: | precision floating-point number.
414: | The input significand `zSig' has its binary point between bits 62
415: | and 61, which is 10 bits to the left of the usual location. This shifted
416: | significand must be normalized or smaller. If `zSig' is not normalized,
417: | `zExp' must be 0; in that case, the result returned is a subnormal number,
418: | and it must not require rounding. In the usual case that `zSig' is
419: | normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
420: | The handling of underflow and overflow follows the IEC/IEEE Standard for
421: | Binary Floating-Point Arithmetic.
422: *----------------------------------------------------------------------------*/
423:
424: static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
425: {
426: int8 roundingMode;
427: flag roundNearestEven;
428: int16 roundIncrement, roundBits;
429: flag isTiny;
430:
431: roundingMode = STATUS(float_rounding_mode);
432: roundNearestEven = ( roundingMode == float_round_nearest_even );
433: roundIncrement = 0x200;
434: if ( ! roundNearestEven ) {
435: if ( roundingMode == float_round_to_zero ) {
436: roundIncrement = 0;
437: }
438: else {
439: roundIncrement = 0x3FF;
440: if ( zSign ) {
441: if ( roundingMode == float_round_up ) roundIncrement = 0;
442: }
443: else {
444: if ( roundingMode == float_round_down ) roundIncrement = 0;
445: }
446: }
447: }
448: roundBits = zSig & 0x3FF;
449: if ( 0x7FD <= (bits16) zExp ) {
450: if ( ( 0x7FD < zExp )
451: || ( ( zExp == 0x7FD )
452: && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
453: ) {
454: float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
455: return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );
456: }
457: if ( zExp < 0 ) {
458: isTiny =
459: ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
460: || ( zExp < -1 )
461: || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
462: shift64RightJamming( zSig, - zExp, &zSig );
463: zExp = 0;
464: roundBits = zSig & 0x3FF;
465: if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
466: }
467: }
468: if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
469: zSig = ( zSig + roundIncrement )>>10;
470: zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
471: if ( zSig == 0 ) zExp = 0;
472: return packFloat64( zSign, zExp, zSig );
473:
474: }
475:
476: /*----------------------------------------------------------------------------
477: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
478: | and significand `zSig', and returns the proper double-precision floating-
479: | point value corresponding to the abstract input. This routine is just like
480: | `roundAndPackFloat64' except that `zSig' does not have to be normalized.
481: | Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
482: | floating-point exponent.
483: *----------------------------------------------------------------------------*/
484:
485: static float64
486: normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
487: {
488: int8 shiftCount;
489:
490: shiftCount = countLeadingZeros64( zSig ) - 1;
491: return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount STATUS_VAR);
492:
493: }
494:
495: #ifdef FLOATX80
496:
497: /*----------------------------------------------------------------------------
498: | Returns the fraction bits of the extended double-precision floating-point
499: | value `a'.
500: *----------------------------------------------------------------------------*/
501:
502: INLINE bits64 extractFloatx80Frac( floatx80 a )
503: {
504:
505: return a.low;
506:
507: }
508:
509: /*----------------------------------------------------------------------------
510: | Returns the exponent bits of the extended double-precision floating-point
511: | value `a'.
512: *----------------------------------------------------------------------------*/
513:
514: INLINE int32 extractFloatx80Exp( floatx80 a )
515: {
516:
517: return a.high & 0x7FFF;
518:
519: }
520:
521: /*----------------------------------------------------------------------------
522: | Returns the sign bit of the extended double-precision floating-point value
523: | `a'.
524: *----------------------------------------------------------------------------*/
525:
526: INLINE flag extractFloatx80Sign( floatx80 a )
527: {
528:
529: return a.high>>15;
530:
531: }
532:
533: /*----------------------------------------------------------------------------
534: | Normalizes the subnormal extended double-precision floating-point value
535: | represented by the denormalized significand `aSig'. The normalized exponent
536: | and significand are stored at the locations pointed to by `zExpPtr' and
537: | `zSigPtr', respectively.
538: *----------------------------------------------------------------------------*/
539:
540: static void
541: normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
542: {
543: int8 shiftCount;
544:
545: shiftCount = countLeadingZeros64( aSig );
546: *zSigPtr = aSig<<shiftCount;
547: *zExpPtr = 1 - shiftCount;
548:
549: }
550:
551: /*----------------------------------------------------------------------------
552: | Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
553: | extended double-precision floating-point value, returning the result.
554: *----------------------------------------------------------------------------*/
555:
556: INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
557: {
558: floatx80 z;
559:
560: z.low = zSig;
561: z.high = ( ( (bits16) zSign )<<15 ) + zExp;
562: return z;
563:
564: }
565:
566: /*----------------------------------------------------------------------------
567: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
568: | and extended significand formed by the concatenation of `zSig0' and `zSig1',
569: | and returns the proper extended double-precision floating-point value
570: | corresponding to the abstract input. Ordinarily, the abstract value is
571: | rounded and packed into the extended double-precision format, with the
572: | inexact exception raised if the abstract input cannot be represented
573: | exactly. However, if the abstract value is too large, the overflow and
574: | inexact exceptions are raised and an infinity or maximal finite value is
575: | returned. If the abstract value is too small, the input value is rounded to
576: | a subnormal number, and the underflow and inexact exceptions are raised if
577: | the abstract input cannot be represented exactly as a subnormal extended
578: | double-precision floating-point number.
579: | If `roundingPrecision' is 32 or 64, the result is rounded to the same
580: | number of bits as single or double precision, respectively. Otherwise, the
581: | result is rounded to the full precision of the extended double-precision
582: | format.
583: | The input significand must be normalized or smaller. If the input
584: | significand is not normalized, `zExp' must be 0; in that case, the result
585: | returned is a subnormal number, and it must not require rounding. The
586: | handling of underflow and overflow follows the IEC/IEEE Standard for Binary
587: | Floating-Point Arithmetic.
588: *----------------------------------------------------------------------------*/
589:
590: static floatx80
591: roundAndPackFloatx80(
592: int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
593: STATUS_PARAM)
594: {
595: int8 roundingMode;
596: flag roundNearestEven, increment, isTiny;
597: int64 roundIncrement, roundMask, roundBits;
598:
599: roundingMode = STATUS(float_rounding_mode);
600: roundNearestEven = ( roundingMode == float_round_nearest_even );
601: if ( roundingPrecision == 80 ) goto precision80;
602: if ( roundingPrecision == 64 ) {
603: roundIncrement = LIT64( 0x0000000000000400 );
604: roundMask = LIT64( 0x00000000000007FF );
605: }
606: else if ( roundingPrecision == 32 ) {
607: roundIncrement = LIT64( 0x0000008000000000 );
608: roundMask = LIT64( 0x000000FFFFFFFFFF );
609: }
610: else {
611: goto precision80;
612: }
613: zSig0 |= ( zSig1 != 0 );
614: if ( ! roundNearestEven ) {
615: if ( roundingMode == float_round_to_zero ) {
616: roundIncrement = 0;
617: }
618: else {
619: roundIncrement = roundMask;
620: if ( zSign ) {
621: if ( roundingMode == float_round_up ) roundIncrement = 0;
622: }
623: else {
624: if ( roundingMode == float_round_down ) roundIncrement = 0;
625: }
626: }
627: }
628: roundBits = zSig0 & roundMask;
629: if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
630: if ( ( 0x7FFE < zExp )
631: || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
632: ) {
633: goto overflow;
634: }
635: if ( zExp <= 0 ) {
636: isTiny =
637: ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
638: || ( zExp < 0 )
639: || ( zSig0 <= zSig0 + roundIncrement );
640: shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
641: zExp = 0;
642: roundBits = zSig0 & roundMask;
643: if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
644: if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
645: zSig0 += roundIncrement;
646: if ( (sbits64) zSig0 < 0 ) zExp = 1;
647: roundIncrement = roundMask + 1;
648: if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
649: roundMask |= roundIncrement;
650: }
651: zSig0 &= ~ roundMask;
652: return packFloatx80( zSign, zExp, zSig0 );
653: }
654: }
655: if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
656: zSig0 += roundIncrement;
657: if ( zSig0 < roundIncrement ) {
658: ++zExp;
659: zSig0 = LIT64( 0x8000000000000000 );
660: }
661: roundIncrement = roundMask + 1;
662: if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
663: roundMask |= roundIncrement;
664: }
665: zSig0 &= ~ roundMask;
666: if ( zSig0 == 0 ) zExp = 0;
667: return packFloatx80( zSign, zExp, zSig0 );
668: precision80:
669: increment = ( (sbits64) zSig1 < 0 );
670: if ( ! roundNearestEven ) {
671: if ( roundingMode == float_round_to_zero ) {
672: increment = 0;
673: }
674: else {
675: if ( zSign ) {
676: increment = ( roundingMode == float_round_down ) && zSig1;
677: }
678: else {
679: increment = ( roundingMode == float_round_up ) && zSig1;
680: }
681: }
682: }
683: if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
684: if ( ( 0x7FFE < zExp )
685: || ( ( zExp == 0x7FFE )
686: && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
687: && increment
688: )
689: ) {
690: roundMask = 0;
691: overflow:
692: float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
693: if ( ( roundingMode == float_round_to_zero )
694: || ( zSign && ( roundingMode == float_round_up ) )
695: || ( ! zSign && ( roundingMode == float_round_down ) )
696: ) {
697: return packFloatx80( zSign, 0x7FFE, ~ roundMask );
698: }
699: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
700: }
701: if ( zExp <= 0 ) {
702: isTiny =
703: ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
704: || ( zExp < 0 )
705: || ! increment
706: || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
707: shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
708: zExp = 0;
709: if ( isTiny && zSig1 ) float_raise( float_flag_underflow STATUS_VAR);
710: if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
711: if ( roundNearestEven ) {
712: increment = ( (sbits64) zSig1 < 0 );
713: }
714: else {
715: if ( zSign ) {
716: increment = ( roundingMode == float_round_down ) && zSig1;
717: }
718: else {
719: increment = ( roundingMode == float_round_up ) && zSig1;
720: }
721: }
722: if ( increment ) {
723: ++zSig0;
724: zSig0 &=
725: ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
726: if ( (sbits64) zSig0 < 0 ) zExp = 1;
727: }
728: return packFloatx80( zSign, zExp, zSig0 );
729: }
730: }
731: if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
732: if ( increment ) {
733: ++zSig0;
734: if ( zSig0 == 0 ) {
735: ++zExp;
736: zSig0 = LIT64( 0x8000000000000000 );
737: }
738: else {
739: zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
740: }
741: }
742: else {
743: if ( zSig0 == 0 ) zExp = 0;
744: }
745: return packFloatx80( zSign, zExp, zSig0 );
746:
747: }
748:
749: /*----------------------------------------------------------------------------
750: | Takes an abstract floating-point value having sign `zSign', exponent
751: | `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
752: | and returns the proper extended double-precision floating-point value
753: | corresponding to the abstract input. This routine is just like
754: | `roundAndPackFloatx80' except that the input significand does not have to be
755: | normalized.
756: *----------------------------------------------------------------------------*/
757:
758: static floatx80
759: normalizeRoundAndPackFloatx80(
760: int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
761: STATUS_PARAM)
762: {
763: int8 shiftCount;
764:
765: if ( zSig0 == 0 ) {
766: zSig0 = zSig1;
767: zSig1 = 0;
768: zExp -= 64;
769: }
770: shiftCount = countLeadingZeros64( zSig0 );
771: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
772: zExp -= shiftCount;
773: return
774: roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 STATUS_VAR);
775:
776: }
777:
778: #endif
779:
780: #ifdef FLOAT128
781:
782: /*----------------------------------------------------------------------------
783: | Returns the least-significant 64 fraction bits of the quadruple-precision
784: | floating-point value `a'.
785: *----------------------------------------------------------------------------*/
786:
787: INLINE bits64 extractFloat128Frac1( float128 a )
788: {
789:
790: return a.low;
791:
792: }
793:
794: /*----------------------------------------------------------------------------
795: | Returns the most-significant 48 fraction bits of the quadruple-precision
796: | floating-point value `a'.
797: *----------------------------------------------------------------------------*/
798:
799: INLINE bits64 extractFloat128Frac0( float128 a )
800: {
801:
802: return a.high & LIT64( 0x0000FFFFFFFFFFFF );
803:
804: }
805:
806: /*----------------------------------------------------------------------------
807: | Returns the exponent bits of the quadruple-precision floating-point value
808: | `a'.
809: *----------------------------------------------------------------------------*/
810:
811: INLINE int32 extractFloat128Exp( float128 a )
812: {
813:
814: return ( a.high>>48 ) & 0x7FFF;
815:
816: }
817:
818: /*----------------------------------------------------------------------------
819: | Returns the sign bit of the quadruple-precision floating-point value `a'.
820: *----------------------------------------------------------------------------*/
821:
822: INLINE flag extractFloat128Sign( float128 a )
823: {
824:
825: return a.high>>63;
826:
827: }
828:
829: /*----------------------------------------------------------------------------
830: | Normalizes the subnormal quadruple-precision floating-point value
831: | represented by the denormalized significand formed by the concatenation of
832: | `aSig0' and `aSig1'. The normalized exponent is stored at the location
833: | pointed to by `zExpPtr'. The most significant 49 bits of the normalized
834: | significand are stored at the location pointed to by `zSig0Ptr', and the
835: | least significant 64 bits of the normalized significand are stored at the
836: | location pointed to by `zSig1Ptr'.
837: *----------------------------------------------------------------------------*/
838:
839: static void
840: normalizeFloat128Subnormal(
841: bits64 aSig0,
842: bits64 aSig1,
843: int32 *zExpPtr,
844: bits64 *zSig0Ptr,
845: bits64 *zSig1Ptr
846: )
847: {
848: int8 shiftCount;
849:
850: if ( aSig0 == 0 ) {
851: shiftCount = countLeadingZeros64( aSig1 ) - 15;
852: if ( shiftCount < 0 ) {
853: *zSig0Ptr = aSig1>>( - shiftCount );
854: *zSig1Ptr = aSig1<<( shiftCount & 63 );
855: }
856: else {
857: *zSig0Ptr = aSig1<<shiftCount;
858: *zSig1Ptr = 0;
859: }
860: *zExpPtr = - shiftCount - 63;
861: }
862: else {
863: shiftCount = countLeadingZeros64( aSig0 ) - 15;
864: shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
865: *zExpPtr = 1 - shiftCount;
866: }
867:
868: }
869:
870: /*----------------------------------------------------------------------------
871: | Packs the sign `zSign', the exponent `zExp', and the significand formed
872: | by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
873: | floating-point value, returning the result. After being shifted into the
874: | proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
875: | added together to form the most significant 32 bits of the result. This
876: | means that any integer portion of `zSig0' will be added into the exponent.
877: | Since a properly normalized significand will have an integer portion equal
878: | to 1, the `zExp' input should be 1 less than the desired result exponent
879: | whenever `zSig0' and `zSig1' concatenated form a complete, normalized
880: | significand.
881: *----------------------------------------------------------------------------*/
882:
883: INLINE float128
884: packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
885: {
886: float128 z;
887:
888: z.low = zSig1;
889: z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
890: return z;
891:
892: }
893:
894: /*----------------------------------------------------------------------------
895: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
896: | and extended significand formed by the concatenation of `zSig0', `zSig1',
897: | and `zSig2', and returns the proper quadruple-precision floating-point value
898: | corresponding to the abstract input. Ordinarily, the abstract value is
899: | simply rounded and packed into the quadruple-precision format, with the
900: | inexact exception raised if the abstract input cannot be represented
901: | exactly. However, if the abstract value is too large, the overflow and
902: | inexact exceptions are raised and an infinity or maximal finite value is
903: | returned. If the abstract value is too small, the input value is rounded to
904: | a subnormal number, and the underflow and inexact exceptions are raised if
905: | the abstract input cannot be represented exactly as a subnormal quadruple-
906: | precision floating-point number.
907: | The input significand must be normalized or smaller. If the input
908: | significand is not normalized, `zExp' must be 0; in that case, the result
909: | returned is a subnormal number, and it must not require rounding. In the
910: | usual case that the input significand is normalized, `zExp' must be 1 less
911: | than the ``true'' floating-point exponent. The handling of underflow and
912: | overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
913: *----------------------------------------------------------------------------*/
914:
915: static float128
916: roundAndPackFloat128(
917: flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 STATUS_PARAM)
918: {
919: int8 roundingMode;
920: flag roundNearestEven, increment, isTiny;
921:
922: roundingMode = STATUS(float_rounding_mode);
923: roundNearestEven = ( roundingMode == float_round_nearest_even );
924: increment = ( (sbits64) zSig2 < 0 );
925: if ( ! roundNearestEven ) {
926: if ( roundingMode == float_round_to_zero ) {
927: increment = 0;
928: }
929: else {
930: if ( zSign ) {
931: increment = ( roundingMode == float_round_down ) && zSig2;
932: }
933: else {
934: increment = ( roundingMode == float_round_up ) && zSig2;
935: }
936: }
937: }
938: if ( 0x7FFD <= (bits32) zExp ) {
939: if ( ( 0x7FFD < zExp )
940: || ( ( zExp == 0x7FFD )
941: && eq128(
942: LIT64( 0x0001FFFFFFFFFFFF ),
943: LIT64( 0xFFFFFFFFFFFFFFFF ),
944: zSig0,
945: zSig1
946: )
947: && increment
948: )
949: ) {
950: float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
951: if ( ( roundingMode == float_round_to_zero )
952: || ( zSign && ( roundingMode == float_round_up ) )
953: || ( ! zSign && ( roundingMode == float_round_down ) )
954: ) {
955: return
956: packFloat128(
957: zSign,
958: 0x7FFE,
959: LIT64( 0x0000FFFFFFFFFFFF ),
960: LIT64( 0xFFFFFFFFFFFFFFFF )
961: );
962: }
963: return packFloat128( zSign, 0x7FFF, 0, 0 );
964: }
965: if ( zExp < 0 ) {
966: isTiny =
967: ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
968: || ( zExp < -1 )
969: || ! increment
970: || lt128(
971: zSig0,
972: zSig1,
973: LIT64( 0x0001FFFFFFFFFFFF ),
974: LIT64( 0xFFFFFFFFFFFFFFFF )
975: );
976: shift128ExtraRightJamming(
977: zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
978: zExp = 0;
979: if ( isTiny && zSig2 ) float_raise( float_flag_underflow STATUS_VAR);
980: if ( roundNearestEven ) {
981: increment = ( (sbits64) zSig2 < 0 );
982: }
983: else {
984: if ( zSign ) {
985: increment = ( roundingMode == float_round_down ) && zSig2;
986: }
987: else {
988: increment = ( roundingMode == float_round_up ) && zSig2;
989: }
990: }
991: }
992: }
993: if ( zSig2 ) STATUS(float_exception_flags) |= float_flag_inexact;
994: if ( increment ) {
995: add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
996: zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
997: }
998: else {
999: if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
1000: }
1001: return packFloat128( zSign, zExp, zSig0, zSig1 );
1002:
1003: }
1004:
1005: /*----------------------------------------------------------------------------
1006: | Takes an abstract floating-point value having sign `zSign', exponent `zExp',
1007: | and significand formed by the concatenation of `zSig0' and `zSig1', and
1008: | returns the proper quadruple-precision floating-point value corresponding
1009: | to the abstract input. This routine is just like `roundAndPackFloat128'
1010: | except that the input significand has fewer bits and does not have to be
1011: | normalized. In all cases, `zExp' must be 1 less than the ``true'' floating-
1012: | point exponent.
1013: *----------------------------------------------------------------------------*/
1014:
1015: static float128
1016: normalizeRoundAndPackFloat128(
1017: flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 STATUS_PARAM)
1018: {
1019: int8 shiftCount;
1020: bits64 zSig2;
1021:
1022: if ( zSig0 == 0 ) {
1023: zSig0 = zSig1;
1024: zSig1 = 0;
1025: zExp -= 64;
1026: }
1027: shiftCount = countLeadingZeros64( zSig0 ) - 15;
1028: if ( 0 <= shiftCount ) {
1029: zSig2 = 0;
1030: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1031: }
1032: else {
1033: shift128ExtraRightJamming(
1034: zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
1035: }
1036: zExp -= shiftCount;
1037: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR);
1038:
1039: }
1040:
1041: #endif
1042:
1043: /*----------------------------------------------------------------------------
1044: | Returns the result of converting the 32-bit two's complement integer `a'
1045: | to the single-precision floating-point format. The conversion is performed
1046: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1047: *----------------------------------------------------------------------------*/
1048:
1049: float32 int32_to_float32( int32 a STATUS_PARAM )
1050: {
1051: flag zSign;
1052:
1053: if ( a == 0 ) return 0;
1054: if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
1055: zSign = ( a < 0 );
1056: return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a STATUS_VAR );
1057:
1058: }
1059:
1060: /*----------------------------------------------------------------------------
1061: | Returns the result of converting the 32-bit two's complement integer `a'
1062: | to the double-precision floating-point format. The conversion is performed
1063: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1064: *----------------------------------------------------------------------------*/
1065:
1066: float64 int32_to_float64( int32 a STATUS_PARAM )
1067: {
1068: flag zSign;
1069: uint32 absA;
1070: int8 shiftCount;
1071: bits64 zSig;
1072:
1073: if ( a == 0 ) return 0;
1074: zSign = ( a < 0 );
1075: absA = zSign ? - a : a;
1076: shiftCount = countLeadingZeros32( absA ) + 21;
1077: zSig = absA;
1078: return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );
1079:
1080: }
1081:
1082: #ifdef FLOATX80
1083:
1084: /*----------------------------------------------------------------------------
1085: | Returns the result of converting the 32-bit two's complement integer `a'
1086: | to the extended double-precision floating-point format. The conversion
1087: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1088: | Arithmetic.
1089: *----------------------------------------------------------------------------*/
1090:
1091: floatx80 int32_to_floatx80( int32 a STATUS_PARAM )
1092: {
1093: flag zSign;
1094: uint32 absA;
1095: int8 shiftCount;
1096: bits64 zSig;
1097:
1098: if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1099: zSign = ( a < 0 );
1100: absA = zSign ? - a : a;
1101: shiftCount = countLeadingZeros32( absA ) + 32;
1102: zSig = absA;
1103: return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
1104:
1105: }
1106:
1107: #endif
1108:
1109: #ifdef FLOAT128
1110:
1111: /*----------------------------------------------------------------------------
1112: | Returns the result of converting the 32-bit two's complement integer `a' to
1113: | the quadruple-precision floating-point format. The conversion is performed
1114: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1115: *----------------------------------------------------------------------------*/
1116:
1117: float128 int32_to_float128( int32 a STATUS_PARAM )
1118: {
1119: flag zSign;
1120: uint32 absA;
1121: int8 shiftCount;
1122: bits64 zSig0;
1123:
1124: if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1125: zSign = ( a < 0 );
1126: absA = zSign ? - a : a;
1127: shiftCount = countLeadingZeros32( absA ) + 17;
1128: zSig0 = absA;
1129: return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );
1130:
1131: }
1132:
1133: #endif
1134:
1135: /*----------------------------------------------------------------------------
1136: | Returns the result of converting the 64-bit two's complement integer `a'
1137: | to the single-precision floating-point format. The conversion is performed
1138: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1139: *----------------------------------------------------------------------------*/
1140:
1141: float32 int64_to_float32( int64 a STATUS_PARAM )
1142: {
1143: flag zSign;
1144: uint64 absA;
1145: int8 shiftCount;
1146:
1147: if ( a == 0 ) return 0;
1148: zSign = ( a < 0 );
1149: absA = zSign ? - a : a;
1150: shiftCount = countLeadingZeros64( absA ) - 40;
1151: if ( 0 <= shiftCount ) {
1152: return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );
1153: }
1154: else {
1155: shiftCount += 7;
1156: if ( shiftCount < 0 ) {
1157: shift64RightJamming( absA, - shiftCount, &absA );
1158: }
1159: else {
1160: absA <<= shiftCount;
1161: }
1162: return roundAndPackFloat32( zSign, 0x9C - shiftCount, absA STATUS_VAR );
1163: }
1164:
1165: }
1166:
1167: /*----------------------------------------------------------------------------
1168: | Returns the result of converting the 64-bit two's complement integer `a'
1169: | to the double-precision floating-point format. The conversion is performed
1170: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1171: *----------------------------------------------------------------------------*/
1172:
1173: float64 int64_to_float64( int64 a STATUS_PARAM )
1174: {
1175: flag zSign;
1176:
1177: if ( a == 0 ) return 0;
1178: if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
1179: return packFloat64( 1, 0x43E, 0 );
1180: }
1181: zSign = ( a < 0 );
1182: return normalizeRoundAndPackFloat64( zSign, 0x43C, zSign ? - a : a STATUS_VAR );
1183:
1184: }
1185:
1186: #ifdef FLOATX80
1187:
1188: /*----------------------------------------------------------------------------
1189: | Returns the result of converting the 64-bit two's complement integer `a'
1190: | to the extended double-precision floating-point format. The conversion
1191: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1192: | Arithmetic.
1193: *----------------------------------------------------------------------------*/
1194:
1195: floatx80 int64_to_floatx80( int64 a STATUS_PARAM )
1196: {
1197: flag zSign;
1198: uint64 absA;
1199: int8 shiftCount;
1200:
1201: if ( a == 0 ) return packFloatx80( 0, 0, 0 );
1202: zSign = ( a < 0 );
1203: absA = zSign ? - a : a;
1204: shiftCount = countLeadingZeros64( absA );
1205: return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
1206:
1207: }
1208:
1209: #endif
1210:
1211: #ifdef FLOAT128
1212:
1213: /*----------------------------------------------------------------------------
1214: | Returns the result of converting the 64-bit two's complement integer `a' to
1215: | the quadruple-precision floating-point format. The conversion is performed
1216: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1217: *----------------------------------------------------------------------------*/
1218:
1219: float128 int64_to_float128( int64 a STATUS_PARAM )
1220: {
1221: flag zSign;
1222: uint64 absA;
1223: int8 shiftCount;
1224: int32 zExp;
1225: bits64 zSig0, zSig1;
1226:
1227: if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
1228: zSign = ( a < 0 );
1229: absA = zSign ? - a : a;
1230: shiftCount = countLeadingZeros64( absA ) + 49;
1231: zExp = 0x406E - shiftCount;
1232: if ( 64 <= shiftCount ) {
1233: zSig1 = 0;
1234: zSig0 = absA;
1235: shiftCount -= 64;
1236: }
1237: else {
1238: zSig1 = absA;
1239: zSig0 = 0;
1240: }
1241: shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
1242: return packFloat128( zSign, zExp, zSig0, zSig1 );
1243:
1244: }
1245:
1246: #endif
1247:
1248: /*----------------------------------------------------------------------------
1249: | Returns the result of converting the single-precision floating-point value
1250: | `a' to the 32-bit two's complement integer format. The conversion is
1251: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1252: | Arithmetic---which means in particular that the conversion is rounded
1253: | according to the current rounding mode. If `a' is a NaN, the largest
1254: | positive integer is returned. Otherwise, if the conversion overflows, the
1255: | largest integer with the same sign as `a' is returned.
1256: *----------------------------------------------------------------------------*/
1257:
1258: int32 float32_to_int32( float32 a STATUS_PARAM )
1259: {
1260: flag aSign;
1261: int16 aExp, shiftCount;
1262: bits32 aSig;
1263: bits64 aSig64;
1264:
1265: aSig = extractFloat32Frac( a );
1266: aExp = extractFloat32Exp( a );
1267: aSign = extractFloat32Sign( a );
1268: if ( ( aExp == 0xFF ) && aSig ) aSign = 0;
1269: if ( aExp ) aSig |= 0x00800000;
1270: shiftCount = 0xAF - aExp;
1271: aSig64 = aSig;
1272: aSig64 <<= 32;
1273: if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );
1274: return roundAndPackInt32( aSign, aSig64 STATUS_VAR );
1275:
1276: }
1277:
1278: /*----------------------------------------------------------------------------
1279: | Returns the result of converting the single-precision floating-point value
1280: | `a' to the 32-bit two's complement integer format. The conversion is
1281: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1282: | Arithmetic, except that the conversion is always rounded toward zero.
1283: | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
1284: | the conversion overflows, the largest integer with the same sign as `a' is
1285: | returned.
1286: *----------------------------------------------------------------------------*/
1287:
1288: int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
1289: {
1290: flag aSign;
1291: int16 aExp, shiftCount;
1292: bits32 aSig;
1293: int32 z;
1294:
1295: aSig = extractFloat32Frac( a );
1296: aExp = extractFloat32Exp( a );
1297: aSign = extractFloat32Sign( a );
1298: shiftCount = aExp - 0x9E;
1299: if ( 0 <= shiftCount ) {
1300: if ( a != 0xCF000000 ) {
1301: float_raise( float_flag_invalid STATUS_VAR);
1302: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
1303: }
1304: return (sbits32) 0x80000000;
1305: }
1306: else if ( aExp <= 0x7E ) {
1307: if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1308: return 0;
1309: }
1310: aSig = ( aSig | 0x00800000 )<<8;
1311: z = aSig>>( - shiftCount );
1312: if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
1313: STATUS(float_exception_flags) |= float_flag_inexact;
1314: }
1315: if ( aSign ) z = - z;
1316: return z;
1317:
1318: }
1319:
1320: /*----------------------------------------------------------------------------
1321: | Returns the result of converting the single-precision floating-point value
1322: | `a' to the 64-bit two's complement integer format. The conversion is
1323: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1324: | Arithmetic---which means in particular that the conversion is rounded
1325: | according to the current rounding mode. If `a' is a NaN, the largest
1326: | positive integer is returned. Otherwise, if the conversion overflows, the
1327: | largest integer with the same sign as `a' is returned.
1328: *----------------------------------------------------------------------------*/
1329:
1330: int64 float32_to_int64( float32 a STATUS_PARAM )
1331: {
1332: flag aSign;
1333: int16 aExp, shiftCount;
1334: bits32 aSig;
1335: bits64 aSig64, aSigExtra;
1336:
1337: aSig = extractFloat32Frac( a );
1338: aExp = extractFloat32Exp( a );
1339: aSign = extractFloat32Sign( a );
1340: shiftCount = 0xBE - aExp;
1341: if ( shiftCount < 0 ) {
1342: float_raise( float_flag_invalid STATUS_VAR);
1343: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1344: return LIT64( 0x7FFFFFFFFFFFFFFF );
1345: }
1346: return (sbits64) LIT64( 0x8000000000000000 );
1347: }
1348: if ( aExp ) aSig |= 0x00800000;
1349: aSig64 = aSig;
1350: aSig64 <<= 40;
1351: shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );
1352: return roundAndPackInt64( aSign, aSig64, aSigExtra STATUS_VAR );
1353:
1354: }
1355:
1356: /*----------------------------------------------------------------------------
1357: | Returns the result of converting the single-precision floating-point value
1358: | `a' to the 64-bit two's complement integer format. The conversion is
1359: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1360: | Arithmetic, except that the conversion is always rounded toward zero. If
1361: | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
1362: | conversion overflows, the largest integer with the same sign as `a' is
1363: | returned.
1364: *----------------------------------------------------------------------------*/
1365:
1366: int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
1367: {
1368: flag aSign;
1369: int16 aExp, shiftCount;
1370: bits32 aSig;
1371: bits64 aSig64;
1372: int64 z;
1373:
1374: aSig = extractFloat32Frac( a );
1375: aExp = extractFloat32Exp( a );
1376: aSign = extractFloat32Sign( a );
1377: shiftCount = aExp - 0xBE;
1378: if ( 0 <= shiftCount ) {
1379: if ( a != 0xDF000000 ) {
1380: float_raise( float_flag_invalid STATUS_VAR);
1381: if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
1382: return LIT64( 0x7FFFFFFFFFFFFFFF );
1383: }
1384: }
1385: return (sbits64) LIT64( 0x8000000000000000 );
1386: }
1387: else if ( aExp <= 0x7E ) {
1388: if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
1389: return 0;
1390: }
1391: aSig64 = aSig | 0x00800000;
1392: aSig64 <<= 40;
1393: z = aSig64>>( - shiftCount );
1394: if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
1395: STATUS(float_exception_flags) |= float_flag_inexact;
1396: }
1397: if ( aSign ) z = - z;
1398: return z;
1399:
1400: }
1401:
1402: /*----------------------------------------------------------------------------
1403: | Returns the result of converting the single-precision floating-point value
1404: | `a' to the double-precision floating-point format. The conversion is
1405: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1406: | Arithmetic.
1407: *----------------------------------------------------------------------------*/
1408:
1409: float64 float32_to_float64( float32 a STATUS_PARAM )
1410: {
1411: flag aSign;
1412: int16 aExp;
1413: bits32 aSig;
1414:
1415: aSig = extractFloat32Frac( a );
1416: aExp = extractFloat32Exp( a );
1417: aSign = extractFloat32Sign( a );
1418: if ( aExp == 0xFF ) {
1419: if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ));
1420: return packFloat64( aSign, 0x7FF, 0 );
1421: }
1422: if ( aExp == 0 ) {
1423: if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
1424: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1425: --aExp;
1426: }
1427: return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
1428:
1429: }
1430:
1431: #ifdef FLOATX80
1432:
1433: /*----------------------------------------------------------------------------
1434: | Returns the result of converting the single-precision floating-point value
1435: | `a' to the extended double-precision floating-point format. The conversion
1436: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
1437: | Arithmetic.
1438: *----------------------------------------------------------------------------*/
1439:
1440: floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
1441: {
1442: flag aSign;
1443: int16 aExp;
1444: bits32 aSig;
1445:
1446: aSig = extractFloat32Frac( a );
1447: aExp = extractFloat32Exp( a );
1448: aSign = extractFloat32Sign( a );
1449: if ( aExp == 0xFF ) {
1450: if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) );
1451: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
1452: }
1453: if ( aExp == 0 ) {
1454: if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
1455: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1456: }
1457: aSig |= 0x00800000;
1458: return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
1459:
1460: }
1461:
1462: #endif
1463:
1464: #ifdef FLOAT128
1465:
1466: /*----------------------------------------------------------------------------
1467: | Returns the result of converting the single-precision floating-point value
1468: | `a' to the double-precision floating-point format. The conversion is
1469: | performed according to the IEC/IEEE Standard for Binary Floating-Point
1470: | Arithmetic.
1471: *----------------------------------------------------------------------------*/
1472:
1473: float128 float32_to_float128( float32 a STATUS_PARAM )
1474: {
1475: flag aSign;
1476: int16 aExp;
1477: bits32 aSig;
1478:
1479: aSig = extractFloat32Frac( a );
1480: aExp = extractFloat32Exp( a );
1481: aSign = extractFloat32Sign( a );
1482: if ( aExp == 0xFF ) {
1483: if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) );
1484: return packFloat128( aSign, 0x7FFF, 0, 0 );
1485: }
1486: if ( aExp == 0 ) {
1487: if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
1488: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1489: --aExp;
1490: }
1491: return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
1492:
1493: }
1494:
1495: #endif
1496:
1497: /*----------------------------------------------------------------------------
1498: | Rounds the single-precision floating-point value `a' to an integer, and
1499: | returns the result as a single-precision floating-point value. The
1500: | operation is performed according to the IEC/IEEE Standard for Binary
1501: | Floating-Point Arithmetic.
1502: *----------------------------------------------------------------------------*/
1503:
1504: float32 float32_round_to_int( float32 a STATUS_PARAM)
1505: {
1506: flag aSign;
1507: int16 aExp;
1508: bits32 lastBitMask, roundBitsMask;
1509: int8 roundingMode;
1510: float32 z;
1511:
1512: aExp = extractFloat32Exp( a );
1513: if ( 0x96 <= aExp ) {
1514: if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
1515: return propagateFloat32NaN( a, a STATUS_VAR );
1516: }
1517: return a;
1518: }
1519: if ( aExp <= 0x7E ) {
1520: if ( (bits32) ( a<<1 ) == 0 ) return a;
1521: STATUS(float_exception_flags) |= float_flag_inexact;
1522: aSign = extractFloat32Sign( a );
1523: switch ( STATUS(float_rounding_mode) ) {
1524: case float_round_nearest_even:
1525: if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
1526: return packFloat32( aSign, 0x7F, 0 );
1527: }
1528: break;
1529: case float_round_down:
1530: return aSign ? 0xBF800000 : 0;
1531: case float_round_up:
1532: return aSign ? 0x80000000 : 0x3F800000;
1533: }
1534: return packFloat32( aSign, 0, 0 );
1535: }
1536: lastBitMask = 1;
1537: lastBitMask <<= 0x96 - aExp;
1538: roundBitsMask = lastBitMask - 1;
1539: z = a;
1540: roundingMode = STATUS(float_rounding_mode);
1541: if ( roundingMode == float_round_nearest_even ) {
1542: z += lastBitMask>>1;
1543: if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
1544: }
1545: else if ( roundingMode != float_round_to_zero ) {
1546: if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
1547: z += roundBitsMask;
1548: }
1549: }
1550: z &= ~ roundBitsMask;
1551: if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
1552: return z;
1553:
1554: }
1555:
1556: /*----------------------------------------------------------------------------
1557: | Returns the result of adding the absolute values of the single-precision
1558: | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
1559: | before being returned. `zSign' is ignored if the result is a NaN.
1560: | The addition is performed according to the IEC/IEEE Standard for Binary
1561: | Floating-Point Arithmetic.
1562: *----------------------------------------------------------------------------*/
1563:
1564: static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1565: {
1566: int16 aExp, bExp, zExp;
1567: bits32 aSig, bSig, zSig;
1568: int16 expDiff;
1569:
1570: aSig = extractFloat32Frac( a );
1571: aExp = extractFloat32Exp( a );
1572: bSig = extractFloat32Frac( b );
1573: bExp = extractFloat32Exp( b );
1574: expDiff = aExp - bExp;
1575: aSig <<= 6;
1576: bSig <<= 6;
1577: if ( 0 < expDiff ) {
1578: if ( aExp == 0xFF ) {
1579: if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1580: return a;
1581: }
1582: if ( bExp == 0 ) {
1583: --expDiff;
1584: }
1585: else {
1586: bSig |= 0x20000000;
1587: }
1588: shift32RightJamming( bSig, expDiff, &bSig );
1589: zExp = aExp;
1590: }
1591: else if ( expDiff < 0 ) {
1592: if ( bExp == 0xFF ) {
1593: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1594: return packFloat32( zSign, 0xFF, 0 );
1595: }
1596: if ( aExp == 0 ) {
1597: ++expDiff;
1598: }
1599: else {
1600: aSig |= 0x20000000;
1601: }
1602: shift32RightJamming( aSig, - expDiff, &aSig );
1603: zExp = bExp;
1604: }
1605: else {
1606: if ( aExp == 0xFF ) {
1607: if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1608: return a;
1609: }
1610: if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
1611: zSig = 0x40000000 + aSig + bSig;
1612: zExp = aExp;
1613: goto roundAndPack;
1614: }
1615: aSig |= 0x20000000;
1616: zSig = ( aSig + bSig )<<1;
1617: --zExp;
1618: if ( (sbits32) zSig < 0 ) {
1619: zSig = aSig + bSig;
1620: ++zExp;
1621: }
1622: roundAndPack:
1623: return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1624:
1625: }
1626:
1627: /*----------------------------------------------------------------------------
1628: | Returns the result of subtracting the absolute values of the single-
1629: | precision floating-point values `a' and `b'. If `zSign' is 1, the
1630: | difference is negated before being returned. `zSign' is ignored if the
1631: | result is a NaN. The subtraction is performed according to the IEC/IEEE
1632: | Standard for Binary Floating-Point Arithmetic.
1633: *----------------------------------------------------------------------------*/
1634:
1635: static float32 subFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
1636: {
1637: int16 aExp, bExp, zExp;
1638: bits32 aSig, bSig, zSig;
1639: int16 expDiff;
1640:
1641: aSig = extractFloat32Frac( a );
1642: aExp = extractFloat32Exp( a );
1643: bSig = extractFloat32Frac( b );
1644: bExp = extractFloat32Exp( b );
1645: expDiff = aExp - bExp;
1646: aSig <<= 7;
1647: bSig <<= 7;
1648: if ( 0 < expDiff ) goto aExpBigger;
1649: if ( expDiff < 0 ) goto bExpBigger;
1650: if ( aExp == 0xFF ) {
1651: if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1652: float_raise( float_flag_invalid STATUS_VAR);
1653: return float32_default_nan;
1654: }
1655: if ( aExp == 0 ) {
1656: aExp = 1;
1657: bExp = 1;
1658: }
1659: if ( bSig < aSig ) goto aBigger;
1660: if ( aSig < bSig ) goto bBigger;
1661: return packFloat32( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
1662: bExpBigger:
1663: if ( bExp == 0xFF ) {
1664: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1665: return packFloat32( zSign ^ 1, 0xFF, 0 );
1666: }
1667: if ( aExp == 0 ) {
1668: ++expDiff;
1669: }
1670: else {
1671: aSig |= 0x40000000;
1672: }
1673: shift32RightJamming( aSig, - expDiff, &aSig );
1674: bSig |= 0x40000000;
1675: bBigger:
1676: zSig = bSig - aSig;
1677: zExp = bExp;
1678: zSign ^= 1;
1679: goto normalizeRoundAndPack;
1680: aExpBigger:
1681: if ( aExp == 0xFF ) {
1682: if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1683: return a;
1684: }
1685: if ( bExp == 0 ) {
1686: --expDiff;
1687: }
1688: else {
1689: bSig |= 0x40000000;
1690: }
1691: shift32RightJamming( bSig, expDiff, &bSig );
1692: aSig |= 0x40000000;
1693: aBigger:
1694: zSig = aSig - bSig;
1695: zExp = aExp;
1696: normalizeRoundAndPack:
1697: --zExp;
1698: return normalizeRoundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1699:
1700: }
1701:
1702: /*----------------------------------------------------------------------------
1703: | Returns the result of adding the single-precision floating-point values `a'
1704: | and `b'. The operation is performed according to the IEC/IEEE Standard for
1705: | Binary Floating-Point Arithmetic.
1706: *----------------------------------------------------------------------------*/
1707:
1708: float32 float32_add( float32 a, float32 b STATUS_PARAM )
1709: {
1710: flag aSign, bSign;
1711:
1712: aSign = extractFloat32Sign( a );
1713: bSign = extractFloat32Sign( b );
1714: if ( aSign == bSign ) {
1715: return addFloat32Sigs( a, b, aSign STATUS_VAR);
1716: }
1717: else {
1718: return subFloat32Sigs( a, b, aSign STATUS_VAR );
1719: }
1720:
1721: }
1722:
1723: /*----------------------------------------------------------------------------
1724: | Returns the result of subtracting the single-precision floating-point values
1725: | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1726: | for Binary Floating-Point Arithmetic.
1727: *----------------------------------------------------------------------------*/
1728:
1729: float32 float32_sub( float32 a, float32 b STATUS_PARAM )
1730: {
1731: flag aSign, bSign;
1732:
1733: aSign = extractFloat32Sign( a );
1734: bSign = extractFloat32Sign( b );
1735: if ( aSign == bSign ) {
1736: return subFloat32Sigs( a, b, aSign STATUS_VAR );
1737: }
1738: else {
1739: return addFloat32Sigs( a, b, aSign STATUS_VAR );
1740: }
1741:
1742: }
1743:
1744: /*----------------------------------------------------------------------------
1745: | Returns the result of multiplying the single-precision floating-point values
1746: | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
1747: | for Binary Floating-Point Arithmetic.
1748: *----------------------------------------------------------------------------*/
1749:
1750: float32 float32_mul( float32 a, float32 b STATUS_PARAM )
1751: {
1752: flag aSign, bSign, zSign;
1753: int16 aExp, bExp, zExp;
1754: bits32 aSig, bSig;
1755: bits64 zSig64;
1756: bits32 zSig;
1757:
1758: aSig = extractFloat32Frac( a );
1759: aExp = extractFloat32Exp( a );
1760: aSign = extractFloat32Sign( a );
1761: bSig = extractFloat32Frac( b );
1762: bExp = extractFloat32Exp( b );
1763: bSign = extractFloat32Sign( b );
1764: zSign = aSign ^ bSign;
1765: if ( aExp == 0xFF ) {
1766: if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1767: return propagateFloat32NaN( a, b STATUS_VAR );
1768: }
1769: if ( ( bExp | bSig ) == 0 ) {
1770: float_raise( float_flag_invalid STATUS_VAR);
1771: return float32_default_nan;
1772: }
1773: return packFloat32( zSign, 0xFF, 0 );
1774: }
1775: if ( bExp == 0xFF ) {
1776: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1777: if ( ( aExp | aSig ) == 0 ) {
1778: float_raise( float_flag_invalid STATUS_VAR);
1779: return float32_default_nan;
1780: }
1781: return packFloat32( zSign, 0xFF, 0 );
1782: }
1783: if ( aExp == 0 ) {
1784: if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1785: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1786: }
1787: if ( bExp == 0 ) {
1788: if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
1789: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1790: }
1791: zExp = aExp + bExp - 0x7F;
1792: aSig = ( aSig | 0x00800000 )<<7;
1793: bSig = ( bSig | 0x00800000 )<<8;
1794: shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
1795: zSig = zSig64;
1796: if ( 0 <= (sbits32) ( zSig<<1 ) ) {
1797: zSig <<= 1;
1798: --zExp;
1799: }
1800: return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1801:
1802: }
1803:
1804: /*----------------------------------------------------------------------------
1805: | Returns the result of dividing the single-precision floating-point value `a'
1806: | by the corresponding value `b'. The operation is performed according to the
1807: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1808: *----------------------------------------------------------------------------*/
1809:
1810: float32 float32_div( float32 a, float32 b STATUS_PARAM )
1811: {
1812: flag aSign, bSign, zSign;
1813: int16 aExp, bExp, zExp;
1814: bits32 aSig, bSig, zSig;
1815:
1816: aSig = extractFloat32Frac( a );
1817: aExp = extractFloat32Exp( a );
1818: aSign = extractFloat32Sign( a );
1819: bSig = extractFloat32Frac( b );
1820: bExp = extractFloat32Exp( b );
1821: bSign = extractFloat32Sign( b );
1822: zSign = aSign ^ bSign;
1823: if ( aExp == 0xFF ) {
1824: if ( aSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1825: if ( bExp == 0xFF ) {
1826: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1827: float_raise( float_flag_invalid STATUS_VAR);
1828: return float32_default_nan;
1829: }
1830: return packFloat32( zSign, 0xFF, 0 );
1831: }
1832: if ( bExp == 0xFF ) {
1833: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1834: return packFloat32( zSign, 0, 0 );
1835: }
1836: if ( bExp == 0 ) {
1837: if ( bSig == 0 ) {
1838: if ( ( aExp | aSig ) == 0 ) {
1839: float_raise( float_flag_invalid STATUS_VAR);
1840: return float32_default_nan;
1841: }
1842: float_raise( float_flag_divbyzero STATUS_VAR);
1843: return packFloat32( zSign, 0xFF, 0 );
1844: }
1845: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1846: }
1847: if ( aExp == 0 ) {
1848: if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1849: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1850: }
1851: zExp = aExp - bExp + 0x7D;
1852: aSig = ( aSig | 0x00800000 )<<7;
1853: bSig = ( bSig | 0x00800000 )<<8;
1854: if ( bSig <= ( aSig + aSig ) ) {
1855: aSig >>= 1;
1856: ++zExp;
1857: }
1858: zSig = ( ( (bits64) aSig )<<32 ) / bSig;
1859: if ( ( zSig & 0x3F ) == 0 ) {
1860: zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
1861: }
1862: return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
1863:
1864: }
1865:
1866: /*----------------------------------------------------------------------------
1867: | Returns the remainder of the single-precision floating-point value `a'
1868: | with respect to the corresponding value `b'. The operation is performed
1869: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1870: *----------------------------------------------------------------------------*/
1871:
1872: float32 float32_rem( float32 a, float32 b STATUS_PARAM )
1873: {
1874: flag aSign, bSign, zSign;
1875: int16 aExp, bExp, expDiff;
1876: bits32 aSig, bSig;
1877: bits32 q;
1878: bits64 aSig64, bSig64, q64;
1879: bits32 alternateASig;
1880: sbits32 sigMean;
1881:
1882: aSig = extractFloat32Frac( a );
1883: aExp = extractFloat32Exp( a );
1884: aSign = extractFloat32Sign( a );
1885: bSig = extractFloat32Frac( b );
1886: bExp = extractFloat32Exp( b );
1887: bSign = extractFloat32Sign( b );
1888: if ( aExp == 0xFF ) {
1889: if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1890: return propagateFloat32NaN( a, b STATUS_VAR );
1891: }
1892: float_raise( float_flag_invalid STATUS_VAR);
1893: return float32_default_nan;
1894: }
1895: if ( bExp == 0xFF ) {
1896: if ( bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
1897: return a;
1898: }
1899: if ( bExp == 0 ) {
1900: if ( bSig == 0 ) {
1901: float_raise( float_flag_invalid STATUS_VAR);
1902: return float32_default_nan;
1903: }
1904: normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1905: }
1906: if ( aExp == 0 ) {
1907: if ( aSig == 0 ) return a;
1908: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1909: }
1910: expDiff = aExp - bExp;
1911: aSig |= 0x00800000;
1912: bSig |= 0x00800000;
1913: if ( expDiff < 32 ) {
1914: aSig <<= 8;
1915: bSig <<= 8;
1916: if ( expDiff < 0 ) {
1917: if ( expDiff < -1 ) return a;
1918: aSig >>= 1;
1919: }
1920: q = ( bSig <= aSig );
1921: if ( q ) aSig -= bSig;
1922: if ( 0 < expDiff ) {
1923: q = ( ( (bits64) aSig )<<32 ) / bSig;
1924: q >>= 32 - expDiff;
1925: bSig >>= 2;
1926: aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
1927: }
1928: else {
1929: aSig >>= 2;
1930: bSig >>= 2;
1931: }
1932: }
1933: else {
1934: if ( bSig <= aSig ) aSig -= bSig;
1935: aSig64 = ( (bits64) aSig )<<40;
1936: bSig64 = ( (bits64) bSig )<<40;
1937: expDiff -= 64;
1938: while ( 0 < expDiff ) {
1939: q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1940: q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1941: aSig64 = - ( ( bSig * q64 )<<38 );
1942: expDiff -= 62;
1943: }
1944: expDiff += 64;
1945: q64 = estimateDiv128To64( aSig64, 0, bSig64 );
1946: q64 = ( 2 < q64 ) ? q64 - 2 : 0;
1947: q = q64>>( 64 - expDiff );
1948: bSig <<= 6;
1949: aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
1950: }
1951: do {
1952: alternateASig = aSig;
1953: ++q;
1954: aSig -= bSig;
1955: } while ( 0 <= (sbits32) aSig );
1956: sigMean = aSig + alternateASig;
1957: if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
1958: aSig = alternateASig;
1959: }
1960: zSign = ( (sbits32) aSig < 0 );
1961: if ( zSign ) aSig = - aSig;
1962: return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig STATUS_VAR );
1963:
1964: }
1965:
1966: /*----------------------------------------------------------------------------
1967: | Returns the square root of the single-precision floating-point value `a'.
1968: | The operation is performed according to the IEC/IEEE Standard for Binary
1969: | Floating-Point Arithmetic.
1970: *----------------------------------------------------------------------------*/
1971:
1972: float32 float32_sqrt( float32 a STATUS_PARAM )
1973: {
1974: flag aSign;
1975: int16 aExp, zExp;
1976: bits32 aSig, zSig;
1977: bits64 rem, term;
1978:
1979: aSig = extractFloat32Frac( a );
1980: aExp = extractFloat32Exp( a );
1981: aSign = extractFloat32Sign( a );
1982: if ( aExp == 0xFF ) {
1983: if ( aSig ) return propagateFloat32NaN( a, 0 STATUS_VAR );
1984: if ( ! aSign ) return a;
1985: float_raise( float_flag_invalid STATUS_VAR);
1986: return float32_default_nan;
1987: }
1988: if ( aSign ) {
1989: if ( ( aExp | aSig ) == 0 ) return a;
1990: float_raise( float_flag_invalid STATUS_VAR);
1991: return float32_default_nan;
1992: }
1993: if ( aExp == 0 ) {
1994: if ( aSig == 0 ) return 0;
1995: normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1996: }
1997: zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
1998: aSig = ( aSig | 0x00800000 )<<8;
1999: zSig = estimateSqrt32( aExp, aSig ) + 2;
2000: if ( ( zSig & 0x7F ) <= 5 ) {
2001: if ( zSig < 2 ) {
2002: zSig = 0x7FFFFFFF;
2003: goto roundAndPack;
2004: }
2005: aSig >>= aExp & 1;
2006: term = ( (bits64) zSig ) * zSig;
2007: rem = ( ( (bits64) aSig )<<32 ) - term;
2008: while ( (sbits64) rem < 0 ) {
2009: --zSig;
2010: rem += ( ( (bits64) zSig )<<1 ) | 1;
2011: }
2012: zSig |= ( rem != 0 );
2013: }
2014: shift32RightJamming( zSig, 1, &zSig );
2015: roundAndPack:
2016: return roundAndPackFloat32( 0, zExp, zSig STATUS_VAR );
2017:
2018: }
2019:
2020: /*----------------------------------------------------------------------------
2021: | Returns 1 if the single-precision floating-point value `a' is equal to
2022: | the corresponding value `b', and 0 otherwise. The comparison is performed
2023: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2024: *----------------------------------------------------------------------------*/
2025:
2026: flag float32_eq( float32 a, float32 b STATUS_PARAM )
2027: {
2028:
2029: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2030: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2031: ) {
2032: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2033: float_raise( float_flag_invalid STATUS_VAR);
2034: }
2035: return 0;
2036: }
2037: return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
2038:
2039: }
2040:
2041: /*----------------------------------------------------------------------------
2042: | Returns 1 if the single-precision floating-point value `a' is less than
2043: | or equal to the corresponding value `b', and 0 otherwise. The comparison
2044: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2045: | Arithmetic.
2046: *----------------------------------------------------------------------------*/
2047:
2048: flag float32_le( float32 a, float32 b STATUS_PARAM )
2049: {
2050: flag aSign, bSign;
2051:
2052: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2053: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2054: ) {
2055: float_raise( float_flag_invalid STATUS_VAR);
2056: return 0;
2057: }
2058: aSign = extractFloat32Sign( a );
2059: bSign = extractFloat32Sign( b );
2060: if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
2061: return ( a == b ) || ( aSign ^ ( a < b ) );
2062:
2063: }
2064:
2065: /*----------------------------------------------------------------------------
2066: | Returns 1 if the single-precision floating-point value `a' is less than
2067: | the corresponding value `b', and 0 otherwise. The comparison is performed
2068: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2069: *----------------------------------------------------------------------------*/
2070:
2071: flag float32_lt( float32 a, float32 b STATUS_PARAM )
2072: {
2073: flag aSign, bSign;
2074:
2075: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2076: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2077: ) {
2078: float_raise( float_flag_invalid STATUS_VAR);
2079: return 0;
2080: }
2081: aSign = extractFloat32Sign( a );
2082: bSign = extractFloat32Sign( b );
2083: if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
2084: return ( a != b ) && ( aSign ^ ( a < b ) );
2085:
2086: }
2087:
2088: /*----------------------------------------------------------------------------
2089: | Returns 1 if the single-precision floating-point value `a' is equal to
2090: | the corresponding value `b', and 0 otherwise. The invalid exception is
2091: | raised if either operand is a NaN. Otherwise, the comparison is performed
2092: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2093: *----------------------------------------------------------------------------*/
2094:
2095: flag float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
2096: {
2097:
2098: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2099: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2100: ) {
2101: float_raise( float_flag_invalid STATUS_VAR);
2102: return 0;
2103: }
2104: return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
2105:
2106: }
2107:
2108: /*----------------------------------------------------------------------------
2109: | Returns 1 if the single-precision floating-point value `a' is less than or
2110: | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
2111: | cause an exception. Otherwise, the comparison is performed according to the
2112: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2113: *----------------------------------------------------------------------------*/
2114:
2115: flag float32_le_quiet( float32 a, float32 b STATUS_PARAM )
2116: {
2117: flag aSign, bSign;
2118:
2119: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2120: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2121: ) {
2122: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2123: float_raise( float_flag_invalid STATUS_VAR);
2124: }
2125: return 0;
2126: }
2127: aSign = extractFloat32Sign( a );
2128: bSign = extractFloat32Sign( b );
2129: if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
2130: return ( a == b ) || ( aSign ^ ( a < b ) );
2131:
2132: }
2133:
2134: /*----------------------------------------------------------------------------
2135: | Returns 1 if the single-precision floating-point value `a' is less than
2136: | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
2137: | exception. Otherwise, the comparison is performed according to the IEC/IEEE
2138: | Standard for Binary Floating-Point Arithmetic.
2139: *----------------------------------------------------------------------------*/
2140:
2141: flag float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
2142: {
2143: flag aSign, bSign;
2144:
2145: if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
2146: || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
2147: ) {
2148: if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
2149: float_raise( float_flag_invalid STATUS_VAR);
2150: }
2151: return 0;
2152: }
2153: aSign = extractFloat32Sign( a );
2154: bSign = extractFloat32Sign( b );
2155: if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
2156: return ( a != b ) && ( aSign ^ ( a < b ) );
2157:
2158: }
2159:
2160: /*----------------------------------------------------------------------------
2161: | Returns the result of converting the double-precision floating-point value
2162: | `a' to the 32-bit two's complement integer format. The conversion is
2163: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2164: | Arithmetic---which means in particular that the conversion is rounded
2165: | according to the current rounding mode. If `a' is a NaN, the largest
2166: | positive integer is returned. Otherwise, if the conversion overflows, the
2167: | largest integer with the same sign as `a' is returned.
2168: *----------------------------------------------------------------------------*/
2169:
2170: int32 float64_to_int32( float64 a STATUS_PARAM )
2171: {
2172: flag aSign;
2173: int16 aExp, shiftCount;
2174: bits64 aSig;
2175:
2176: aSig = extractFloat64Frac( a );
2177: aExp = extractFloat64Exp( a );
2178: aSign = extractFloat64Sign( a );
2179: if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2180: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2181: shiftCount = 0x42C - aExp;
2182: if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
2183: return roundAndPackInt32( aSign, aSig STATUS_VAR );
2184:
2185: }
2186:
2187: /*----------------------------------------------------------------------------
2188: | Returns the result of converting the double-precision floating-point value
2189: | `a' to the 32-bit two's complement integer format. The conversion is
2190: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2191: | Arithmetic, except that the conversion is always rounded toward zero.
2192: | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2193: | the conversion overflows, the largest integer with the same sign as `a' is
2194: | returned.
2195: *----------------------------------------------------------------------------*/
2196:
2197: int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
2198: {
2199: flag aSign;
2200: int16 aExp, shiftCount;
2201: bits64 aSig, savedASig;
2202: int32 z;
2203:
2204: aSig = extractFloat64Frac( a );
2205: aExp = extractFloat64Exp( a );
2206: aSign = extractFloat64Sign( a );
2207: if ( 0x41E < aExp ) {
2208: if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
2209: goto invalid;
2210: }
2211: else if ( aExp < 0x3FF ) {
2212: if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2213: return 0;
2214: }
2215: aSig |= LIT64( 0x0010000000000000 );
2216: shiftCount = 0x433 - aExp;
2217: savedASig = aSig;
2218: aSig >>= shiftCount;
2219: z = aSig;
2220: if ( aSign ) z = - z;
2221: if ( ( z < 0 ) ^ aSign ) {
2222: invalid:
2223: float_raise( float_flag_invalid STATUS_VAR);
2224: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
2225: }
2226: if ( ( aSig<<shiftCount ) != savedASig ) {
2227: STATUS(float_exception_flags) |= float_flag_inexact;
2228: }
2229: return z;
2230:
2231: }
2232:
2233: /*----------------------------------------------------------------------------
2234: | Returns the result of converting the double-precision floating-point value
2235: | `a' to the 64-bit two's complement integer format. The conversion is
2236: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2237: | Arithmetic---which means in particular that the conversion is rounded
2238: | according to the current rounding mode. If `a' is a NaN, the largest
2239: | positive integer is returned. Otherwise, if the conversion overflows, the
2240: | largest integer with the same sign as `a' is returned.
2241: *----------------------------------------------------------------------------*/
2242:
2243: int64 float64_to_int64( float64 a STATUS_PARAM )
2244: {
2245: flag aSign;
2246: int16 aExp, shiftCount;
2247: bits64 aSig, aSigExtra;
2248:
2249: aSig = extractFloat64Frac( a );
2250: aExp = extractFloat64Exp( a );
2251: aSign = extractFloat64Sign( a );
2252: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2253: shiftCount = 0x433 - aExp;
2254: if ( shiftCount <= 0 ) {
2255: if ( 0x43E < aExp ) {
2256: float_raise( float_flag_invalid STATUS_VAR);
2257: if ( ! aSign
2258: || ( ( aExp == 0x7FF )
2259: && ( aSig != LIT64( 0x0010000000000000 ) ) )
2260: ) {
2261: return LIT64( 0x7FFFFFFFFFFFFFFF );
2262: }
2263: return (sbits64) LIT64( 0x8000000000000000 );
2264: }
2265: aSigExtra = 0;
2266: aSig <<= - shiftCount;
2267: }
2268: else {
2269: shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
2270: }
2271: return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
2272:
2273: }
2274:
2275: /*----------------------------------------------------------------------------
2276: | Returns the result of converting the double-precision floating-point value
2277: | `a' to the 64-bit two's complement integer format. The conversion is
2278: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2279: | Arithmetic, except that the conversion is always rounded toward zero.
2280: | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
2281: | the conversion overflows, the largest integer with the same sign as `a' is
2282: | returned.
2283: *----------------------------------------------------------------------------*/
2284:
2285: int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
2286: {
2287: flag aSign;
2288: int16 aExp, shiftCount;
2289: bits64 aSig;
2290: int64 z;
2291:
2292: aSig = extractFloat64Frac( a );
2293: aExp = extractFloat64Exp( a );
2294: aSign = extractFloat64Sign( a );
2295: if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
2296: shiftCount = aExp - 0x433;
2297: if ( 0 <= shiftCount ) {
2298: if ( 0x43E <= aExp ) {
2299: if ( a != LIT64( 0xC3E0000000000000 ) ) {
2300: float_raise( float_flag_invalid STATUS_VAR);
2301: if ( ! aSign
2302: || ( ( aExp == 0x7FF )
2303: && ( aSig != LIT64( 0x0010000000000000 ) ) )
2304: ) {
2305: return LIT64( 0x7FFFFFFFFFFFFFFF );
2306: }
2307: }
2308: return (sbits64) LIT64( 0x8000000000000000 );
2309: }
2310: z = aSig<<shiftCount;
2311: }
2312: else {
2313: if ( aExp < 0x3FE ) {
2314: if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
2315: return 0;
2316: }
2317: z = aSig>>( - shiftCount );
2318: if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
2319: STATUS(float_exception_flags) |= float_flag_inexact;
2320: }
2321: }
2322: if ( aSign ) z = - z;
2323: return z;
2324:
2325: }
2326:
2327: /*----------------------------------------------------------------------------
2328: | Returns the result of converting the double-precision floating-point value
2329: | `a' to the single-precision floating-point format. The conversion is
2330: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2331: | Arithmetic.
2332: *----------------------------------------------------------------------------*/
2333:
2334: float32 float64_to_float32( float64 a STATUS_PARAM )
2335: {
2336: flag aSign;
2337: int16 aExp;
2338: bits64 aSig;
2339: bits32 zSig;
2340:
2341: aSig = extractFloat64Frac( a );
2342: aExp = extractFloat64Exp( a );
2343: aSign = extractFloat64Sign( a );
2344: if ( aExp == 0x7FF ) {
2345: if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR ) );
2346: return packFloat32( aSign, 0xFF, 0 );
2347: }
2348: shift64RightJamming( aSig, 22, &aSig );
2349: zSig = aSig;
2350: if ( aExp || zSig ) {
2351: zSig |= 0x40000000;
2352: aExp -= 0x381;
2353: }
2354: return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
2355:
2356: }
2357:
2358: #ifdef FLOATX80
2359:
2360: /*----------------------------------------------------------------------------
2361: | Returns the result of converting the double-precision floating-point value
2362: | `a' to the extended double-precision floating-point format. The conversion
2363: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
2364: | Arithmetic.
2365: *----------------------------------------------------------------------------*/
2366:
2367: floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
2368: {
2369: flag aSign;
2370: int16 aExp;
2371: bits64 aSig;
2372:
2373: aSig = extractFloat64Frac( a );
2374: aExp = extractFloat64Exp( a );
2375: aSign = extractFloat64Sign( a );
2376: if ( aExp == 0x7FF ) {
2377: if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) );
2378: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
2379: }
2380: if ( aExp == 0 ) {
2381: if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
2382: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2383: }
2384: return
2385: packFloatx80(
2386: aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
2387:
2388: }
2389:
2390: #endif
2391:
2392: #ifdef FLOAT128
2393:
2394: /*----------------------------------------------------------------------------
2395: | Returns the result of converting the double-precision floating-point value
2396: | `a' to the quadruple-precision floating-point format. The conversion is
2397: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2398: | Arithmetic.
2399: *----------------------------------------------------------------------------*/
2400:
2401: float128 float64_to_float128( float64 a STATUS_PARAM )
2402: {
2403: flag aSign;
2404: int16 aExp;
2405: bits64 aSig, zSig0, zSig1;
2406:
2407: aSig = extractFloat64Frac( a );
2408: aExp = extractFloat64Exp( a );
2409: aSign = extractFloat64Sign( a );
2410: if ( aExp == 0x7FF ) {
2411: if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR ) );
2412: return packFloat128( aSign, 0x7FFF, 0, 0 );
2413: }
2414: if ( aExp == 0 ) {
2415: if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
2416: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2417: --aExp;
2418: }
2419: shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
2420: return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
2421:
2422: }
2423:
2424: #endif
2425:
2426: /*----------------------------------------------------------------------------
2427: | Rounds the double-precision floating-point value `a' to an integer, and
2428: | returns the result as a double-precision floating-point value. The
2429: | operation is performed according to the IEC/IEEE Standard for Binary
2430: | Floating-Point Arithmetic.
2431: *----------------------------------------------------------------------------*/
2432:
2433: float64 float64_round_to_int( float64 a STATUS_PARAM )
2434: {
2435: flag aSign;
2436: int16 aExp;
2437: bits64 lastBitMask, roundBitsMask;
2438: int8 roundingMode;
2439: float64 z;
2440:
2441: aExp = extractFloat64Exp( a );
2442: if ( 0x433 <= aExp ) {
2443: if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
2444: return propagateFloat64NaN( a, a STATUS_VAR );
2445: }
2446: return a;
2447: }
2448: if ( aExp < 0x3FF ) {
2449: if ( (bits64) ( a<<1 ) == 0 ) return a;
2450: STATUS(float_exception_flags) |= float_flag_inexact;
2451: aSign = extractFloat64Sign( a );
2452: switch ( STATUS(float_rounding_mode) ) {
2453: case float_round_nearest_even:
2454: if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
2455: return packFloat64( aSign, 0x3FF, 0 );
2456: }
2457: break;
2458: case float_round_down:
2459: return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
2460: case float_round_up:
2461: return
2462: aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
2463: }
2464: return packFloat64( aSign, 0, 0 );
2465: }
2466: lastBitMask = 1;
2467: lastBitMask <<= 0x433 - aExp;
2468: roundBitsMask = lastBitMask - 1;
2469: z = a;
2470: roundingMode = STATUS(float_rounding_mode);
2471: if ( roundingMode == float_round_nearest_even ) {
2472: z += lastBitMask>>1;
2473: if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
2474: }
2475: else if ( roundingMode != float_round_to_zero ) {
2476: if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
2477: z += roundBitsMask;
2478: }
2479: }
2480: z &= ~ roundBitsMask;
2481: if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
2482: return z;
2483:
2484: }
2485:
2486: /*----------------------------------------------------------------------------
2487: | Returns the result of adding the absolute values of the double-precision
2488: | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
2489: | before being returned. `zSign' is ignored if the result is a NaN.
2490: | The addition is performed according to the IEC/IEEE Standard for Binary
2491: | Floating-Point Arithmetic.
2492: *----------------------------------------------------------------------------*/
2493:
2494: static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2495: {
2496: int16 aExp, bExp, zExp;
2497: bits64 aSig, bSig, zSig;
2498: int16 expDiff;
2499:
2500: aSig = extractFloat64Frac( a );
2501: aExp = extractFloat64Exp( a );
2502: bSig = extractFloat64Frac( b );
2503: bExp = extractFloat64Exp( b );
2504: expDiff = aExp - bExp;
2505: aSig <<= 9;
2506: bSig <<= 9;
2507: if ( 0 < expDiff ) {
2508: if ( aExp == 0x7FF ) {
2509: if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2510: return a;
2511: }
2512: if ( bExp == 0 ) {
2513: --expDiff;
2514: }
2515: else {
2516: bSig |= LIT64( 0x2000000000000000 );
2517: }
2518: shift64RightJamming( bSig, expDiff, &bSig );
2519: zExp = aExp;
2520: }
2521: else if ( expDiff < 0 ) {
2522: if ( bExp == 0x7FF ) {
2523: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2524: return packFloat64( zSign, 0x7FF, 0 );
2525: }
2526: if ( aExp == 0 ) {
2527: ++expDiff;
2528: }
2529: else {
2530: aSig |= LIT64( 0x2000000000000000 );
2531: }
2532: shift64RightJamming( aSig, - expDiff, &aSig );
2533: zExp = bExp;
2534: }
2535: else {
2536: if ( aExp == 0x7FF ) {
2537: if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2538: return a;
2539: }
2540: if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
2541: zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
2542: zExp = aExp;
2543: goto roundAndPack;
2544: }
2545: aSig |= LIT64( 0x2000000000000000 );
2546: zSig = ( aSig + bSig )<<1;
2547: --zExp;
2548: if ( (sbits64) zSig < 0 ) {
2549: zSig = aSig + bSig;
2550: ++zExp;
2551: }
2552: roundAndPack:
2553: return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2554:
2555: }
2556:
2557: /*----------------------------------------------------------------------------
2558: | Returns the result of subtracting the absolute values of the double-
2559: | precision floating-point values `a' and `b'. If `zSign' is 1, the
2560: | difference is negated before being returned. `zSign' is ignored if the
2561: | result is a NaN. The subtraction is performed according to the IEC/IEEE
2562: | Standard for Binary Floating-Point Arithmetic.
2563: *----------------------------------------------------------------------------*/
2564:
2565: static float64 subFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
2566: {
2567: int16 aExp, bExp, zExp;
2568: bits64 aSig, bSig, zSig;
2569: int16 expDiff;
2570:
2571: aSig = extractFloat64Frac( a );
2572: aExp = extractFloat64Exp( a );
2573: bSig = extractFloat64Frac( b );
2574: bExp = extractFloat64Exp( b );
2575: expDiff = aExp - bExp;
2576: aSig <<= 10;
2577: bSig <<= 10;
2578: if ( 0 < expDiff ) goto aExpBigger;
2579: if ( expDiff < 0 ) goto bExpBigger;
2580: if ( aExp == 0x7FF ) {
2581: if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2582: float_raise( float_flag_invalid STATUS_VAR);
2583: return float64_default_nan;
2584: }
2585: if ( aExp == 0 ) {
2586: aExp = 1;
2587: bExp = 1;
2588: }
2589: if ( bSig < aSig ) goto aBigger;
2590: if ( aSig < bSig ) goto bBigger;
2591: return packFloat64( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
2592: bExpBigger:
2593: if ( bExp == 0x7FF ) {
2594: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2595: return packFloat64( zSign ^ 1, 0x7FF, 0 );
2596: }
2597: if ( aExp == 0 ) {
2598: ++expDiff;
2599: }
2600: else {
2601: aSig |= LIT64( 0x4000000000000000 );
2602: }
2603: shift64RightJamming( aSig, - expDiff, &aSig );
2604: bSig |= LIT64( 0x4000000000000000 );
2605: bBigger:
2606: zSig = bSig - aSig;
2607: zExp = bExp;
2608: zSign ^= 1;
2609: goto normalizeRoundAndPack;
2610: aExpBigger:
2611: if ( aExp == 0x7FF ) {
2612: if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2613: return a;
2614: }
2615: if ( bExp == 0 ) {
2616: --expDiff;
2617: }
2618: else {
2619: bSig |= LIT64( 0x4000000000000000 );
2620: }
2621: shift64RightJamming( bSig, expDiff, &bSig );
2622: aSig |= LIT64( 0x4000000000000000 );
2623: aBigger:
2624: zSig = aSig - bSig;
2625: zExp = aExp;
2626: normalizeRoundAndPack:
2627: --zExp;
2628: return normalizeRoundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2629:
2630: }
2631:
2632: /*----------------------------------------------------------------------------
2633: | Returns the result of adding the double-precision floating-point values `a'
2634: | and `b'. The operation is performed according to the IEC/IEEE Standard for
2635: | Binary Floating-Point Arithmetic.
2636: *----------------------------------------------------------------------------*/
2637:
2638: float64 float64_add( float64 a, float64 b STATUS_PARAM )
2639: {
2640: flag aSign, bSign;
2641:
2642: aSign = extractFloat64Sign( a );
2643: bSign = extractFloat64Sign( b );
2644: if ( aSign == bSign ) {
2645: return addFloat64Sigs( a, b, aSign STATUS_VAR );
2646: }
2647: else {
2648: return subFloat64Sigs( a, b, aSign STATUS_VAR );
2649: }
2650:
2651: }
2652:
2653: /*----------------------------------------------------------------------------
2654: | Returns the result of subtracting the double-precision floating-point values
2655: | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2656: | for Binary Floating-Point Arithmetic.
2657: *----------------------------------------------------------------------------*/
2658:
2659: float64 float64_sub( float64 a, float64 b STATUS_PARAM )
2660: {
2661: flag aSign, bSign;
2662:
2663: aSign = extractFloat64Sign( a );
2664: bSign = extractFloat64Sign( b );
2665: if ( aSign == bSign ) {
2666: return subFloat64Sigs( a, b, aSign STATUS_VAR );
2667: }
2668: else {
2669: return addFloat64Sigs( a, b, aSign STATUS_VAR );
2670: }
2671:
2672: }
2673:
2674: /*----------------------------------------------------------------------------
2675: | Returns the result of multiplying the double-precision floating-point values
2676: | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
2677: | for Binary Floating-Point Arithmetic.
2678: *----------------------------------------------------------------------------*/
2679:
2680: float64 float64_mul( float64 a, float64 b STATUS_PARAM )
2681: {
2682: flag aSign, bSign, zSign;
2683: int16 aExp, bExp, zExp;
2684: bits64 aSig, bSig, zSig0, zSig1;
2685:
2686: aSig = extractFloat64Frac( a );
2687: aExp = extractFloat64Exp( a );
2688: aSign = extractFloat64Sign( a );
2689: bSig = extractFloat64Frac( b );
2690: bExp = extractFloat64Exp( b );
2691: bSign = extractFloat64Sign( b );
2692: zSign = aSign ^ bSign;
2693: if ( aExp == 0x7FF ) {
2694: if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2695: return propagateFloat64NaN( a, b STATUS_VAR );
2696: }
2697: if ( ( bExp | bSig ) == 0 ) {
2698: float_raise( float_flag_invalid STATUS_VAR);
2699: return float64_default_nan;
2700: }
2701: return packFloat64( zSign, 0x7FF, 0 );
2702: }
2703: if ( bExp == 0x7FF ) {
2704: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2705: if ( ( aExp | aSig ) == 0 ) {
2706: float_raise( float_flag_invalid STATUS_VAR);
2707: return float64_default_nan;
2708: }
2709: return packFloat64( zSign, 0x7FF, 0 );
2710: }
2711: if ( aExp == 0 ) {
2712: if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2713: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2714: }
2715: if ( bExp == 0 ) {
2716: if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
2717: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2718: }
2719: zExp = aExp + bExp - 0x3FF;
2720: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2721: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2722: mul64To128( aSig, bSig, &zSig0, &zSig1 );
2723: zSig0 |= ( zSig1 != 0 );
2724: if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
2725: zSig0 <<= 1;
2726: --zExp;
2727: }
2728: return roundAndPackFloat64( zSign, zExp, zSig0 STATUS_VAR );
2729:
2730: }
2731:
2732: /*----------------------------------------------------------------------------
2733: | Returns the result of dividing the double-precision floating-point value `a'
2734: | by the corresponding value `b'. The operation is performed according to
2735: | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2736: *----------------------------------------------------------------------------*/
2737:
2738: float64 float64_div( float64 a, float64 b STATUS_PARAM )
2739: {
2740: flag aSign, bSign, zSign;
2741: int16 aExp, bExp, zExp;
2742: bits64 aSig, bSig, zSig;
2743: bits64 rem0, rem1;
2744: bits64 term0, term1;
2745:
2746: aSig = extractFloat64Frac( a );
2747: aExp = extractFloat64Exp( a );
2748: aSign = extractFloat64Sign( a );
2749: bSig = extractFloat64Frac( b );
2750: bExp = extractFloat64Exp( b );
2751: bSign = extractFloat64Sign( b );
2752: zSign = aSign ^ bSign;
2753: if ( aExp == 0x7FF ) {
2754: if ( aSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2755: if ( bExp == 0x7FF ) {
2756: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2757: float_raise( float_flag_invalid STATUS_VAR);
2758: return float64_default_nan;
2759: }
2760: return packFloat64( zSign, 0x7FF, 0 );
2761: }
2762: if ( bExp == 0x7FF ) {
2763: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2764: return packFloat64( zSign, 0, 0 );
2765: }
2766: if ( bExp == 0 ) {
2767: if ( bSig == 0 ) {
2768: if ( ( aExp | aSig ) == 0 ) {
2769: float_raise( float_flag_invalid STATUS_VAR);
2770: return float64_default_nan;
2771: }
2772: float_raise( float_flag_divbyzero STATUS_VAR);
2773: return packFloat64( zSign, 0x7FF, 0 );
2774: }
2775: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2776: }
2777: if ( aExp == 0 ) {
2778: if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
2779: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2780: }
2781: zExp = aExp - bExp + 0x3FD;
2782: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
2783: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2784: if ( bSig <= ( aSig + aSig ) ) {
2785: aSig >>= 1;
2786: ++zExp;
2787: }
2788: zSig = estimateDiv128To64( aSig, 0, bSig );
2789: if ( ( zSig & 0x1FF ) <= 2 ) {
2790: mul64To128( bSig, zSig, &term0, &term1 );
2791: sub128( aSig, 0, term0, term1, &rem0, &rem1 );
2792: while ( (sbits64) rem0 < 0 ) {
2793: --zSig;
2794: add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
2795: }
2796: zSig |= ( rem1 != 0 );
2797: }
2798: return roundAndPackFloat64( zSign, zExp, zSig STATUS_VAR );
2799:
2800: }
2801:
2802: /*----------------------------------------------------------------------------
2803: | Returns the remainder of the double-precision floating-point value `a'
2804: | with respect to the corresponding value `b'. The operation is performed
2805: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2806: *----------------------------------------------------------------------------*/
2807:
2808: float64 float64_rem( float64 a, float64 b STATUS_PARAM )
2809: {
2810: flag aSign, bSign, zSign;
2811: int16 aExp, bExp, expDiff;
2812: bits64 aSig, bSig;
2813: bits64 q, alternateASig;
2814: sbits64 sigMean;
2815:
2816: aSig = extractFloat64Frac( a );
2817: aExp = extractFloat64Exp( a );
2818: aSign = extractFloat64Sign( a );
2819: bSig = extractFloat64Frac( b );
2820: bExp = extractFloat64Exp( b );
2821: bSign = extractFloat64Sign( b );
2822: if ( aExp == 0x7FF ) {
2823: if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
2824: return propagateFloat64NaN( a, b STATUS_VAR );
2825: }
2826: float_raise( float_flag_invalid STATUS_VAR);
2827: return float64_default_nan;
2828: }
2829: if ( bExp == 0x7FF ) {
2830: if ( bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
2831: return a;
2832: }
2833: if ( bExp == 0 ) {
2834: if ( bSig == 0 ) {
2835: float_raise( float_flag_invalid STATUS_VAR);
2836: return float64_default_nan;
2837: }
2838: normalizeFloat64Subnormal( bSig, &bExp, &bSig );
2839: }
2840: if ( aExp == 0 ) {
2841: if ( aSig == 0 ) return a;
2842: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2843: }
2844: expDiff = aExp - bExp;
2845: aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
2846: bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
2847: if ( expDiff < 0 ) {
2848: if ( expDiff < -1 ) return a;
2849: aSig >>= 1;
2850: }
2851: q = ( bSig <= aSig );
2852: if ( q ) aSig -= bSig;
2853: expDiff -= 64;
2854: while ( 0 < expDiff ) {
2855: q = estimateDiv128To64( aSig, 0, bSig );
2856: q = ( 2 < q ) ? q - 2 : 0;
2857: aSig = - ( ( bSig>>2 ) * q );
2858: expDiff -= 62;
2859: }
2860: expDiff += 64;
2861: if ( 0 < expDiff ) {
2862: q = estimateDiv128To64( aSig, 0, bSig );
2863: q = ( 2 < q ) ? q - 2 : 0;
2864: q >>= 64 - expDiff;
2865: bSig >>= 2;
2866: aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
2867: }
2868: else {
2869: aSig >>= 2;
2870: bSig >>= 2;
2871: }
2872: do {
2873: alternateASig = aSig;
2874: ++q;
2875: aSig -= bSig;
2876: } while ( 0 <= (sbits64) aSig );
2877: sigMean = aSig + alternateASig;
2878: if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
2879: aSig = alternateASig;
2880: }
2881: zSign = ( (sbits64) aSig < 0 );
2882: if ( zSign ) aSig = - aSig;
2883: return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig STATUS_VAR );
2884:
2885: }
2886:
2887: /*----------------------------------------------------------------------------
2888: | Returns the square root of the double-precision floating-point value `a'.
2889: | The operation is performed according to the IEC/IEEE Standard for Binary
2890: | Floating-Point Arithmetic.
2891: *----------------------------------------------------------------------------*/
2892:
2893: float64 float64_sqrt( float64 a STATUS_PARAM )
2894: {
2895: flag aSign;
2896: int16 aExp, zExp;
2897: bits64 aSig, zSig, doubleZSig;
2898: bits64 rem0, rem1, term0, term1;
2899:
2900: aSig = extractFloat64Frac( a );
2901: aExp = extractFloat64Exp( a );
2902: aSign = extractFloat64Sign( a );
2903: if ( aExp == 0x7FF ) {
2904: if ( aSig ) return propagateFloat64NaN( a, a STATUS_VAR );
2905: if ( ! aSign ) return a;
2906: float_raise( float_flag_invalid STATUS_VAR);
2907: return float64_default_nan;
2908: }
2909: if ( aSign ) {
2910: if ( ( aExp | aSig ) == 0 ) return a;
2911: float_raise( float_flag_invalid STATUS_VAR);
2912: return float64_default_nan;
2913: }
2914: if ( aExp == 0 ) {
2915: if ( aSig == 0 ) return 0;
2916: normalizeFloat64Subnormal( aSig, &aExp, &aSig );
2917: }
2918: zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
2919: aSig |= LIT64( 0x0010000000000000 );
2920: zSig = estimateSqrt32( aExp, aSig>>21 );
2921: aSig <<= 9 - ( aExp & 1 );
2922: zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
2923: if ( ( zSig & 0x1FF ) <= 5 ) {
2924: doubleZSig = zSig<<1;
2925: mul64To128( zSig, zSig, &term0, &term1 );
2926: sub128( aSig, 0, term0, term1, &rem0, &rem1 );
2927: while ( (sbits64) rem0 < 0 ) {
2928: --zSig;
2929: doubleZSig -= 2;
2930: add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
2931: }
2932: zSig |= ( ( rem0 | rem1 ) != 0 );
2933: }
2934: return roundAndPackFloat64( 0, zExp, zSig STATUS_VAR );
2935:
2936: }
2937:
2938: /*----------------------------------------------------------------------------
2939: | Returns 1 if the double-precision floating-point value `a' is equal to the
2940: | corresponding value `b', and 0 otherwise. The comparison is performed
2941: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2942: *----------------------------------------------------------------------------*/
2943:
2944: flag float64_eq( float64 a, float64 b STATUS_PARAM )
2945: {
2946:
2947: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
2948: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
2949: ) {
2950: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
2951: float_raise( float_flag_invalid STATUS_VAR);
2952: }
2953: return 0;
2954: }
2955: return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
2956:
2957: }
2958:
2959: /*----------------------------------------------------------------------------
2960: | Returns 1 if the double-precision floating-point value `a' is less than or
2961: | equal to the corresponding value `b', and 0 otherwise. The comparison is
2962: | performed according to the IEC/IEEE Standard for Binary Floating-Point
2963: | Arithmetic.
2964: *----------------------------------------------------------------------------*/
2965:
2966: flag float64_le( float64 a, float64 b STATUS_PARAM )
2967: {
2968: flag aSign, bSign;
2969:
2970: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
2971: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
2972: ) {
2973: float_raise( float_flag_invalid STATUS_VAR);
2974: return 0;
2975: }
2976: aSign = extractFloat64Sign( a );
2977: bSign = extractFloat64Sign( b );
2978: if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
2979: return ( a == b ) || ( aSign ^ ( a < b ) );
2980:
2981: }
2982:
2983: /*----------------------------------------------------------------------------
2984: | Returns 1 if the double-precision floating-point value `a' is less than
2985: | the corresponding value `b', and 0 otherwise. The comparison is performed
2986: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2987: *----------------------------------------------------------------------------*/
2988:
2989: flag float64_lt( float64 a, float64 b STATUS_PARAM )
2990: {
2991: flag aSign, bSign;
2992:
2993: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
2994: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
2995: ) {
2996: float_raise( float_flag_invalid STATUS_VAR);
2997: return 0;
2998: }
2999: aSign = extractFloat64Sign( a );
3000: bSign = extractFloat64Sign( b );
3001: if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
3002: return ( a != b ) && ( aSign ^ ( a < b ) );
3003:
3004: }
3005:
3006: /*----------------------------------------------------------------------------
3007: | Returns 1 if the double-precision floating-point value `a' is equal to the
3008: | corresponding value `b', and 0 otherwise. The invalid exception is raised
3009: | if either operand is a NaN. Otherwise, the comparison is performed
3010: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3011: *----------------------------------------------------------------------------*/
3012:
3013: flag float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
3014: {
3015:
3016: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3017: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3018: ) {
3019: float_raise( float_flag_invalid STATUS_VAR);
3020: return 0;
3021: }
3022: return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
3023:
3024: }
3025:
3026: /*----------------------------------------------------------------------------
3027: | Returns 1 if the double-precision floating-point value `a' is less than or
3028: | equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
3029: | cause an exception. Otherwise, the comparison is performed according to the
3030: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3031: *----------------------------------------------------------------------------*/
3032:
3033: flag float64_le_quiet( float64 a, float64 b STATUS_PARAM )
3034: {
3035: flag aSign, bSign;
3036:
3037: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3038: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3039: ) {
3040: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3041: float_raise( float_flag_invalid STATUS_VAR);
3042: }
3043: return 0;
3044: }
3045: aSign = extractFloat64Sign( a );
3046: bSign = extractFloat64Sign( b );
3047: if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
3048: return ( a == b ) || ( aSign ^ ( a < b ) );
3049:
3050: }
3051:
3052: /*----------------------------------------------------------------------------
3053: | Returns 1 if the double-precision floating-point value `a' is less than
3054: | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
3055: | exception. Otherwise, the comparison is performed according to the IEC/IEEE
3056: | Standard for Binary Floating-Point Arithmetic.
3057: *----------------------------------------------------------------------------*/
3058:
3059: flag float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
3060: {
3061: flag aSign, bSign;
3062:
3063: if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
3064: || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
3065: ) {
3066: if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
3067: float_raise( float_flag_invalid STATUS_VAR);
3068: }
3069: return 0;
3070: }
3071: aSign = extractFloat64Sign( a );
3072: bSign = extractFloat64Sign( b );
3073: if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
3074: return ( a != b ) && ( aSign ^ ( a < b ) );
3075:
3076: }
3077:
3078: #ifdef FLOATX80
3079:
3080: /*----------------------------------------------------------------------------
3081: | Returns the result of converting the extended double-precision floating-
3082: | point value `a' to the 32-bit two's complement integer format. The
3083: | conversion is performed according to the IEC/IEEE Standard for Binary
3084: | Floating-Point Arithmetic---which means in particular that the conversion
3085: | is rounded according to the current rounding mode. If `a' is a NaN, the
3086: | largest positive integer is returned. Otherwise, if the conversion
3087: | overflows, the largest integer with the same sign as `a' is returned.
3088: *----------------------------------------------------------------------------*/
3089:
3090: int32 floatx80_to_int32( floatx80 a STATUS_PARAM )
3091: {
3092: flag aSign;
3093: int32 aExp, shiftCount;
3094: bits64 aSig;
3095:
3096: aSig = extractFloatx80Frac( a );
3097: aExp = extractFloatx80Exp( a );
3098: aSign = extractFloatx80Sign( a );
3099: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3100: shiftCount = 0x4037 - aExp;
3101: if ( shiftCount <= 0 ) shiftCount = 1;
3102: shift64RightJamming( aSig, shiftCount, &aSig );
3103: return roundAndPackInt32( aSign, aSig STATUS_VAR );
3104:
3105: }
3106:
3107: /*----------------------------------------------------------------------------
3108: | Returns the result of converting the extended double-precision floating-
3109: | point value `a' to the 32-bit two's complement integer format. The
3110: | conversion is performed according to the IEC/IEEE Standard for Binary
3111: | Floating-Point Arithmetic, except that the conversion is always rounded
3112: | toward zero. If `a' is a NaN, the largest positive integer is returned.
3113: | Otherwise, if the conversion overflows, the largest integer with the same
3114: | sign as `a' is returned.
3115: *----------------------------------------------------------------------------*/
3116:
3117: int32 floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM )
3118: {
3119: flag aSign;
3120: int32 aExp, shiftCount;
3121: bits64 aSig, savedASig;
3122: int32 z;
3123:
3124: aSig = extractFloatx80Frac( a );
3125: aExp = extractFloatx80Exp( a );
3126: aSign = extractFloatx80Sign( a );
3127: if ( 0x401E < aExp ) {
3128: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
3129: goto invalid;
3130: }
3131: else if ( aExp < 0x3FFF ) {
3132: if ( aExp || aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3133: return 0;
3134: }
3135: shiftCount = 0x403E - aExp;
3136: savedASig = aSig;
3137: aSig >>= shiftCount;
3138: z = aSig;
3139: if ( aSign ) z = - z;
3140: if ( ( z < 0 ) ^ aSign ) {
3141: invalid:
3142: float_raise( float_flag_invalid STATUS_VAR);
3143: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
3144: }
3145: if ( ( aSig<<shiftCount ) != savedASig ) {
3146: STATUS(float_exception_flags) |= float_flag_inexact;
3147: }
3148: return z;
3149:
3150: }
3151:
3152: /*----------------------------------------------------------------------------
3153: | Returns the result of converting the extended double-precision floating-
3154: | point value `a' to the 64-bit two's complement integer format. The
3155: | conversion is performed according to the IEC/IEEE Standard for Binary
3156: | Floating-Point Arithmetic---which means in particular that the conversion
3157: | is rounded according to the current rounding mode. If `a' is a NaN,
3158: | the largest positive integer is returned. Otherwise, if the conversion
3159: | overflows, the largest integer with the same sign as `a' is returned.
3160: *----------------------------------------------------------------------------*/
3161:
3162: int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
3163: {
3164: flag aSign;
3165: int32 aExp, shiftCount;
3166: bits64 aSig, aSigExtra;
3167:
3168: aSig = extractFloatx80Frac( a );
3169: aExp = extractFloatx80Exp( a );
3170: aSign = extractFloatx80Sign( a );
3171: shiftCount = 0x403E - aExp;
3172: if ( shiftCount <= 0 ) {
3173: if ( shiftCount ) {
3174: float_raise( float_flag_invalid STATUS_VAR);
3175: if ( ! aSign
3176: || ( ( aExp == 0x7FFF )
3177: && ( aSig != LIT64( 0x8000000000000000 ) ) )
3178: ) {
3179: return LIT64( 0x7FFFFFFFFFFFFFFF );
3180: }
3181: return (sbits64) LIT64( 0x8000000000000000 );
3182: }
3183: aSigExtra = 0;
3184: }
3185: else {
3186: shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
3187: }
3188: return roundAndPackInt64( aSign, aSig, aSigExtra STATUS_VAR );
3189:
3190: }
3191:
3192: /*----------------------------------------------------------------------------
3193: | Returns the result of converting the extended double-precision floating-
3194: | point value `a' to the 64-bit two's complement integer format. The
3195: | conversion is performed according to the IEC/IEEE Standard for Binary
3196: | Floating-Point Arithmetic, except that the conversion is always rounded
3197: | toward zero. If `a' is a NaN, the largest positive integer is returned.
3198: | Otherwise, if the conversion overflows, the largest integer with the same
3199: | sign as `a' is returned.
3200: *----------------------------------------------------------------------------*/
3201:
3202: int64 floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM )
3203: {
3204: flag aSign;
3205: int32 aExp, shiftCount;
3206: bits64 aSig;
3207: int64 z;
3208:
3209: aSig = extractFloatx80Frac( a );
3210: aExp = extractFloatx80Exp( a );
3211: aSign = extractFloatx80Sign( a );
3212: shiftCount = aExp - 0x403E;
3213: if ( 0 <= shiftCount ) {
3214: aSig &= LIT64( 0x7FFFFFFFFFFFFFFF );
3215: if ( ( a.high != 0xC03E ) || aSig ) {
3216: float_raise( float_flag_invalid STATUS_VAR);
3217: if ( ! aSign || ( ( aExp == 0x7FFF ) && aSig ) ) {
3218: return LIT64( 0x7FFFFFFFFFFFFFFF );
3219: }
3220: }
3221: return (sbits64) LIT64( 0x8000000000000000 );
3222: }
3223: else if ( aExp < 0x3FFF ) {
3224: if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
3225: return 0;
3226: }
3227: z = aSig>>( - shiftCount );
3228: if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
3229: STATUS(float_exception_flags) |= float_flag_inexact;
3230: }
3231: if ( aSign ) z = - z;
3232: return z;
3233:
3234: }
3235:
3236: /*----------------------------------------------------------------------------
3237: | Returns the result of converting the extended double-precision floating-
3238: | point value `a' to the single-precision floating-point format. The
3239: | conversion is performed according to the IEC/IEEE Standard for Binary
3240: | Floating-Point Arithmetic.
3241: *----------------------------------------------------------------------------*/
3242:
3243: float32 floatx80_to_float32( floatx80 a STATUS_PARAM )
3244: {
3245: flag aSign;
3246: int32 aExp;
3247: bits64 aSig;
3248:
3249: aSig = extractFloatx80Frac( a );
3250: aExp = extractFloatx80Exp( a );
3251: aSign = extractFloatx80Sign( a );
3252: if ( aExp == 0x7FFF ) {
3253: if ( (bits64) ( aSig<<1 ) ) {
3254: return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) );
3255: }
3256: return packFloat32( aSign, 0xFF, 0 );
3257: }
3258: shift64RightJamming( aSig, 33, &aSig );
3259: if ( aExp || aSig ) aExp -= 0x3F81;
3260: return roundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
3261:
3262: }
3263:
3264: /*----------------------------------------------------------------------------
3265: | Returns the result of converting the extended double-precision floating-
3266: | point value `a' to the double-precision floating-point format. The
3267: | conversion is performed according to the IEC/IEEE Standard for Binary
3268: | Floating-Point Arithmetic.
3269: *----------------------------------------------------------------------------*/
3270:
3271: float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
3272: {
3273: flag aSign;
3274: int32 aExp;
3275: bits64 aSig, zSig;
3276:
3277: aSig = extractFloatx80Frac( a );
3278: aExp = extractFloatx80Exp( a );
3279: aSign = extractFloatx80Sign( a );
3280: if ( aExp == 0x7FFF ) {
3281: if ( (bits64) ( aSig<<1 ) ) {
3282: return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) );
3283: }
3284: return packFloat64( aSign, 0x7FF, 0 );
3285: }
3286: shift64RightJamming( aSig, 1, &zSig );
3287: if ( aExp || aSig ) aExp -= 0x3C01;
3288: return roundAndPackFloat64( aSign, aExp, zSig STATUS_VAR );
3289:
3290: }
3291:
3292: #ifdef FLOAT128
3293:
3294: /*----------------------------------------------------------------------------
3295: | Returns the result of converting the extended double-precision floating-
3296: | point value `a' to the quadruple-precision floating-point format. The
3297: | conversion is performed according to the IEC/IEEE Standard for Binary
3298: | Floating-Point Arithmetic.
3299: *----------------------------------------------------------------------------*/
3300:
3301: float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
3302: {
3303: flag aSign;
3304: int16 aExp;
3305: bits64 aSig, zSig0, zSig1;
3306:
3307: aSig = extractFloatx80Frac( a );
3308: aExp = extractFloatx80Exp( a );
3309: aSign = extractFloatx80Sign( a );
3310: if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
3311: return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) );
3312: }
3313: shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
3314: return packFloat128( aSign, aExp, zSig0, zSig1 );
3315:
3316: }
3317:
3318: #endif
3319:
3320: /*----------------------------------------------------------------------------
3321: | Rounds the extended double-precision floating-point value `a' to an integer,
3322: | and returns the result as an extended quadruple-precision floating-point
3323: | value. The operation is performed according to the IEC/IEEE Standard for
3324: | Binary Floating-Point Arithmetic.
3325: *----------------------------------------------------------------------------*/
3326:
3327: floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM )
3328: {
3329: flag aSign;
3330: int32 aExp;
3331: bits64 lastBitMask, roundBitsMask;
3332: int8 roundingMode;
3333: floatx80 z;
3334:
3335: aExp = extractFloatx80Exp( a );
3336: if ( 0x403E <= aExp ) {
3337: if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
3338: return propagateFloatx80NaN( a, a STATUS_VAR );
3339: }
3340: return a;
3341: }
3342: if ( aExp < 0x3FFF ) {
3343: if ( ( aExp == 0 )
3344: && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
3345: return a;
3346: }
3347: STATUS(float_exception_flags) |= float_flag_inexact;
3348: aSign = extractFloatx80Sign( a );
3349: switch ( STATUS(float_rounding_mode) ) {
3350: case float_round_nearest_even:
3351: if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
3352: ) {
3353: return
3354: packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
3355: }
3356: break;
3357: case float_round_down:
3358: return
3359: aSign ?
3360: packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
3361: : packFloatx80( 0, 0, 0 );
3362: case float_round_up:
3363: return
3364: aSign ? packFloatx80( 1, 0, 0 )
3365: : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
3366: }
3367: return packFloatx80( aSign, 0, 0 );
3368: }
3369: lastBitMask = 1;
3370: lastBitMask <<= 0x403E - aExp;
3371: roundBitsMask = lastBitMask - 1;
3372: z = a;
3373: roundingMode = STATUS(float_rounding_mode);
3374: if ( roundingMode == float_round_nearest_even ) {
3375: z.low += lastBitMask>>1;
3376: if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
3377: }
3378: else if ( roundingMode != float_round_to_zero ) {
3379: if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
3380: z.low += roundBitsMask;
3381: }
3382: }
3383: z.low &= ~ roundBitsMask;
3384: if ( z.low == 0 ) {
3385: ++z.high;
3386: z.low = LIT64( 0x8000000000000000 );
3387: }
3388: if ( z.low != a.low ) STATUS(float_exception_flags) |= float_flag_inexact;
3389: return z;
3390:
3391: }
3392:
3393: /*----------------------------------------------------------------------------
3394: | Returns the result of adding the absolute values of the extended double-
3395: | precision floating-point values `a' and `b'. If `zSign' is 1, the sum is
3396: | negated before being returned. `zSign' is ignored if the result is a NaN.
3397: | The addition is performed according to the IEC/IEEE Standard for Binary
3398: | Floating-Point Arithmetic.
3399: *----------------------------------------------------------------------------*/
3400:
3401: static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM)
3402: {
3403: int32 aExp, bExp, zExp;
3404: bits64 aSig, bSig, zSig0, zSig1;
3405: int32 expDiff;
3406:
3407: aSig = extractFloatx80Frac( a );
3408: aExp = extractFloatx80Exp( a );
3409: bSig = extractFloatx80Frac( b );
3410: bExp = extractFloatx80Exp( b );
3411: expDiff = aExp - bExp;
3412: if ( 0 < expDiff ) {
3413: if ( aExp == 0x7FFF ) {
3414: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3415: return a;
3416: }
3417: if ( bExp == 0 ) --expDiff;
3418: shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3419: zExp = aExp;
3420: }
3421: else if ( expDiff < 0 ) {
3422: if ( bExp == 0x7FFF ) {
3423: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3424: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3425: }
3426: if ( aExp == 0 ) ++expDiff;
3427: shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3428: zExp = bExp;
3429: }
3430: else {
3431: if ( aExp == 0x7FFF ) {
3432: if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3433: return propagateFloatx80NaN( a, b STATUS_VAR );
3434: }
3435: return a;
3436: }
3437: zSig1 = 0;
3438: zSig0 = aSig + bSig;
3439: if ( aExp == 0 ) {
3440: normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
3441: goto roundAndPack;
3442: }
3443: zExp = aExp;
3444: goto shiftRight1;
3445: }
3446: zSig0 = aSig + bSig;
3447: if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
3448: shiftRight1:
3449: shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
3450: zSig0 |= LIT64( 0x8000000000000000 );
3451: ++zExp;
3452: roundAndPack:
3453: return
3454: roundAndPackFloatx80(
3455: STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3456:
3457: }
3458:
3459: /*----------------------------------------------------------------------------
3460: | Returns the result of subtracting the absolute values of the extended
3461: | double-precision floating-point values `a' and `b'. If `zSign' is 1, the
3462: | difference is negated before being returned. `zSign' is ignored if the
3463: | result is a NaN. The subtraction is performed according to the IEC/IEEE
3464: | Standard for Binary Floating-Point Arithmetic.
3465: *----------------------------------------------------------------------------*/
3466:
3467: static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM )
3468: {
3469: int32 aExp, bExp, zExp;
3470: bits64 aSig, bSig, zSig0, zSig1;
3471: int32 expDiff;
3472: floatx80 z;
3473:
3474: aSig = extractFloatx80Frac( a );
3475: aExp = extractFloatx80Exp( a );
3476: bSig = extractFloatx80Frac( b );
3477: bExp = extractFloatx80Exp( b );
3478: expDiff = aExp - bExp;
3479: if ( 0 < expDiff ) goto aExpBigger;
3480: if ( expDiff < 0 ) goto bExpBigger;
3481: if ( aExp == 0x7FFF ) {
3482: if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
3483: return propagateFloatx80NaN( a, b STATUS_VAR );
3484: }
3485: float_raise( float_flag_invalid STATUS_VAR);
3486: z.low = floatx80_default_nan_low;
3487: z.high = floatx80_default_nan_high;
3488: return z;
3489: }
3490: if ( aExp == 0 ) {
3491: aExp = 1;
3492: bExp = 1;
3493: }
3494: zSig1 = 0;
3495: if ( bSig < aSig ) goto aBigger;
3496: if ( aSig < bSig ) goto bBigger;
3497: return packFloatx80( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
3498: bExpBigger:
3499: if ( bExp == 0x7FFF ) {
3500: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3501: return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
3502: }
3503: if ( aExp == 0 ) ++expDiff;
3504: shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
3505: bBigger:
3506: sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
3507: zExp = bExp;
3508: zSign ^= 1;
3509: goto normalizeRoundAndPack;
3510: aExpBigger:
3511: if ( aExp == 0x7FFF ) {
3512: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3513: return a;
3514: }
3515: if ( bExp == 0 ) --expDiff;
3516: shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
3517: aBigger:
3518: sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
3519: zExp = aExp;
3520: normalizeRoundAndPack:
3521: return
3522: normalizeRoundAndPackFloatx80(
3523: STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3524:
3525: }
3526:
3527: /*----------------------------------------------------------------------------
3528: | Returns the result of adding the extended double-precision floating-point
3529: | values `a' and `b'. The operation is performed according to the IEC/IEEE
3530: | Standard for Binary Floating-Point Arithmetic.
3531: *----------------------------------------------------------------------------*/
3532:
3533: floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM )
3534: {
3535: flag aSign, bSign;
3536:
3537: aSign = extractFloatx80Sign( a );
3538: bSign = extractFloatx80Sign( b );
3539: if ( aSign == bSign ) {
3540: return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3541: }
3542: else {
3543: return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3544: }
3545:
3546: }
3547:
3548: /*----------------------------------------------------------------------------
3549: | Returns the result of subtracting the extended double-precision floating-
3550: | point values `a' and `b'. The operation is performed according to the
3551: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3552: *----------------------------------------------------------------------------*/
3553:
3554: floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM )
3555: {
3556: flag aSign, bSign;
3557:
3558: aSign = extractFloatx80Sign( a );
3559: bSign = extractFloatx80Sign( b );
3560: if ( aSign == bSign ) {
3561: return subFloatx80Sigs( a, b, aSign STATUS_VAR );
3562: }
3563: else {
3564: return addFloatx80Sigs( a, b, aSign STATUS_VAR );
3565: }
3566:
3567: }
3568:
3569: /*----------------------------------------------------------------------------
3570: | Returns the result of multiplying the extended double-precision floating-
3571: | point values `a' and `b'. The operation is performed according to the
3572: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3573: *----------------------------------------------------------------------------*/
3574:
3575: floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
3576: {
3577: flag aSign, bSign, zSign;
3578: int32 aExp, bExp, zExp;
3579: bits64 aSig, bSig, zSig0, zSig1;
3580: floatx80 z;
3581:
3582: aSig = extractFloatx80Frac( a );
3583: aExp = extractFloatx80Exp( a );
3584: aSign = extractFloatx80Sign( a );
3585: bSig = extractFloatx80Frac( b );
3586: bExp = extractFloatx80Exp( b );
3587: bSign = extractFloatx80Sign( b );
3588: zSign = aSign ^ bSign;
3589: if ( aExp == 0x7FFF ) {
3590: if ( (bits64) ( aSig<<1 )
3591: || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3592: return propagateFloatx80NaN( a, b STATUS_VAR );
3593: }
3594: if ( ( bExp | bSig ) == 0 ) goto invalid;
3595: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3596: }
3597: if ( bExp == 0x7FFF ) {
3598: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3599: if ( ( aExp | aSig ) == 0 ) {
3600: invalid:
3601: float_raise( float_flag_invalid STATUS_VAR);
3602: z.low = floatx80_default_nan_low;
3603: z.high = floatx80_default_nan_high;
3604: return z;
3605: }
3606: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3607: }
3608: if ( aExp == 0 ) {
3609: if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3610: normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3611: }
3612: if ( bExp == 0 ) {
3613: if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
3614: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3615: }
3616: zExp = aExp + bExp - 0x3FFE;
3617: mul64To128( aSig, bSig, &zSig0, &zSig1 );
3618: if ( 0 < (sbits64) zSig0 ) {
3619: shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
3620: --zExp;
3621: }
3622: return
3623: roundAndPackFloatx80(
3624: STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3625:
3626: }
3627:
3628: /*----------------------------------------------------------------------------
3629: | Returns the result of dividing the extended double-precision floating-point
3630: | value `a' by the corresponding value `b'. The operation is performed
3631: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3632: *----------------------------------------------------------------------------*/
3633:
3634: floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
3635: {
3636: flag aSign, bSign, zSign;
3637: int32 aExp, bExp, zExp;
3638: bits64 aSig, bSig, zSig0, zSig1;
3639: bits64 rem0, rem1, rem2, term0, term1, term2;
3640: floatx80 z;
3641:
3642: aSig = extractFloatx80Frac( a );
3643: aExp = extractFloatx80Exp( a );
3644: aSign = extractFloatx80Sign( a );
3645: bSig = extractFloatx80Frac( b );
3646: bExp = extractFloatx80Exp( b );
3647: bSign = extractFloatx80Sign( b );
3648: zSign = aSign ^ bSign;
3649: if ( aExp == 0x7FFF ) {
3650: if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3651: if ( bExp == 0x7FFF ) {
3652: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3653: goto invalid;
3654: }
3655: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3656: }
3657: if ( bExp == 0x7FFF ) {
3658: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3659: return packFloatx80( zSign, 0, 0 );
3660: }
3661: if ( bExp == 0 ) {
3662: if ( bSig == 0 ) {
3663: if ( ( aExp | aSig ) == 0 ) {
3664: invalid:
3665: float_raise( float_flag_invalid STATUS_VAR);
3666: z.low = floatx80_default_nan_low;
3667: z.high = floatx80_default_nan_high;
3668: return z;
3669: }
3670: float_raise( float_flag_divbyzero STATUS_VAR);
3671: return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
3672: }
3673: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3674: }
3675: if ( aExp == 0 ) {
3676: if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
3677: normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
3678: }
3679: zExp = aExp - bExp + 0x3FFE;
3680: rem1 = 0;
3681: if ( bSig <= aSig ) {
3682: shift128Right( aSig, 0, 1, &aSig, &rem1 );
3683: ++zExp;
3684: }
3685: zSig0 = estimateDiv128To64( aSig, rem1, bSig );
3686: mul64To128( bSig, zSig0, &term0, &term1 );
3687: sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
3688: while ( (sbits64) rem0 < 0 ) {
3689: --zSig0;
3690: add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
3691: }
3692: zSig1 = estimateDiv128To64( rem1, 0, bSig );
3693: if ( (bits64) ( zSig1<<1 ) <= 8 ) {
3694: mul64To128( bSig, zSig1, &term1, &term2 );
3695: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3696: while ( (sbits64) rem1 < 0 ) {
3697: --zSig1;
3698: add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
3699: }
3700: zSig1 |= ( ( rem1 | rem2 ) != 0 );
3701: }
3702: return
3703: roundAndPackFloatx80(
3704: STATUS(floatx80_rounding_precision), zSign, zExp, zSig0, zSig1 STATUS_VAR );
3705:
3706: }
3707:
3708: /*----------------------------------------------------------------------------
3709: | Returns the remainder of the extended double-precision floating-point value
3710: | `a' with respect to the corresponding value `b'. The operation is performed
3711: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3712: *----------------------------------------------------------------------------*/
3713:
3714: floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
3715: {
3716: flag aSign, bSign, zSign;
3717: int32 aExp, bExp, expDiff;
3718: bits64 aSig0, aSig1, bSig;
3719: bits64 q, term0, term1, alternateASig0, alternateASig1;
3720: floatx80 z;
3721:
3722: aSig0 = extractFloatx80Frac( a );
3723: aExp = extractFloatx80Exp( a );
3724: aSign = extractFloatx80Sign( a );
3725: bSig = extractFloatx80Frac( b );
3726: bExp = extractFloatx80Exp( b );
3727: bSign = extractFloatx80Sign( b );
3728: if ( aExp == 0x7FFF ) {
3729: if ( (bits64) ( aSig0<<1 )
3730: || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
3731: return propagateFloatx80NaN( a, b STATUS_VAR );
3732: }
3733: goto invalid;
3734: }
3735: if ( bExp == 0x7FFF ) {
3736: if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
3737: return a;
3738: }
3739: if ( bExp == 0 ) {
3740: if ( bSig == 0 ) {
3741: invalid:
3742: float_raise( float_flag_invalid STATUS_VAR);
3743: z.low = floatx80_default_nan_low;
3744: z.high = floatx80_default_nan_high;
3745: return z;
3746: }
3747: normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
3748: }
3749: if ( aExp == 0 ) {
3750: if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
3751: normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3752: }
3753: bSig |= LIT64( 0x8000000000000000 );
3754: zSign = aSign;
3755: expDiff = aExp - bExp;
3756: aSig1 = 0;
3757: if ( expDiff < 0 ) {
3758: if ( expDiff < -1 ) return a;
3759: shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
3760: expDiff = 0;
3761: }
3762: q = ( bSig <= aSig0 );
3763: if ( q ) aSig0 -= bSig;
3764: expDiff -= 64;
3765: while ( 0 < expDiff ) {
3766: q = estimateDiv128To64( aSig0, aSig1, bSig );
3767: q = ( 2 < q ) ? q - 2 : 0;
3768: mul64To128( bSig, q, &term0, &term1 );
3769: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3770: shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
3771: expDiff -= 62;
3772: }
3773: expDiff += 64;
3774: if ( 0 < expDiff ) {
3775: q = estimateDiv128To64( aSig0, aSig1, bSig );
3776: q = ( 2 < q ) ? q - 2 : 0;
3777: q >>= 64 - expDiff;
3778: mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
3779: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3780: shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
3781: while ( le128( term0, term1, aSig0, aSig1 ) ) {
3782: ++q;
3783: sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
3784: }
3785: }
3786: else {
3787: term1 = 0;
3788: term0 = bSig;
3789: }
3790: sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
3791: if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
3792: || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
3793: && ( q & 1 ) )
3794: ) {
3795: aSig0 = alternateASig0;
3796: aSig1 = alternateASig1;
3797: zSign = ! zSign;
3798: }
3799: return
3800: normalizeRoundAndPackFloatx80(
3801: 80, zSign, bExp + expDiff, aSig0, aSig1 STATUS_VAR );
3802:
3803: }
3804:
3805: /*----------------------------------------------------------------------------
3806: | Returns the square root of the extended double-precision floating-point
3807: | value `a'. The operation is performed according to the IEC/IEEE Standard
3808: | for Binary Floating-Point Arithmetic.
3809: *----------------------------------------------------------------------------*/
3810:
3811: floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
3812: {
3813: flag aSign;
3814: int32 aExp, zExp;
3815: bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
3816: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
3817: floatx80 z;
3818:
3819: aSig0 = extractFloatx80Frac( a );
3820: aExp = extractFloatx80Exp( a );
3821: aSign = extractFloatx80Sign( a );
3822: if ( aExp == 0x7FFF ) {
3823: if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a STATUS_VAR );
3824: if ( ! aSign ) return a;
3825: goto invalid;
3826: }
3827: if ( aSign ) {
3828: if ( ( aExp | aSig0 ) == 0 ) return a;
3829: invalid:
3830: float_raise( float_flag_invalid STATUS_VAR);
3831: z.low = floatx80_default_nan_low;
3832: z.high = floatx80_default_nan_high;
3833: return z;
3834: }
3835: if ( aExp == 0 ) {
3836: if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
3837: normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
3838: }
3839: zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
3840: zSig0 = estimateSqrt32( aExp, aSig0>>32 );
3841: shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
3842: zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
3843: doubleZSig0 = zSig0<<1;
3844: mul64To128( zSig0, zSig0, &term0, &term1 );
3845: sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
3846: while ( (sbits64) rem0 < 0 ) {
3847: --zSig0;
3848: doubleZSig0 -= 2;
3849: add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
3850: }
3851: zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
3852: if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
3853: if ( zSig1 == 0 ) zSig1 = 1;
3854: mul64To128( doubleZSig0, zSig1, &term1, &term2 );
3855: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
3856: mul64To128( zSig1, zSig1, &term2, &term3 );
3857: sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
3858: while ( (sbits64) rem1 < 0 ) {
3859: --zSig1;
3860: shortShift128Left( 0, zSig1, 1, &term2, &term3 );
3861: term3 |= 1;
3862: term2 |= doubleZSig0;
3863: add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
3864: }
3865: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
3866: }
3867: shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
3868: zSig0 |= doubleZSig0;
3869: return
3870: roundAndPackFloatx80(
3871: STATUS(floatx80_rounding_precision), 0, zExp, zSig0, zSig1 STATUS_VAR );
3872:
3873: }
3874:
3875: /*----------------------------------------------------------------------------
3876: | Returns 1 if the extended double-precision floating-point value `a' is
3877: | equal to the corresponding value `b', and 0 otherwise. The comparison is
3878: | performed according to the IEC/IEEE Standard for Binary Floating-Point
3879: | Arithmetic.
3880: *----------------------------------------------------------------------------*/
3881:
3882: flag floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
3883: {
3884:
3885: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3886: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3887: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3888: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3889: ) {
3890: if ( floatx80_is_signaling_nan( a )
3891: || floatx80_is_signaling_nan( b ) ) {
3892: float_raise( float_flag_invalid STATUS_VAR);
3893: }
3894: return 0;
3895: }
3896: return
3897: ( a.low == b.low )
3898: && ( ( a.high == b.high )
3899: || ( ( a.low == 0 )
3900: && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
3901: );
3902:
3903: }
3904:
3905: /*----------------------------------------------------------------------------
3906: | Returns 1 if the extended double-precision floating-point value `a' is
3907: | less than or equal to the corresponding value `b', and 0 otherwise. The
3908: | comparison is performed according to the IEC/IEEE Standard for Binary
3909: | Floating-Point Arithmetic.
3910: *----------------------------------------------------------------------------*/
3911:
3912: flag floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
3913: {
3914: flag aSign, bSign;
3915:
3916: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3917: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3918: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3919: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3920: ) {
3921: float_raise( float_flag_invalid STATUS_VAR);
3922: return 0;
3923: }
3924: aSign = extractFloatx80Sign( a );
3925: bSign = extractFloatx80Sign( b );
3926: if ( aSign != bSign ) {
3927: return
3928: aSign
3929: || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3930: == 0 );
3931: }
3932: return
3933: aSign ? le128( b.high, b.low, a.high, a.low )
3934: : le128( a.high, a.low, b.high, b.low );
3935:
3936: }
3937:
3938: /*----------------------------------------------------------------------------
3939: | Returns 1 if the extended double-precision floating-point value `a' is
3940: | less than the corresponding value `b', and 0 otherwise. The comparison
3941: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
3942: | Arithmetic.
3943: *----------------------------------------------------------------------------*/
3944:
3945: flag floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
3946: {
3947: flag aSign, bSign;
3948:
3949: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3950: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3951: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3952: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3953: ) {
3954: float_raise( float_flag_invalid STATUS_VAR);
3955: return 0;
3956: }
3957: aSign = extractFloatx80Sign( a );
3958: bSign = extractFloatx80Sign( b );
3959: if ( aSign != bSign ) {
3960: return
3961: aSign
3962: && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
3963: != 0 );
3964: }
3965: return
3966: aSign ? lt128( b.high, b.low, a.high, a.low )
3967: : lt128( a.high, a.low, b.high, b.low );
3968:
3969: }
3970:
3971: /*----------------------------------------------------------------------------
3972: | Returns 1 if the extended double-precision floating-point value `a' is equal
3973: | to the corresponding value `b', and 0 otherwise. The invalid exception is
3974: | raised if either operand is a NaN. Otherwise, the comparison is performed
3975: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
3976: *----------------------------------------------------------------------------*/
3977:
3978: flag floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
3979: {
3980:
3981: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
3982: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
3983: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
3984: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
3985: ) {
3986: float_raise( float_flag_invalid STATUS_VAR);
3987: return 0;
3988: }
3989: return
3990: ( a.low == b.low )
3991: && ( ( a.high == b.high )
3992: || ( ( a.low == 0 )
3993: && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
3994: );
3995:
3996: }
3997:
3998: /*----------------------------------------------------------------------------
3999: | Returns 1 if the extended double-precision floating-point value `a' is less
4000: | than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs
4001: | do not cause an exception. Otherwise, the comparison is performed according
4002: | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4003: *----------------------------------------------------------------------------*/
4004:
4005: flag floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4006: {
4007: flag aSign, bSign;
4008:
4009: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4010: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4011: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4012: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4013: ) {
4014: if ( floatx80_is_signaling_nan( a )
4015: || floatx80_is_signaling_nan( b ) ) {
4016: float_raise( float_flag_invalid STATUS_VAR);
4017: }
4018: return 0;
4019: }
4020: aSign = extractFloatx80Sign( a );
4021: bSign = extractFloatx80Sign( b );
4022: if ( aSign != bSign ) {
4023: return
4024: aSign
4025: || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4026: == 0 );
4027: }
4028: return
4029: aSign ? le128( b.high, b.low, a.high, a.low )
4030: : le128( a.high, a.low, b.high, b.low );
4031:
4032: }
4033:
4034: /*----------------------------------------------------------------------------
4035: | Returns 1 if the extended double-precision floating-point value `a' is less
4036: | than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause
4037: | an exception. Otherwise, the comparison is performed according to the
4038: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4039: *----------------------------------------------------------------------------*/
4040:
4041: flag floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
4042: {
4043: flag aSign, bSign;
4044:
4045: if ( ( ( extractFloatx80Exp( a ) == 0x7FFF )
4046: && (bits64) ( extractFloatx80Frac( a )<<1 ) )
4047: || ( ( extractFloatx80Exp( b ) == 0x7FFF )
4048: && (bits64) ( extractFloatx80Frac( b )<<1 ) )
4049: ) {
4050: if ( floatx80_is_signaling_nan( a )
4051: || floatx80_is_signaling_nan( b ) ) {
4052: float_raise( float_flag_invalid STATUS_VAR);
4053: }
4054: return 0;
4055: }
4056: aSign = extractFloatx80Sign( a );
4057: bSign = extractFloatx80Sign( b );
4058: if ( aSign != bSign ) {
4059: return
4060: aSign
4061: && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
4062: != 0 );
4063: }
4064: return
4065: aSign ? lt128( b.high, b.low, a.high, a.low )
4066: : lt128( a.high, a.low, b.high, b.low );
4067:
4068: }
4069:
4070: #endif
4071:
4072: #ifdef FLOAT128
4073:
4074: /*----------------------------------------------------------------------------
4075: | Returns the result of converting the quadruple-precision floating-point
4076: | value `a' to the 32-bit two's complement integer format. The conversion
4077: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4078: | Arithmetic---which means in particular that the conversion is rounded
4079: | according to the current rounding mode. If `a' is a NaN, the largest
4080: | positive integer is returned. Otherwise, if the conversion overflows, the
4081: | largest integer with the same sign as `a' is returned.
4082: *----------------------------------------------------------------------------*/
4083:
4084: int32 float128_to_int32( float128 a STATUS_PARAM )
4085: {
4086: flag aSign;
4087: int32 aExp, shiftCount;
4088: bits64 aSig0, aSig1;
4089:
4090: aSig1 = extractFloat128Frac1( a );
4091: aSig0 = extractFloat128Frac0( a );
4092: aExp = extractFloat128Exp( a );
4093: aSign = extractFloat128Sign( a );
4094: if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
4095: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4096: aSig0 |= ( aSig1 != 0 );
4097: shiftCount = 0x4028 - aExp;
4098: if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
4099: return roundAndPackInt32( aSign, aSig0 STATUS_VAR );
4100:
4101: }
4102:
4103: /*----------------------------------------------------------------------------
4104: | Returns the result of converting the quadruple-precision floating-point
4105: | value `a' to the 32-bit two's complement integer format. The conversion
4106: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4107: | Arithmetic, except that the conversion is always rounded toward zero. If
4108: | `a' is a NaN, the largest positive integer is returned. Otherwise, if the
4109: | conversion overflows, the largest integer with the same sign as `a' is
4110: | returned.
4111: *----------------------------------------------------------------------------*/
4112:
4113: int32 float128_to_int32_round_to_zero( float128 a STATUS_PARAM )
4114: {
4115: flag aSign;
4116: int32 aExp, shiftCount;
4117: bits64 aSig0, aSig1, savedASig;
4118: int32 z;
4119:
4120: aSig1 = extractFloat128Frac1( a );
4121: aSig0 = extractFloat128Frac0( a );
4122: aExp = extractFloat128Exp( a );
4123: aSign = extractFloat128Sign( a );
4124: aSig0 |= ( aSig1 != 0 );
4125: if ( 0x401E < aExp ) {
4126: if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
4127: goto invalid;
4128: }
4129: else if ( aExp < 0x3FFF ) {
4130: if ( aExp || aSig0 ) STATUS(float_exception_flags) |= float_flag_inexact;
4131: return 0;
4132: }
4133: aSig0 |= LIT64( 0x0001000000000000 );
4134: shiftCount = 0x402F - aExp;
4135: savedASig = aSig0;
4136: aSig0 >>= shiftCount;
4137: z = aSig0;
4138: if ( aSign ) z = - z;
4139: if ( ( z < 0 ) ^ aSign ) {
4140: invalid:
4141: float_raise( float_flag_invalid STATUS_VAR);
4142: return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
4143: }
4144: if ( ( aSig0<<shiftCount ) != savedASig ) {
4145: STATUS(float_exception_flags) |= float_flag_inexact;
4146: }
4147: return z;
4148:
4149: }
4150:
4151: /*----------------------------------------------------------------------------
4152: | Returns the result of converting the quadruple-precision floating-point
4153: | value `a' to the 64-bit two's complement integer format. The conversion
4154: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4155: | Arithmetic---which means in particular that the conversion is rounded
4156: | according to the current rounding mode. If `a' is a NaN, the largest
4157: | positive integer is returned. Otherwise, if the conversion overflows, the
4158: | largest integer with the same sign as `a' is returned.
4159: *----------------------------------------------------------------------------*/
4160:
4161: int64 float128_to_int64( float128 a STATUS_PARAM )
4162: {
4163: flag aSign;
4164: int32 aExp, shiftCount;
4165: bits64 aSig0, aSig1;
4166:
4167: aSig1 = extractFloat128Frac1( a );
4168: aSig0 = extractFloat128Frac0( a );
4169: aExp = extractFloat128Exp( a );
4170: aSign = extractFloat128Sign( a );
4171: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4172: shiftCount = 0x402F - aExp;
4173: if ( shiftCount <= 0 ) {
4174: if ( 0x403E < aExp ) {
4175: float_raise( float_flag_invalid STATUS_VAR);
4176: if ( ! aSign
4177: || ( ( aExp == 0x7FFF )
4178: && ( aSig1 || ( aSig0 != LIT64( 0x0001000000000000 ) ) )
4179: )
4180: ) {
4181: return LIT64( 0x7FFFFFFFFFFFFFFF );
4182: }
4183: return (sbits64) LIT64( 0x8000000000000000 );
4184: }
4185: shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
4186: }
4187: else {
4188: shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
4189: }
4190: return roundAndPackInt64( aSign, aSig0, aSig1 STATUS_VAR );
4191:
4192: }
4193:
4194: /*----------------------------------------------------------------------------
4195: | Returns the result of converting the quadruple-precision floating-point
4196: | value `a' to the 64-bit two's complement integer format. The conversion
4197: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4198: | Arithmetic, except that the conversion is always rounded toward zero.
4199: | If `a' is a NaN, the largest positive integer is returned. Otherwise, if
4200: | the conversion overflows, the largest integer with the same sign as `a' is
4201: | returned.
4202: *----------------------------------------------------------------------------*/
4203:
4204: int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
4205: {
4206: flag aSign;
4207: int32 aExp, shiftCount;
4208: bits64 aSig0, aSig1;
4209: int64 z;
4210:
4211: aSig1 = extractFloat128Frac1( a );
4212: aSig0 = extractFloat128Frac0( a );
4213: aExp = extractFloat128Exp( a );
4214: aSign = extractFloat128Sign( a );
4215: if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
4216: shiftCount = aExp - 0x402F;
4217: if ( 0 < shiftCount ) {
4218: if ( 0x403E <= aExp ) {
4219: aSig0 &= LIT64( 0x0000FFFFFFFFFFFF );
4220: if ( ( a.high == LIT64( 0xC03E000000000000 ) )
4221: && ( aSig1 < LIT64( 0x0002000000000000 ) ) ) {
4222: if ( aSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
4223: }
4224: else {
4225: float_raise( float_flag_invalid STATUS_VAR);
4226: if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
4227: return LIT64( 0x7FFFFFFFFFFFFFFF );
4228: }
4229: }
4230: return (sbits64) LIT64( 0x8000000000000000 );
4231: }
4232: z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
4233: if ( (bits64) ( aSig1<<shiftCount ) ) {
4234: STATUS(float_exception_flags) |= float_flag_inexact;
4235: }
4236: }
4237: else {
4238: if ( aExp < 0x3FFF ) {
4239: if ( aExp | aSig0 | aSig1 ) {
4240: STATUS(float_exception_flags) |= float_flag_inexact;
4241: }
4242: return 0;
4243: }
4244: z = aSig0>>( - shiftCount );
4245: if ( aSig1
4246: || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
4247: STATUS(float_exception_flags) |= float_flag_inexact;
4248: }
4249: }
4250: if ( aSign ) z = - z;
4251: return z;
4252:
4253: }
4254:
4255: /*----------------------------------------------------------------------------
4256: | Returns the result of converting the quadruple-precision floating-point
4257: | value `a' to the single-precision floating-point format. The conversion
4258: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4259: | Arithmetic.
4260: *----------------------------------------------------------------------------*/
4261:
4262: float32 float128_to_float32( float128 a STATUS_PARAM )
4263: {
4264: flag aSign;
4265: int32 aExp;
4266: bits64 aSig0, aSig1;
4267: bits32 zSig;
4268:
4269: aSig1 = extractFloat128Frac1( a );
4270: aSig0 = extractFloat128Frac0( a );
4271: aExp = extractFloat128Exp( a );
4272: aSign = extractFloat128Sign( a );
4273: if ( aExp == 0x7FFF ) {
4274: if ( aSig0 | aSig1 ) {
4275: return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR ) );
4276: }
4277: return packFloat32( aSign, 0xFF, 0 );
4278: }
4279: aSig0 |= ( aSig1 != 0 );
4280: shift64RightJamming( aSig0, 18, &aSig0 );
4281: zSig = aSig0;
4282: if ( aExp || zSig ) {
4283: zSig |= 0x40000000;
4284: aExp -= 0x3F81;
4285: }
4286: return roundAndPackFloat32( aSign, aExp, zSig STATUS_VAR );
4287:
4288: }
4289:
4290: /*----------------------------------------------------------------------------
4291: | Returns the result of converting the quadruple-precision floating-point
4292: | value `a' to the double-precision floating-point format. The conversion
4293: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
4294: | Arithmetic.
4295: *----------------------------------------------------------------------------*/
4296:
4297: float64 float128_to_float64( float128 a STATUS_PARAM )
4298: {
4299: flag aSign;
4300: int32 aExp;
4301: bits64 aSig0, aSig1;
4302:
4303: aSig1 = extractFloat128Frac1( a );
4304: aSig0 = extractFloat128Frac0( a );
4305: aExp = extractFloat128Exp( a );
4306: aSign = extractFloat128Sign( a );
4307: if ( aExp == 0x7FFF ) {
4308: if ( aSig0 | aSig1 ) {
4309: return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR ) );
4310: }
4311: return packFloat64( aSign, 0x7FF, 0 );
4312: }
4313: shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4314: aSig0 |= ( aSig1 != 0 );
4315: if ( aExp || aSig0 ) {
4316: aSig0 |= LIT64( 0x4000000000000000 );
4317: aExp -= 0x3C01;
4318: }
4319: return roundAndPackFloat64( aSign, aExp, aSig0 STATUS_VAR );
4320:
4321: }
4322:
4323: #ifdef FLOATX80
4324:
4325: /*----------------------------------------------------------------------------
4326: | Returns the result of converting the quadruple-precision floating-point
4327: | value `a' to the extended double-precision floating-point format. The
4328: | conversion is performed according to the IEC/IEEE Standard for Binary
4329: | Floating-Point Arithmetic.
4330: *----------------------------------------------------------------------------*/
4331:
4332: floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
4333: {
4334: flag aSign;
4335: int32 aExp;
4336: bits64 aSig0, aSig1;
4337:
4338: aSig1 = extractFloat128Frac1( a );
4339: aSig0 = extractFloat128Frac0( a );
4340: aExp = extractFloat128Exp( a );
4341: aSign = extractFloat128Sign( a );
4342: if ( aExp == 0x7FFF ) {
4343: if ( aSig0 | aSig1 ) {
4344: return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR ) );
4345: }
4346: return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
4347: }
4348: if ( aExp == 0 ) {
4349: if ( ( aSig0 | aSig1 ) == 0 ) return packFloatx80( aSign, 0, 0 );
4350: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4351: }
4352: else {
4353: aSig0 |= LIT64( 0x0001000000000000 );
4354: }
4355: shortShift128Left( aSig0, aSig1, 15, &aSig0, &aSig1 );
4356: return roundAndPackFloatx80( 80, aSign, aExp, aSig0, aSig1 STATUS_VAR );
4357:
4358: }
4359:
4360: #endif
4361:
4362: /*----------------------------------------------------------------------------
4363: | Rounds the quadruple-precision floating-point value `a' to an integer, and
4364: | returns the result as a quadruple-precision floating-point value. The
4365: | operation is performed according to the IEC/IEEE Standard for Binary
4366: | Floating-Point Arithmetic.
4367: *----------------------------------------------------------------------------*/
4368:
4369: float128 float128_round_to_int( float128 a STATUS_PARAM )
4370: {
4371: flag aSign;
4372: int32 aExp;
4373: bits64 lastBitMask, roundBitsMask;
4374: int8 roundingMode;
4375: float128 z;
4376:
4377: aExp = extractFloat128Exp( a );
4378: if ( 0x402F <= aExp ) {
4379: if ( 0x406F <= aExp ) {
4380: if ( ( aExp == 0x7FFF )
4381: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
4382: ) {
4383: return propagateFloat128NaN( a, a STATUS_VAR );
4384: }
4385: return a;
4386: }
4387: lastBitMask = 1;
4388: lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
4389: roundBitsMask = lastBitMask - 1;
4390: z = a;
4391: roundingMode = STATUS(float_rounding_mode);
4392: if ( roundingMode == float_round_nearest_even ) {
4393: if ( lastBitMask ) {
4394: add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
4395: if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
4396: }
4397: else {
4398: if ( (sbits64) z.low < 0 ) {
4399: ++z.high;
4400: if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
4401: }
4402: }
4403: }
4404: else if ( roundingMode != float_round_to_zero ) {
4405: if ( extractFloat128Sign( z )
4406: ^ ( roundingMode == float_round_up ) ) {
4407: add128( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
4408: }
4409: }
4410: z.low &= ~ roundBitsMask;
4411: }
4412: else {
4413: if ( aExp < 0x3FFF ) {
4414: if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
4415: STATUS(float_exception_flags) |= float_flag_inexact;
4416: aSign = extractFloat128Sign( a );
4417: switch ( STATUS(float_rounding_mode) ) {
4418: case float_round_nearest_even:
4419: if ( ( aExp == 0x3FFE )
4420: && ( extractFloat128Frac0( a )
4421: | extractFloat128Frac1( a ) )
4422: ) {
4423: return packFloat128( aSign, 0x3FFF, 0, 0 );
4424: }
4425: break;
4426: case float_round_down:
4427: return
4428: aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
4429: : packFloat128( 0, 0, 0, 0 );
4430: case float_round_up:
4431: return
4432: aSign ? packFloat128( 1, 0, 0, 0 )
4433: : packFloat128( 0, 0x3FFF, 0, 0 );
4434: }
4435: return packFloat128( aSign, 0, 0, 0 );
4436: }
4437: lastBitMask = 1;
4438: lastBitMask <<= 0x402F - aExp;
4439: roundBitsMask = lastBitMask - 1;
4440: z.low = 0;
4441: z.high = a.high;
4442: roundingMode = STATUS(float_rounding_mode);
4443: if ( roundingMode == float_round_nearest_even ) {
4444: z.high += lastBitMask>>1;
4445: if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
4446: z.high &= ~ lastBitMask;
4447: }
4448: }
4449: else if ( roundingMode != float_round_to_zero ) {
4450: if ( extractFloat128Sign( z )
4451: ^ ( roundingMode == float_round_up ) ) {
4452: z.high |= ( a.low != 0 );
4453: z.high += roundBitsMask;
4454: }
4455: }
4456: z.high &= ~ roundBitsMask;
4457: }
4458: if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
4459: STATUS(float_exception_flags) |= float_flag_inexact;
4460: }
4461: return z;
4462:
4463: }
4464:
4465: /*----------------------------------------------------------------------------
4466: | Returns the result of adding the absolute values of the quadruple-precision
4467: | floating-point values `a' and `b'. If `zSign' is 1, the sum is negated
4468: | before being returned. `zSign' is ignored if the result is a NaN.
4469: | The addition is performed according to the IEC/IEEE Standard for Binary
4470: | Floating-Point Arithmetic.
4471: *----------------------------------------------------------------------------*/
4472:
4473: static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4474: {
4475: int32 aExp, bExp, zExp;
4476: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4477: int32 expDiff;
4478:
4479: aSig1 = extractFloat128Frac1( a );
4480: aSig0 = extractFloat128Frac0( a );
4481: aExp = extractFloat128Exp( a );
4482: bSig1 = extractFloat128Frac1( b );
4483: bSig0 = extractFloat128Frac0( b );
4484: bExp = extractFloat128Exp( b );
4485: expDiff = aExp - bExp;
4486: if ( 0 < expDiff ) {
4487: if ( aExp == 0x7FFF ) {
4488: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4489: return a;
4490: }
4491: if ( bExp == 0 ) {
4492: --expDiff;
4493: }
4494: else {
4495: bSig0 |= LIT64( 0x0001000000000000 );
4496: }
4497: shift128ExtraRightJamming(
4498: bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
4499: zExp = aExp;
4500: }
4501: else if ( expDiff < 0 ) {
4502: if ( bExp == 0x7FFF ) {
4503: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4504: return packFloat128( zSign, 0x7FFF, 0, 0 );
4505: }
4506: if ( aExp == 0 ) {
4507: ++expDiff;
4508: }
4509: else {
4510: aSig0 |= LIT64( 0x0001000000000000 );
4511: }
4512: shift128ExtraRightJamming(
4513: aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
4514: zExp = bExp;
4515: }
4516: else {
4517: if ( aExp == 0x7FFF ) {
4518: if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4519: return propagateFloat128NaN( a, b STATUS_VAR );
4520: }
4521: return a;
4522: }
4523: add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4524: if ( aExp == 0 ) return packFloat128( zSign, 0, zSig0, zSig1 );
4525: zSig2 = 0;
4526: zSig0 |= LIT64( 0x0002000000000000 );
4527: zExp = aExp;
4528: goto shiftRight1;
4529: }
4530: aSig0 |= LIT64( 0x0001000000000000 );
4531: add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4532: --zExp;
4533: if ( zSig0 < LIT64( 0x0002000000000000 ) ) goto roundAndPack;
4534: ++zExp;
4535: shiftRight1:
4536: shift128ExtraRightJamming(
4537: zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4538: roundAndPack:
4539: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4540:
4541: }
4542:
4543: /*----------------------------------------------------------------------------
4544: | Returns the result of subtracting the absolute values of the quadruple-
4545: | precision floating-point values `a' and `b'. If `zSign' is 1, the
4546: | difference is negated before being returned. `zSign' is ignored if the
4547: | result is a NaN. The subtraction is performed according to the IEC/IEEE
4548: | Standard for Binary Floating-Point Arithmetic.
4549: *----------------------------------------------------------------------------*/
4550:
4551: static float128 subFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
4552: {
4553: int32 aExp, bExp, zExp;
4554: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
4555: int32 expDiff;
4556: float128 z;
4557:
4558: aSig1 = extractFloat128Frac1( a );
4559: aSig0 = extractFloat128Frac0( a );
4560: aExp = extractFloat128Exp( a );
4561: bSig1 = extractFloat128Frac1( b );
4562: bSig0 = extractFloat128Frac0( b );
4563: bExp = extractFloat128Exp( b );
4564: expDiff = aExp - bExp;
4565: shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
4566: shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
4567: if ( 0 < expDiff ) goto aExpBigger;
4568: if ( expDiff < 0 ) goto bExpBigger;
4569: if ( aExp == 0x7FFF ) {
4570: if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
4571: return propagateFloat128NaN( a, b STATUS_VAR );
4572: }
4573: float_raise( float_flag_invalid STATUS_VAR);
4574: z.low = float128_default_nan_low;
4575: z.high = float128_default_nan_high;
4576: return z;
4577: }
4578: if ( aExp == 0 ) {
4579: aExp = 1;
4580: bExp = 1;
4581: }
4582: if ( bSig0 < aSig0 ) goto aBigger;
4583: if ( aSig0 < bSig0 ) goto bBigger;
4584: if ( bSig1 < aSig1 ) goto aBigger;
4585: if ( aSig1 < bSig1 ) goto bBigger;
4586: return packFloat128( STATUS(float_rounding_mode) == float_round_down, 0, 0, 0 );
4587: bExpBigger:
4588: if ( bExp == 0x7FFF ) {
4589: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4590: return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
4591: }
4592: if ( aExp == 0 ) {
4593: ++expDiff;
4594: }
4595: else {
4596: aSig0 |= LIT64( 0x4000000000000000 );
4597: }
4598: shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
4599: bSig0 |= LIT64( 0x4000000000000000 );
4600: bBigger:
4601: sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
4602: zExp = bExp;
4603: zSign ^= 1;
4604: goto normalizeRoundAndPack;
4605: aExpBigger:
4606: if ( aExp == 0x7FFF ) {
4607: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4608: return a;
4609: }
4610: if ( bExp == 0 ) {
4611: --expDiff;
4612: }
4613: else {
4614: bSig0 |= LIT64( 0x4000000000000000 );
4615: }
4616: shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
4617: aSig0 |= LIT64( 0x4000000000000000 );
4618: aBigger:
4619: sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
4620: zExp = aExp;
4621: normalizeRoundAndPack:
4622: --zExp;
4623: return normalizeRoundAndPackFloat128( zSign, zExp - 14, zSig0, zSig1 STATUS_VAR );
4624:
4625: }
4626:
4627: /*----------------------------------------------------------------------------
4628: | Returns the result of adding the quadruple-precision floating-point values
4629: | `a' and `b'. The operation is performed according to the IEC/IEEE Standard
4630: | for Binary Floating-Point Arithmetic.
4631: *----------------------------------------------------------------------------*/
4632:
4633: float128 float128_add( float128 a, float128 b STATUS_PARAM )
4634: {
4635: flag aSign, bSign;
4636:
4637: aSign = extractFloat128Sign( a );
4638: bSign = extractFloat128Sign( b );
4639: if ( aSign == bSign ) {
4640: return addFloat128Sigs( a, b, aSign STATUS_VAR );
4641: }
4642: else {
4643: return subFloat128Sigs( a, b, aSign STATUS_VAR );
4644: }
4645:
4646: }
4647:
4648: /*----------------------------------------------------------------------------
4649: | Returns the result of subtracting the quadruple-precision floating-point
4650: | values `a' and `b'. The operation is performed according to the IEC/IEEE
4651: | Standard for Binary Floating-Point Arithmetic.
4652: *----------------------------------------------------------------------------*/
4653:
4654: float128 float128_sub( float128 a, float128 b STATUS_PARAM )
4655: {
4656: flag aSign, bSign;
4657:
4658: aSign = extractFloat128Sign( a );
4659: bSign = extractFloat128Sign( b );
4660: if ( aSign == bSign ) {
4661: return subFloat128Sigs( a, b, aSign STATUS_VAR );
4662: }
4663: else {
4664: return addFloat128Sigs( a, b, aSign STATUS_VAR );
4665: }
4666:
4667: }
4668:
4669: /*----------------------------------------------------------------------------
4670: | Returns the result of multiplying the quadruple-precision floating-point
4671: | values `a' and `b'. The operation is performed according to the IEC/IEEE
4672: | Standard for Binary Floating-Point Arithmetic.
4673: *----------------------------------------------------------------------------*/
4674:
4675: float128 float128_mul( float128 a, float128 b STATUS_PARAM )
4676: {
4677: flag aSign, bSign, zSign;
4678: int32 aExp, bExp, zExp;
4679: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
4680: float128 z;
4681:
4682: aSig1 = extractFloat128Frac1( a );
4683: aSig0 = extractFloat128Frac0( a );
4684: aExp = extractFloat128Exp( a );
4685: aSign = extractFloat128Sign( a );
4686: bSig1 = extractFloat128Frac1( b );
4687: bSig0 = extractFloat128Frac0( b );
4688: bExp = extractFloat128Exp( b );
4689: bSign = extractFloat128Sign( b );
4690: zSign = aSign ^ bSign;
4691: if ( aExp == 0x7FFF ) {
4692: if ( ( aSig0 | aSig1 )
4693: || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
4694: return propagateFloat128NaN( a, b STATUS_VAR );
4695: }
4696: if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
4697: return packFloat128( zSign, 0x7FFF, 0, 0 );
4698: }
4699: if ( bExp == 0x7FFF ) {
4700: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4701: if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4702: invalid:
4703: float_raise( float_flag_invalid STATUS_VAR);
4704: z.low = float128_default_nan_low;
4705: z.high = float128_default_nan_high;
4706: return z;
4707: }
4708: return packFloat128( zSign, 0x7FFF, 0, 0 );
4709: }
4710: if ( aExp == 0 ) {
4711: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4712: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4713: }
4714: if ( bExp == 0 ) {
4715: if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4716: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4717: }
4718: zExp = aExp + bExp - 0x4000;
4719: aSig0 |= LIT64( 0x0001000000000000 );
4720: shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
4721: mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
4722: add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
4723: zSig2 |= ( zSig3 != 0 );
4724: if ( LIT64( 0x0002000000000000 ) <= zSig0 ) {
4725: shift128ExtraRightJamming(
4726: zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
4727: ++zExp;
4728: }
4729: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4730:
4731: }
4732:
4733: /*----------------------------------------------------------------------------
4734: | Returns the result of dividing the quadruple-precision floating-point value
4735: | `a' by the corresponding value `b'. The operation is performed according to
4736: | the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4737: *----------------------------------------------------------------------------*/
4738:
4739: float128 float128_div( float128 a, float128 b STATUS_PARAM )
4740: {
4741: flag aSign, bSign, zSign;
4742: int32 aExp, bExp, zExp;
4743: bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
4744: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
4745: float128 z;
4746:
4747: aSig1 = extractFloat128Frac1( a );
4748: aSig0 = extractFloat128Frac0( a );
4749: aExp = extractFloat128Exp( a );
4750: aSign = extractFloat128Sign( a );
4751: bSig1 = extractFloat128Frac1( b );
4752: bSig0 = extractFloat128Frac0( b );
4753: bExp = extractFloat128Exp( b );
4754: bSign = extractFloat128Sign( b );
4755: zSign = aSign ^ bSign;
4756: if ( aExp == 0x7FFF ) {
4757: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4758: if ( bExp == 0x7FFF ) {
4759: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4760: goto invalid;
4761: }
4762: return packFloat128( zSign, 0x7FFF, 0, 0 );
4763: }
4764: if ( bExp == 0x7FFF ) {
4765: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4766: return packFloat128( zSign, 0, 0, 0 );
4767: }
4768: if ( bExp == 0 ) {
4769: if ( ( bSig0 | bSig1 ) == 0 ) {
4770: if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
4771: invalid:
4772: float_raise( float_flag_invalid STATUS_VAR);
4773: z.low = float128_default_nan_low;
4774: z.high = float128_default_nan_high;
4775: return z;
4776: }
4777: float_raise( float_flag_divbyzero STATUS_VAR);
4778: return packFloat128( zSign, 0x7FFF, 0, 0 );
4779: }
4780: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4781: }
4782: if ( aExp == 0 ) {
4783: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
4784: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4785: }
4786: zExp = aExp - bExp + 0x3FFD;
4787: shortShift128Left(
4788: aSig0 | LIT64( 0x0001000000000000 ), aSig1, 15, &aSig0, &aSig1 );
4789: shortShift128Left(
4790: bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
4791: if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
4792: shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
4793: ++zExp;
4794: }
4795: zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
4796: mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
4797: sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
4798: while ( (sbits64) rem0 < 0 ) {
4799: --zSig0;
4800: add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
4801: }
4802: zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
4803: if ( ( zSig1 & 0x3FFF ) <= 4 ) {
4804: mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
4805: sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
4806: while ( (sbits64) rem1 < 0 ) {
4807: --zSig1;
4808: add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
4809: }
4810: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
4811: }
4812: shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
4813: return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4814:
4815: }
4816:
4817: /*----------------------------------------------------------------------------
4818: | Returns the remainder of the quadruple-precision floating-point value `a'
4819: | with respect to the corresponding value `b'. The operation is performed
4820: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
4821: *----------------------------------------------------------------------------*/
4822:
4823: float128 float128_rem( float128 a, float128 b STATUS_PARAM )
4824: {
4825: flag aSign, bSign, zSign;
4826: int32 aExp, bExp, expDiff;
4827: bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
4828: bits64 allZero, alternateASig0, alternateASig1, sigMean1;
4829: sbits64 sigMean0;
4830: float128 z;
4831:
4832: aSig1 = extractFloat128Frac1( a );
4833: aSig0 = extractFloat128Frac0( a );
4834: aExp = extractFloat128Exp( a );
4835: aSign = extractFloat128Sign( a );
4836: bSig1 = extractFloat128Frac1( b );
4837: bSig0 = extractFloat128Frac0( b );
4838: bExp = extractFloat128Exp( b );
4839: bSign = extractFloat128Sign( b );
4840: if ( aExp == 0x7FFF ) {
4841: if ( ( aSig0 | aSig1 )
4842: || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
4843: return propagateFloat128NaN( a, b STATUS_VAR );
4844: }
4845: goto invalid;
4846: }
4847: if ( bExp == 0x7FFF ) {
4848: if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b STATUS_VAR );
4849: return a;
4850: }
4851: if ( bExp == 0 ) {
4852: if ( ( bSig0 | bSig1 ) == 0 ) {
4853: invalid:
4854: float_raise( float_flag_invalid STATUS_VAR);
4855: z.low = float128_default_nan_low;
4856: z.high = float128_default_nan_high;
4857: return z;
4858: }
4859: normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
4860: }
4861: if ( aExp == 0 ) {
4862: if ( ( aSig0 | aSig1 ) == 0 ) return a;
4863: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4864: }
4865: expDiff = aExp - bExp;
4866: if ( expDiff < -1 ) return a;
4867: shortShift128Left(
4868: aSig0 | LIT64( 0x0001000000000000 ),
4869: aSig1,
4870: 15 - ( expDiff < 0 ),
4871: &aSig0,
4872: &aSig1
4873: );
4874: shortShift128Left(
4875: bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
4876: q = le128( bSig0, bSig1, aSig0, aSig1 );
4877: if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
4878: expDiff -= 64;
4879: while ( 0 < expDiff ) {
4880: q = estimateDiv128To64( aSig0, aSig1, bSig0 );
4881: q = ( 4 < q ) ? q - 4 : 0;
4882: mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
4883: shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZero );
4884: shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero );
4885: sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 );
4886: expDiff -= 61;
4887: }
4888: if ( -64 < expDiff ) {
4889: q = estimateDiv128To64( aSig0, aSig1, bSig0 );
4890: q = ( 4 < q ) ? q - 4 : 0;
4891: q >>= - expDiff;
4892: shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
4893: expDiff += 52;
4894: if ( expDiff < 0 ) {
4895: shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
4896: }
4897: else {
4898: shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
4899: }
4900: mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
4901: sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
4902: }
4903: else {
4904: shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 );
4905: shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
4906: }
4907: do {
4908: alternateASig0 = aSig0;
4909: alternateASig1 = aSig1;
4910: ++q;
4911: sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
4912: } while ( 0 <= (sbits64) aSig0 );
4913: add128(
4914: aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
4915: if ( ( sigMean0 < 0 )
4916: || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
4917: aSig0 = alternateASig0;
4918: aSig1 = alternateASig1;
4919: }
4920: zSign = ( (sbits64) aSig0 < 0 );
4921: if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
4922: return
4923: normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 STATUS_VAR );
4924:
4925: }
4926:
4927: /*----------------------------------------------------------------------------
4928: | Returns the square root of the quadruple-precision floating-point value `a'.
4929: | The operation is performed according to the IEC/IEEE Standard for Binary
4930: | Floating-Point Arithmetic.
4931: *----------------------------------------------------------------------------*/
4932:
4933: float128 float128_sqrt( float128 a STATUS_PARAM )
4934: {
4935: flag aSign;
4936: int32 aExp, zExp;
4937: bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
4938: bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
4939: float128 z;
4940:
4941: aSig1 = extractFloat128Frac1( a );
4942: aSig0 = extractFloat128Frac0( a );
4943: aExp = extractFloat128Exp( a );
4944: aSign = extractFloat128Sign( a );
4945: if ( aExp == 0x7FFF ) {
4946: if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, a STATUS_VAR );
4947: if ( ! aSign ) return a;
4948: goto invalid;
4949: }
4950: if ( aSign ) {
4951: if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
4952: invalid:
4953: float_raise( float_flag_invalid STATUS_VAR);
4954: z.low = float128_default_nan_low;
4955: z.high = float128_default_nan_high;
4956: return z;
4957: }
4958: if ( aExp == 0 ) {
4959: if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( 0, 0, 0, 0 );
4960: normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
4961: }
4962: zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFE;
4963: aSig0 |= LIT64( 0x0001000000000000 );
4964: zSig0 = estimateSqrt32( aExp, aSig0>>17 );
4965: shortShift128Left( aSig0, aSig1, 13 - ( aExp & 1 ), &aSig0, &aSig1 );
4966: zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
4967: doubleZSig0 = zSig0<<1;
4968: mul64To128( zSig0, zSig0, &term0, &term1 );
4969: sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
4970: while ( (sbits64) rem0 < 0 ) {
4971: --zSig0;
4972: doubleZSig0 -= 2;
4973: add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
4974: }
4975: zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
4976: if ( ( zSig1 & 0x1FFF ) <= 5 ) {
4977: if ( zSig1 == 0 ) zSig1 = 1;
4978: mul64To128( doubleZSig0, zSig1, &term1, &term2 );
4979: sub128( rem1, 0, term1, term2, &rem1, &rem2 );
4980: mul64To128( zSig1, zSig1, &term2, &term3 );
4981: sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
4982: while ( (sbits64) rem1 < 0 ) {
4983: --zSig1;
4984: shortShift128Left( 0, zSig1, 1, &term2, &term3 );
4985: term3 |= 1;
4986: term2 |= doubleZSig0;
4987: add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
4988: }
4989: zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
4990: }
4991: shift128ExtraRightJamming( zSig0, zSig1, 0, 14, &zSig0, &zSig1, &zSig2 );
4992: return roundAndPackFloat128( 0, zExp, zSig0, zSig1, zSig2 STATUS_VAR );
4993:
4994: }
4995:
4996: /*----------------------------------------------------------------------------
4997: | Returns 1 if the quadruple-precision floating-point value `a' is equal to
4998: | the corresponding value `b', and 0 otherwise. The comparison is performed
4999: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5000: *----------------------------------------------------------------------------*/
5001:
5002: flag float128_eq( float128 a, float128 b STATUS_PARAM )
5003: {
5004:
5005: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5006: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5007: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5008: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5009: ) {
5010: if ( float128_is_signaling_nan( a )
5011: || float128_is_signaling_nan( b ) ) {
5012: float_raise( float_flag_invalid STATUS_VAR);
5013: }
5014: return 0;
5015: }
5016: return
5017: ( a.low == b.low )
5018: && ( ( a.high == b.high )
5019: || ( ( a.low == 0 )
5020: && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5021: );
5022:
5023: }
5024:
5025: /*----------------------------------------------------------------------------
5026: | Returns 1 if the quadruple-precision floating-point value `a' is less than
5027: | or equal to the corresponding value `b', and 0 otherwise. The comparison
5028: | is performed according to the IEC/IEEE Standard for Binary Floating-Point
5029: | Arithmetic.
5030: *----------------------------------------------------------------------------*/
5031:
5032: flag float128_le( float128 a, float128 b STATUS_PARAM )
5033: {
5034: flag aSign, bSign;
5035:
5036: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5037: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5038: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5039: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5040: ) {
5041: float_raise( float_flag_invalid STATUS_VAR);
5042: return 0;
5043: }
5044: aSign = extractFloat128Sign( a );
5045: bSign = extractFloat128Sign( b );
5046: if ( aSign != bSign ) {
5047: return
5048: aSign
5049: || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5050: == 0 );
5051: }
5052: return
5053: aSign ? le128( b.high, b.low, a.high, a.low )
5054: : le128( a.high, a.low, b.high, b.low );
5055:
5056: }
5057:
5058: /*----------------------------------------------------------------------------
5059: | Returns 1 if the quadruple-precision floating-point value `a' is less than
5060: | the corresponding value `b', and 0 otherwise. The comparison is performed
5061: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5062: *----------------------------------------------------------------------------*/
5063:
5064: flag float128_lt( float128 a, float128 b STATUS_PARAM )
5065: {
5066: flag aSign, bSign;
5067:
5068: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5069: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5070: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5071: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5072: ) {
5073: float_raise( float_flag_invalid STATUS_VAR);
5074: return 0;
5075: }
5076: aSign = extractFloat128Sign( a );
5077: bSign = extractFloat128Sign( b );
5078: if ( aSign != bSign ) {
5079: return
5080: aSign
5081: && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5082: != 0 );
5083: }
5084: return
5085: aSign ? lt128( b.high, b.low, a.high, a.low )
5086: : lt128( a.high, a.low, b.high, b.low );
5087:
5088: }
5089:
5090: /*----------------------------------------------------------------------------
5091: | Returns 1 if the quadruple-precision floating-point value `a' is equal to
5092: | the corresponding value `b', and 0 otherwise. The invalid exception is
5093: | raised if either operand is a NaN. Otherwise, the comparison is performed
5094: | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5095: *----------------------------------------------------------------------------*/
5096:
5097: flag float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
5098: {
5099:
5100: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5101: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5102: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5103: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5104: ) {
5105: float_raise( float_flag_invalid STATUS_VAR);
5106: return 0;
5107: }
5108: return
5109: ( a.low == b.low )
5110: && ( ( a.high == b.high )
5111: || ( ( a.low == 0 )
5112: && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
5113: );
5114:
5115: }
5116:
5117: /*----------------------------------------------------------------------------
5118: | Returns 1 if the quadruple-precision floating-point value `a' is less than
5119: | or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
5120: | cause an exception. Otherwise, the comparison is performed according to the
5121: | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
5122: *----------------------------------------------------------------------------*/
5123:
5124: flag float128_le_quiet( float128 a, float128 b STATUS_PARAM )
5125: {
5126: flag aSign, bSign;
5127:
5128: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5129: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5130: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5131: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5132: ) {
5133: if ( float128_is_signaling_nan( a )
5134: || float128_is_signaling_nan( b ) ) {
5135: float_raise( float_flag_invalid STATUS_VAR);
5136: }
5137: return 0;
5138: }
5139: aSign = extractFloat128Sign( a );
5140: bSign = extractFloat128Sign( b );
5141: if ( aSign != bSign ) {
5142: return
5143: aSign
5144: || ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5145: == 0 );
5146: }
5147: return
5148: aSign ? le128( b.high, b.low, a.high, a.low )
5149: : le128( a.high, a.low, b.high, b.low );
5150:
5151: }
5152:
5153: /*----------------------------------------------------------------------------
5154: | Returns 1 if the quadruple-precision floating-point value `a' is less than
5155: | the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an
5156: | exception. Otherwise, the comparison is performed according to the IEC/IEEE
5157: | Standard for Binary Floating-Point Arithmetic.
5158: *----------------------------------------------------------------------------*/
5159:
5160: flag float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
5161: {
5162: flag aSign, bSign;
5163:
5164: if ( ( ( extractFloat128Exp( a ) == 0x7FFF )
5165: && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
5166: || ( ( extractFloat128Exp( b ) == 0x7FFF )
5167: && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
5168: ) {
5169: if ( float128_is_signaling_nan( a )
5170: || float128_is_signaling_nan( b ) ) {
5171: float_raise( float_flag_invalid STATUS_VAR);
5172: }
5173: return 0;
5174: }
5175: aSign = extractFloat128Sign( a );
5176: bSign = extractFloat128Sign( b );
5177: if ( aSign != bSign ) {
5178: return
5179: aSign
5180: && ( ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
5181: != 0 );
5182: }
5183: return
5184: aSign ? lt128( b.high, b.low, a.high, a.low )
5185: : lt128( a.high, a.low, b.high, b.low );
5186:
5187: }
5188:
5189: #endif
5190:
5191: /* misc functions */
5192: float32 uint32_to_float32( unsigned int a STATUS_PARAM )
5193: {
5194: return int64_to_float32(a STATUS_VAR);
5195: }
5196:
5197: float64 uint32_to_float64( unsigned int a STATUS_PARAM )
5198: {
5199: return int64_to_float64(a STATUS_VAR);
5200: }
5201:
5202: unsigned int float32_to_uint32( float32 a STATUS_PARAM )
5203: {
5204: int64_t v;
5205: unsigned int res;
5206:
5207: v = float32_to_int64(a STATUS_VAR);
5208: if (v < 0) {
5209: res = 0;
5210: float_raise( float_flag_invalid STATUS_VAR);
5211: } else if (v > 0xffffffff) {
5212: res = 0xffffffff;
5213: float_raise( float_flag_invalid STATUS_VAR);
5214: } else {
5215: res = v;
5216: }
5217: return res;
5218: }
5219:
5220: unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
5221: {
5222: int64_t v;
5223: unsigned int res;
5224:
5225: v = float32_to_int64_round_to_zero(a STATUS_VAR);
5226: if (v < 0) {
5227: res = 0;
5228: float_raise( float_flag_invalid STATUS_VAR);
5229: } else if (v > 0xffffffff) {
5230: res = 0xffffffff;
5231: float_raise( float_flag_invalid STATUS_VAR);
5232: } else {
5233: res = v;
5234: }
5235: return res;
5236: }
5237:
5238: unsigned int float64_to_uint32( float64 a STATUS_PARAM )
5239: {
5240: int64_t v;
5241: unsigned int res;
5242:
5243: v = float64_to_int64(a STATUS_VAR);
5244: if (v < 0) {
5245: res = 0;
5246: float_raise( float_flag_invalid STATUS_VAR);
5247: } else if (v > 0xffffffff) {
5248: res = 0xffffffff;
5249: float_raise( float_flag_invalid STATUS_VAR);
5250: } else {
5251: res = v;
5252: }
5253: return res;
5254: }
5255:
5256: unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
5257: {
5258: int64_t v;
5259: unsigned int res;
5260:
5261: v = float64_to_int64_round_to_zero(a STATUS_VAR);
5262: if (v < 0) {
5263: res = 0;
5264: float_raise( float_flag_invalid STATUS_VAR);
5265: } else if (v > 0xffffffff) {
5266: res = 0xffffffff;
5267: float_raise( float_flag_invalid STATUS_VAR);
5268: } else {
5269: res = v;
5270: }
5271: return res;
5272: }
5273:
5274: #define COMPARE(s, nan_exp) \
5275: INLINE char float ## s ## _compare_internal( float ## s a, float ## s b, \
5276: int is_quiet STATUS_PARAM ) \
5277: { \
5278: flag aSign, bSign; \
5279: \
5280: if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) && \
5281: extractFloat ## s ## Frac( a ) ) || \
5282: ( ( extractFloat ## s ## Exp( b ) == nan_exp ) && \
5283: extractFloat ## s ## Frac( b ) )) { \
5284: if (!is_quiet || \
5285: float ## s ## _is_signaling_nan( a ) || \
5286: float ## s ## _is_signaling_nan( b ) ) { \
5287: float_raise( float_flag_invalid STATUS_VAR); \
5288: } \
5289: return float_relation_unordered; \
5290: } \
5291: aSign = extractFloat ## s ## Sign( a ); \
5292: bSign = extractFloat ## s ## Sign( b ); \
5293: if ( aSign != bSign ) { \
5294: if ( (bits ## s) ( ( a | b )<<1 ) == 0 ) { \
5295: /* zero case */ \
5296: return float_relation_equal; \
5297: } else { \
5298: return 1 - (2 * aSign); \
5299: } \
5300: } else { \
5301: if (a == b) { \
5302: return float_relation_equal; \
5303: } else { \
5304: return 1 - 2 * (aSign ^ ( a < b )); \
5305: } \
5306: } \
5307: } \
5308: \
5309: char float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM ) \
5310: { \
5311: return float ## s ## _compare_internal(a, b, 0 STATUS_VAR); \
5312: } \
5313: \
5314: char float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
5315: { \
5316: return float ## s ## _compare_internal(a, b, 1 STATUS_VAR); \
5317: }
5318:
5319: COMPARE(32, 0xff)
5320: COMPARE(64, 0x7ff)
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.