Note that orbiter's current physics engine already takes some measures to improve numerical precision over large dynamic distance ranges compared to a naive floating point implementation.
The main problem is the loss of accuracy when adding floating point numbers of widely different magnitude. The way this is implemented in the processor is
1. step one: equalize the exponents of the operands by shifting the digits in the mantissa.
2. step two: add the numbers by adding their mantissae.
Step one is the one that can lead to a loss of precision, since it may shift any number of significant digits out of the mantissa. In the extreme case, if the exponents differ by more than the number of digits stored in the mantissa, you will shift out all digits and end up with zero. You can find out this precision limit by finding an eps > 0 such that
1+eps = 1
Orbiter tries to overcome the worst of the loss of precision by storing position variables in two separate numbers:
- a base position
- an offset position
At each time step, the position change is added to the offset part of the position variable. At regular time intervals, the offset is flushed into the base, and the offset is reset to zero. The update intervals are chosen such that there is a compromise between the loss of precision when adding a single step to the offset, and the loss of precision when adding the offset to the base. Effectively, this scheme simulates a floating point representation with a longer mantissa.
If higher precisions are required, this scheme could be augmented by adding more stages. This would make calculations slightly less efficient, since each time you need the actual position of an object, you first have to add the contributions from all stages.
A better solution would probably be a reference frame that is moving with the observer, rather than one that is fixed at the solar system barycentre, as suggested above. This would however make a lot of computations more awkward than they are now.