sys/kern/subr_clockcalib.c

   1 /*-
   2  * Copyright (c) 2022 Colin Percival
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  */
  26
  27 #include <sys/param.h>
  28 #include <sys/systm.h>
  29 #include <sys/timetc.h>
  30 #include <sys/tslog.h>
  31 #include <machine/cpu.h>
  32
  33 /**
  34  * clockcalib(clk, clkname):
  35  * Return the frequency of the provided timer, as calibrated against the
  36  * current best-available timecounter.
  37  */
  38 uint64_t
  39 clockcalib(uint64_t (*clk)(void), const char *clkname)
  40 {
  41         struct timecounter *tc = atomic_load_ptr(&timecounter);
  42         uint64_t clk0, clk1, clk_delay, n, passes = 0;
  43         uint64_t t0, t1, tadj, tlast;
  44         double mu_clk = 0;
  45         double mu_t = 0;
  46         double va_clk = 0;
  47         double va_t = 0;
  48         double cva = 0;
  49         double d1, d2;
  50         double inv_n;
  51         uint64_t freq;
  52
  53         TSENTER();
  54         /*-
  55          * The idea here is to compute a best-fit linear regression between
  56          * the clock we're calibrating and the reference clock; the slope of
  57          * that line multiplied by the frequency of the reference clock gives
  58          * us the frequency we're looking for.
  59          *
  60          * To do this, we calculate the
  61          * (a) mean of the target clock measurements,
  62          * (b) variance of the target clock measurements,
  63          * (c) mean of the reference clock measurements,
  64          * (d) variance of the reference clock measurements, and
  65          * (e) covariance of the target clock and reference clock measurements
  66          * on an ongoing basis, updating all five values after each new data
  67          * point arrives, stopping when we're confident that we've accurately
  68          * measured the target clock frequency.
  69          *
  70          * Given those five values, the important formulas to remember from
  71          * introductory statistics are:
  72          * 1. slope of regression line = covariance(x, y) / variance(x)
  73          * 2. (relative uncertainty in slope)^2 =
  74          *    (variance(x) * variance(y) - covariance(x, y)^2)
  75          *    ------------------------------------------------
  76          *              covariance(x, y)^2 * (N - 2)
  77          *
  78          * We adjust the second formula slightly, adding a term to each of
  79          * the variance values to reflect the measurement quantization.
  80          *
  81          * Finally, we need to determine when to stop gathering data.  We
  82          * can't simply stop as soon as the computed uncertainty estimate
  83          * is below our threshold; this would make us overconfident since it
  84          * would introduce a multiple-comparisons problem (cf. sequential
  85          * analysis in clinical trials).  Instead, we stop with N data points
  86          * if the estimated uncertainty of the first k data points meets our
  87          * target for all N/2 < k <= N; this is not theoretically optimal,
  88          * but in practice works well enough.
  89          */
  90
  91         /*
  92          * Initial values for clocks; we'll subtract these off from values
  93          * we measure later in order to reduce floating-point rounding errors.
  94          * We keep track of an adjustment for values read from the reference
  95          * timecounter, since it can wrap.
  96          */
  97         clk0 = clk();
  98         t0 = tc->tc_get_timecount(tc) & tc->tc_counter_mask;
  99         tadj = 0;
 100         tlast = t0;
 101
 102         /* Loop until we give up or decide that we're calibrated. */
 103         for (n = 1; ; n++) {
 104                 /* Get a new data point. */
 105                 clk1 = clk() - clk0;
 106                 t1 = tc->tc_get_timecount(tc) & tc->tc_counter_mask;
 107                 while (t1 + tadj < tlast)
 108                         tadj += (uint64_t)tc->tc_counter_mask + 1;
 109                 tlast = t1 + tadj;
 110                 t1 += tadj - t0;
 111
 112                 /* If we spent too long, bail. */
 113                 if (t1 > tc->tc_frequency) {
 114                         printf("Statistical %s calibration failed!  "
 115                             "Clocks might be ticking at variable rates.\n",
 116                              clkname);
 117                         printf("Falling back to slow %s calibration.\n",
 118                             clkname);
 119                         freq = (double)(tc->tc_frequency) * clk1 / t1;
 120                         break;
 121                 }
 122
 123                 /* Precompute to save on divisions later. */
 124                 inv_n = 1.0 / n;
 125
 126                 /* Update mean and variance of recorded TSC values. */
 127                 d1 = clk1 - mu_clk;
 128                 mu_clk += d1 * inv_n;
 129                 d2 = d1 * (clk1 - mu_clk);
 130                 va_clk += (d2 - va_clk) * inv_n;
 131
 132                 /* Update mean and variance of recorded time values. */
 133                 d1 = t1 - mu_t;
 134                 mu_t += d1 * inv_n;
 135                 d2 = d1 * (t1 - mu_t);
 136                 va_t += (d2 - va_t) * inv_n;
 137
 138                 /* Update covariance. */
 139                 d2 = d1 * (clk1 - mu_clk);
 140                 cva += (d2 - cva) * inv_n;
 141
 142                 /*
 143                  * Count low-uncertainty iterations.  This is a rearrangement
 144                  * of "relative uncertainty < 1 PPM" avoiding division.
 145                  */
 146 #define TSC_PPM_UNCERTAINTY     1
 147 #define TSC_UNCERTAINTY         TSC_PPM_UNCERTAINTY * 0.000001
 148 #define TSC_UNCERTAINTY_SQR     TSC_UNCERTAINTY * TSC_UNCERTAINTY
 149                 if (TSC_UNCERTAINTY_SQR * (n - 2) * cva * cva >
 150                     (va_t + 4) * (va_clk + 4) - cva * cva)
 151                         passes++;
 152                 else
 153                         passes = 0;
 154
 155                 /* Break if we're consistently certain. */
 156                 if (passes * 2 > n) {
 157                         freq = (double)(tc->tc_frequency) * cva / va_t;
 158                         if (bootverbose)
 159                                 printf("Statistical %s calibration took"
 160                                     " %lu us and %lu data points\n",
 161                                     clkname, (unsigned long)(t1 *
 162                                         1000000.0 / tc->tc_frequency),
 163                                     (unsigned long)n);
 164                         break;
 165                 }
 166
 167                 /*
 168                  * Add variable delay to avoid theoretical risk of aliasing
 169                  * resulting from this loop synchronizing with the frequency
 170                  * of the reference clock.  On the nth iteration, we spend
 171                  * O(1 / n) time here -- long enough to avoid aliasing, but
 172                  * short enough to be insignificant as n grows.
 173                  */
 174                 clk_delay = clk() + (clk() - clk0) / (n * n);
 175                 while (clk() < clk_delay)
 176                         cpu_spinwait(); /* Do nothing. */
 177         }
 178         TSEXIT();
 179         return (freq);
 180 }