inaccuracy: Complex ^ Int

[bsxfan]

@aswart found another complex function that causes problems for complex-step differentiation (see also #2889).

When raising a complex number with negative real part and small imaginary part to an integer power (an everyday occurrence with complex step diferentiation), the imaginary part is wrong.

Example:

This gives an accurate result :

julia> imag(complex(-2,1e-20)|x->x*x)
-4.0e-20

But, when using ^(z::Complex, n::Integer) = z^complex(n)

julia> imag(complex(-2,1e-20)|x->x^2)
-9.797174393178826e-16 

then we get something entirely different. (MATLAB and Python give -4.0e-20).

One way to fix this would be to simply defer this case to Base.power_by_squaring(A,p), but the more general case ^(Complex,Complex) would then probably still need some attention.

julia> imag(complex(-2,1e-20)|x->Base.power_by_squaring(x,2))
-4.0e-20

CC: @jiahao

[JeffBezanson]

Probably limited by the precision of pi in exp(p*log(z)). Using power_by_squaring is probably the best solution.

[aswart]

It looks like the difference is the result of the accuracy of sin(x) for values of x close to n*pi. For small x, sin(x) \approx x. The sine function seems to handle this case, but not for sin(pi - x), which is also \approx x for small x.

Compare

julia> sin(1e-20)
1.0e-20
julia> sin(pi - 1e-20)
1.2246467991473532e-16

For negative real value and small complex component, atan2 returns an angle close to pi, resulting in the precision loss at complex.jl#L253.

[bsxfan]

julia> sin(0)==sin(pi)
false

[aswart]

Ah, yes the precision loss is due to the fact that atan2(1e-20,-1) == pi, not because of the sine function.

[bsxfan]

Sorry, I'm going off on a tangent, but has the following been reported somewhere?:

julia> sin(1e6pi)
2.785584461607201e94

[simonbyrne]

How about this as a crazy solution:

import  Base.convert, Base.show, Base.sin, Base.cos

# define a type that is equivalent to angle + quad * pi/2
type AngleQuad
    angle::Float64
    quad::Int
end
convert(::Type{Float64},a::AngleQuad) = a.angle + a.quad*pi*0.5
show(io::IO,a::AngleQuad) = print(io,a.angle," + ",a.quad,"π/2")

# atan2 for AngleQuads
function atan2q(y,x)
    if abs(x) > abs(y) 
        AngleQuad(atan2(copysign(y,x),abs(x)), 2*signbit(x))
    else
        AngleQuad(atan2(-copysign(x,y),abs(y)), 2*signbit(y)-1)
    end
end

# more accurate sin/cos
function sin(a::AngleQuad)
    r = mod(a.quad,4)
    if r == 0
        return sin(a.angle)
    elseif r == 1
        return cos(a.angle)
    elseif r == 2
        return -sin(a.angle)
    else
        return -cos(a.angle)
    end
end

function cos(a::AngleQuad)
    r = mod(a.quad,4)
    if r == 0
        return cos(a.angle)
    elseif r == 1
        return -sin(a.angle)
    elseif r == 2
        return -cos(a.angle)
    else
        return sin(a.angle)
    end
end

Then we have:

julia> x = atan2q(1e-20,-2)
-5.0e-21 + 2π/2

julia> sin(x)
5.0e-21

[bsxfan]

Could we re-open this issue, because 9470bb2 has introduced some new problems, for example:

julia> complex(1//2,1//3)^2
ERROR: no method power_by_squaring(Complex{Float64},Float64)
 in ^ at complex.jl:349

julia> complex(2.0,2.0)^(-2)
ERROR: DomainError()
 in power_by_squaring at intfuncs.jl:84

In the previous build, these examples worked.

Then, there are also some type stability issues:

julia> complex(2,2)^2
 0.0 + 8.0im

Here one would probably want Complex{Int} as return type, analogously to 4==2^2.

As far as I understand the new complex ^, there are three available solutions now:

1. Power by squaring: works for non-negative, real integer exponents. This can be extended for use with negative exponents, by pre or post inversion, but again for type stability, one may not want to allow Complex{Integer} ^ Int for negative powers. Inversion of Complex{FloatingPoint} and Complex{Rational} do preserve their types, so they can be raised to negative integer powers via inversion.
2. log(p*exp(z))
3. an algorithm involving atan2, sin and cos.

It seems some care is needed here to decide which argument types to send to which algorithm, but I think it makes sense just to immediately send Complex^Integer to the power_by squaring solution. I tried the following:

^{T<:FloatingPoint}(z::Complex{T}, n::Bool) = n ? z : one(z)  # to suppress ambiguity warning
^{T<:Rational}(z::Complex{T}, n::Bool) = n ? z : one(z)           # to suppress ambiguity warning
^{T<:Integer}(z::Complex{T}, n::Bool) = n ? z : one(z)             # to suppress ambiguity warning

^{T<:FloatingPoint}(z::Complex{T}, n::Integer) = n>=0?power_by_squaring(z,n):power_by_squaring(inv(z),-n)
^{T<:Rational}(z::Complex{T}, n::Integer) = n>=0?power_by_squaring(z,n):power_by_squaring(inv(z),-n)
^{T<:Integer}(z::Complex{T}, n::Integer) = power_by_squaring(z,n) # DomainError for n<0

This fixes the three examples above:

julia> complex(1//2,1//3)^2
5//36 + 1//3im

julia> complex(2.0,2.0)^(-2)
0.0 - 0.125im

julia> complex(2,2)^2
0 + 8im

[simonbyrne]

Also, I noticed

julia> complex(-2.0,1e-20)^complex(2.0,1e-20)
4.0 - 9.797174393178826e-16im

We might be able to avoid that with some atan2 trickery.

[JeffBezanson]

Matlab gives the same answer we do for complex(-2.0,1e-20)^complex(2.0,1e-20), which seems to indicate it uses the same integer special case trick.

[bsxfan]

Yes, that is my guess too. I have observed that in both languages we also get the same bad result when using the exp(p*log(z)) algorithm:

imag(exp(2*log(complex(-2,1e-20)))) == -9.797e-16

[bsxfan]

So if we adopt @simonbyrne's clever AngleQuad solution, maybe Julia can be more accurate than MATLAB for such cases ...

[JeffBezanson]

It is cool, but we have to consider how it integrates with the rest of the system. It's essentially a full new numeric type.