ntp: add server time smoothing

Time smoothing determines an offset that needs to be applied to the
cooked time to make it smooth for external observers. Observed offset
and frequency change slowly and there are no discontinuities. This can
be used on an NTP server to make it easier for the clients to track the
time and keep their clocks close together even when large offset or
frequency corrections are applied to the server's clock (e.g. after
being offline for longer time).

Accumulated offset and frequency are smoothed out in three stages. In
the first stage, the frequency is changed at a constant rate (wander) up
to a maximum, in the second stage the frequency stays at the maximum for
as long as needed and in the third stage the frequency is brought back
to zero.

Time smoothing is configured by the smoothtime directive. It takes two
arguments, maximum frequency offset and maximum wander. It's disabled by
default.

Makefile.in

@@ -40,7 +40,7 @@ HASH_OBJ = @HASH_OBJ@
 
 OBJS = array.o cmdparse.o conf.o local.o logging.o main.o memory.o mkdirpp.o \
        reference.o regress.o rtc.o sched.o sources.o sourcestats.o stubs.o \
-       sys.o tempcomp.o util.o $(HASH_OBJ)
+       sys.o smooth.o tempcomp.o util.o $(HASH_OBJ)
 
 EXTRA_OBJS=@EXTRA_OBJECTS@
 

conf.c

@@ -73,6 +73,7 @@ static void parse_peer(char *);
 static void parse_pool(char *);
 static void parse_refclock(char *);
 static void parse_server(char *);
+static void parse_smoothtime(char *);
 static void parse_tempcomp(char *);
 
 /* ================================================== */
@@ -185,6 +186,10 @@ static IPAddr bind_cmd_address4, bind_cmd_address6;
  * chronyds being started. */
 static char *pidfile;
 
+/* Smoothing constants */
+static double smooth_max_freq = 0.0; /* in ppm */
+static double smooth_max_wander = 0.0; /* in ppm/s */
+
 /* Temperature sensor, update interval and compensation coefficients */
 static char *tempcomp_sensor_file = NULL;
 static char *tempcomp_point_file = NULL;
@@ -512,6 +517,8 @@ CNF_ParseLine(const char *filename, int number, char *line)
     parse_int(p, &sched_priority);
   } else if (!strcasecmp(command, "server")) {
     parse_server(p);
+  } else if (!strcasecmp(command, "smoothtime")) {
+    parse_smoothtime(p);
   } else if (!strcasecmp(command, "stratumweight")) {
     parse_double(p, &stratum_weight);
   } else if (!strcasecmp(command, "tempcomp")) {
@@ -1164,6 +1171,17 @@ parse_broadcast(char *line)
 
 /* ================================================== */
 
+static void
+parse_smoothtime(char *line)
+{
+  check_number_of_args(line, 2);
+  if (sscanf(line, "%lf %lf", &smooth_max_freq, &smooth_max_wander) != 2) {
+    smooth_max_freq = 0.0;
+    command_parse_error();
+  }
+}
+
+/* ================================================== */
 static void
 parse_tempcomp(char *line)
 {
@@ -1719,6 +1737,15 @@ CNF_GetLockMemory(void)
 
 /* ================================================== */
 
+void
+CNF_GetSmooth(double *max_freq, double *max_wander)
+{
+  *max_freq = smooth_max_freq;
+  *max_wander = smooth_max_wander;
+}
+
+/* ================================================== */
+
 void
 CNF_GetTempComp(char **file, double *interval, char **point_file, double *T0, double *k0, double *k1, double *k2)
 {

conf.h

@@ -96,6 +96,7 @@ extern void CNF_SetupAccessRestrictions(void);
 extern int CNF_GetSchedPriority(void);
 extern int CNF_GetLockMemory(void);
 
+extern void CNF_GetSmooth(double *max_freq, double *max_wander);
 extern void CNF_GetTempComp(char **file, double *interval, char **point_file, double *T0, double *k0, double *k1, double *k2);
 
 extern char *CNF_GetUser(void);

local.c

@@ -36,6 +36,7 @@
 #include "local.h"
 #include "localp.h"
 #include "memory.h"
+#include "smooth.h"
 #include "util.h"
 #include "logging.h"
 
@@ -488,6 +489,9 @@ LCL_ApplyStepOffset(double offset)
 
   (*drv_apply_step_offset)(offset);
 
+  /* Reset smoothing on all clock steps */
+  SMT_Reset(&cooked);
+
   /* Dispatch to all handlers */
   invoke_parameter_change_handlers(&raw, &cooked, 0.0, offset, LCL_ChangeStep);
 }

logging.h

@@ -100,7 +100,8 @@ typedef enum {
   LOGF_SysWinnt,
   LOGF_TempComp,
   LOGF_RtcLinux,
-  LOGF_Refclock
+  LOGF_Refclock,
+  LOGF_Smooth,
 } LOG_Facility;
 
 /* Init function */

main.c

@@ -49,6 +49,7 @@
 #include "refclock.h"
 #include "clientlog.h"
 #include "nameserv.h"
+#include "smooth.h"
 #include "tempcomp.h"
 
 /* ================================================== */
@@ -88,6 +89,7 @@ MAI_CleanupAndExit(void)
   /* Don't update clock when removing sources */
   REF_SetMode(REF_ModeIgnore);
 
+  SMT_Finalise();
   TMC_Finalise();
   MNL_Finalise();
   CLG_Finalise();
@@ -495,6 +497,7 @@ int main
   CLG_Initialise();
   MNL_Initialise();
   TMC_Initialise();
+  SMT_Initialise();
 
   /* From now on, it is safe to do finalisation on exit */
   initialised = 1;

ntp_core.c

@@ -36,6 +36,7 @@
 #include "sched.h"
 #include "reference.h"
 #include "local.h"
+#include "smooth.h"
 #include "sources.h"
 #include "util.h"
 #include "conf.h"
@@ -757,14 +758,14 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
 {
   NTP_Packet message;
   int leap, auth_len, length, ret;
-  struct timeval local_transmit;
+  struct timeval local_receive, local_transmit;
 
   /* Parameters read from reference module */
-  int are_we_synchronised, our_stratum;
+  int are_we_synchronised, our_stratum, smooth_time;
   NTP_Leap leap_status;
   uint32_t our_ref_id, ts_fuzz;
   struct timeval our_ref_time;
-  double our_root_delay, our_root_dispersion;
+  double our_root_delay, our_root_dispersion, smooth_offset;
 
   /* Don't reply with version higher than ours */
   if (version > NTP_VERSION) {
@@ -781,6 +782,15 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
                          &our_ref_id, &our_ref_time,
                          &our_root_delay, &our_root_dispersion);
 
+  /* Get current smoothing offset when sending packet to a client */
+  if (SMT_IsEnabled() && (my_mode == MODE_SERVER || my_mode == MODE_BROADCAST)) {
+    smooth_time = 1;
+    smooth_offset = SMT_GetOffset(&local_transmit);
+  } else {
+    smooth_time = 0;
+    smooth_offset = 0.0;
+  }
+
   if (are_we_synchronised) {
     leap = (int) leap_status;
   } else {
@@ -807,6 +817,14 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
   message.reference_id = htonl((NTP_int32) our_ref_id);
 
   /* Now fill in timestamps */
+
+  if (smooth_time) {
+    UTI_AddDoubleToTimeval(&our_ref_time, smooth_offset, &our_ref_time);
+    UTI_AddDoubleToTimeval(local_rx, smooth_offset, &local_receive);
+  } else {
+    local_receive = *local_rx;
+  }
+
   UTI_TimevalToInt64(&our_ref_time, &message.reference_ts, 0);
 
   /* Originate - this comes from the last packet the source sent us */
@@ -816,7 +834,7 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
      This timestamp will have been adjusted so that it will now look to
      the source like we have been running on our latest estimate of
      frequency all along */
-  UTI_TimevalToInt64(local_rx, &message.receive_ts, 0);
+  UTI_TimevalToInt64(&local_receive, &message.receive_ts, 0);
 
   /* Prepare random bits which will be added to the transmit timestamp. */
   ts_fuzz = UTI_GetNTPTsFuzz(message.precision);
@@ -826,6 +844,9 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
      from the source we're sending to now. */
   LCL_ReadCookedTime(&local_transmit, NULL);
 
+  if (smooth_time)
+    UTI_AddDoubleToTimeval(&local_transmit, smooth_offset, &local_transmit);
+
   length = NTP_NORMAL_PACKET_LENGTH;
 
   /* Authenticate */

reference.c

@@ -1241,6 +1241,14 @@ REF_GetOurStratum(void)
 
 /* ================================================== */
 
+double
+REF_GetSkew(void)
+{
+  return our_skew;
+}
+
+/* ================================================== */
+
 void
 REF_ModifyMaxupdateskew(double new_max_update_skew)
 {

reference.h

@@ -162,6 +162,9 @@ REF_SetUnsynchronised(void);
    synchronised */
 extern int REF_GetOurStratum(void);
 
+/* Return the current skew */
+extern double REF_GetSkew(void);
+
 /* Modify the setting for the maximum skew we are prepared to allow updates on (in ppm). */
 extern void REF_ModifyMaxupdateskew(double new_max_update_skew);
 

smooth.c

@@ -0,0 +1,267 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2015
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ * 
+ **********************************************************************
+
+  =======================================================================
+
+  Routines implementing time smoothing.
+
+  */
+
+#include "config.h"
+
+#include "sysincl.h"
+
+#include "conf.h"
+#include "local.h"
+#include "logging.h"
+#include "reference.h"
+#include "smooth.h"
+#include "util.h"
+
+/*
+  Time smoothing determines an offset that needs to be applied to the cooked
+  time to make it smooth for external observers.  Observed offset and frequency
+  change slowly and there are no discontinuities.  This can be used on an NTP
+  server to make it easier for the clients to track the time and keep their
+  clocks close together even when large offset or frequency corrections are
+  applied to the server's clock (e.g. after being offline for longer time).
+
+  Accumulated offset and frequency are smoothed out in three stages.  In the
+  first stage, the frequency is changed at a constant rate (wander) up to a
+  maximum, in the second stage the frequency stays at the maximum for as long
+  as needed and in the third stage the frequency is brought back to zero.
+
+              |
+    max_freq  +-------/--------\-------------
+              |      /|        |\
+        freq  |     / |        | \
+              |    /  |        |  \
+              |   /   |        |   \
+           0  +--/----+--------+----\--------
+              | /     |        |    |    time
+              |/      |        |    |
+
+        stage     1       2      3
+
+  Integral of this function is the smoothed out offset.  It's a continuous
+  piecewise polynomial with two quadratic parts and one linear.
+*/
+
+struct stage {
+  double wander;
+  double length;
+};
+
+#define NUM_STAGES 3
+
+static struct stage stages[NUM_STAGES];
+
+/* Enabled/disabled smoothing */
+static int enabled;
+
+/* Maximum skew/max_wander ratio to start updating offset and frequency */
+#define UNLOCK_SKEW_WANDER_RATIO 10000
+
+static int locked;
+
+/* Maximum wander and frequency offset */
+static double max_wander;
+static double max_freq;
+
+/* Frequency offset, time offset and the time of the last smoothing update */
+static double smooth_freq;
+static double smooth_offset;
+static struct timeval last_update;
+
+
+static void
+get_offset_freq(struct timeval *now, double *offset, double *freq)
+{
+  double elapsed, length;
+  int i;
+
+  UTI_DiffTimevalsToDouble(&elapsed, now, &last_update);
+
+  *offset = smooth_offset;
+  *freq = smooth_freq;
+
+  for (i = 0; i < NUM_STAGES; i++) {
+    if (elapsed <= 0.0)
+      break;
+
+    length = stages[i].length;
+    if (length >= elapsed)
+      length = elapsed;
+
+    *offset -= length * (2.0 * *freq + stages[i].wander * length) / 2.0;
+    *freq += stages[i].wander * length;
+    elapsed -= length;
+  }
+
+  if (elapsed > 0.0)
+    *offset -= elapsed * *freq;
+}
+
+static void
+update_stages(void)
+{
+  double s1, s2, s, l1, l2, l3, lc, f, f2;
+  int i, dir;
+
+  /* Prepare the three stages so that the integral of the frequency offset
+     is equal to the offset that should be smoothed out */
+
+  s1 = smooth_offset / max_wander;
+  s2 = smooth_freq * smooth_freq / (2.0 * max_wander * max_wander);
+  
+  l1 = l2 = l3 = 0.0;
+
+  /* Calculate the lengths of the 1st and 3rd stage assuming there is no
+     frequency limit.  If length of the 1st stage comes out negative, switch
+     its direction. */
+  for (dir = -1; dir <= 1; dir += 2) {
+    s = dir * s1 + s2;
+    if (s >= 0.0) {
+      l3 = sqrt(s);
+      l1 = l3 - dir * smooth_freq / max_wander;
+      if (l1 >= 0.0)
+        break;
+    }
+  }
+
+  assert(dir <= 1 && l1 >= 0.0 && l3 >= 0.0);
+
+  /* If the limit was reached, shorten 1st+3rd stages and set a 2nd stage */
+  f = dir * smooth_freq + l1 * max_wander - max_freq;
+  if (f > 0.0) {
+    lc = f / max_wander;
+
+    /* No 1st stage if the frequency is already above the maximum */
+    if (lc > l1) {
+      lc = l1;
+      f2 = dir * smooth_freq;
+    } else {
+      f2 = max_freq;
+    }
+
+    l2 = lc * (2.0 + f / f2);
+    l1 -= lc;
+    l3 -= lc;
+  }
+
+  stages[0].wander = dir * max_wander;
+  stages[0].length = l1;
+  stages[1].wander = 0.0;
+  stages[1].length = l2;
+  stages[2].wander = -dir * max_wander;
+  stages[2].length = l3;
+
+  for (i = 0; i < NUM_STAGES; i++) {
+    DEBUG_LOG(LOGF_Smooth, "Smooth stage %d wander %e length %f",
+              i + 1, stages[i].wander, stages[i].length);
+  }
+}
+
+static void
+update_smoothing(struct timeval *now, double offset, double freq)
+{
+  /* Don't accept offset/frequency until the clock has stabilized */
+  if (locked) {
+    if (REF_GetSkew() / max_wander < UNLOCK_SKEW_WANDER_RATIO) {
+      LOG(LOGS_INFO, LOGF_Smooth, "Time smoothing activated");
+      locked = 0;
+    }
+    return;
+  }
+
+  get_offset_freq(now, &smooth_offset, &smooth_freq);
+  smooth_offset += offset;
+  smooth_freq = (smooth_freq - freq) / (1.0 - freq);
+  last_update = *now;
+
+  update_stages();
+
+  DEBUG_LOG(LOGF_Smooth, "Smooth offset %e freq %e", smooth_offset, smooth_freq);
+}
+
+static void
+handle_slew(struct timeval *raw, struct timeval *cooked, double dfreq,
+            double doffset, LCL_ChangeType change_type, void *anything)
+{
+  double delta;
+
+  if (change_type == LCL_ChangeAdjust)
+    update_smoothing(cooked, doffset, dfreq);
+
+  UTI_AdjustTimeval(&last_update, cooked, &last_update, &delta, dfreq, doffset);
+}
+
+void SMT_Initialise(void)
+{
+  CNF_GetSmooth(&max_freq, &max_wander);
+  if (max_freq <= 0.0 || max_wander <= 0.0) {
+      enabled = 0;
+      return;
+  }
+
+  enabled = 1;
+  locked = 1;
+
+  /* Convert from ppm */
+  max_freq *= 1e-6;
+  max_wander *= 1e-6;
+
+  LCL_AddParameterChangeHandler(handle_slew, NULL);
+}
+
+void SMT_Finalise(void)
+{
+}
+
+int SMT_IsEnabled(void)
+{
+  return enabled;
+}
+
+double
+SMT_GetOffset(struct timeval *now)
+{
+  double offset, freq;
+
+  if (!enabled)
+    return 0.0;
+  
+  get_offset_freq(now, &offset, &freq);
+
+  return offset;
+}
+
+void
+SMT_Reset(struct timeval *now)
+{
+  if (!enabled)
+    return;
+
+  locked = 1;
+  smooth_offset = 0.0;
+  smooth_freq = 0.0;
+  last_update = *now;
+}

smooth.h

@@ -0,0 +1,40 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2015
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ * 
+ **********************************************************************
+
+  =======================================================================
+
+  This module implements time smoothing.
+  */
+
+#ifndef GOT_SMOOTH_H
+#define GOT_SMOOTH_H
+
+extern void SMT_Initialise(void);
+
+extern void SMT_Finalise(void);
+
+extern int SMT_IsEnabled(void);
+
+extern double SMT_GetOffset(struct timeval *now);
+
+extern void SMT_Reset(struct timeval *now);
+
+#endif