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 }
 276
 277 int SMT_IsEnabled(void)
 278 {
 279   return enabled;
 280 }
 281
 282 double
 283 SMT_GetOffset(struct timespec *now)
 284 {
 285   double offset, freq;
 286
 287   if (!enabled)
 288     return 0.0;
 289
 290   get_smoothing(now, &offset, &freq, NULL);
 291
 292   return offset;
 293 }
 294
 295 void
 296 SMT_Activate(struct timespec *now)
 297 {
 298   if (!enabled || !locked)
 299     return;
 300
 301   LOG(LOGS_INFO, "Time smoothing activated%s", leap_only_mode ?
 302       " (leap seconds only)" : "");
 303   locked = 0;
 304   last_update = *now;
 305 }
 306
 307 void
 308 SMT_Reset(struct timespec *now)
 309 {
 310   int i;
 311
 312   if (!enabled)
 313     return;
 314
 315   smooth_offset = 0.0;
 316   smooth_freq = 0.0;
 317   last_update = *now;
 318
 319   for (i = 0; i < NUM_STAGES; i++)
 320     stages[i].wander = stages[i].length = 0.0;
 321 }
 322
 323 void
 324 SMT_Leap(struct timespec *now, int leap)
 325 {
 326   /* When the leap-only mode is disabled, the leap second will be accumulated
 327      in handle_slew() as a normal offset */
 328   if (!enabled || !leap_only_mode)
 329     return;
 330
 331   update_smoothing(now, leap, 0.0);
 332 }
 333
 334 int
 335 SMT_GetSmoothingReport(RPT_SmoothingReport *report, struct timespec *now)
 336 {
 337   double length, elapsed;
 338   int i;
 339
 340   if (!enabled)
 341     return 0;
 342
 343   report->active = !locked;
 344   report->leap_only = leap_only_mode;
 345
 346   get_smoothing(now, &report->offset, &report->freq_ppm, &report->wander_ppm);
 347
 348   /* Convert to ppm and negate (positive values mean faster/speeding up) */
 349   report->freq_ppm *= -1.0e6;
 350   report->wander_ppm *= -1.0e6;
 351
 352   elapsed = UTI_DiffTimespecsToDouble(now, &last_update);
 353   if (!locked && elapsed >= 0.0) {
 354     for (i = 0, length = 0.0; i < NUM_STAGES; i++)
 355       length += stages[i].length;
 356     report->last_update_ago = elapsed;
 357     report->remaining_time = elapsed < length ? length - elapsed : 0.0;
 358   } else {
 359     report->last_update_ago = 0.0;
 360     report->remaining_time = 0.0;
 361   }
 362
 363   return 1;
 364 }