smooth.c

   1 /*
   2   chronyd/chronyc - Programs for keeping computer clocks accurate.
   3
   4  **********************************************************************
   5  * Copyright (C) Miroslav Lichvar  2015
   6  *
   7  * This program is free software; you can redistribute it and/or modify
   8  * it under the terms of version 2 of the GNU General Public License as
   9  * published by the Free Software Foundation.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License along
  17  * with this program; if not, write to the Free Software Foundation, Inc.,
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  **********************************************************************
  21
  22   =======================================================================
  23
  24   Routines implementing time smoothing.
  25
  26   */
  27
  28 #include "config.h"
  29
  30 #include "sysincl.h"
  31
  32 #include "conf.h"
  33 #include "local.h"
  34 #include "logging.h"
  35 #include "reference.h"
  36 #include "smooth.h"
  37 #include "util.h"
  38
  39 /*
  40   Time smoothing determines an offset that needs to be applied to the cooked
  41   time to make it smooth for external observers.  Observed offset and frequency
  42   change slowly and there are no discontinuities.  This can be used on an NTP
  43   server to make it easier for the clients to track the time and keep their
  44   clocks close together even when large offset or frequency corrections are
  45   applied to the server's clock (e.g. after being offline for longer time).
  46
  47   Accumulated offset and frequency are smoothed out in three stages.  In the
  48   first stage, the frequency is changed at a constant rate (wander) up to a
  49   maximum, in the second stage the frequency stays at the maximum for as long
  50   as needed and in the third stage the frequency is brought back to zero.
  51
  52               |
  53     max_freq  +-------/--------\-------------
  54               |      /|        |\
  55         freq  |     / |        | \
  56               |    /  |        |  \
  57               |   /   |        |   \
  58            0  +--/----+--------+----\--------
  59               | /     |        |    |    time
  60               |/      |        |    |
  61
  62         stage     1       2      3
  63
  64   Integral of this function is the smoothed out offset.  It's a continuous
  65   piecewise polynomial with two quadratic parts and one linear.
  66 */
  67
  68 struct stage {
  69   double wander;
  70   double length;
  71 };
  72
  73 #define NUM_STAGES 3
  74
  75 static struct stage stages[NUM_STAGES];
  76
  77 /* Enabled/disabled smoothing */
  78 static int enabled;
  79
  80 /* Enabled/disabled mode where only leap seconds are smoothed out and normal
  81    offset/frequency changes are ignored */
  82 static int leap_only_mode;
  83
  84 /* Maximum skew/max_wander ratio to start updating offset and frequency */
  85 #define UNLOCK_SKEW_WANDER_RATIO 10000
  86
  87 static int locked;
  88
  89 /* Maximum wander and frequency offset */
  90 static double max_wander;
  91 static double max_freq;
  92
  93 /* Frequency offset, time offset and the time of the last smoothing update */
  94 static double smooth_freq;
  95 static double smooth_offset;
  96 static struct timespec last_update;
  97
  98
  99 static void
 100 get_smoothing(struct timespec *now, double *poffset, double *pfreq,
 101               double *pwander)
 102 {
 103   double elapsed, length, offset, freq, wander;
 104   int i;
 105
 106   elapsed = UTI_DiffTimespecsToDouble(now, &last_update);
 107
 108   offset = smooth_offset;
 109   freq = smooth_freq;
 110   wander = 0.0;
 111
 112   for (i = 0; i < NUM_STAGES; i++) {
 113     if (elapsed <= 0.0)
 114       break;
 115
 116     length = stages[i].length;
 117     if (length >= elapsed)
 118       length = elapsed;
 119
 120     wander = stages[i].wander;
 121     offset -= length * (2.0 * freq + wander * length) / 2.0;
 122     freq += wander * length;
 123     elapsed -= length;
 124   }
 125
 126   if (elapsed > 0.0) {
 127     wander = 0.0;
 128     offset -= elapsed * freq;
 129   }
 130
 131   *poffset = offset;
 132   *pfreq = freq;
 133   if (pwander)
 134     *pwander = wander;
 135 }
 136
 137 static void
 138 update_stages(void)
 139 {
 140   double s1, s2, s, l1, l2, l3, lc, f, f2, l1t[2], l3t[2], err[2];
 141   int i, dir;
 142
 143   /* Prepare the three stages so that the integral of the frequency offset
 144      is equal to the offset that should be smoothed out */
 145
 146   s1 = smooth_offset / max_wander;
 147   s2 = SQUARE(smooth_freq) / (2.0 * SQUARE(max_wander));
 148
 149   /* Calculate the lengths of the 1st and 3rd stage assuming there is no
 150      frequency limit.  The direction of the 1st stage is selected so that
 151      the lengths will not be negative.  With extremely small offsets both
 152      directions may give a negative length due to numerical errors, so select
 153      the one which gives a smaller error. */
 154
 155   for (i = 0, dir = -1; i <= 1; i++, dir += 2) {
 156     err[i] = 0.0;
 157     s = dir * s1 + s2;
 158
 159     if (s < 0.0) {
 160       err[i] += -s;
 161       s = 0.0;
 162     }
 163
 164     l3t[i] = sqrt(s);
 165     l1t[i] = l3t[i] - dir * smooth_freq / max_wander;
 166
 167     if (l1t[i] < 0.0) {
 168       err[i] += l1t[i] * l1t[i];
 169       l1t[i] = 0.0;
 170     }
 171   }
 172
 173   if (err[0] < err[1]) {
 174     l1 = l1t[0];
 175     l3 = l3t[0];
 176     dir = -1;
 177   } else {
 178     l1 = l1t[1];
 179     l3 = l3t[1];
 180     dir = 1;
 181   }
 182
 183   l2 = 0.0;
 184
 185   /* If the limit was reached, shorten 1st+3rd stages and set a 2nd stage */
 186   f = dir * smooth_freq + l1 * max_wander - max_freq;
 187   if (f > 0.0) {
 188     lc = f / max_wander;
 189
 190     /* No 1st stage if the frequency is already above the maximum */
 191     if (lc > l1) {
 192       lc = l1;
 193       f2 = dir * smooth_freq;
 194     } else {
 195       f2 = max_freq;
 196     }
 197
 198     l2 = lc * (2.0 + f / f2);
 199     l1 -= lc;
 200     l3 -= lc;
 201   }
 202
 203   stages[0].wander = dir * max_wander;
 204   stages[0].length = l1;
 205   stages[1].wander = 0.0;
 206   stages[1].length = l2;
 207   stages[2].wander = -dir * max_wander;
 208   stages[2].length = l3;
 209
 210   for (i = 0; i < NUM_STAGES; i++) {
 211     DEBUG_LOG("Smooth stage %d wander %e length %f",
 212               i + 1, stages[i].wander, stages[i].length);
 213   }
 214 }
 215
 216 static void
 217 update_smoothing(struct timespec *now, double offset, double freq)
 218 {
 219   /* Don't accept offset/frequency until the clock has stabilized */
 220   if (locked) {
 221     if (REF_GetSkew() / max_wander < UNLOCK_SKEW_WANDER_RATIO || leap_only_mode)
 222       SMT_Activate(now);
 223     return;
 224   }
 225
 226   get_smoothing(now, &smooth_offset, &smooth_freq, NULL);
 227   smooth_offset += offset;
 228   smooth_freq = (smooth_freq - freq) / (1.0 - freq);
 229   last_update = *now;
 230
 231   update_stages();
 232
 233   DEBUG_LOG("Smooth offset %e freq %e", smooth_offset, smooth_freq);
 234 }
 235
 236 static void
 237 handle_slew(struct timespec *raw, struct timespec *cooked, double dfreq,
 238             double doffset, LCL_ChangeType change_type, void *anything)
 239 {
 240   double delta;
 241
 242   if (change_type == LCL_ChangeAdjust) {
 243     if (leap_only_mode)
 244       update_smoothing(cooked, 0.0, 0.0);
 245     else
 246       update_smoothing(cooked, doffset, dfreq);
 247   }
 248
 249   if (!UTI_IsZeroTimespec(&last_update))
 250     UTI_AdjustTimespec(&last_update, cooked, &last_update, &delta, dfreq, doffset);
 251 }
 252
 253 void SMT_Initialise(void)
 254 {
 255   CNF_GetSmooth(&max_freq, &max_wander, &leap_only_mode);
 256   if (max_freq <= 0.0 || max_wander <= 0.0) {
 257       enabled = 0;
 258       return;
 259   }
 260
 261   enabled = 1;
 262   locked = 1;
 263
 264   /* Convert from ppm */
 265   max_freq *= 1e-6;
 266   max_wander *= 1e-6;
 267
 268   UTI_ZeroTimespec(&last_update);
 269
 270   LCL_AddParameterChangeHandler(handle_slew, NULL);
 271 }
 272
 273 void SMT_Finalise(void)
 274 {
 275   if (!enabled)
 276     return;
 277
 278   LCL_RemoveParameterChangeHandler(handle_slew, NULL);
 279 }
 280
 281 int SMT_IsEnabled(void)
 282 {
 283   return enabled;
 284 }
 285
 286 double
 287 SMT_GetOffset(struct timespec *now)
 288 {
 289   double offset, freq;
 290
 291   if (!enabled)
 292     return 0.0;
 293
 294   get_smoothing(now, &offset, &freq, NULL);
 295
 296   return offset;
 297 }
 298
 299 void
 300 SMT_Activate(struct timespec *now)
 301 {
 302   if (!enabled || !locked)
 303     return;
 304
 305   LOG(LOGS_INFO, "Time smoothing activated%s", leap_only_mode ?
 306       " (leap seconds only)" : "");
 307   locked = 0;
 308   last_update = *now;
 309 }
 310
 311 void
 312 SMT_Reset(struct timespec *now)
 313 {
 314   int i;
 315
 316   if (!enabled)
 317     return;
 318
 319   smooth_offset = 0.0;
 320   smooth_freq = 0.0;
 321   last_update = *now;
 322
 323   for (i = 0; i < NUM_STAGES; i++)
 324     stages[i].wander = stages[i].length = 0.0;
 325 }
 326
 327 void
 328 SMT_Leap(struct timespec *now, int leap)
 329 {
 330   /* When the leap-only mode is disabled, the leap second will be accumulated
 331      in handle_slew() as a normal offset */
 332   if (!enabled || !leap_only_mode)
 333     return;
 334
 335   update_smoothing(now, leap, 0.0);
 336 }
 337
 338 int
 339 SMT_GetSmoothingReport(RPT_SmoothingReport *report, struct timespec *now)
 340 {
 341   double length, elapsed;
 342   int i;
 343
 344   if (!enabled)
 345     return 0;
 346
 347   report->active = !locked;
 348   report->leap_only = leap_only_mode;
 349
 350   get_smoothing(now, &report->offset, &report->freq_ppm, &report->wander_ppm);
 351
 352   /* Convert to ppm and negate (positive values mean faster/speeding up) */
 353   report->freq_ppm *= -1.0e6;
 354   report->wander_ppm *= -1.0e6;
 355
 356   elapsed = UTI_DiffTimespecsToDouble(now, &last_update);
 357   if (!locked && elapsed >= 0.0) {
 358     for (i = 0, length = 0.0; i < NUM_STAGES; i++)
 359       length += stages[i].length;
 360     report->last_update_ago = elapsed;
 361     report->remaining_time = elapsed < length ? length - elapsed : 0.0;
 362   } else {
 363     report->last_update_ago = 0.0;
 364     report->remaining_time = 0.0;
 365   }
 366
 367   return 1;
 368 }