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 */ 38uint64_t 39clockcalib(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} 181