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