/arch/arm/nwfpe/softfloat.c
C | 1426 lines | 922 code | 118 blank | 386 comment | 350 complexity | 7565592ccc1ed64672923ca7b357ff5f MD5 | raw file
- /*
- ===============================================================================
- This C source file is part of the SoftFloat IEC/IEEE Floating-point
- Arithmetic Package, Release 2.
- Written by John R. Hauser. This work was made possible in part by the
- International Computer Science Institute, located at Suite 600, 1947 Center
- Street, Berkeley, California 94704. Funding was partially provided by the
- National Science Foundation under grant MIP-9311980. The original version
- of this code was written as part of a project to build a fixed-point vector
- processor in collaboration with the University of California at Berkeley,
- overseen by Profs. Nelson Morgan and John Wawrzynek. More information
- is available through the web page
- http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
- THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort
- has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
- TIMES RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO
- PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
- AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
- Derivative works are acceptable, even for commercial purposes, so long as
- (1) they include prominent notice that the work is derivative, and (2) they
- include prominent notice akin to these three paragraphs for those parts of
- this code that are retained.
- ===============================================================================
- */
- #include <asm/div64.h>
- #include "fpa11.h"
- //#include "milieu.h"
- //#include "softfloat.h"
- /*
- -------------------------------------------------------------------------------
- Primitive arithmetic functions, including multi-word arithmetic, and
- division and square root approximations. (Can be specialized to target if
- desired.)
- -------------------------------------------------------------------------------
- */
- #include "softfloat-macros"
- /*
- -------------------------------------------------------------------------------
- Functions and definitions to determine: (1) whether tininess for underflow
- is detected before or after rounding by default, (2) what (if anything)
- happens when exceptions are raised, (3) how signaling NaNs are distinguished
- from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
- are propagated from function inputs to output. These details are target-
- specific.
- -------------------------------------------------------------------------------
- */
- #include "softfloat-specialize"
- /*
- -------------------------------------------------------------------------------
- Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
- and 7, and returns the properly rounded 32-bit integer corresponding to the
- input. If `zSign' is nonzero, the input is negated before being converted
- to an integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point
- input is simply rounded to an integer, with the inexact exception raised if
- the input cannot be represented exactly as an integer. If the fixed-point
- input is too large, however, the invalid exception is raised and the largest
- positive or negative integer is returned.
- -------------------------------------------------------------------------------
- */
- static int32 roundAndPackInt32( struct roundingData *roundData, flag zSign, bits64 absZ )
- {
- int8 roundingMode;
- flag roundNearestEven;
- int8 roundIncrement, roundBits;
- int32 z;
- roundingMode = roundData->mode;
- roundNearestEven = ( roundingMode == float_round_nearest_even );
- roundIncrement = 0x40;
- if ( ! roundNearestEven ) {
- if ( roundingMode == float_round_to_zero ) {
- roundIncrement = 0;
- }
- else {
- roundIncrement = 0x7F;
- if ( zSign ) {
- if ( roundingMode == float_round_up ) roundIncrement = 0;
- }
- else {
- if ( roundingMode == float_round_down ) roundIncrement = 0;
- }
- }
- }
- roundBits = absZ & 0x7F;
- absZ = ( absZ + roundIncrement )>>7;
- absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
- z = absZ;
- if ( zSign ) z = - z;
- if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
- roundData->exception |= float_flag_invalid;
- return zSign ? 0x80000000 : 0x7FFFFFFF;
- }
- if ( roundBits ) roundData->exception |= float_flag_inexact;
- return z;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the fraction bits of the single-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE bits32 extractFloat32Frac( float32 a )
- {
- return a & 0x007FFFFF;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the exponent bits of the single-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE int16 extractFloat32Exp( float32 a )
- {
- return ( a>>23 ) & 0xFF;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the sign bit of the single-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- #if 0 /* in softfloat.h */
- INLINE flag extractFloat32Sign( float32 a )
- {
- return a>>31;
- }
- #endif
- /*
- -------------------------------------------------------------------------------
- Normalizes the subnormal single-precision floating-point value represented
- by the denormalized significand `aSig'. The normalized exponent and
- significand are stored at the locations pointed to by `zExpPtr' and
- `zSigPtr', respectively.
- -------------------------------------------------------------------------------
- */
- static void
- normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
- {
- int8 shiftCount;
- shiftCount = countLeadingZeros32( aSig ) - 8;
- *zSigPtr = aSig<<shiftCount;
- *zExpPtr = 1 - shiftCount;
- }
- /*
- -------------------------------------------------------------------------------
- Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
- single-precision floating-point value, returning the result. After being
- shifted into the proper positions, the three fields are simply added
- together to form the result. This means that any integer portion of `zSig'
- will be added into the exponent. Since a properly normalized significand
- will have an integer portion equal to 1, the `zExp' input should be 1 less
- than the desired result exponent whenever `zSig' is a complete, normalized
- significand.
- -------------------------------------------------------------------------------
- */
- INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
- {
- #if 0
- float32 f;
- __asm__("@ packFloat32 \n\
- mov %0, %1, asl #31 \n\
- orr %0, %2, asl #23 \n\
- orr %0, %3"
- : /* no outputs */
- : "g" (f), "g" (zSign), "g" (zExp), "g" (zSig)
- : "cc");
- return f;
- #else
- return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
- #endif
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent `zExp',
- and significand `zSig', and returns the proper single-precision floating-
- point value corresponding to the abstract input. Ordinarily, the abstract
- value is simply rounded and packed into the single-precision format, with
- the inexact exception raised if the abstract input cannot be represented
- exactly. If the abstract value is too large, however, the overflow and
- inexact exceptions are raised and an infinity or maximal finite value is
- returned. If the abstract value is too small, the input value is rounded to
- a subnormal number, and the underflow and inexact exceptions are raised if
- the abstract input cannot be represented exactly as a subnormal single-
- precision floating-point number.
- The input significand `zSig' has its binary point between bits 30
- and 29, which is 7 bits to the left of the usual location. This shifted
- significand must be normalized or smaller. If `zSig' is not normalized,
- `zExp' must be 0; in that case, the result returned is a subnormal number,
- and it must not require rounding. In the usual case that `zSig' is
- normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
- The handling of underflow and overflow follows the IEC/IEEE Standard for
- Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- static float32 roundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig )
- {
- int8 roundingMode;
- flag roundNearestEven;
- int8 roundIncrement, roundBits;
- flag isTiny;
- roundingMode = roundData->mode;
- roundNearestEven = ( roundingMode == float_round_nearest_even );
- roundIncrement = 0x40;
- if ( ! roundNearestEven ) {
- if ( roundingMode == float_round_to_zero ) {
- roundIncrement = 0;
- }
- else {
- roundIncrement = 0x7F;
- if ( zSign ) {
- if ( roundingMode == float_round_up ) roundIncrement = 0;
- }
- else {
- if ( roundingMode == float_round_down ) roundIncrement = 0;
- }
- }
- }
- roundBits = zSig & 0x7F;
- if ( 0xFD <= (bits16) zExp ) {
- if ( ( 0xFD < zExp )
- || ( ( zExp == 0xFD )
- && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
- ) {
- roundData->exception |= float_flag_overflow | float_flag_inexact;
- return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
- }
- if ( zExp < 0 ) {
- isTiny =
- ( float_detect_tininess == float_tininess_before_rounding )
- || ( zExp < -1 )
- || ( zSig + roundIncrement < 0x80000000 );
- shift32RightJamming( zSig, - zExp, &zSig );
- zExp = 0;
- roundBits = zSig & 0x7F;
- if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
- }
- }
- if ( roundBits ) roundData->exception |= float_flag_inexact;
- zSig = ( zSig + roundIncrement )>>7;
- zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
- if ( zSig == 0 ) zExp = 0;
- return packFloat32( zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent `zExp',
- and significand `zSig', and returns the proper single-precision floating-
- point value corresponding to the abstract input. This routine is just like
- `roundAndPackFloat32' except that `zSig' does not have to be normalized in
- any way. In all cases, `zExp' must be 1 less than the ``true'' floating-
- point exponent.
- -------------------------------------------------------------------------------
- */
- static float32
- normalizeRoundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig )
- {
- int8 shiftCount;
- shiftCount = countLeadingZeros32( zSig ) - 1;
- return roundAndPackFloat32( roundData, zSign, zExp - shiftCount, zSig<<shiftCount );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the fraction bits of the double-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE bits64 extractFloat64Frac( float64 a )
- {
- return a & LIT64( 0x000FFFFFFFFFFFFF );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the exponent bits of the double-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE int16 extractFloat64Exp( float64 a )
- {
- return ( a>>52 ) & 0x7FF;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the sign bit of the double-precision floating-point value `a'.
- -------------------------------------------------------------------------------
- */
- #if 0 /* in softfloat.h */
- INLINE flag extractFloat64Sign( float64 a )
- {
- return a>>63;
- }
- #endif
- /*
- -------------------------------------------------------------------------------
- Normalizes the subnormal double-precision floating-point value represented
- by the denormalized significand `aSig'. The normalized exponent and
- significand are stored at the locations pointed to by `zExpPtr' and
- `zSigPtr', respectively.
- -------------------------------------------------------------------------------
- */
- static void
- normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
- {
- int8 shiftCount;
- shiftCount = countLeadingZeros64( aSig ) - 11;
- *zSigPtr = aSig<<shiftCount;
- *zExpPtr = 1 - shiftCount;
- }
- /*
- -------------------------------------------------------------------------------
- Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
- double-precision floating-point value, returning the result. After being
- shifted into the proper positions, the three fields are simply added
- together to form the result. This means that any integer portion of `zSig'
- will be added into the exponent. Since a properly normalized significand
- will have an integer portion equal to 1, the `zExp' input should be 1 less
- than the desired result exponent whenever `zSig' is a complete, normalized
- significand.
- -------------------------------------------------------------------------------
- */
- INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
- {
- return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent `zExp',
- and significand `zSig', and returns the proper double-precision floating-
- point value corresponding to the abstract input. Ordinarily, the abstract
- value is simply rounded and packed into the double-precision format, with
- the inexact exception raised if the abstract input cannot be represented
- exactly. If the abstract value is too large, however, the overflow and
- inexact exceptions are raised and an infinity or maximal finite value is
- returned. If the abstract value is too small, the input value is rounded to
- a subnormal number, and the underflow and inexact exceptions are raised if
- the abstract input cannot be represented exactly as a subnormal double-
- precision floating-point number.
- The input significand `zSig' has its binary point between bits 62
- and 61, which is 10 bits to the left of the usual location. This shifted
- significand must be normalized or smaller. If `zSig' is not normalized,
- `zExp' must be 0; in that case, the result returned is a subnormal number,
- and it must not require rounding. In the usual case that `zSig' is
- normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
- The handling of underflow and overflow follows the IEC/IEEE Standard for
- Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- static float64 roundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig )
- {
- int8 roundingMode;
- flag roundNearestEven;
- int16 roundIncrement, roundBits;
- flag isTiny;
- roundingMode = roundData->mode;
- roundNearestEven = ( roundingMode == float_round_nearest_even );
- roundIncrement = 0x200;
- if ( ! roundNearestEven ) {
- if ( roundingMode == float_round_to_zero ) {
- roundIncrement = 0;
- }
- else {
- roundIncrement = 0x3FF;
- if ( zSign ) {
- if ( roundingMode == float_round_up ) roundIncrement = 0;
- }
- else {
- if ( roundingMode == float_round_down ) roundIncrement = 0;
- }
- }
- }
- roundBits = zSig & 0x3FF;
- if ( 0x7FD <= (bits16) zExp ) {
- if ( ( 0x7FD < zExp )
- || ( ( zExp == 0x7FD )
- && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
- ) {
- //register int lr = __builtin_return_address(0);
- //printk("roundAndPackFloat64 called from 0x%08x\n",lr);
- roundData->exception |= float_flag_overflow | float_flag_inexact;
- return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );
- }
- if ( zExp < 0 ) {
- isTiny =
- ( float_detect_tininess == float_tininess_before_rounding )
- || ( zExp < -1 )
- || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
- shift64RightJamming( zSig, - zExp, &zSig );
- zExp = 0;
- roundBits = zSig & 0x3FF;
- if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
- }
- }
- if ( roundBits ) roundData->exception |= float_flag_inexact;
- zSig = ( zSig + roundIncrement )>>10;
- zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
- if ( zSig == 0 ) zExp = 0;
- return packFloat64( zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent `zExp',
- and significand `zSig', and returns the proper double-precision floating-
- point value corresponding to the abstract input. This routine is just like
- `roundAndPackFloat64' except that `zSig' does not have to be normalized in
- any way. In all cases, `zExp' must be 1 less than the ``true'' floating-
- point exponent.
- -------------------------------------------------------------------------------
- */
- static float64
- normalizeRoundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig )
- {
- int8 shiftCount;
- shiftCount = countLeadingZeros64( zSig ) - 1;
- return roundAndPackFloat64( roundData, zSign, zExp - shiftCount, zSig<<shiftCount );
- }
- #ifdef FLOATX80
- /*
- -------------------------------------------------------------------------------
- Returns the fraction bits of the extended double-precision floating-point
- value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE bits64 extractFloatx80Frac( floatx80 a )
- {
- return a.low;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the exponent bits of the extended double-precision floating-point
- value `a'.
- -------------------------------------------------------------------------------
- */
- INLINE int32 extractFloatx80Exp( floatx80 a )
- {
- return a.high & 0x7FFF;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the sign bit of the extended double-precision floating-point value
- `a'.
- -------------------------------------------------------------------------------
- */
- INLINE flag extractFloatx80Sign( floatx80 a )
- {
- return a.high>>15;
- }
- /*
- -------------------------------------------------------------------------------
- Normalizes the subnormal extended double-precision floating-point value
- represented by the denormalized significand `aSig'. The normalized exponent
- and significand are stored at the locations pointed to by `zExpPtr' and
- `zSigPtr', respectively.
- -------------------------------------------------------------------------------
- */
- static void
- normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
- {
- int8 shiftCount;
- shiftCount = countLeadingZeros64( aSig );
- *zSigPtr = aSig<<shiftCount;
- *zExpPtr = 1 - shiftCount;
- }
- /*
- -------------------------------------------------------------------------------
- Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
- extended double-precision floating-point value, returning the result.
- -------------------------------------------------------------------------------
- */
- INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
- {
- floatx80 z;
- z.low = zSig;
- z.high = ( ( (bits16) zSign )<<15 ) + zExp;
- z.__padding = 0;
- return z;
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent `zExp',
- and extended significand formed by the concatenation of `zSig0' and `zSig1',
- and returns the proper extended double-precision floating-point value
- corresponding to the abstract input. Ordinarily, the abstract value is
- rounded and packed into the extended double-precision format, with the
- inexact exception raised if the abstract input cannot be represented
- exactly. If the abstract value is too large, however, the overflow and
- inexact exceptions are raised and an infinity or maximal finite value is
- returned. If the abstract value is too small, the input value is rounded to
- a subnormal number, and the underflow and inexact exceptions are raised if
- the abstract input cannot be represented exactly as a subnormal extended
- double-precision floating-point number.
- If `roundingPrecision' is 32 or 64, the result is rounded to the same
- number of bits as single or double precision, respectively. Otherwise, the
- result is rounded to the full precision of the extended double-precision
- format.
- The input significand must be normalized or smaller. If the input
- significand is not normalized, `zExp' must be 0; in that case, the result
- returned is a subnormal number, and it must not require rounding. The
- handling of underflow and overflow follows the IEC/IEEE Standard for Binary
- Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- static floatx80
- roundAndPackFloatx80(
- struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
- )
- {
- int8 roundingMode, roundingPrecision;
- flag roundNearestEven, increment, isTiny;
- int64 roundIncrement, roundMask, roundBits;
- roundingMode = roundData->mode;
- roundingPrecision = roundData->precision;
- roundNearestEven = ( roundingMode == float_round_nearest_even );
- if ( roundingPrecision == 80 ) goto precision80;
- if ( roundingPrecision == 64 ) {
- roundIncrement = LIT64( 0x0000000000000400 );
- roundMask = LIT64( 0x00000000000007FF );
- }
- else if ( roundingPrecision == 32 ) {
- roundIncrement = LIT64( 0x0000008000000000 );
- roundMask = LIT64( 0x000000FFFFFFFFFF );
- }
- else {
- goto precision80;
- }
- zSig0 |= ( zSig1 != 0 );
- if ( ! roundNearestEven ) {
- if ( roundingMode == float_round_to_zero ) {
- roundIncrement = 0;
- }
- else {
- roundIncrement = roundMask;
- if ( zSign ) {
- if ( roundingMode == float_round_up ) roundIncrement = 0;
- }
- else {
- if ( roundingMode == float_round_down ) roundIncrement = 0;
- }
- }
- }
- roundBits = zSig0 & roundMask;
- if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
- if ( ( 0x7FFE < zExp )
- || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
- ) {
- goto overflow;
- }
- if ( zExp <= 0 ) {
- isTiny =
- ( float_detect_tininess == float_tininess_before_rounding )
- || ( zExp < 0 )
- || ( zSig0 <= zSig0 + roundIncrement );
- shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
- zExp = 0;
- roundBits = zSig0 & roundMask;
- if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
- if ( roundBits ) roundData->exception |= float_flag_inexact;
- zSig0 += roundIncrement;
- if ( (sbits64) zSig0 < 0 ) zExp = 1;
- roundIncrement = roundMask + 1;
- if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
- roundMask |= roundIncrement;
- }
- zSig0 &= ~ roundMask;
- return packFloatx80( zSign, zExp, zSig0 );
- }
- }
- if ( roundBits ) roundData->exception |= float_flag_inexact;
- zSig0 += roundIncrement;
- if ( zSig0 < roundIncrement ) {
- ++zExp;
- zSig0 = LIT64( 0x8000000000000000 );
- }
- roundIncrement = roundMask + 1;
- if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
- roundMask |= roundIncrement;
- }
- zSig0 &= ~ roundMask;
- if ( zSig0 == 0 ) zExp = 0;
- return packFloatx80( zSign, zExp, zSig0 );
- precision80:
- increment = ( (sbits64) zSig1 < 0 );
- if ( ! roundNearestEven ) {
- if ( roundingMode == float_round_to_zero ) {
- increment = 0;
- }
- else {
- if ( zSign ) {
- increment = ( roundingMode == float_round_down ) && zSig1;
- }
- else {
- increment = ( roundingMode == float_round_up ) && zSig1;
- }
- }
- }
- if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
- if ( ( 0x7FFE < zExp )
- || ( ( zExp == 0x7FFE )
- && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
- && increment
- )
- ) {
- roundMask = 0;
- overflow:
- roundData->exception |= float_flag_overflow | float_flag_inexact;
- if ( ( roundingMode == float_round_to_zero )
- || ( zSign && ( roundingMode == float_round_up ) )
- || ( ! zSign && ( roundingMode == float_round_down ) )
- ) {
- return packFloatx80( zSign, 0x7FFE, ~ roundMask );
- }
- return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
- }
- if ( zExp <= 0 ) {
- isTiny =
- ( float_detect_tininess == float_tininess_before_rounding )
- || ( zExp < 0 )
- || ! increment
- || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
- shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
- zExp = 0;
- if ( isTiny && zSig1 ) roundData->exception |= float_flag_underflow;
- if ( zSig1 ) roundData->exception |= float_flag_inexact;
- if ( roundNearestEven ) {
- increment = ( (sbits64) zSig1 < 0 );
- }
- else {
- if ( zSign ) {
- increment = ( roundingMode == float_round_down ) && zSig1;
- }
- else {
- increment = ( roundingMode == float_round_up ) && zSig1;
- }
- }
- if ( increment ) {
- ++zSig0;
- zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );
- if ( (sbits64) zSig0 < 0 ) zExp = 1;
- }
- return packFloatx80( zSign, zExp, zSig0 );
- }
- }
- if ( zSig1 ) roundData->exception |= float_flag_inexact;
- if ( increment ) {
- ++zSig0;
- if ( zSig0 == 0 ) {
- ++zExp;
- zSig0 = LIT64( 0x8000000000000000 );
- }
- else {
- zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );
- }
- }
- else {
- if ( zSig0 == 0 ) zExp = 0;
- }
-
- return packFloatx80( zSign, zExp, zSig0 );
- }
- /*
- -------------------------------------------------------------------------------
- Takes an abstract floating-point value having sign `zSign', exponent
- `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
- and returns the proper extended double-precision floating-point value
- corresponding to the abstract input. This routine is just like
- `roundAndPackFloatx80' except that the input significand does not have to be
- normalized.
- -------------------------------------------------------------------------------
- */
- static floatx80
- normalizeRoundAndPackFloatx80(
- struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
- )
- {
- int8 shiftCount;
- if ( zSig0 == 0 ) {
- zSig0 = zSig1;
- zSig1 = 0;
- zExp -= 64;
- }
- shiftCount = countLeadingZeros64( zSig0 );
- shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
- zExp -= shiftCount;
- return
- roundAndPackFloatx80( roundData, zSign, zExp, zSig0, zSig1 );
- }
- #endif
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the 32-bit two's complement integer `a' to
- the single-precision floating-point format. The conversion is performed
- according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 int32_to_float32(struct roundingData *roundData, int32 a)
- {
- flag zSign;
- if ( a == 0 ) return 0;
- if ( a == 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
- zSign = ( a < 0 );
- return normalizeRoundAndPackFloat32( roundData, zSign, 0x9C, zSign ? - a : a );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the 32-bit two's complement integer `a' to
- the double-precision floating-point format. The conversion is performed
- according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float64 int32_to_float64( int32 a )
- {
- flag aSign;
- uint32 absA;
- int8 shiftCount;
- bits64 zSig;
- if ( a == 0 ) return 0;
- aSign = ( a < 0 );
- absA = aSign ? - a : a;
- shiftCount = countLeadingZeros32( absA ) + 21;
- zSig = absA;
- return packFloat64( aSign, 0x432 - shiftCount, zSig<<shiftCount );
- }
- #ifdef FLOATX80
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the 32-bit two's complement integer `a'
- to the extended double-precision floating-point format. The conversion
- is performed according to the IEC/IEEE Standard for Binary Floating-point
- Arithmetic.
- -------------------------------------------------------------------------------
- */
- floatx80 int32_to_floatx80( int32 a )
- {
- flag zSign;
- uint32 absA;
- int8 shiftCount;
- bits64 zSig;
- if ( a == 0 ) return packFloatx80( 0, 0, 0 );
- zSign = ( a < 0 );
- absA = zSign ? - a : a;
- shiftCount = countLeadingZeros32( absA ) + 32;
- zSig = absA;
- return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
- }
- #endif
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the single-precision floating-point value
- `a' to the 32-bit two's complement integer format. The conversion is
- performed according to the IEC/IEEE Standard for Binary Floating-point
- Arithmetic---which means in particular that the conversion is rounded
- according to the current rounding mode. If `a' is a NaN, the largest
- positive integer is returned. Otherwise, if the conversion overflows, the
- largest integer with the same sign as `a' is returned.
- -------------------------------------------------------------------------------
- */
- int32 float32_to_int32( struct roundingData *roundData, float32 a )
- {
- flag aSign;
- int16 aExp, shiftCount;
- bits32 aSig;
- bits64 zSig;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
- if ( aExp ) aSig |= 0x00800000;
- shiftCount = 0xAF - aExp;
- zSig = aSig;
- zSig <<= 32;
- if ( 0 < shiftCount ) shift64RightJamming( zSig, shiftCount, &zSig );
- return roundAndPackInt32( roundData, aSign, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the single-precision floating-point value
- `a' to the 32-bit two's complement integer format. The conversion is
- performed according to the IEC/IEEE Standard for Binary Floating-point
- Arithmetic, except that the conversion is always rounded toward zero. If
- `a' is a NaN, the largest positive integer is returned. Otherwise, if the
- conversion overflows, the largest integer with the same sign as `a' is
- returned.
- -------------------------------------------------------------------------------
- */
- int32 float32_to_int32_round_to_zero( float32 a )
- {
- flag aSign;
- int16 aExp, shiftCount;
- bits32 aSig;
- int32 z;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- shiftCount = aExp - 0x9E;
- if ( 0 <= shiftCount ) {
- if ( a == 0xCF000000 ) return 0x80000000;
- float_raise( float_flag_invalid );
- if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
- return 0x80000000;
- }
- else if ( aExp <= 0x7E ) {
- if ( aExp | aSig ) float_raise( float_flag_inexact );
- return 0;
- }
- aSig = ( aSig | 0x00800000 )<<8;
- z = aSig>>( - shiftCount );
- if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
- float_raise( float_flag_inexact );
- }
- return aSign ? - z : z;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the single-precision floating-point value
- `a' to the double-precision floating-point format. The conversion is
- performed according to the IEC/IEEE Standard for Binary Floating-point
- Arithmetic.
- -------------------------------------------------------------------------------
- */
- float64 float32_to_float64( float32 a )
- {
- flag aSign;
- int16 aExp;
- bits32 aSig;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- if ( aExp == 0xFF ) {
- if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) );
- return packFloat64( aSign, 0x7FF, 0 );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- --aExp;
- }
- return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
- }
- #ifdef FLOATX80
- /*
- -------------------------------------------------------------------------------
- Returns the result of converting the single-precision floating-point value
- `a' to the extended double-precision floating-point format. The conversion
- is performed according to the IEC/IEEE Standard for Binary Floating-point
- Arithmetic.
- -------------------------------------------------------------------------------
- */
- floatx80 float32_to_floatx80( float32 a )
- {
- flag aSign;
- int16 aExp;
- bits32 aSig;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- if ( aExp == 0xFF ) {
- if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a ) );
- return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- }
- aSig |= 0x00800000;
- return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
- }
- #endif
- /*
- -------------------------------------------------------------------------------
- Rounds the single-precision floating-point value `a' to an integer, and
- returns the result as a single-precision floating-point value. The
- operation is performed according to the IEC/IEEE Standard for Binary
- Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_round_to_int( struct roundingData *roundData, float32 a )
- {
- flag aSign;
- int16 aExp;
- bits32 lastBitMask, roundBitsMask;
- int8 roundingMode;
- float32 z;
- aExp = extractFloat32Exp( a );
- if ( 0x96 <= aExp ) {
- if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
- return propagateFloat32NaN( a, a );
- }
- return a;
- }
- roundingMode = roundData->mode;
- if ( aExp <= 0x7E ) {
- if ( (bits32) ( a<<1 ) == 0 ) return a;
- roundData->exception |= float_flag_inexact;
- aSign = extractFloat32Sign( a );
- switch ( roundingMode ) {
- case float_round_nearest_even:
- if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
- return packFloat32( aSign, 0x7F, 0 );
- }
- break;
- case float_round_down:
- return aSign ? 0xBF800000 : 0;
- case float_round_up:
- return aSign ? 0x80000000 : 0x3F800000;
- }
- return packFloat32( aSign, 0, 0 );
- }
- lastBitMask = 1;
- lastBitMask <<= 0x96 - aExp;
- roundBitsMask = lastBitMask - 1;
- z = a;
- if ( roundingMode == float_round_nearest_even ) {
- z += lastBitMask>>1;
- if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
- }
- else if ( roundingMode != float_round_to_zero ) {
- if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
- z += roundBitsMask;
- }
- }
- z &= ~ roundBitsMask;
- if ( z != a ) roundData->exception |= float_flag_inexact;
- return z;
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of adding the absolute values of the single-precision
- floating-point values `a' and `b'. If `zSign' is true, the sum is negated
- before being returned. `zSign' is ignored if the result is a NaN. The
- addition is performed according to the IEC/IEEE Standard for Binary
- Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- static float32 addFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign )
- {
- int16 aExp, bExp, zExp;
- bits32 aSig, bSig, zSig;
- int16 expDiff;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- bSig = extractFloat32Frac( b );
- bExp = extractFloat32Exp( b );
- expDiff = aExp - bExp;
- aSig <<= 6;
- bSig <<= 6;
- if ( 0 < expDiff ) {
- if ( aExp == 0xFF ) {
- if ( aSig ) return propagateFloat32NaN( a, b );
- return a;
- }
- if ( bExp == 0 ) {
- --expDiff;
- }
- else {
- bSig |= 0x20000000;
- }
- shift32RightJamming( bSig, expDiff, &bSig );
- zExp = aExp;
- }
- else if ( expDiff < 0 ) {
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- return packFloat32( zSign, 0xFF, 0 );
- }
- if ( aExp == 0 ) {
- ++expDiff;
- }
- else {
- aSig |= 0x20000000;
- }
- shift32RightJamming( aSig, - expDiff, &aSig );
- zExp = bExp;
- }
- else {
- if ( aExp == 0xFF ) {
- if ( aSig | bSig ) return propagateFloat32NaN( a, b );
- return a;
- }
- if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
- zSig = 0x40000000 + aSig + bSig;
- zExp = aExp;
- goto roundAndPack;
- }
- aSig |= 0x20000000;
- zSig = ( aSig + bSig )<<1;
- --zExp;
- if ( (sbits32) zSig < 0 ) {
- zSig = aSig + bSig;
- ++zExp;
- }
- roundAndPack:
- return roundAndPackFloat32( roundData, zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of subtracting the absolute values of the single-
- precision floating-point values `a' and `b'. If `zSign' is true, the
- difference is negated before being returned. `zSign' is ignored if the
- result is a NaN. The subtraction is performed according to the IEC/IEEE
- Standard for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- static float32 subFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign )
- {
- int16 aExp, bExp, zExp;
- bits32 aSig, bSig, zSig;
- int16 expDiff;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- bSig = extractFloat32Frac( b );
- bExp = extractFloat32Exp( b );
- expDiff = aExp - bExp;
- aSig <<= 7;
- bSig <<= 7;
- if ( 0 < expDiff ) goto aExpBigger;
- if ( expDiff < 0 ) goto bExpBigger;
- if ( aExp == 0xFF ) {
- if ( aSig | bSig ) return propagateFloat32NaN( a, b );
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- if ( aExp == 0 ) {
- aExp = 1;
- bExp = 1;
- }
- if ( bSig < aSig ) goto aBigger;
- if ( aSig < bSig ) goto bBigger;
- return packFloat32( roundData->mode == float_round_down, 0, 0 );
- bExpBigger:
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- return packFloat32( zSign ^ 1, 0xFF, 0 );
- }
- if ( aExp == 0 ) {
- ++expDiff;
- }
- else {
- aSig |= 0x40000000;
- }
- shift32RightJamming( aSig, - expDiff, &aSig );
- bSig |= 0x40000000;
- bBigger:
- zSig = bSig - aSig;
- zExp = bExp;
- zSign ^= 1;
- goto normalizeRoundAndPack;
- aExpBigger:
- if ( aExp == 0xFF ) {
- if ( aSig ) return propagateFloat32NaN( a, b );
- return a;
- }
- if ( bExp == 0 ) {
- --expDiff;
- }
- else {
- bSig |= 0x40000000;
- }
- shift32RightJamming( bSig, expDiff, &bSig );
- aSig |= 0x40000000;
- aBigger:
- zSig = aSig - bSig;
- zExp = aExp;
- normalizeRoundAndPack:
- --zExp;
- return normalizeRoundAndPackFloat32( roundData, zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of adding the single-precision floating-point values `a'
- and `b'. The operation is performed according to the IEC/IEEE Standard for
- Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_add( struct roundingData *roundData, float32 a, float32 b )
- {
- flag aSign, bSign;
- aSign = extractFloat32Sign( a );
- bSign = extractFloat32Sign( b );
- if ( aSign == bSign ) {
- return addFloat32Sigs( roundData, a, b, aSign );
- }
- else {
- return subFloat32Sigs( roundData, a, b, aSign );
- }
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of subtracting the single-precision floating-point values
- `a' and `b'. The operation is performed according to the IEC/IEEE Standard
- for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_sub( struct roundingData *roundData, float32 a, float32 b )
- {
- flag aSign, bSign;
- aSign = extractFloat32Sign( a );
- bSign = extractFloat32Sign( b );
- if ( aSign == bSign ) {
- return subFloat32Sigs( roundData, a, b, aSign );
- }
- else {
- return addFloat32Sigs( roundData, a, b, aSign );
- }
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of multiplying the single-precision floating-point values
- `a' and `b'. The operation is performed according to the IEC/IEEE Standard
- for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_mul( struct roundingData *roundData, float32 a, float32 b )
- {
- flag aSign, bSign, zSign;
- int16 aExp, bExp, zExp;
- bits32 aSig, bSig;
- bits64 zSig64;
- bits32 zSig;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- bSig = extractFloat32Frac( b );
- bExp = extractFloat32Exp( b );
- bSign = extractFloat32Sign( b );
- zSign = aSign ^ bSign;
- if ( aExp == 0xFF ) {
- if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
- return propagateFloat32NaN( a, b );
- }
- if ( ( bExp | bSig ) == 0 ) {
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- return packFloat32( zSign, 0xFF, 0 );
- }
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- if ( ( aExp | aSig ) == 0 ) {
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- return packFloat32( zSign, 0xFF, 0 );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- }
- if ( bExp == 0 ) {
- if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
- normalizeFloat32Subnormal( bSig, &bExp, &bSig );
- }
- zExp = aExp + bExp - 0x7F;
- aSig = ( aSig | 0x00800000 )<<7;
- bSig = ( bSig | 0x00800000 )<<8;
- shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
- zSig = zSig64;
- if ( 0 <= (sbits32) ( zSig<<1 ) ) {
- zSig <<= 1;
- --zExp;
- }
- return roundAndPackFloat32( roundData, zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the result of dividing the single-precision floating-point value `a'
- by the corresponding value `b'. The operation is performed according to the
- IEC/IEEE Standard for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_div( struct roundingData *roundData, float32 a, float32 b )
- {
- flag aSign, bSign, zSign;
- int16 aExp, bExp, zExp;
- bits32 aSig, bSig, zSig;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- bSig = extractFloat32Frac( b );
- bExp = extractFloat32Exp( b );
- bSign = extractFloat32Sign( b );
- zSign = aSign ^ bSign;
- if ( aExp == 0xFF ) {
- if ( aSig ) return propagateFloat32NaN( a, b );
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- return packFloat32( zSign, 0xFF, 0 );
- }
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- return packFloat32( zSign, 0, 0 );
- }
- if ( bExp == 0 ) {
- if ( bSig == 0 ) {
- if ( ( aExp | aSig ) == 0 ) {
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- roundData->exception |= float_flag_divbyzero;
- return packFloat32( zSign, 0xFF, 0 );
- }
- normalizeFloat32Subnormal( bSig, &bExp, &bSig );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- }
- zExp = aExp - bExp + 0x7D;
- aSig = ( aSig | 0x00800000 )<<7;
- bSig = ( bSig | 0x00800000 )<<8;
- if ( bSig <= ( aSig + aSig ) ) {
- aSig >>= 1;
- ++zExp;
- }
- {
- bits64 tmp = ( (bits64) aSig )<<32;
- do_div( tmp, bSig );
- zSig = tmp;
- }
- if ( ( zSig & 0x3F ) == 0 ) {
- zSig |= ( ( (bits64) bSig ) * zSig != ( (bits64) aSig )<<32 );
- }
- return roundAndPackFloat32( roundData, zSign, zExp, zSig );
- }
- /*
- -------------------------------------------------------------------------------
- Returns the remainder of the single-precision floating-point value `a'
- with respect to the corresponding value `b'. The operation is performed
- according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
- -------------------------------------------------------------------------------
- */
- float32 float32_rem( struct roundingData *roundData, float32 a, float32 b )
- {
- flag aSign, bSign, zSign;
- int16 aExp, bExp, expDiff;
- bits32 aSig, bSig;
- bits32 q;
- bits64 aSig64, bSig64, q64;
- bits32 alternateASig;
- sbits32 sigMean;
- aSig = extractFloat32Frac( a );
- aExp = extractFloat32Exp( a );
- aSign = extractFloat32Sign( a );
- bSig = extractFloat32Frac( b );
- bExp = extractFloat32Exp( b );
- bSign = extractFloat32Sign( b );
- if ( aExp == 0xFF ) {
- if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
- return propagateFloat32NaN( a, b );
- }
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- if ( bExp == 0xFF ) {
- if ( bSig ) return propagateFloat32NaN( a, b );
- return a;
- }
- if ( bExp == 0 ) {
- if ( bSig == 0 ) {
- roundData->exception |= float_flag_invalid;
- return float32_default_nan;
- }
- normalizeFloat32Subnormal( bSig, &bExp, &bSig );
- }
- if ( aExp == 0 ) {
- if ( aSig == 0 ) return a;
- normalizeFloat32Subnormal( aSig, &aExp, &aSig );
- }
- expDiff = aExp - bExp;
- aSig |= 0x00800000;
- bSig |= 0x00800000;
- if ( expDiff < 32 ) {
- aSig <<= 8;
- bSig <<= 8;
- if ( expDiff < 0 ) {
- if ( expDiff < -1 ) return a;
- aSig >>= 1;
- }
- q = ( bSig <= aSig );
- if ( q ) aSig -= bSig;
- if ( 0 < expDiff ) {
- bits64 tmp = ( (bits64) aSig )<<32;
- do_div( tmp, bSig );
- q = tmp;
- q >>= 32 - expDiff;
- bSig >>= 2;
- aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
- }
- else {
- aSig >>= 2;
- bSig >>= 2;
- }
- }
- else {
- if ( bSig <= aSig ) aSig -= bSig;
- aSig64 = ( (bits64) aSig )<<40;
- bSig64 = ( (bits64) bSig )<<40;
- expDiff -= 64;
- while ( 0 < expDiff ) {
- q64 = estimateDiv128To64( aSig64, 0, bSig64 );
- q64 = ( 2 < q64 ) ? q64 - 2 : 0;
- aSig64 = - ( ( bSig * q64 )<<38 );
- expDiff -= 62;
- }
- expDiff += 64;
- q64 = estimateDiv128To64( aSig64, 0, bSig64 );
- q64 = ( 2 < q64 ) ? q64 - 2 : 0;
- q = q64>>( 64 - expDiff );
- bSig <<= 6;
- aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) )