|
@@ -2333,3 +2333,163 @@ float64 float64_sqrt( struct roundingData *roundData, float64 a )
|
|
|
flag aSign;
|
|
|
int16 aExp, zExp;
|
|
|
bits64 aSig, zSig;
|
|
|
+ bits64 rem0, rem1, term0, term1; //, shiftedRem;
|
|
|
+ //float64 z;
|
|
|
+
|
|
|
+ aSig = extractFloat64Frac( a );
|
|
|
+ aExp = extractFloat64Exp( a );
|
|
|
+ aSign = extractFloat64Sign( a );
|
|
|
+ if ( aExp == 0x7FF ) {
|
|
|
+ if ( aSig ) return propagateFloat64NaN( a, a );
|
|
|
+ if ( ! aSign ) return a;
|
|
|
+ roundData->exception |= float_flag_invalid;
|
|
|
+ return float64_default_nan;
|
|
|
+ }
|
|
|
+ if ( aSign ) {
|
|
|
+ if ( ( aExp | aSig ) == 0 ) return a;
|
|
|
+ roundData->exception |= float_flag_invalid;
|
|
|
+ return float64_default_nan;
|
|
|
+ }
|
|
|
+ if ( aExp == 0 ) {
|
|
|
+ if ( aSig == 0 ) return 0;
|
|
|
+ normalizeFloat64Subnormal( aSig, &aExp, &aSig );
|
|
|
+ }
|
|
|
+ zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
|
|
|
+ aSig |= LIT64( 0x0010000000000000 );
|
|
|
+ zSig = estimateSqrt32( aExp, aSig>>21 );
|
|
|
+ zSig <<= 31;
|
|
|
+ aSig <<= 9 - ( aExp & 1 );
|
|
|
+ zSig = estimateDiv128To64( aSig, 0, zSig ) + zSig + 2;
|
|
|
+ if ( ( zSig & 0x3FF ) <= 5 ) {
|
|
|
+ if ( zSig < 2 ) {
|
|
|
+ zSig = LIT64( 0xFFFFFFFFFFFFFFFF );
|
|
|
+ }
|
|
|
+ else {
|
|
|
+ aSig <<= 2;
|
|
|
+ mul64To128( zSig, zSig, &term0, &term1 );
|
|
|
+ sub128( aSig, 0, term0, term1, &rem0, &rem1 );
|
|
|
+ while ( (sbits64) rem0 < 0 ) {
|
|
|
+ --zSig;
|
|
|
+ shortShift128Left( 0, zSig, 1, &term0, &term1 );
|
|
|
+ term1 |= 1;
|
|
|
+ add128( rem0, rem1, term0, term1, &rem0, &rem1 );
|
|
|
+ }
|
|
|
+ zSig |= ( ( rem0 | rem1 ) != 0 );
|
|
|
+ }
|
|
|
+ }
|
|
|
+ shift64RightJamming( zSig, 1, &zSig );
|
|
|
+ return roundAndPackFloat64( roundData, 0, zExp, zSig );
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+/*
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+Returns 1 if the double-precision floating-point value `a' is equal to the
|
|
|
+corresponding value `b', and 0 otherwise. The comparison is performed
|
|
|
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+*/
|
|
|
+flag float64_eq( float64 a, float64 b )
|
|
|
+{
|
|
|
+
|
|
|
+ if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
|
|
|
+ || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
|
|
|
+ ) {
|
|
|
+ if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
|
|
|
+ float_raise( float_flag_invalid );
|
|
|
+ }
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+ return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+/*
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+Returns 1 if the double-precision floating-point value `a' is less than or
|
|
|
+equal to the corresponding value `b', and 0 otherwise. The comparison is
|
|
|
+performed according to the IEC/IEEE Standard for Binary Floating-point
|
|
|
+Arithmetic.
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+*/
|
|
|
+flag float64_le( float64 a, float64 b )
|
|
|
+{
|
|
|
+ flag aSign, bSign;
|
|
|
+
|
|
|
+ if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
|
|
|
+ || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
|
|
|
+ ) {
|
|
|
+ float_raise( float_flag_invalid );
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+ aSign = extractFloat64Sign( a );
|
|
|
+ bSign = extractFloat64Sign( b );
|
|
|
+ if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
|
|
|
+ return ( a == b ) || ( aSign ^ ( a < b ) );
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+/*
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+Returns 1 if the double-precision floating-point value `a' is less than
|
|
|
+the corresponding value `b', and 0 otherwise. The comparison is performed
|
|
|
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+*/
|
|
|
+flag float64_lt( float64 a, float64 b )
|
|
|
+{
|
|
|
+ flag aSign, bSign;
|
|
|
+
|
|
|
+ if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
|
|
|
+ || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
|
|
|
+ ) {
|
|
|
+ float_raise( float_flag_invalid );
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+ aSign = extractFloat64Sign( a );
|
|
|
+ bSign = extractFloat64Sign( b );
|
|
|
+ if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
|
|
|
+ return ( a != b ) && ( aSign ^ ( a < b ) );
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+/*
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+Returns 1 if the double-precision floating-point value `a' is equal to the
|
|
|
+corresponding value `b', and 0 otherwise. The invalid exception is raised
|
|
|
+if either operand is a NaN. Otherwise, the comparison is performed
|
|
|
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+*/
|
|
|
+flag float64_eq_signaling( float64 a, float64 b )
|
|
|
+{
|
|
|
+
|
|
|
+ if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
|
|
|
+ || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
|
|
|
+ ) {
|
|
|
+ float_raise( float_flag_invalid );
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+ return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
|
|
|
+
|
|
|
+}
|
|
|
+
|
|
|
+/*
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+Returns 1 if the double-precision floating-point value `a' is less than or
|
|
|
+equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not
|
|
|
+cause an exception. Otherwise, the comparison is performed according to the
|
|
|
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
|
|
|
+-------------------------------------------------------------------------------
|
|
|
+*/
|
|
|
+flag float64_le_quiet( float64 a, float64 b )
|
|
|
+{
|
|
|
+ flag aSign, bSign;
|
|
|
+ //int16 aExp, bExp;
|
|
|
+
|
|
|
+ if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
|
|
|
+ || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
|
|
|
+ ) {
|
|
|
+ /* Do nothing, even if NaN as we're quiet */
|
|
|
+ return 0;
|
|
|
+ }
|