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