1 /*-
2  * SPDX-License-Identifier: BSD-3-Clause
3  *
4  * Copyright (c) 1985, 1993
5  *	The Regents of the University of California.  All rights reserved.
6  *
7  * Redistribution and use in source and binary forms, with or without
8  * modification, are permitted provided that the following conditions
9  * are met:
10  * 1. Redistributions of source code must retain the above copyright
11  *    notice, this list of conditions and the following disclaimer.
12  * 2. Redistributions in binary form must reproduce the above copyright
13  *    notice, this list of conditions and the following disclaimer in the
14  *    documentation and/or other materials provided with the distribution.
15  * 3. Neither the name of the University nor the names of its contributors
16  *    may be used to endorse or promote products derived from this software
17  *    without specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
20  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
23  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29  * SUCH DAMAGE.
30  */
31 
32 #ifndef lint
33 #if 0
34 static char sccsid[] = "@(#)networkdelta.c	8.1 (Berkeley) 6/6/93";
35 #endif
36 static const char rcsid[] =
37   "$FreeBSD$";
38 #endif /* not lint */
39 
40 #include "globals.h"
41 
42 static long median(float, float *, long *, long *, unsigned int);
43 
44 /*
45  * Compute a corrected date.
46  *	Compute the median of the reasonable differences.  First compute
47  *	the median of all authorized differences, and then compute the
48  *	median of all differences that are reasonably close to the first
49  *	median.
50  *
51  * This differs from the original BSD implementation, which looked for
52  *	the largest group of machines with essentially the same date.
53  *	That assumed that machines with bad clocks would be uniformly
54  *	distributed.  Unfortunately, in real life networks, the distribution
55  *	of machines is not uniform among models of machines, and the
56  *	distribution of errors in clocks tends to be quite consistent
57  *	for a given model.  In other words, all model VI Supre Servres
58  *	from GoFast Inc. tend to have about the same error.
59  *	The original BSD implementation would chose the clock of the
60  *	most common model, and discard all others.
61  *
62  *	Therefore, get best we can do is to try to average over all
63  *	of the machines in the network, while discarding "obviously"
64  *	bad values.
65  */
66 long
networkdelta(void)67 networkdelta(void)
68 {
69 	struct hosttbl *htp;
70 	long med;
71 	long lodelta, hidelta;
72 	long logood, higood;
73 	long x[NHOSTS];
74 	long *xp;
75 	int numdelta;
76 	float eps;
77 
78 	/*
79 	 * compute the median of the good values
80 	 */
81 	med = 0;
82 	numdelta = 1;
83 	xp = &x[0];
84 	*xp = 0;			/* account for ourself */
85 	for (htp = self.l_fwd; htp != &self; htp = htp->l_fwd) {
86 		if (htp->good
87 		    && htp->noanswer == 0
88 		    && htp->delta != HOSTDOWN) {
89 			med += htp->delta;
90 			numdelta++;
91 			*++xp = htp->delta;
92 		}
93 	}
94 
95 	/*
96 	 * If we are the only trusted time keeper, then do not change our
97 	 * clock.  There may be another time keeping service active.
98 	 */
99 	if (numdelta == 1)
100 		return 0;
101 
102 	med /= numdelta;
103 	eps = med - x[0];
104 	if (trace)
105 		fprintf(fd, "median of %d values starting at %ld is about ",
106 			numdelta, med);
107 	med = median(med, &eps, &x[0], xp+1, VALID_RANGE);
108 
109 	/*
110 	 * compute the median of all values near the good median
111 	 */
112 	hidelta = med + GOOD_RANGE;
113 	lodelta = med - GOOD_RANGE;
114 	higood = med + VGOOD_RANGE;
115 	logood = med - VGOOD_RANGE;
116 	xp = &x[0];
117 	htp = &self;
118 	do {
119 		if (htp->noanswer == 0
120 		    && htp->delta >= lodelta
121 		    && htp->delta <= hidelta
122 		    && (htp->good
123 			|| (htp->delta >= logood
124 			    && htp->delta <= higood))) {
125 			*xp++ = htp->delta;
126 		}
127 	} while (&self != (htp = htp->l_fwd));
128 
129 	if (xp == &x[0]) {
130 		if (trace)
131 			fprintf(fd, "nothing close to median %ld\n", med);
132 		return med;
133 	}
134 
135 	if (xp == &x[1]) {
136 		if (trace)
137 			fprintf(fd, "only value near median is %ld\n", x[0]);
138 		return x[0];
139 	}
140 
141 	if (trace)
142 		fprintf(fd, "median of %td values starting at %ld is ",
143 		        xp-&x[0], med);
144 	return median(med, &eps, &x[0], xp, 1);
145 }
146 
147 
148 /*
149  * compute the median of an array of signed integers, using the idea
150  *	in <<Numerical Recipes>>.
151  */
152 static long
median(float a,float * eps_ptr,long * x,long * xlim,unsigned int gnuf)153 median(float a, float *eps_ptr, long *x, long *xlim, unsigned int gnuf)
154 	/* float a; */			/* initial guess for the median */
155 	/* float *eps_ptr; */		/* spacing near the median */
156 	/* long *x, *xlim; */		/* the data */
157 	/* unsigned int gnuf; */	/* good enough estimate */
158 {
159 	long *xptr;
160 	float ap = LONG_MAX;		/* bounds on the median */
161 	float am = -LONG_MAX;
162 	float aa;
163 	int npts;			/* # of points above & below guess */
164 	float xp;			/* closet point above the guess */
165 	float xm;			/* closet point below the guess */
166 	float eps;
167 	float dum, sum, sumx;
168 	int pass;
169 #define AMP	1.5			/* smoothing constants */
170 #define AFAC	1.5
171 
172 	eps = *eps_ptr;
173 	if (eps < 1.0) {
174 		eps = -eps;
175 		if (eps < 1.0)
176 			eps = 1.0;
177 	}
178 
179 	for (pass = 1; ; pass++) {	/* loop over the data */
180 		sum = 0.0;
181 		sumx = 0.0;
182 		npts = 0;
183 		xp = LONG_MAX;
184 		xm = -LONG_MAX;
185 
186 		for (xptr = x; xptr != xlim; xptr++) {
187 			float xx = *xptr;
188 
189 			dum = xx - a;
190 			if (dum != 0.0) {	/* avoid dividing by 0 */
191 				if (dum > 0.0) {
192 					npts++;
193 					if (xx < xp)
194 						xp = xx;
195 				} else {
196 					npts--;
197 					if (xx > xm)
198 						xm = xx;
199 					dum = -dum;
200 				}
201 				dum = 1.0/(eps + dum);
202 				sum += dum;
203 				sumx += xx * dum;
204 			}
205 		}
206 
207 		if (ap-am < gnuf || sum == 0) {
208 			if (trace)
209 				fprintf(fd,
210 			           "%ld in %d passes; early out balance=%d\n",
211 				        (long)a, pass, npts);
212 			return a;	/* guess was good enough */
213 		}
214 
215 		aa = (sumx/sum-a)*AMP;
216 		if (npts >= 2) {	/* guess was too low */
217 			am = a;
218 			aa = xp + max(0.0, aa);
219 			if (aa > ap)
220 				aa = (a + ap)/2;
221 
222 		} else if (npts <= -2) {  /* guess was two high */
223 			ap = a;
224 			aa = xm + min(0.0, aa);
225 			if (aa < am)
226 				aa = (a + am)/2;
227 
228 		} else {
229 			break;		/* got it */
230 		}
231 
232 		if (a == aa) {
233 			if (trace)
234 				fprintf(fd,
235 				  "%ld in %d passes; force out balance=%d\n",
236 				        (long)a, pass, npts);
237 			return a;
238 		}
239 		eps = AFAC*abs(aa - a);
240 		*eps_ptr = eps;
241 		a = aa;
242 	}
243 
244 	if (((x - xlim) % 2) != 0) {    /* even number of points? */
245 		if (npts == 0)		/* yes, return an average */
246 			a = (xp+xm)/2;
247 		else if (npts > 0)
248 			a =  (a+xp)/2;
249 		else
250 			a = (xm+a)/2;
251 
252 	} else 	if (npts != 0) {	/* odd number of points */
253 		if (npts > 0)
254 			a = xp;
255 		else
256 			a = xm;
257 	}
258 
259 	if (trace)
260 		fprintf(fd, "%ld in %d passes\n", (long)a, pass);
261 	return a;
262 }
263