musl erf/erfc port

[musm]

musl port

There is also the possibility of a sleef port and to benchmark the two implementations.

Edit:
These implementations are just as fast as Base.

[simonbyrne]

This can be done more efficiently by

function GET_HIGH_WORD(d::Float64)
    u = reinterpret(UInt64, d)
    (u >> 32) % UInt32
end

[simonbyrne]

Aside from the notes by @alyst and myself, this actually looks rather reasonable. Try some profiling to see where the remaining bottlenecks are.

[alyst]

SET_LOW_WORD(d::Float64, lo::UInt32) = reinterpret(Float64, reinterpret(UInt64, d) & 0xFFFFFFFF00000000 | lo)

should be faster.

Also note that in general to get faster code you should avoid changing the type of the variable throughout the function as it makes Julia generating type-agnostic (much slower) code:

julia> @code_warntype SET_LOW_WORD(1.0, UInt(5))
Variables:
  #self#::#SET_LOW_WORD
  d@_2::Float64
  lo@_3::UInt64
  db::String
  lb::String
  d@_6::Any
  lo@_7::Any

Body:
  begin 
      lo@_7::Any = lo@_3::UInt64
      d@_6::Any = d@_2::Float64
      db::String = $(Expr(:invoke, LambdaInfo for bin(::UInt64, ::Int64, ::Bool), :(Base.bin), :((Base.box)(UInt64,d@_6::Float64)), 64, false)) # line 3:
      lo@_7::Any = (Base.box)(UInt32,(Base.checked_trunc_uint)(UInt32,lo@_7::UInt64)) # line 4:
      lb::String = $(Expr(:invoke, LambdaInfo for bin(::UInt32, ::Int64, ::Bool), :(Base.bin), :((Base.box)(UInt32,lo@_7::UInt32)), 32, false)) # line 5:
      SSAValue(0) = $(Expr(:invoke, LambdaInfo for getindex(::String, ::UnitRange{Int64}), :(Main.getindex), :(db), :($(Expr(:new, UnitRange{Int64}, 1, :((Base.select_value)((Base.sle_int)(1,32)::Bool,32,(Base.box)(Int64,(Base.sub_int)(1,1)))::Int64))))))
      d@_6::Any = $(Expr(:invoke, LambdaInfo for string(::String, ::String, ::Vararg{String,N}), :(Base.string), SSAValue(0), :(lb))) # line 6:
      d@_6::Any = $(Expr(:invoke, LambdaInfo for parse(::Type{Int64}, ::String, ::Int64), :(Main.parse), :(Main.Int), :(d@_6::String), 2)) # line 7:
      d@_6::Any = (Base.box)(Float64,d@_6::Int64) # line 8:
      return d@_6::Float64
  end::Float64

(::Any signifies type inference problem)

[musm]

Thank you @simonbyrne & @alyst

[tkelman]

should list this kind of thing in the top-level license file

[codecov-io]

Current coverage is 100% (diff: 100%)

No coverage report found for master at e4c2c15.

Powered by Codecov. Last update e4c2c15...8c35271

[tkelman]

add a blank line before this or it'll be highlighted as if it's part of the quote

[musm]

what's going on with codecov, it looks like it has my wrong username on there?

[ararslan]

In general, we should probably try to stick to the convention of constants in uppercase, types in upper camel case, and functions and variables in lowercase.

[simonbyrne]

This is perhaps also a good candidate as a utility function to add in a different file.

[tkelman]

I wonder whether we shouldn't take a lesson-learned from openlibm and try to make this an automated translation from the upstream musl source code, if that would be possible? That way future updates to pull in bug fixes would be simpler and less dependent on when we happen to do individual ports.

[musm]

Perhaps, but most of the bug reports are not related to math functions, which are infrequent.

[simonbyrne]

Do we go with this implementation what about sleef ?

I don't think it's a problem to have multiple implementations in here: indeed, I would encourage it. I suspect there will be some tradeoff between accuracy, performance and vectorisation between different implementations that would be useful to compare.

[alyst]

For me this code and Base.Math.@horner s ra0 ra1 ra2 ra3 ra4 ra5 ra6 ra7 generate the same @code_native, but @horner approach looks clearer. Yet maybe on other machines @horner would generate more efficient machine code.

[ararslan]

Thinking more about this, would it make more sense to have error functions as part of SpecialFunctions.jl?

[ararslan]

Why remove the spaces? IMO it's much more readable with them in there.

[ararslan]

probably easiest to turn on appveyor

@simonbyrne Would you be willing to put this repo on your AppVeyor account?

[simonbyrne]

This should be import, not using

[tkelman]

you're not extending it, so actually using is better if you're asking for a specific function

[simonbyrne]

Probably needs a copyright notice here: i.e. based on musl code, copyright whoever is on that code. Also would be good to add which git commit you used to make it easy to check for updates.

[musm]

With the latest changes both erf and erfc are just as fast as Base!
I have to thank alyst and simonbyrne for their help

[tkelman]

needs a copyright notice here: i.e. based on musl code, copyright whoever is on that code. Also would be good to add which git commit you used to make it easy to check for updates.

And please clean your commit messages up.

[musm]

Does that look ok now?

[ViralBShah]

For the first few functions, it would be great to see if a port from sleef gives us code that can benefit from SIMD. Would be good to compare and see how that goes.

+1 to @tkelman's comment about trying to do source-to-source wherever possible so that we can benefit from upstream as well. This is of course, good to have.

[simonbyrne]

Great, thanks for doing this, and your patience with the reviews.