Specialized math functions for use with Duals?

[dbeach24]

I was examining some trig code with Dual values using @code_llvm and @code_native, and noticed that the generated code included redundant calls to sin(x) and cos(x); this derivative could be coded more optimally.

In my case, the function evaluates both sin and cos of several input variables. However, I can take advantage of this fact when using Duals, by creating a specialized sincos function:

sincos(x) = (sin(x), cos(x))

function sincos(x::Dual)
    sx = sin(value(x))
    cx = cos(value(x))
    return (
        Dual(sx, cx * partials(x)),
        Dual(cx, -sx * partials(x)),
    )
end

I understand ForwardDiff already implements a similar optimization in the case of exp, for example. I'm not really pushing this as a must-have addition (though not opposed if others think it's good). What I'm really curious to know is if this sort of specialized math function + optimization has other similar use cases and whether it makes sense for other users.

Footnote: I didn't notice a significant runtime improvement from implementing this is my code, so, in my case the idea may be a dud. (Or I may have too much overhead from other sources.)

[jrevels]

More specialized methods for Duals are always welcome!

It would be nice to have a sincos method in Base that leveraged FSINCOS - then we could overload that for Dual numbers. It's possible, though, that LLVM already optimizes away the extra function call anyway. EDIT: Just read your comment where you said that that optimization wasn't happening.

[mlubin]

See also JuliaLang/julia#10442

[dbeach24]

Agree that it would be better if it could be supplied as an overload to a function in Base.

I remember looking into that a little bit. As I understand it, FSINCOS is an 387 instruction, but is slower than separately computing sin and cos using SSE instructions.

[dbeach24]

Based on what I saw in JuliaLang/julia#10442, I feel I must update my draft implementation:

@inline sincos(x) = (sin(x), cos(x))

@inline function sincos(x::Dual)
    sx, cx = sincos(value(x))
    return (
        Dual(sx, cx * partials(x)),
        Dual(cx, -sx * partials(x)),
    )
end

Now, if sincos is ever implemented in Base, presumably we would just specialize that implementation for Dual numbers. This should also recurse better in the case of nested Duals.