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