Indicate problem with MJD initialization - #9328

Why isn't day_frac, or anything in astropy/time/utils.py in fact, in the online documents? Or at least, where is it? Neither I nor Google are having much luck finding it anywhere outside the source code, in spite of there being reasonable docstrings.

Thanks for working on this! If we decide to go ahead with #9017, it would be good to add the unit test as an actual hypothesis test with explicit @example for the confirmed problematic cases.

Let's ask @taldcroft and @mhvk about the lack of docs - whether it's deliberate if the functions are not meant to be public, or an oversight in the API docs?

Deliberate: these are utilities for internal use which were thought unlikely to be generally useful. Though I readily admit I have wanted to use it myself (and ended up reimplementing day_frac as I found a small buglet - #8763).
While I can see the point of making the utilities public, it has the quite big disadvantage that we're getting stuck with their behaviour. My sense is that it may be preferable to keep day_frac as an implementation detail (for this PR, that probably means one should test only Time itself); if not, then arguably it should be moved out of astropy.time altogether.

Well, there is a test that uses Time only, but in the interest of narrowing things down I used Decimals to show that two_sum does not appear to do what it claims to. That is, the sum coming out differs from the sum going in by considerably more than one part in 10**32. In fact it's probably worth writing a whole wodge of Time precision tests using Decimals as a check, though I think hypothesis is the way to make them reasonably thorough.

Really I think we're better off using something like libqd (a compiled double-double and quad-double library with a reasonable license) to implement Time; with a little cleverness double-doubles and quad-doubles could be made available as numpy dtypes when the appropriate module is imported. This is partially done here: https://github.com/spacetelescope/spherical_geometry

@aarchiba - thanks very much for uncovering these problems. I'll try to have a look. I definitely agree with you that we should move towards just using quad-precision, especially once numpy has proper support for (real) float128.

@aarchiba - thanks very much for uncovering these problems. I'll try to have a look. I definitely agree with you that we should move towards just using quad-precision, especially once numpy has proper support for (real) float128.

I've been thinking about that (and occasionally chipping in to the sporadic discussions). Numpy actually already has real quad precision on platforms where longdouble is 128-bit floats - this includes arm64 (which can run on Raspberry Pi 3 and 4, albeit with an unusual OS! also the pine64 single-board computers and maybe your telephone too) and PowerPC (we have some monster GPU PPC server here). These are probably the platforms that have hardware support for quad precision, though how fast exactly that hardware support is may vary.

If you don't have hardware quad precision, a robust compiled double-double type, especially if you get quad-double for free, is probably preferable for most applications: it will run faster and carry almost as many digits of accuracy (with loads more in the quad-double for a backstop check). And this can be done by importing a portable compiled C library with a reasonable license. So I think that's the way to go. The C library could ultimately become part of numpy, if integration there were wanted.

Really I think we're better off using something like libqd (a compiled double-double and quad-double library with a reasonable license) to implement Time.

Here when you say "... to implement Time", I'm not sure quite sure of the scope. We can't lose sight of the fact that all the (hard) computations for Time are done with ERFA, so I think we are fundamentally bound to stick with their implementation of time represented as two 64-bit floats jd1 and jd2.

Of course taking advantage of more precision in the astropy interface to ERFA might reduce some of the pain and bugs.

@taldcroft makes a good point - we are stuck with erfa/sofa, which works with two doubles, and their implementations mean that it is by far the safest to have the larger double be an exact (fraction of a) day.
Of course, none of this distracts from the fact that we need to ensure our code works well...

OK, found the problem (why can I not simply copy a prescription?!?!?). Will push to this branch.

p.s. I'm now totally swayed about testing two_sum and day_frac directly. After all, if we ever change implementation, we can remove the tests.

@bsipocz - I'm wondering how to schedule this - is there still a 2.0 buxfix release, or are we done?

@mhvk - I think we're done. I've push 2.0.15 to pypi yesterday and not planning on adding any more to that series (with the exception of serious security issues).

I haven't got to finishing 3.2.2 as there are still some CI failures, so feel free to milestone stuff that is expected to be merged very soon (~today/tomorrow) to that one.

@bsipocz - OK, I milestoned for 3.2.2. Though this really is quite bad, and perhaps worth another 2.0 release.

OK, I've added a label for it. However the timelines just won't work out. I'm definitely a 👎 on pushing another release today before getting 3.2.2 released, and if this is being included in 3.2.2 then the whole timeline will be messed up to be also included in a last after last 2.0.x. 🤷‍♀

Well, if it matters, PINT will have to include its own implementation of two_sum and day_frac and just hope that its Time objects aren't broken otherwise unless this makes it into a python2-usable form. We just can't give up on python2 until key related software (e.g. the entire Fermi toolchain, PRESTO, PSRCHIVE) is moved out of the Stone Age.

FWIW: I can cherry-pick at least the solution (218ad7f) cleanly on 2.0.x

@mhvk - could you please approve and merge? I would like to get on with the 3.2.2 release.

Based on discussion between the release managers, the decision that there will be one more LTS release to include this (and other strong cases if there will be any). We'll likely make it alongside 4.0rc1 in 1-1.5 months time.