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/cdefs.h>
  28 #include <sys/param.h>
  29 #include <sys/systm.h>
  30 #include <sys/timetc.h>
  31 #include <sys/tslog.h>
  32 #include <machine/cpu.h>
  33
  34 /**
  35  * clockcalib(clk, clkname):
  36  * Return the frequency of the provided timer, as calibrated against the
  37  * current best-available timecounter.
  38  */
  39 uint64_t
  40 clockcalib(uint64_t (*clk)(void), const char *clkname)
  41 {
  42         struct timecounter *tc = atomic_load_ptr(&timecounter);
  43         uint64_t clk0, clk1, clk_delay, n, passes = 0;
  44         uint64_t t0, t1, tadj, tlast;
  45         double mu_clk = 0;
  46         double mu_t = 0;
  47         double va_clk = 0;
  48         double va_t = 0;
  49         double cva = 0;
  50         double d1, d2;
  51         double inv_n;
  52         uint64_t freq;
  53
  54         TSENTER();
  55         /*-
  56          * The idea here is to compute a best-fit linear regression between
  57          * the clock we're calibrating and the reference clock; the slope of
  58          * that line multiplied by the frequency of the reference clock gives
  59          * us the frequency we're looking for.
  60          *
  61          * To do this, we calculate the
  62          * (a) mean of the target clock measurements,
  63          * (b) variance of the target clock measurements,
  64          * (c) mean of the reference clock measurements,
  65          * (d) variance of the reference clock measurements, and
  66          * (e) covariance of the target clock and reference clock measurements
  67          * on an ongoing basis, updating all five values after each new data
  68          * point arrives, stopping when we're confident that we've accurately
  69          * measured the target clock frequency.
  70          *
  71          * Given those five values, the important formulas to remember from
  72          * introductory statistics are:
  73          * 1. slope of regression line = covariance(x, y) / variance(x)
  74          * 2. (relative uncertainty in slope)^2 =
  75          *    (variance(x) * variance(y) - covariance(x, y)^2)
  76          *    ------------------------------------------------
  77          *              covariance(x, y)^2 * (N - 2)
  78          *
  79          * We adjust the second formula slightly, adding a term to each of
  80          * the variance values to reflect the measurement quantization.
  81          *
  82          * Finally, we need to determine when to stop gathering data.  We
  83          * can't simply stop as soon as the computed uncertainty estimate
  84          * is below our threshold; this would make us overconfident since it
  85          * would introduce a multiple-comparisons problem (cf. sequential
  86          * analysis in clinical trials).  Instead, we stop with N data points
  87          * if the estimated uncertainty of the first k data points meets our
  88          * target for all N/2 < k <= N; this is not theoretically optimal,
  89          * but in practice works well enough.
  90          */
  91
  92         /*
  93          * Initial values for clocks; we'll subtract these off from values
  94          * we measure later in order to reduce floating-point rounding errors.
  95          * We keep track of an adjustment for values read from the reference
  96          * timecounter, since it can wrap.
  97          */
  98         clk0 = clk();
  99         t0 = tc->tc_get_timecount(tc) & tc->tc_counter_mask;
 100         tadj = 0;
 101         tlast = t0;
 102
 103         /* Loop until we give up or decide that we're calibrated. */
 104         for (n = 1; ; n++) {
 105                 /* Get a new data point. */
 106                 clk1 = clk() - clk0;
 107                 t1 = tc->tc_get_timecount(tc) & tc->tc_counter_mask;
 108                 while (t1 + tadj < tlast)
 109                         tadj += (uint64_t)tc->tc_counter_mask + 1;
 110                 tlast = t1 + tadj;
 111                 t1 += tadj - t0;
 112
 113                 /* If we spent too long, bail. */
 114                 if (t1 > tc->tc_frequency) {
 115                         printf("Statistical %s calibration failed!  "
 116                             "Clocks might be ticking at variable rates.\n",
 117                              clkname);
 118                         printf("Falling back to slow %s calibration.\n",
 119                             clkname);
 120                         freq = (double)(tc->tc_frequency) * clk1 / t1;
 121                         break;
 122                 }
 123
 124                 /* Precompute to save on divisions later. */
 125                 inv_n = 1.0 / n;
 126
 127                 /* Update mean and variance of recorded TSC values. */
 128                 d1 = clk1 - mu_clk;
 129                 mu_clk += d1 * inv_n;
 130                 d2 = d1 * (clk1 - mu_clk);
 131                 va_clk += (d2 - va_clk) * inv_n;
 132
 133                 /* Update mean and variance of recorded time values. */
 134                 d1 = t1 - mu_t;
 135                 mu_t += d1 * inv_n;
 136                 d2 = d1 * (t1 - mu_t);
 137                 va_t += (d2 - va_t) * inv_n;
 138
 139                 /* Update covariance. */
 140                 d2 = d1 * (clk1 - mu_clk);
 141                 cva += (d2 - cva) * inv_n;
 142
 143                 /*
 144                  * Count low-uncertainty iterations.  This is a rearrangement
 145                  * of "relative uncertainty < 1 PPM" avoiding division.
 146                  */
 147 #define TSC_PPM_UNCERTAINTY     1
 148 #define TSC_UNCERTAINTY         TSC_PPM_UNCERTAINTY * 0.000001
 149 #define TSC_UNCERTAINTY_SQR     TSC_UNCERTAINTY * TSC_UNCERTAINTY
 150                 if (TSC_UNCERTAINTY_SQR * (n - 2) * cva * cva >
 151                     (va_t + 4) * (va_clk + 4) - cva * cva)
 152                         passes++;
 153                 else
 154                         passes = 0;
 155
 156                 /* Break if we're consistently certain. */
 157                 if (passes * 2 > n) {
 158                         freq = (double)(tc->tc_frequency) * cva / va_t;
 159                         if (bootverbose)
 160                                 printf("Statistical %s calibration took"
 161                                     " %lu us and %lu data points\n",
 162                                     clkname, (unsigned long)(t1 *
 163                                         1000000.0 / tc->tc_frequency),
 164                                     (unsigned long)n);
 165                         break;
 166                 }
 167
 168                 /*
 169                  * Add variable delay to avoid theoretical risk of aliasing
 170                  * resulting from this loop synchronizing with the frequency
 171                  * of the reference clock.  On the nth iteration, we spend
 172                  * O(1 / n) time here -- long enough to avoid aliasing, but
 173                  * short enough to be insignificant as n grows.
 174                  */
 175                 clk_delay = clk() + (clk() - clk0) / (n * n);
 176                 while (clk() < clk_delay)
 177                         cpu_spinwait(); /* Do nothing. */
 178         }
 179         TSEXIT();
 180         return (freq);
 181 }