astropy/astropy — pull request #1189 by mhvk (2013-07-09 19:17:38+00:00)
Title: Rounding errors in Time and TimeDelta and internal representation of times

This is partly a placeholder, to be sure we address that currently the arithmetic in `Time` suffers from rounding errors, since there is no carrying between the two parts of the JD. This follows on discussions on astropy-dev [1] and as part of #1082 (also see #366). Copying the specific example of a problem that might easily catch someone by surprise:
```
Time('2012-01-01', scale='utc')+TimeDelta(1.e-6/24./3600., format='jd')
```
Here, one looses the microsecond (but adding `TimeDelta(0., 1.e-6/24/3600, format='jd')` does not).

A complication is that the `sofa_time` internals do not properly carry between the two either, but rather just add to the smaller of the two if what is added is small, or have constructions like `t = ((date1 - DJ00) + date2) / DJC`, which are safe only when `date1-DJ00` can be represented exactly.

Two possible solutions are
1) Use a double-double technique [2,3] to mimic `float128` with two `float64`'s (cannot do `np.float128` directly since it is not supported on all platforms; this was tried for #1082 and #366)
2) Always split the JD by (half-)integer and fractional parts.

The advantage of (1) is better precision (1e-19 s for current time, so good to ns even at the big bang), but care would have to be taken in `sofa` calls. Disadvantage is that double-double only works if the intermediate results are strictly stored as `float64` (*not* the default for 64-bit intel/amd); currently, easy to do, but might break with if, e.g., the python interpreter starts to try to "optimize" math (pypy?).
Advantage of (2) is that `sofa_time` is OK (numbers subtracted are always (half-)integer, so one has no rounding errors), the format is easy to understand (and is in fact what one gets already when entering ISO dates+times), and the span and precision seem "good enough" ( (2**52 day = 1e13 yr; 1day/2**52 ~ 10 ps).

So far, the discussion has been leaning towards (2); an implementation may soon be added here, but in the meantime, it would be good to collect test cases.

[1] https://groups.google.com/forum/#!topic/astropy-dev/JVifTBx8w50
[2] http://en.wikipedia.org/wiki/Quadruple_precision#Double-double_arithmetic
[3] http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf