Cos(pi/2) does not return 0

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?

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

You can use

julia> cospi(0.5)
0.0
6 Likes

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

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

4 Likes

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…

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

Alright! Thanks for your answer.