# 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