1 /*-
2  * Copyright (c) 1985, 1993
3  *	The Regents of the University of California.  All rights reserved.
4  *
5  * %sccs.include.redist.c%
6  */
7 
8 #ifndef lint
9 static char sccsid[] = "@(#)networkdelta.c	8.3 (Berkeley) 04/27/95";
10 #endif /* not lint */
11 
12 #ifdef sgi
13 #ident "$Revision: 1.4 $"
14 #endif
15 
16 #include "globals.h"
17 
18 static long median __P((double, float *, long *, long *, unsigned int));
19 
20 /*
21  * Compute a corrected date.
22  *	Compute the median of the reasonable differences.  First compute
23  *	the median of all authorized differences, and then compute the
24  *	median of all differences that are reasonably close to the first
25  *	median.
26  *
27  * This differs from the original BSD implementation, which looked for
28  *	the largest group of machines with essentially the same date.
29  *	That assumed that machines with bad clocks would be uniformly
30  *	distributed.  Unfortunately, in real life networks, the distribution
31  *	of machines is not uniform among models of machines, and the
32  *	distribution of errors in clocks tends to be quite consistent
33  *	for a given model.  In other words, all model VI Supre Servres
34  *	from GoFast Inc. tend to have about the same error.
35  *	The original BSD implementation would chose the clock of the
36  *	most common model, and discard all others.
37  *
38  *	Therefore, get best we can do is to try to average over all
39  *	of the machines in the network, while discarding "obviously"
40  *	bad values.
41  */
42 long
43 networkdelta()
44 {
45 	struct hosttbl *htp;
46 	long med;
47 	long lodelta, hidelta;
48 	long logood, higood;
49 	long x[NHOSTS];
50 	long *xp;
51 	int numdelta;
52 	float eps;
53 
54 	/*
55 	 * compute the median of the good values
56 	 */
57 	med = 0;
58 	numdelta = 1;
59 	xp = &x[0];
60 	*xp = 0;			/* account for ourself */
61 	for (htp = self.l_fwd; htp != &self; htp = htp->l_fwd) {
62 		if (htp->good
63 		    && htp->noanswer == 0
64 		    && htp->delta != HOSTDOWN) {
65 			med += htp->delta;
66 			numdelta++;
67 			*++xp = htp->delta;
68 		}
69 	}
70 
71 	/*
72 	 * If we are the only trusted time keeper, then do not change our
73 	 * clock.  There may be another time keeping service active.
74 	 */
75 	if (numdelta == 1)
76 		return 0;
77 
78 	/* get average of trusted values */
79 	med /= numdelta;
80 
81 	if (trace)
82 		fprintf(fd, "median of %d values starting at %ld is about ",
83 			numdelta, med);
84 	/* get median of all trusted values, using average as initial guess */
85 	eps = med - x[0];
86 	med = median(med, &eps, &x[0], xp+1, VALID_RANGE);
87 
88 	/* Compute the median of all good values.
89 	 * Good values are those of all clocks, including untrusted clocks,
90 	 * that are
91 	 *	- trusted and somewhat close to the median of the
92 	 *		trusted clocks
93 	 *	- trusted or untrusted and very close to the median of the
94 	 *		trusted clocks
95 	 */
96 	hidelta = med + GOOD_RANGE;
97 	lodelta = med - GOOD_RANGE;
98 	higood = med + VGOOD_RANGE;
99 	logood = med - VGOOD_RANGE;
100 	xp = &x[0];
101 	htp = &self;
102 	do {
103 		if (htp->noanswer == 0
104 		    && htp->delta >= lodelta
105 		    && htp->delta <= hidelta
106 		    && (htp->good
107 			|| (htp->delta >= logood
108 			    && htp->delta <= higood))) {
109 			*xp++ = htp->delta;
110 		}
111 	} while (&self != (htp = htp->l_fwd));
112 
113 	if (xp == &x[0]) {
114 		if (trace)
115 			fprintf(fd, "nothing close to median %ld\n", med);
116 		return med;
117 	}
118 
119 	if (xp == &x[1]) {
120 		if (trace)
121 			fprintf(fd, "only value near median is %ld\n", x[0]);
122 		return x[0];
123 	}
124 
125 	if (trace)
126 		fprintf(fd, "median of %d values starting at %ld is ",
127 		        xp-&x[0], med);
128 	return median(med, &eps, &x[0], xp, 1);
129 }
130 
131 
132 /*
133  * compute the median of an array of signed integers, using the idea
134  *	in <<Numerical Recipes>>.
135  */
136 static long
137 median(a0, eps_ptr, x, xlim, gnuf)
138 	double a0;			/* initial guess for the median */
139 	float *eps_ptr;			/* spacing near the median */
140 	long *x, *xlim;			/* the data */
141 	unsigned int gnuf;		/* good enough estimate */
142 {
143 	long *xptr;
144 	float a = a0;
145 	float ap = LONG_MAX;		/* bounds on the median */
146 	float am = -LONG_MAX;
147 	float aa;
148 	int npts;			/* # of points above & below guess */
149 	float xp;			/* closet point above the guess */
150 	float xm;			/* closet point below the guess */
151 	float eps;
152 	float dum, sum, sumx;
153 	int pass;
154 #define AMP	1.5			/* smoothing constants */
155 #define AFAC	1.5
156 
157 	eps = *eps_ptr;
158 	if (eps < 1.0) {
159 		eps = -eps;
160 		if (eps < 1.0)
161 			eps = 1.0;
162 	}
163 
164 	for (pass = 1; ; pass++) {	/* loop over the data */
165 		sum = 0.0;
166 		sumx = 0.0;
167 		npts = 0;
168 		xp = LONG_MAX;
169 		xm = -LONG_MAX;
170 
171 		for (xptr = x; xptr != xlim; xptr++) {
172 			float xx = *xptr;
173 
174 			dum = xx - a;
175 			if (dum != 0.0) {   /* avoid dividing by 0 */
176 				if (dum > 0.0) {
177 					npts++;
178 					if (xx < xp)
179 						xp = xx;
180 				} else {
181 					npts--;
182 					if (xx > xm)
183 						xm = xx;
184 					dum = -dum;
185 				}
186 				dum = 1.0/(eps + dum);
187 				sum += dum;
188 				sumx += xx * dum;
189 			}
190 		}
191 
192 		if (ap-am < gnuf || sum == 0) {
193 			if (trace)
194 				fprintf(fd,
195 					"%ld in %d passes;"
196 					" early out balance=%d\n",
197 				        (long)a, pass, npts);
198 			return a;	/* guess was good enough */
199 		}
200 
201 		aa = (sumx/sum-a)*AMP;
202 		if (npts >= 2) {	/* guess was too low */
203 			am = a;
204 			aa = xp + max(0.0, aa);
205 			if (aa >= ap)
206 				aa = (a + ap)/2;
207 
208 		} else if (npts <= -2) {    /* guess was two high */
209 			ap = a;
210 			aa = xm + min(0.0, aa);
211 			if (aa <= am)
212 				aa = (a + am)/2;
213 
214 		} else {
215 			break;		/* got it */
216 		}
217 
218 		if (a == aa) {
219 			if (trace)
220 				fprintf(fd, "%ld in %d passes;"
221 					" force out balance=%d\n",
222 				        (long)a, pass, npts);
223 			return a;
224 		}
225 		eps = AFAC*abs(aa - a);
226 		*eps_ptr = eps;
227 		a = aa;
228 	}
229 
230 	if (((x - xlim) % 2) != 0) {    /* even number of points? */
231 		if (npts == 0)		/* yes, return an average */
232 			a = (xp+xm)/2;
233 		else if (npts > 0)
234 			a =  (a+xp)/2;
235 		else
236 			a = (xm+a)/2;
237 
238 	} else if (npts != 0) {		/* odd number of points */
239 		if (npts > 0)
240 			a = xp;
241 		else
242 			a = xm;
243 	}
244 
245 	if (trace)
246 		fprintf(fd, "%ld in %d passes\n", (long)a, pass);
247 	return a;
248 }
249