+++ /dev/null
-\r
-#include <stdbool.h>\r
-#include <stdint.h>\r
-#include "platform.h"\r
-#include "primitives.h"\r
-#include "internals.h"\r
-#include "specialize.h"\r
-#include "softfloat.h"\r
-\r
-float32_t f32_sqrt( float32_t a )\r
-{\r
- union ui32_f32 uA;\r
- uint_fast32_t uiA;\r
- bool signA;\r
- int_fast16_t expA;\r
- uint_fast32_t sigA, uiZ;\r
- struct exp16_sig32 normExpSig;\r
- int_fast16_t expZ;\r
- uint_fast32_t sigZ;\r
- uint_fast64_t term, rem;\r
- union ui32_f32 uZ;\r
-\r
- uA.f = a;\r
- uiA = uA.ui;\r
- signA = signF32UI( uiA );\r
- expA = expF32UI( uiA );\r
- sigA = fracF32UI( uiA );\r
- if ( expA == 0xFF ) {\r
- if ( sigA ) {\r
- uiZ = softfloat_propagateNaNF32UI( uiA, 0 );\r
- goto uiZ;\r
- }\r
- if ( ! signA ) return a;\r
- goto invalid;\r
- }\r
- if ( signA ) {\r
- if ( ! ( expA | sigA ) ) return a;\r
- goto invalid;\r
- }\r
- if ( ! expA ) {\r
- if ( ! sigA ) return a;\r
- normExpSig = softfloat_normSubnormalF32Sig( sigA );\r
- expA = normExpSig.exp;\r
- sigA = normExpSig.sig;\r
- }\r
- expZ = ( ( expA - 0x7F )>>1 ) + 0x7E;\r
- sigA = ( sigA | 0x00800000 )<<8;\r
- sigZ = softfloat_estimateSqrt32( expA, sigA ) + 2;\r
- if ( ( sigZ & 0x7F ) <= 5 ) {\r
- if ( sigZ < 2 ) {\r
- sigZ = 0x7FFFFFFF;\r
- goto roundPack;\r
- }\r
- sigA >>= expA & 1;\r
- term = (uint_fast64_t) sigZ * sigZ;\r
- rem = ( (uint_fast64_t) sigA<<32 ) - term;\r
- while ( UINT64_C( 0x8000000000000000 ) <= rem ) {\r
- --sigZ;\r
- rem += ( (uint_fast64_t) sigZ<<1 ) | 1;\r
- }\r
- sigZ |= ( rem != 0 );\r
- }\r
- sigZ = softfloat_shortShift32Right1Jam( sigZ );\r
- roundPack:\r
- return softfloat_roundPackToF32( 0, expZ, sigZ );\r
- invalid:\r
- softfloat_raiseFlags( softfloat_flag_invalid );\r
- uiZ = defaultNaNF32UI;\r
- uiZ:\r
- uZ.ui = uiZ;\r
- return uZ.f;\r
-\r
-}\r
-\r