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