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
32static char sccsid[] = "@(#)networkdelta.c	8.1 (Berkeley) 6/6/93";
33#endif
34static const char rcsid[] =
35  "$FreeBSD: stable/11/usr.sbin/timed/timed/networkdelta.c 360338 2020-04-26 15:50:32Z dim $";
36#endif /* not lint */
37
38#include "globals.h"
39
40static 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 */
64long
65networkdelta(void)
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 */
150static long
151median(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 = (float)LONG_MAX;		/* bounds on the median */
159	float am = -(float)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 = (float)LONG_MAX;
182		xm = -(float)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