--- /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