IEEE 754 Double Precision uses (effectively) 53 bits of mantissa, able to reflect relative changes down to 2^(-53) . If you are using values in the approximate range of 1.0, then 2^(-53) of that is 1.1102230246251565404236316680908203125 * 10^(-16) which is not quite as good as 10^(-16) . In the area of 1.0, a single bit difference translates to about 1E-16 change in output.
It is a truism in floating point calculation that the same calculation performed two different ways that are "algebraically equivalent" can give you different results due to differences in rounding. http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F
Also, I seem to recall that some of the hardware instructions are permitted to have an accuracy less than 1 ULP (Unit In The Last Place); the details of the standard would have to be checked for that kind of detail.
If you need higher precision you will need to use a higher precision package such as the Symbolic Toolbox
