/*===============================================================================This C source file is part of the SoftFloat IEC/IEEE Floating-pointArithmetic Package, Release 2.Written by John R. Hauser.  This work was made possible in part by theInternational Computer Science Institute, located at Suite 600, 1947 CenterStreet, Berkeley, California 94704.  Funding was partially provided by theNational Science Foundation under grant MIP-9311980.  The original versionof this code was written as part of a project to build a fixed-point vectorprocessor in collaboration with the University of California at Berkeley,overseen by Profs. Nelson Morgan and John Wawrzynek.  More informationis available through the web page `http://HTTP.CS.Berkeley.EDU/~jhauser/arithmetic/softfloat.html'.THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable efforthas been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL ATTIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TOPERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANYAND 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) theyinclude prominent notice akin to these three paragraphs for those parts ofthis code that are retained.===============================================================================*/#include "fpa11.h"#include "milieu.h"#include "softfloat.h"/*-------------------------------------------------------------------------------Floating-point rounding mode, extended double-precision rounding precision,and exception flags.-------------------------------------------------------------------------------*/int8 float_rounding_mode = float_round_nearest_even;int8 floatx80_rounding_precision = 80;int8 float_exception_flags;/*-------------------------------------------------------------------------------Primitive arithmetic functions, including multi-word arithmetic, anddivision and square root approximations.  (Can be specialized to target ifdesired.)-------------------------------------------------------------------------------*/#include "softfloat-macros"/*-------------------------------------------------------------------------------Functions and definitions to determine:  (1) whether tininess for underflowis detected before or after rounding by default, (2) what (if anything)happens when exceptions are raised, (3) how signaling NaNs are distinguishedfrom quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNsare 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 6and 7, and returns the properly rounded 32-bit integer corresponding to theinput.  If `zSign' is nonzero, the input is negated before being convertedto an integer.  Bit 63 of `absZ' must be zero.  Ordinarily, the fixed-pointinput is simply rounded to an integer, with the inexact exception raised ifthe input cannot be represented exactly as an integer.  If the fixed-pointinput is too large, however, the invalid exception is raised and the largestpositive or negative integer is returned.-------------------------------------------------------------------------------*/static int32 roundAndPackInt32( flag zSign, bits64 absZ ){    int8 roundingMode;    flag roundNearestEven;    int8 roundIncrement, roundBits;    int32 z;    roundingMode = float_rounding_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 ) ) ) {        float_exception_flags |= float_flag_invalid;        return zSign ? 0x80000000 : 0x7FFFFFFF;    }    if ( roundBits ) float_exception_flags |= 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'.-------------------------------------------------------------------------------*/INLINE flag extractFloat32Sign( float32 a ){    return a>>31;}/*-------------------------------------------------------------------------------Normalizes the subnormal single-precision floating-point value representedby the denormalized significand `aSig'.  The normalized exponent andsignificand 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 asingle-precision floating-point value, returning the result.  After beingshifted into the proper positions, the three fields are simply addedtogether to form the result.  This means that any integer portion of `zSig'will be added into the exponent.  Since a properly normalized significandwill have an integer portion equal to 1, the `zExp' input should be 1 lessthan the desired result exponent whenever `zSig' is a complete, normalizedsignificand.-------------------------------------------------------------------------------*/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 abstractvalue is simply rounded and packed into the single-precision format, withthe inexact exception raised if the abstract input cannot be representedexactly.  If the abstract value is too large, however, the overflow andinexact exceptions are raised and an infinity or maximal finite value isreturned.  If the abstract value is too small, the input value is rounded toa subnormal number, and the underflow and inexact exceptions are raised ifthe abstract input cannot be represented exactly as a subnormal single-precision floating-point number.    The input significand `zSig' has its binary point between bits 30and 29, which is 7 bits to the left of the usual location.  This shiftedsignificand 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' isnormalized, `zExp' must be 1 less than the ``true'' floating-point exponent.The handling of underflow and overflow follows the IEC/IEEE Standard forBinary Floating-point Arithmetic.-------------------------------------------------------------------------------*/static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig ){    int8 roundingMode;    flag roundNearestEven;    int8 roundIncrement, roundBits;    flag isTiny;    roundingMode = float_rounding_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 ) )           ) {            float_raise( 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 ) float_raise( float_flag_underflow );        }    }    if ( roundBits ) float_exception_flags |= 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 inany way.  In all cases, `zExp' must be 1 less than the ``true'' floating-point exponent.-------------------------------------------------------------------------------*/static float32 normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig ){    int8 shiftCount;    shiftCount = countLeadingZeros32( zSig ) - 1;    return roundAndPackFloat32( 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'.-------------------------------------------------------------------------------*/INLINE flag extractFloat64Sign( float64 a ){    return a>>63;}/*-------------------------------------------------------------------------------Normalizes the subnormal double-precision floating-point value representedby the denormalized significand `aSig'.  The normalized exponent andsignificand 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 adouble-precision floating-point value, returning the result.  After beingshifted into the proper positions, the three fields are simply addedtogether to form the result.  This means that any integer portion of `zSig'will be added into the exponent.  Since a properly normalized significandwill have an integer portion equal to 1, the `zExp' input should be 1 lessthan the desired result exponent whenever `zSig' is a complete, normalizedsignificand.-------------------------------------------------------------------------------*/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 abstractvalue is simply rounded and packed into the double-precision format, withthe inexact exception raised if the abstract input cannot be representedexactly.  If the abstract value is too large, however, the overflow andinexact exceptions are raised and an infinity or maximal finite value isreturned.  If the abstract value is too small, the input value is rounded to