Where does the floating point error come from? (Finite difference using matrix multiplication versus shifts and adding.)

In Julia it appears that one picks up some error terms when doing finite differences using matrix multiplication versus shifts and addition.

julia> n = 1000
1000

julia> hessM = circshift(eye(n),-1) + circshift(eye(n),1) - 2* eye(n);

julia> function hess(x::Vector)
         return circshift(x,-1) + circshift(x,1) - 2*x
       end
hess (generic function with 1 method)

julia> x = rand(n);

julia> maximum(hessM * x - hess(x))
0.0

julia> x = rand(n)/300;

julia> maximum(hessM * x - hess(x))
8.673617379884035e-19

I understand that this is floating point error, since if I do x = rand(n) / 256 (or any power of 2 for that matter) the error goes away, but if we divide by something as simple as 3 the error crops up.

But where exactly is the floating point error coming in when running the above code? And why is it sensitive to the division?

---

Edit: in fact, I have localized the difference to

julia> x = rand(n)/3;

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
1.1102230246251565e-16

julia> x = rand(n);

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
0.0

julia> x = rand(n) * 3;

julia> maximum(hessM * x + 2 * eye(n) * x - circshift(eye(n),-1) * x - circshift(eye(n),1)*x)
8.881784197001252e-16

What gives?