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
clockcalib(uint64_t (* clk)(void),const char * clkname)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 }
181