fix subnormal sqrt and cbrt

[t-bltg]

Fix #166 (tests added).

Also remove .DS_Store files and add them to .gitignore.

• benchmark performance

PR

julia> using BenchmarkTools, DoubleFloats
julia> x = Double64(1e-10);
julia> @btime sqrt($x)
  230.792 ns (0 allocations: 0 bytes)

master

julia> using BenchmarkTools, DoubleFloats
julia> x = Double64(1e-10)
julia> @btime sqrt($x)
  215.868 ns (0 allocations: 0 bytes)

[codecov]

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.07 🎉

Comparison is base (ef689cc) 49.12% compared to head (5b99418) 49.19%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #167      +/-   ##
==========================================
+ Coverage   49.12%   49.19%   +0.07%     
==========================================
  Files          63       63              
  Lines        3422     3425       +3     
==========================================
+ Hits         1681     1685       +4     
+ Misses       1741     1740       -1     

Impacted Files	Coverage Δ
src/math/ops/op_dd_dd.jl	91.01% <100.00%> (+0.31%)	⬆️
src/type/predicates.jl	65.51% <0.00%> (+3.44%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[JeffreySarnoff]

for sqrt using

HI(x) < 5.562684646268013e-309 && return DoubleFloat{T}(sqrt(HI(x)), zero(T))

for cbrt using

HI(x) < 7.458340731200208e-155 && return DoubleFloat{T}(cbrt(HI(x)), zero(T))

will be faster (and shrinks the special case space).

[t-bltg]

How are these subnormal numbers determined ?

[JeffreySarnoff]

By hand. For example, at the end of the cbrt search

julia> x=prevfloat(y,3_545_284); cbrt(Double64(x))
NaN

julia> x=prevfloat(y,3_545_283); cbrt(Double64(x))
4.209340649576657e-52

julia> x
7.458340731200208e-155

For cbrt, the value is not subnormal.

[t-bltg]

I see, thanks.
Are we sure that we want to hardcode these numbers in the algorithm itself ?

Especially since these numbers are determined from the output of algorithm itself, any change to the aforementioned algorithm is likely to change these values and thus breaking the code again: I find it error prone and brittle.

[JeffreySarnoff]

For the algorithm as designed, they are the precise values -- changing the design (esp for cbrt) would change that, but they are empirical. The alternative would be to check the result for NaN and then apply the approximation. That would be more general at the cost of slowdown (checking all inputs for small inputs and checking results for small inputs to see if they are NaN).

…

On Wed, Mar 8, 2023 at 8:33 AM t-bltg ***@***.***> wrote: I see, thanks. Are we sure that we want to hardcode these numbers in the algorithm itself ? Especially since these numbers are determined from the output of algorithm itself, any change to the aforementioned algorithm is likely to change these values and thus breaking the code again: I find it error prone and brittle. — Reply to this email directly, view it on GitHub <#167 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAM2VRSVGSDWNL6ZVHOIGLDW3CDCRANCNFSM6AAAAAAVTWICMY> . You are receiving this because you commented.Message ID: ***@***.***>

[t-bltg]

I've pushed an alternative in 5b99418, although I must admit that I find it less readable.

[t-bltg]

By hand. For example, at the end of the cbrt search

What is y in your example ?

[JeffreySarnoff]

y is an approximation to the final value x, 3.351952e153 was found by informed search.

julia> yy = (floatmin(Float64) * 3.351952e153)
7.458340770170931e-155

julia> xx=prevfloat(yy, 3_545_283); cbrt(Double64(xx))
4.209340655803547e-52