Add support for Float32 and Float16

[mfherbst]

This PR adds transparent support for Float32 and Float16 datatypes. Closes #32.

• I experimented a bit with timings. The differences between using MPFR and doing the computation in the respective larger IEEE type and than rounding the result are very small. Happy to discuss alternative solutions.
• I also did some small cleanup (removed duplicate definitions, some tabs)

[dpsanders]

Great idea, thanks!

However, there's definitely a huge performance difference between using Float64 and BigFloat as the intermediate step (as I expected there must be):

function cos1(x::Float32, r::RoundingMode)
    return Float32(CRlibm.cos(Float64(x), r), r)
end

function cos2(x::Float32, r::RoundingMode)
           y = setprecision(BigFloat, 24) do
                   setrounding(BigFloat, r) do
                       cos(BigFloat(x))
                   end
           end

           return Float32(y, r)
       end
end

using BenchmarkTools

julia> @btime cos1($x, RoundDown)
  20.465 ns (0 allocations: 0 bytes)
0.9950041f0

julia> @btime cos2($x, RoundDown)
  648.976 ns (7 allocations: 176 bytes)
0.9950041f0

[mfherbst]

Ok, thanks for checking that. I found that a little odd, too.

I'll change the code accordingly (i.e. only use standard FP types when possible). Any reason you use Float32(CRlibm.cos(Float64(x), r), r) in cos1? I would expect Float32(cos(Float64(x), r)) to give the same answer (and be even faster).

[dpsanders]

Maybe you should convert straight to Float64 here?

[mfherbst]

I changed it for Float64, but I doubt it really matters ... but I'm also not really an expert in FP arithmetic either. Do you have any reason why you think Float64 could be better?

[dpsanders]

Because if you use Float32, it will just convert it to Float64 in the following line ;) Although the compiler will probably compile that away, so feel free to leave it.

[mfherbst]

I don't think this is what happens here (unless I'm wrong). At least I try to avoid that by calling e.g. Base.cos and not just cos. I quickly verified using a print-statement: So there is no dispatch to CRlibm.cos in line 136 or 137.

[dpsanders]

Ah, there's a subtlety here: CRlibm doesn't import cos from Base, so "cos" within the module refers to the internal cos.

We may want to write it explicitly as CRlibm.cos to make this clear.

[mfherbst]

Yes, sure, but Base.cos should not refer to CRlibm.cos, right? That would be pretty counter-intuitive. So the @eval in line 136 should generate calls to Base.cos only and not back to CRlibm.cos, no?

To clarify: If I am not mistaken, this means that a call to CRlibm.cos(Float16(1.0)) will, in the current version of the code, widen to Float64 and then dispatch to Base.cos(::Float64) and never get passed on to the CRlib C code.

[dpsanders]

There is no use of Base.cos within the package.
Any use of cos within the CRlibm package refers to the functions defined inside the package. Base.cos never comes into the picture at all.

However, in IntervalArithmetic.jl we define a new method Base.cos(x, r), with a RoundingMode r, to use CRlibm.cos by default. (This is, strictly speaking, type piracy.)

You should be able to check this using the debugger in Juno or using Cthulhu.jl to check that the correct methods are being called.

[mfherbst]

Thanks for pointing out Cthulhu.jl. That's quite a neat package for the @code_warntype mess.

Ah I see, I think we might be misunderstanding each other here ... I am talking about the state of CRlibm after this PR is merged and this PR does introduce an explicit call to functions such as Base.cos, Base.sin and so on. If this is not what you want, please let me know, I'm happy to change this.

For example, one line where this happens is line 137. In the debugger a call to CRlibm.cos with a Float32 type dispatches to Base.cos (with the data widened to Float64):

julia> using CRlibm; using Debugger; @breakpoint Base.cos(1.0)
cos(x::T) where T<:Union{Float32, Float64} in Base.Math at special/trig.jl:100

julia> @run CRlibm.cos(Float32(1.))
Hit breakpoint:
[ Info: tracking Base
In cos(x) at special/trig.jl:100
  99  function cos(x::T) where T<:Union{Float32, Float64}
>100      absx = abs(x)
 101      if absx < T(pi)/4
 102          if absx < sqrt(eps(T)/T(2.0))
 103              return T(1.0)
 104          end

About to run: Core.NewvarNode(:(_4))
1|debug> bt
[1] cos(x) at special/trig.jl:100
  | x::Float64 = 1.0
  | T::DataType = Float64
[2] cos(x, r) at CRlibm.jl/src/CRlibm.jl:137
  | x::Float32 = 1.0f0
  | r::RoundingMode{:Nearest} = RoundingMode{:Nearest}()
[3] cos(x) at CRlibm.jl/src/CRlibm.jl:132
  | x::Float32 = 1.0f0

[dpsanders]

OK I had not noticed that change, sorry.

But Base.cos does not do correct rounding, so can't be used here. You need to replace it with a rounded CRlibm.cos call, as I did in my benchmark code.

[mfherbst]

All right ... did that. All Base.cos is gone again.

[dpsanders]

cos(x, r), i.e. Base.cos(x, r), is not defined in CRlibm -- i.e. CRlibm (deliberately) does not export the functions it defines. Base.cos(x, r) is only defined in IntervalArithmetic.jl.

[dpsanders]

We/you should probably update the Travis and Appveyor files to test on Julia 1.0, 1.3 and nightly, and allow failures on nightly.

[mfherbst]

Sorry about the confusion. I talked about Base.cos when I wrote cos in my post. I can take care of updating Julia.

[dpsanders]

At least part of the test failures seems to be due to the fact that .travis.yml is too old.
You can try copying the one from e.g. IntervalArithmetic.jl and modifying it. Or I will try to get round to it soon.

[dpsanders]

Thanks, LGTM! Is this ready to merge?

[mfherbst]

Yes, sorry. It is from my point of view 😄.

[dpsanders]

Many thanks for this contribution!