bug#26632: TAI<->UTC conversion botches 1961 to 1971

[Mark H Weaver]

Hi Zefram,

Zefram <address@hidden> writes:
> The SRFI-19 library gets TAI<->UTC conversions badly wrong in the years
> 1961 to 1971 (inclusive).

Indeed.

> This has to be examined somewhat indirectly, because SRFI-19 doesn't offer
> any way to display a TAI time in its conventional form as a date-like
> structure, nor to input a TAI time from such a structure.

FWIW, here are two procedures I hacked up to support TAI dates:

--8<---------------cut here---------------start------------->8---
(define (time-tai->date-tai t . tz-offset)
  (apply time-tai->date
         (time-utc->time-tai! (make-time time-utc
                                         (time-nanosecond t)
                                         (time-second t)))
         tz-offset))

(define (date-tai->time-tai d)
  (let ((t (time-tai->time-utc (date->time-tai d))))
    (set-time-type! t time-tai)
    t))
--8<---------------cut here---------------end--------------->8---

> SRFI-19's date structure, as implemented, is always interpreted
> according to UTC.

Indeed, I discovered this as well, and found it surprising.  The text of
SRFI-19 fails to mention it, but the reference implementation makes it
quite clear.

> First I'll consider an ordinary day in 1967:
>
> scheme@(guile-user)> (use-modules (srfi srfi-19))
> scheme@(guile-user)> (time-difference
> ... (time-utc->time-tai (date->time-utc (make-date 0 0 0 0 15 3 1967 0)))
> ... (time-utc->time-tai (date->time-utc (make-date 0 0 0 0 14 3 1967 0))))
> $1 = #<time type: time-duration nanosecond: 0 second: 86400>
>
> This takes the start and end of 1967-03-14, as judged by UTC, converts
> both of these times to TAI, and asks for the duration of that TAI
> interval.  It's asking how many TAI seconds long that UTC day was.
> As described in <http://maia.usno.navy.mil/ser7/tai-utc.dat>, there
> was no UTC leap on that day, but throughout 1967 UTC had a frequency
> offset from TAI such that each UTC second lasted exactly 1.00000003 TAI
> seconds.  The correct answer to the above question is therefore exactly
> 86400.002592 s.  The answer shown above, of 86400.000000 s, is incorrect.
>
> If time-tai->time-utc is applied to the times in the above example,
> it accurately inverts what time-utc->time-tai did.  It is good that the
> conversions are mutually consistent, but in this case it means they are
> both wrong.
>
> Second, I'll consider a less ordinary day:
>
> scheme@(guile-user)> (time-difference
> ... (time-utc->time-tai (date->time-utc (make-date 0 0 0 12 1 2 1968 0)))
> ... (time-utc->time-tai (date->time-utc (make-date 0 0 0 12 31 1 1968 0))))
> $2 = #<time type: time-duration nanosecond: 0 second: 86400>
>
> This time the period considered is from noon 1968-01-31 to noon
> 1968-02-01.  The same frequency offset described above applies throughout
> this period.  The additional complication here is that at the end of
> 1968-01-31 there was a leap of -0.1 (TAI) seconds.  The true duration of
> this day is therefore exactly 86399.902592 s.  The answer shown above,
> of 86400.000000 s, is incorrect in two ways, accounting for neither the
> frequency offset nor the leap.

I've attached two patches to fix this bug.  The first lays the
groundwork by adding support for non-integer TAI-UTC deltas.  The second
patch adds the TAI-UTC tables for 1961-1971 and uses them to implement
TAI<->UTC conversions over that time range with nanosecond accuracy.

Although the code is now written, I'm unsure whether we should add it to
Guile.  I'm vaguely concerned about violating widely-held assumptions,
e.g. that UTC runs at the same rate as TAI (except when leap seconds are
introduced), which might cause some code on top of Guile to misbehave if
the system clock is set pre-1972, although admittedly such a scenario
seems unlikely.

I'm curious to hear opinions on this.

Anyway, here are the patches.

      Mark


>From 3a67fbfd441b39630ff3c3201d2a731b51b1a8ee Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 23 Oct 2018 00:39:30 -0400
Subject: [PATCH 1/2] SRFI-19: TAI<->UTC conversions support non-integer
 deltas.

Previously, 'utc->tai' and 'tai->utc' were unary procedures and returned
a second value, mapping seconds to seconds, based on the assumption that
the nanoseconds field would never be changed by these maps.  To support
non-integer values of TAI-UTC, here we change 'utc->tai' and 'tai->utc'
to accept two arguments and return two values: nanoseconds and seconds.

* module/srfi/srfi-19 (utc->tai, tai->utc): Add an additional 'ns'
argument, and return it as an additional value.
(current-time-tai, priv:time-tai->time-utc!, priv:time-utc->time-tai!)
(time-tai->julian-day, time-monotonic->julian-day): Adapt accordingly.
---
 module/srfi/srfi-19.scm | 56 +++++++++++++++++++++++------------------
 1 file changed, 31 insertions(+), 25 deletions(-)

diff --git a/module/srfi/srfi-19.scm b/module/srfi/srfi-19.scm
index 46de91a7e..f653b1bc3 100644
--- a/module/srfi/srfi-19.scm
+++ b/module/srfi/srfi-19.scm
@@ -237,21 +237,21 @@
   (set! leap-second-table (read-tai-utc-data filename)))
 
 
-(define (utc->tai utc-seconds)
+(define (utc->tai ns utc-seconds)
   (let loop ((table leap-second-table))
     (cond ((null? table)
-           utc-seconds)
+           (values ns utc-seconds))
           ((>= utc-seconds (caar table))
-           (+ utc-seconds (cdar table)))
+           (values ns (+ utc-seconds (cdar table))))
           (else
            (loop (cdr table))))))
 
-(define (tai->utc tai-seconds)
+(define (tai->utc ns tai-seconds)
   (let loop ((table leap-second-table))
     (cond ((null? table)
-           tai-seconds)
+           (values ns tai-seconds))
           ((>= tai-seconds (+ (caar table) (cdar table)))
-           (- tai-seconds (cdar table)))
+           (values ns (- tai-seconds (cdar table))))
           (else
            (loop (cdr table))))))
 
@@ -309,9 +309,8 @@
   (let* ((tod (gettimeofday))
          (sec (car tod))
          (usec (cdr tod)))
-    (make-time time-tai
-               (* usec 1000)
-               (utc->tai sec))))
+    (receive (ns s) (utc->tai (* usec 1000) sec)
+      (make-time time-tai ns s))))
 
 ;;(define (current-time-ms-time time-type proc)
 ;;  (let ((current-ms (proc)))
@@ -460,10 +459,12 @@
 (define (priv:time-tai->time-utc! time-in time-out caller)
   (if (not (eq? (time-type time-in) time-tai))
       (time-error caller 'incompatible-time-types time-in))
-  (set-time-type! time-out time-utc)
-  (set-time-nanosecond! time-out (time-nanosecond time-in))
-  (set-time-second!     time-out (tai->utc (time-second time-in)))
-  time-out)
+  (receive (ns s) (tai->utc (time-nanosecond time-in)
+                            (time-second time-in))
+    (set-time-type! time-out time-utc)
+    (set-time-nanosecond! time-out ns)
+    (set-time-second!     time-out s)
+    time-out))
 
 (define (time-tai->time-utc time-in)
   (priv:time-tai->time-utc! time-in (make-time-unnormalized #f #f #f) 
'time-tai->time-utc))
@@ -475,10 +476,13 @@
 (define (priv:time-utc->time-tai! time-in time-out caller)
   (if (not (eq? (time-type time-in) time-utc))
       (time-error caller 'incompatible-time-types time-in))
-  (set-time-type! time-out time-tai)
-  (set-time-nanosecond! time-out (time-nanosecond time-in))
-  (set-time-second!     time-out (utc->tai (time-second time-in)))
-  time-out)
+  (receive (ns s)
+      (utc->tai (time-nanosecond time-in)
+                (time-second time-in))
+    (set-time-type! time-out time-tai)
+    (set-time-nanosecond! time-out ns)
+    (set-time-second!     time-out s)
+    time-out))
 
 (define (time-utc->time-tai time-in)
   (priv:time-utc->time-tai! time-in (make-time-unnormalized #f #f #f) 
'time-utc->time-tai))
@@ -807,10 +811,11 @@
 (define (time-tai->julian-day time)
   (if (not (eq? (time-type time) time-tai))
       (time-error 'time-tai->julian-day 'incompatible-time-types  time))
-  (+ (/ (+ (tai->utc (time-second time))
-           (/ (time-nanosecond time) nano))
-        sid)
-     tai-epoch-in-jd))
+  (receive (ns s) (tai->utc (time-nanosecond time)
+                            (time-second time))
+    (+ (/ (+ s (/ ns nano))
+          sid)
+       tai-epoch-in-jd)))
 
 (define (time-tai->modified-julian-day time)
   (- (time-tai->julian-day time)
@@ -820,10 +825,11 @@
 (define (time-monotonic->julian-day time)
   (if (not (eq? (time-type time) time-monotonic))
       (time-error 'time-monotonic->julian-day 'incompatible-time-types  time))
-  (+ (/ (+ (tai->utc (time-second time))
-           (/ (time-nanosecond time) nano))
-        sid)
-     tai-epoch-in-jd))
+  (receive (ns s) (tai->utc (time-nanosecond time)
+                            (time-second time))
+    (+ (/ (+ s (/ ns nano))
+          sid)
+       tai-epoch-in-jd)))
 
 (define (time-monotonic->modified-julian-day time)
   (- (time-monotonic->julian-day time)
-- 
2.19.1

>From adfcb5ee624d457a8dae33e8e122cbff6d4a1428 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 23 Oct 2018 00:42:12 -0400
Subject: [PATCH 2/2] SRFI-19: Implement correct TAI<->UTC conversions for
 1961-1971.

* module/srfi/srfi-19.scm (tai-utc-1961-1971-table): New variable.
(pre-1972-utc->tai, pre-1972-tai->utc): New procedures.
(utc->tai, tai->utc): Use the new procedures.
---
 module/srfi/srfi-19.scm | 58 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 56 insertions(+), 2 deletions(-)

diff --git a/module/srfi/srfi-19.scm b/module/srfi/srfi-19.scm
index f653b1bc3..d4afa1931 100644
--- a/module/srfi/srfi-19.scm
+++ b/module/srfi/srfi-19.scm
@@ -240,7 +240,7 @@
 (define (utc->tai ns utc-seconds)
   (let loop ((table leap-second-table))
     (cond ((null? table)
-           (values ns utc-seconds))
+           (pre-1972-utc->tai ns utc-seconds))
           ((>= utc-seconds (caar table))
            (values ns (+ utc-seconds (cdar table))))
           (else
@@ -249,13 +249,67 @@
 (define (tai->utc ns tai-seconds)
   (let loop ((table leap-second-table))
     (cond ((null? table)
-           (values ns tai-seconds))
+           (pre-1972-tai->utc ns tai-seconds))
           ((>= tai-seconds (+ (caar table) (cdar table)))
            (values ns (- tai-seconds (cdar table))))
           (else
            (loop (cdr table))))))
 
 
+(define tai-utc-1961-1971-table
+  ;;      UTC start       address@hidden  slope
+  ;;        (ns)             (ns)      (ns/s)
+  '(( -60480000000000000  6185682000    30)
+    (-126230400000000000  4313170000    30)
+    (-136771200000000000  4155058000    15)
+    (-142128000000000000  3974706000    15)
+    (-152668800000000000  3716594000    15)
+    (-157766400000000000  3540130000    15)
+    (-168307200000000000  3282018000    15)
+    (-181526400000000000  2983730000    15)
+    (-189388800000000000  2765794000    15)
+    (-194659200000000000  2697278800    13)
+    (-252460800000000000  1845858000    13)
+    (-265680000000000000  1647570000    15)
+    (-283996800000000000  1422818000    15)))
+
+(define (pre-1972-utc->tai ns utc-seconds)
+  (let ((utc-ns (+ ns (* nano utc-seconds))))
+    (let loop ((table tai-utc-1961-1971-table))
+      (cond ((null? table)
+             (values ns utc-seconds))
+            ((>= utc-ns (caar table))
+             (apply (lambda (utc-start address@hidden slope)
+                      (let ((tai-ns (+ utc-ns
+                                       address@hidden
+                                       (round-quotient (* slope (- utc-ns
+                                                                   utc-start))
+                                                       nano))))
+                        (values (truncate-remainder tai-ns nano)
+                                (truncate-quotient tai-ns nano))))
+                    (car table)))
+            (else
+             (loop (cdr table)))))))
+
+(define (pre-1972-tai->utc ns tai-seconds)
+  (let ((tai-ns (+ ns (* nano tai-seconds))))
+    (let loop ((table tai-utc-1961-1971-table))
+      (cond ((null? table)
+             (values ns tai-seconds))
+            ((>= tai-ns (+ (caar table) (cadar table)))
+             (apply (lambda (utc-start address@hidden slope)
+                      (let ((utc-ns (+ utc-start
+                                       (round-quotient (* nano (- tai-ns
+                                                                  utc-start
+                                                                  
address@hidden))
+                                                       (+ nano slope)))))
+                        (values (truncate-remainder utc-ns nano)
+                                (truncate-quotient utc-ns nano))))
+                    (car table)))
+            (else
+             (loop (cdr table)))))))
+
+
 ;;; the TIME structure; creates the accessors, too.
 
 (define-record-type time
-- 
2.19.1