Add sincosd function

[ronisbr]

@StefanKarpinski and @KristofferC

This is the first PR as I was instructed. It only adds the sincosd function. After a lot of tests, it turns out that calling sind and cosd separately in an @inline function yielded a much better performance than considering the cases separately.

My first attempt was:

function sincosd(x::Real)
    if isinf(x)
        return throw(DomainError(x, "sincosd(x) is only defined for finite `x`."))
    elseif isnan(x)
        return ( oftype(x,NaN), oftype(x,NaN) )
    end

    rx = copysign(float(rem(x,360)), x)
    arx = abs(rx)

    if rx == oftype(rx, 0)
        return ( float(0), float(1) )

    elseif arx < oftype(rx, 45)
        r = deg2rad_ext(rx)
        return sincos_kernel(deg2rad_ext(rx))

    elseif arx <= oftype(rx, 135)
        y = deg2rad_ext(oftype(rx,90) - arx)
        s,c = sincos_kernel(y)
        return ( copysign(c, rx), s )

    elseif arx == oftype(rx, 180)
        return ( copysign(zero(rx), rx), float(-1) )

    elseif arx < oftype(rx, 225)
        y = deg2rad_ext( (oftype(rx,180) - arx)*sign(rx) )
        s,c = sincos_kernel(y)
        return ( s, -c )

    elseif arx < oftype(rx, 315)
        y = deg2rad_ext( arx - oftype(rx,270) )
        s,c = sincos_kernel(y)
        return ( -copysign(c, rx), s )

    else
        y = deg2rad_ext( rx - copysign(oftype(rx,360), rx) )
        return sincos_kernel(y)

    end
end

which was almost 10x slower than calling sind and cosd for small angles (< 640).

The current implementation has a very good performance:

julia> @btime sincosd(19.6)
  1.382 ns (0 allocations: 0 bytes)
(0.33545156975025503, 0.9420574527872967)

julia> @btime sincos(19.6*pi/180)
  1.649 ns (0 allocations: 0 bytes)
(0.33545156975025503, 0.9420574527872967)

julia> @btime sincosd(1986.1835)
  1.386 ns (0 allocations: 0 bytes)
(-0.10771305630655274, -0.9941820243301029)

julia> @btime sincos(1986.1835*pi/180)
  12.493 ns (0 allocations: 0 bytes)
(-0.1077130563065501, -0.9941820243301033)

and, by design, gives the same answer as sind and cosd.

[KristofferC]

If this just calls sind and cosd, is having this function at all needed?

[ronisbr]

@KristofferC well, IMHO, it keeps the consistency, since we have sind, cosd, tand, atand, asind, etc.

[simonbyrne]

I'm fine with this, though we probably don't need too many more degree functions.

[KristofferC]

If there are no inputs where sincosd(x) is beneficial over sind(x), cosd(x) then it just seems better to me to save the code and documentation.

[ronisbr]

@KristofferC

I think this is good when you want to compute the sine and cosine of a angle that is computed by a lot of terms, like:

sincosd(theta + omega + alpha/2 + gamma - beta - 180)

Otherwise, or you use deg2rad and sincos or you need to store the angle in an auxiliary variable to call sind and cosd. In this case, I think it is good to have the degree version of sincosd.

@simonbyrne

I use degree functions a lot, and I think this is the only trigonometric function without a “degree” version. That’s why I proposed it, to keep consistency. Btw, I applied all of your suggestions. Thanks for the review.

[bramtayl]

Instead of adding another _d function, is it too late to have

struct Degrees end
const ° = Degrees()
sincos(°, angle)

Or whatever clever unitful way there is to go about this?

[KristofferC]

I think this is good when you want to compute the sine and cosine of a angle that is computed by a lot of terms, like:..

Yeah, I guess that is an ok reason.

[oxinabox]

Instead of adding another _d function, is it too late to have

struct Degrees end
const ° = Degrees()
sincos(°, angle)

Or whatever clever unitful way there is to go about this?

Unitful.jl (and suppliments in UnitfulAngles.jl) already does basically this.
After this PR is merged, it will be trival to add sincos here:
https://github.com/ajkeller34/Unitful.jl/blob/20ba423b2587b6dd9772a7cf7904b4dbc6c19bf9/src/pkgdefaults.jl#L64-L69

[ronisbr]

Hi guys! Should I need to do something more here?

[ronisbr]

Bump! :)

[musm]

This is pretty surprising, any idea why? Part of the point of sincos, is it avoid the double cost of calling sin and cos on the same argument.

[ronisbr]

I think the compiler is clever enough to remove those duplicated executions.

[ronisbr]

In fact, the sincos_kernel use the same approach and it has a comment saying that inlining takes care of duplicated computations:

julia/base/special/trig.jl

Line 208 in 8221c87

# There's no need to write specialized kernels, as inlining takes care of remo-

[JeffBezanson]

which was almost 10x slower than calling sind and cosd for small angles (< 640).

That sounds suspicious --- are we sure the function is type stable etc.?

[StefanKarpinski]

Triage approves the addition of this function for the sake of consistency alone. Even if it's currently no faster than calling (sind, cosd)(x) we can have it and people can use it and then if someone figures out a clever work-saving speedup in the future, callers will get the speedup for free.

[ronisbr]

That sounds suspicious --- are we sure the function is type stable etc.?

Me too. I checked everything I could (given my knowledge of how julia works, including type-stability) and I was not capable to come even close to what I got when calling those functions separated.

But since you also think this can be improved, I will try more and then I submit a new PR, is it OK?

[fredrikekre]

Isn't these conditions checked in sind and cosd? Why check here?

[ronisbr]

Because without it, then the user will see an error on the function sind if Inf is passed. I though it could lead to confusion.

With respect of NaN, it can be removed. Do you want me to remove it? This was added when I was trying to make a function faster than calling both separated.

[fredrikekre]

I don't think we do error checking like that elsewhere though.

[ronisbr]

With respect to "compat annotation" what should I do? Is it also necessary to add a new commit with the information in NEWS.md?

[ronisbr]

Hi guys!

Done :) Sorry for the delay, I am still learning...

After some help (thanks @ararslan and @narendrakpatel), I could rebase everything, add NEWS and COMPAT annotation. Is everything fine now?

[ronisbr]

I have no idea why the build failed in some cases... it does not seem related to this PR.

[ronisbr]

Bump! :)

[ronisbr]

Hi guys

Is it possible to add this to 1.3.0? I saw that alpha was released.

[StefanKarpinski]

Kicked off CI again. Triage says let's merge if green.

[ronisbr]

@StefanKarpinski I have absolutely no idea why it failed on Windows. It does not seem related to this PR. Am I right?

[StefanKarpinski]

Can you rebase the PR and push it again? This seems like it might be unrelated and just happening because of some kind of old cruft on the version of master that this was based on.

[ronisbr]

Hi @StefanKarpinski

Sorry, I just had access to Github now. I saw this was already merged. Should I need to do something?

[JeffBezanson]

No, it's fine.