Yes, the mesh.geometry.x is more accurate when the nodes of the mesh coincide with the degrees of freedom, which is only the case in iso-parametric FEM methods.
As I have already stated, there will always be floating error in tabulation (usually of somewhere around sqrt(floating_precision)), which you cannot avoid. It rarely causes any problems as long as you use double-precision floats.
Yes and no, this depends on your function space.
What happens is that the function is evaluated at the interpolation points of the finite element (which is the dof coordinates for Lagrange, but not for Nedelec and RT), and these are then pushed forward to the physical element for evaluation.
For Nedelec/RT etc then a linear combination and a pull-back routine is used to account for the global orientation vs local orientation.