sin and cos inaccurate at arguments with small imaginary components

[bsxfan]

Julia version: 0.2.0

Examples of problem:

• julia> imag(sin(complex(1,1e-20))) gives
  0.0
• julia> imag(cos(complex(1,1e-20))) gives
  0.0

Reason: complex.jl does as a first step of its sin(z) and cos(z) implementations: u = exp(imag(z)). This becomes exactly u==1.0 for small imag(z), while the imaginary part of z is not referenced again in the rest of the implementation.

Comparison: Python and MATLAB give non-zero values.

Why this is important:
The 'complex-step differentiation' trick is a powerful and very convenient way of implementing 'forward-mode automatic differentiation'. It can be used to get extremely accurate derivatives of complicated functions, with no programmer effort. For some function f(x), we can get a very accurate derivative at x by doing something like:

• 1e20*imag(f(complex(x,1e-20)))
  A Taylor series expansion will show why. Currectly this trick fails with Julia's complex sin and cos implementations. (It does work for log and exp).

[jiahao]

For now, perhaps you could use the complex sine and cosine in libgsl, currently wrapped by GSL.jl:

julia> using GSL
julia> z = sf_complex_sin_e(1.0, 1e-20)

(gsl_sf_result(0.8414709848078965,3.736881847550928e-16),gsl_sf_result(5.4030230586813976e-21,2.3994242409314903e-36))

julia> z[2].val #imaginary part
5.4030230586813976e-21

julia> z[2].err #error in imaginary part
2.3994242409314903e-36

[ViralBShah]

The easier solution would be to use the complex versions of these functions in openlibm. I am not sure if I have imported these yet from FreeBSD's msun, but that is easy to do.

[jiahao]

Similar problem with tan. Will submit fix.

[bsxfan]

Hi @jiahao. I agree that the old tan had the same problem, but I think your new one is also broken. Try for example:

julia> tan(complex(2))
-0.32254371161915346 + 0.0im

julia> tan(2)
-2.185039863261519

julia> sin(complex(2))/cos(complex(2))
-2.185039863261519 + 0.0im

[JeffBezanson]

Seems like sinh, cosh, and tanh also have the problem.

[bsxfan]

Here is a demo of a tool that can help to check the accuracy of complex step differentiation (and therefore also of functions of complex arguments). I spent the day implementing a new Number type, namely a 'dual number' (see Wikipedia). It is very similar to a complex number, except that the 2nd component is a multiple of 'du', the infinitessimally small 'differential unit', having the property du*du=0.

julia> g(x)=-2sin(x/2)_cos(x^2)+log(2^x)-exp(x)
julia> x=3
3
julia> d=x+du (d is a dual number and du is the new differential unit)
3 | 1 (I show() the 2 components separated by | )
julia> du_du
0 | 0
julia> d^2
9 | 6 (this is x^2 and the derivative 2x)
julia> c=complex(d) (convenience function to facilitate complex step differentiation)
3.0 + 1.0e-20im
julia> g(x) (ordinary function value)
-16.188399644761425
julia> g(d) (with dual number differentiation)
-16.188399644761425 | -14.394905462561129
julia> g(c) (with complex step differentiation)
-16.188399644761425 - 1.4394905462561128e-19im
julia> dual(g(c)) (dual() conveniently converts back to dual for display)
-16.188399644761425 | -14.394905462561129

So for the new sin() and cos() and these simple arithmetic operations, everything works and agrees.

[JeffBezanson]

Thanks @bsxfan . All the complex trig functions should be ok now --- could you apply your test to confirm?

[kmsquire]

@bsxfan, if you haven't already, please take a look at and consider contributing to https://github.com/scidom/AutoDiff.jl.

cc: @Scidom, @tshort, @mlubin

[papamarkou]

I wrote this script

https://github.com/scidom/AutoDiff.jl/blob/master/src/dual.jl

with the Dual type for representing dual numbers. I also wrote a test in the test folder of the prospective AutoDiff package to test forward automatic differentiation on a univariate function.

The dual.jl file has been written irrespectively of AD, abiding to the paradigm of base/complex.jl. It may be good to ship it with Julia as a predefined type in order to represent dual numbers, along with complex numbers.