cd910107c74c046d65c4369d8f41ff6b747303b6
5 #include "primitives.h"
7 #include "specialize.h"
10 float64_t
f64_sqrt( float64_t a
)
16 uint_fast64_t sigA
, uiZ
;
17 struct exp16_sig64 normExpSig
;
21 struct uint128 term
, rem
;
26 signA
= signF64UI( uiA
);
27 expA
= expF64UI( uiA
);
28 sigA
= fracF64UI( uiA
);
29 if ( expA
== 0x7FF ) {
31 uiZ
= softfloat_propagateNaNF64UI( uiA
, 0 );
34 if ( ! signA
) return a
;
38 if ( ! ( expA
| sigA
) ) return a
;
42 if ( ! sigA
) return a
;
43 normExpSig
= softfloat_normSubnormalF64Sig( sigA
);
44 expA
= normExpSig
.exp
;
45 sigA
= normExpSig
.sig
;
47 expZ
= ( ( expA
- 0x3FF )>>1 ) + 0x3FE;
48 sigA
|= UINT64_C( 0x0010000000000000 );
49 sigZ32
= softfloat_estimateSqrt32( expA
, sigA
>>21 );
50 sigA
<<= 9 - ( expA
& 1 );
52 softfloat_estimateDiv128To64( sigA
, 0, (uint_fast64_t) sigZ32
<<32 )
53 + ( (uint_fast64_t) sigZ32
<<30 );
54 if ( ( sigZ
& 0x1FF ) <= 5 ) {
55 term
= softfloat_mul64To128( sigZ
, sigZ
);
56 rem
= softfloat_sub128( sigA
, 0, term
.v64
, term
.v0
);
57 while ( UINT64_C( 0x8000000000000000 ) <= rem
.v64
) {
61 rem
.v64
, rem
.v0
, sigZ
>>63, (uint64_t) ( sigZ
<<1 ) );
63 sigZ
|= ( ( rem
.v64
| rem
.v0
) != 0 );
65 return softfloat_roundPackToF64( 0, expZ
, sigZ
);
67 softfloat_raiseFlags( softfloat_flag_invalid
);
68 uiZ
= defaultNaNF64UI
;