Hello!

While implementing some helper functions to convert SQLite's datetime values, I've gotten slightly stuck trying to understand a precision issue.

Parsing the datetime string 2022-05-17 13:56:12.569Z using julianday gives me 2459717.08070102. (strftime('%J' ...) gives a similar number).

select julianday('2022-05-17 13:56:12.569Z'); 
-- Output: 2459717.08070102

Plugging this value back into strftime:

select strftime('%Y-%m-%d %H:%M:%fZ', 2459717.08070102) 
-- Output: 2022-05-17 13:56:12.568Z

The second strftime call seems to be correct, so my assumption is that something odd happens when trying to parse the original datetime string using julianday.

Unsurprisingly, storing to text works as expected:

select strftime('%Y-%m-%d %H:%M:%fZ', '2022-05-17 13:56:12.569Z');
-- Output: 2022-05-17 13:56:12.569Z

A runnable example may be found here.

While I'm not depending on using the julianday format, I'm still curious as to why this happens?

Marcus

The usual issue of the limited precision of floating point numbers, I expect.

As an aside, this is exactly the reason that the modified-julianday and the reduced-julianday were "invented" -- to allow more bits to be allocated to guard duty rather than result duty.

Right.

I guess that SQLite is already used in software monitoring and controlling telescopes, satellite orbits, etc. The extra precision from using modified julian dates would be of good use.

How about an extension enabling to easily use MJD? For example, it could have a mjulianday function and override the SQLite date and time functions enabling an additional modifier, mjulianday or so. Someone who has already created several very useful extensions might want to have a look, pse?

Alternatively, the built-in functions could of course be updated in a future version.

That would only gain you about five or six bits of extra precision. My guess is that those who need better than millisecond precision will often want microsecond or even nanosecond precision. Also, for the applications you have in mind, one probably needs to take leap seconds into account, meaning it is better to work in TAI than in UTC. I imagine nanoseconds since some suitable (fairly recent) epoch, stored in a 64 bit integer, is more suited for the purpose.

Well, yes, but things should not become too complicated for a modest SQLite extension. Instead of or in addition to MJD(= JD − 2400000.5) also J2000(= JD - 2451545.0) could be enabled. The precision would be enough for milliseconds and guarantee that contemporary times in string format could be converted forth and back to double precision reals nearly lossless.

Regarding some more "advanced" stuff, the clocks by the various time keeping institutes in the world differ between each other by up to a few microseconds. Microseconds is roughly the accuracy that our timing technology presently has. Standard clocks in computers are typically accurate to a few milliseconds when synchronized with the NTP protocol. Storing absolute time stamps with nanosecond accuracy can of course be done, but presently is not terribly useful for the real world (relative timing and its precision is a different story).

TAI is also not useful for controlling telescopes and satellite orbits etc. It just has a fixed offset from UTC. Rather TT (terrestrial time) or UT1 (former GMT) is needed. Clocks in computers, phones etc are synchronized to UTC (with the NTP protocol), while the calculated positions of celestial objects and satellites are often given for UT1, e.g. (https://astro.ukho.gov.uk/data/aph/APh_2021.pdf). Proper conversion of timestamps between UTC and UT1 involves consulting monthly bulletins from the IERS, International Earth Rotation and Reference Systems Service (https://www.iers.org/IERS/EN/Publications/Bulletins/bulletins.html) to obtain the daily UTC-UT1 differences. I would not use SQLite for this.

I have done an examination of this and determined that there is indeed an error somewhere in the calculation of the julian day number. The error pattern is different from a precision loss error due to insufficient bits.

Here are but a few of the (relatively dense) errors:

sqlite> select epsilon(2440587.50001159);
┌───────────────────────────┐
│ epsilon(2440587.50001159) │
├───────────────────────────┤
│ 4.65661287307739e-10      │
└───────────────────────────┘
VM-steps: 9
Run Time: real 0.004 user 0.000000 sys 0.000000
sqlite> select *, (jdn-sjd) from results limit 10;
┌───────┬───────────────────────────┬──────────────────┬──────────────────┬───────────────────────────┬───────────────────────────┬───────────────────────────┬──────────────────────┐
│  uen  │            odt            │       jdn        │       sjd        │            sd1            │            sd2            │            sd3            │      (jdn-sjd)       │
├───────┼───────────────────────────┼──────────────────┼──────────────────┼───────────────────────────┼───────────────────────────┼───────────────────────────┼──────────────────────┤
│ 1.001 │ 1970-01-01 00:00:01.001 Z │ 2440587.50001159 │ 2440587.50001157 │ 1970-01-01 00:00:01.001 Z │ 1970-01-01 00:00:01.001 Z │ 1970-01-01 00:00:01.000 Z │ 1.16415321826935e-08 │
│ 1.003 │ 1970-01-01 00:00:01.003 Z │ 2440587.50001161 │ 2440587.5000116  │ 1970-01-01 00:00:01.003 Z │ 1970-01-01 00:00:01.003 Z │ 1970-01-01 00:00:01.002 Z │ 1.16415321826935e-08 │
│ 1.005 │ 1970-01-01 00:00:01.005 Z │ 2440587.50001163 │ 2440587.50001162 │ 1970-01-01 00:00:01.005 Z │ 1970-01-01 00:00:01.005 Z │ 1970-01-01 00:00:01.004 Z │ 1.11758708953857e-08 │
│ 1.007 │ 1970-01-01 00:00:01.007 Z │ 2440587.50001166 │ 2440587.50001164 │ 1970-01-01 00:00:01.007 Z │ 1970-01-01 00:00:01.007 Z │ 1970-01-01 00:00:01.006 Z │ 1.16415321826935e-08 │
│ 1.009 │ 1970-01-01 00:00:01.009 Z │ 2440587.50001168 │ 2440587.50001167 │ 1970-01-01 00:00:01.009 Z │ 1970-01-01 00:00:01.009 Z │ 1970-01-01 00:00:01.008 Z │ 1.16415321826935e-08 │
│ 1.011 │ 1970-01-01 00:00:01.011 Z │ 2440587.5000117  │ 2440587.50001169 │ 1970-01-01 00:00:01.011 Z │ 1970-01-01 00:00:01.011 Z │ 1970-01-01 00:00:01.010 Z │ 1.16415321826935e-08 │
│ 1.013 │ 1970-01-01 00:00:01.013 Z │ 2440587.50001172 │ 2440587.50001171 │ 1970-01-01 00:00:01.013 Z │ 1970-01-01 00:00:01.013 Z │ 1970-01-01 00:00:01.012 Z │ 1.16415321826935e-08 │
│ 1.015 │ 1970-01-01 00:00:01.015 Z │ 2440587.50001175 │ 2440587.50001174 │ 1970-01-01 00:00:01.015 Z │ 1970-01-01 00:00:01.015 Z │ 1970-01-01 00:00:01.014 Z │ 1.16415321826935e-08 │
│ 1.017 │ 1970-01-01 00:00:01.017 Z │ 2440587.50001177 │ 2440587.50001176 │ 1970-01-01 00:00:01.017 Z │ 1970-01-01 00:00:01.017 Z │ 1970-01-01 00:00:01.016 Z │ 1.16415321826935e-08 │
│ 1.019 │ 1970-01-01 00:00:01.019 Z │ 2440587.50001179 │ 2440587.50001178 │ 1970-01-01 00:00:01.019 Z │ 1970-01-01 00:00:01.019 Z │ 1970-01-01 00:00:01.018 Z │ 1.11758708953857e-08 │
└───────┴───────────────────────────┴──────────────────┴──────────────────┴───────────────────────────┴───────────────────────────┴───────────────────────────┴──────────────────────┘

uen - unixepoch number
odt - datetime string corresponding to uen
jdn - julianday 2440587.5 + (uen / 86400.0)
sjd - julianday(odt) returned by sqlite3
sd1 - datetime returned by sqlite3 from uen
sd2 - datetime returned by sqlite3 from jdn
sd3 - datetime returned by sqlite3 from sjd

You will note that the error is within (by several decimal orders of magnitude) the epsilon limit of the magnitude of the number which indicates that the error likely arises as a result of improper rounding of intermediates.

Note that the error is 25 ULP.

Found the problem, and it is due to naughty rounding of intermediates (or in this particular case, truncating instead of rounding).

In function computeJD in date.c the following line incorrectly truncates an inexact value. The operation should be rounded in order to get around the fact that trunc(p->s*1000.0) will not give the correct answer if p->s is inexact.

p->s*1000 needs to be "properly rounded" before conversion to an integer.

The easiest way to achieve that is to change:

p->iJD += p->h*3600000 + p->m*60000 + (sqlite3_int64)(p->s*1000);

to

p->iJD += p->h*3600000 + p->m*60000 + (sqlite3_int64)(p->s*1000+.5);

This eliminates all the errors.

Thanks for chasing this down. This improvement was just checked in.

Milliseconds should be preserved over the entire datetime range.

with dt(dt) as
     (
      values ('-4713-11-24 12:00:00.000 Z'),
              ('0000-01-01 00:00:00.000 Z'),
              ('1970-01-01 00:00:00.000 Z'),
              ('9999-12-31 23:59:59.999 Z')
     )
select dt as DateTime,
       epsilon(julianday(dt)) * 86400 as JDPrecisionS,
       epsilon(unixepoch(dt)) as UEPrecisionS
  from dt
;

┌──────────────────────────────┬──────────────────────┬──────────────────────┐
│           DateTime           │     JDPrecisionS     │     UEPrecisionS     │
├──────────────────────────────┼──────────────────────┼──────────────────────┤
│ '-4713-11-24 12:00:00.000 Z' │ 9.59232693276135e-12 │ 3.0517578125e-05     │
│ '0000-01-01 00:00:00.000 Z'  │ 2.01165676116943e-05 │ 7.62939453125e-06    │
│ '1970-01-01 00:00:00.000 Z'  │ 4.02331352233887e-05 │ 1.11022302462516e-16 │
│ '9999-12-31 23:59:59.999 Z'  │ 8.04662704467773e-05 │ 3.0517578125e-05     │
└──────────────────────────────┴──────────────────────┴──────────────────────┘

With this improvement, it appears they are.

Thanks for the fix and the insight into the SQLite internals!