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