Why are physical simulations more accurate with floating point numbers "closer to one"? I've often heard it said that any sort of "dimensional" (involving length, time, mass, charge, etc.) calculation should be put in a dimensionless form for two reasons


*

*Getting a value ~ 1e14 when all of your inputs are near 1 is rather obvious (and unlikely) making it easier to pinpoint algorithmic errors.

*Conventional IEEE floating point numbers are denser near zero so dimensionless calculations are more accurate.


I have a pretty intuitive grasp of the first point (if all my lengths ≈ 1 and my masses ≈ 1 and my times ≈ 1 then I can expect all of my forces will reasonably ≈ 1) but I'm having a more difficult time digesting the second. I totally understand why the epsilon between floats depends on your location on the real line, but I don't see how this fact can be exploited for better accuracy. For instance, if you're using the SI speed of light (c = 3e8) and add to it a small epsilon (eps = 1e-8), then c + eps == c and you've "lost" the epsilon to the floating-point void. I don't see how using a dimensionless speed of light (c = 1) diminishes this sort of issue. Why are "dimensionless" simulations said to be more accurate?
 A: So over here what matters isn't as much addition exclusively as it is multiplication/scaling.
Let's consider what it means for floating-point numbers to be "more dense" near zero. The way floating-point works is... imagine the exponent to define which interval of the form [2^{n-1}, 2^n] we are working in - as you can see, as the exponent part grows, your intervals grow larger: [0.5,1], [1,2], [2,4], [4,8] and so on. However, the mantissa - the part of the number that's fixed in size (2^54 for double precision) is, well, fixed. And this number can be viewed as an offset into the interval. (What happens when your number hits the largest possible offset? Well, we reset the offset and increment our exponent by 1. Hence "floating" point, since the binary decimal point floats depending on the magnitude of the number you are trying to represent).Let's look at the interval [1,2]: what's the spacing between successive floats in it? A: (2-1)/2^54 = 2^-54. Now, what about [4096,8192]? A: (8192-4096)/2^54 = 2^-42. And so on. So by "more dense" we literally mean that we can represent smaller differences between numbers in that interval, and so improve accuracy.
Another reason is actually the reason you mentioned - that sufficiently small numbers, those less than "machine epsilon", if added to any floating-point number, will have no effect. Now, machine epsilon for double precision is smaller than you think: it turns out for any binary precision b, your machine epsilon is 2^{-b+1}, so for doubles we have 2^-53 ≈ 2.2204460492503131e-16. Now, consider the large number of physics formulae that have a factor of c^2 ≈ 8.98755179e16 in the denominator. It turns out this becomes ~1.11265e-17: so if whatever else you are multiplying by it is less than ~8, your overall expression is... < DBL_EPSILON and thus eated.
The solution to these issues is to scale intermediate expressions so they indeed land in the densest possible interval for floating-point arithmetic. The easiest way to start is... indeed, to use natural/dimensionless units! It's both physically sound and actually even makes the formulas simpler to write down, and you can always scale back from dimensionless units to dimension-ful ones at the very end if desired. So we do just that. :)
