Cos(pi/2) does not return 0

General Usage
1

Hi!
While trying out some stuff in Julia, I came across the fact that cos(pi/2) does not give exact 0, but some type dependent output (6e-17 for Float64, -4e-8 for Float32, 5e-4 for Float16). Is this intended? And what is the best way to avoid this?

2

Not intended, but hard to avoid because of floating point error. See

You can use

julia> cospi(0.5)
0.0
7 Likes
3

If what you are wondering is why Julia does not know the cosine involving the Irrational object π exactly, this is because it is converted to a float when performing the operation π/2. Base Julia does not implement AffineFunctionOfIrrational or anything like that. Someone might have done so in a package.

3 Likes
4

Something similar to UniformScaling for irrationals could maybe help fix that issue.

4 Likes
5

I see, thanks so much for your responses!

I’m stuck with one more question though, but that might be because I’m starting out with Julia:
I’ve had a look in the source code at /opt/julia/share/julia/base/special/trig.jl 

function cos(x::T) where T<:Union{Float32, Float64}
    absx = abs(x)
    if absx < T(pi)/4
        if absx < sqrt(eps(T)/T(2.0))
            return T(1.0)
        end
        return cos_kernel(x)
    elseif isnan(x)
        return T(NaN)
    elseif isinf(x)
        cos_domain_error(x)
    else
        n, y = rem_pio2_kernel(x)
        n = n&3
        if n == 0
            return cos_kernel(y)
        elseif n == 1
            return -sin_kernel(y)
        elseif n == 2
            return -cos_kernel(y)
        else
            return sin_kernel(y)
        end
    end
end

that is apparently using approximations for cos(x), -sin(x), -cos(x) or sin(x) depending on the interval x is in.
(I must admit I haven’t had a closer look at the rem_pio2_kernel function, but I admit that it returns the div and rem value from the division by π/2)

I defined my own function based on the same idea of evaluating the trigonometric functions only in the interval [-π/4, π/4] around either their extremum (cosine) or root (sine) as suggested by cos_kernel

function mycos(x)
  if x < 0
    return mycos(-x)
  else
    j = div(x+pi/4,pi/2)%4
    newx = x%2pi - j*pi/2
    if j == 0 return cos(newx)
    elseif j == 1 return -sin(newx)
    elseif j == 2 return -cos(newx)
    elseif j == 3 return sin(newx)
    end
    #if that should (in whatever case) not work, return normal value
    return cos(x)
  end
end

This does return the values at multiple of π/2 correctly (not talking about performance here). I don’t see too may differences to the function in trig.jl and wonder why that one is not capable of doing that…

6

Probably because it’s not in fact correct (e.g. the cosine of the closest Float64 representation of pi/2 is not zero) and you get worse results everywhere. For example

julia> x = 1.571
1.571

julia> y1 = mycos(x)
-0.00020367320369523912

julia> y2 = cos(x)
-0.00020367320369517786

julia> y3 = cos(BigFloat(x))
-2.036732036951778709106180717800557926547406375863370495134343067126837115362263e-04

julia> abs(y1 - y3)
6.124442755166129148251053819710905768303105297194328731628846377367955634449768e-17

julia> abs(y2 - y3)
1.299519376970903891054505204264901252558765305671268371153622632044365550231594e-20

If you are specifically concerned about the cosine of multiples of pi, use the cospi function. It’s there for exactly that purpose.

2 Likes
7

Alright! Thanks for your answer.