1/*	$NetBSD$	*/
2
3/* chutest.c,v 3.1 1993/07/06 01:05:21 jbj Exp
4 * chutest - test the CHU clock
5 */
6
7#include <stdio.h>
8#include <sys/types.h>
9#include <sys/socket.h>
10#include <netinet/in.h>
11#include <sys/ioctl.h>
12#include <sys/time.h>
13#include <sys/file.h>
14#include <sgtty.h>
15
16#include "../include/ntp_fp.h"
17#include "../include/ntp.h"
18#include "../include/ntp_unixtime.h"
19
20#ifdef CHULDISC
21#ifdef STREAM
22# ifdef HAVE_SYS_CHUDEFS_H
23#include <sys/chudefs.h>
24#endif
25#include <stropts.h>
26#endif
27#endif
28
29#ifdef CHULDISC
30# ifdef HAVE_SYS_CHUDEFS_H
31#include <sys/chudefs.h>
32#endif
33#endif
34
35#ifndef CHULDISC
36#ifndef STREAM
37#define	NCHUCHARS	(10)
38
39struct chucode {
40	u_char codechars[NCHUCHARS];	/* code characters */
41	u_char ncodechars;		/* number of code characters */
42	u_char chustatus;		/* not used currently */
43	struct timeval codetimes[NCHUCHARS];	/* arrival times */
44};
45#endif
46#endif
47
48#define	STREQ(a, b)	(*(a) == *(b) && strcmp((a), (b)) == 0)
49
50char *progname;
51int debug;
52
53int dofilter = 0;	/* set to 1 when we should run filter algorithm */
54int showtimes = 0;	/* set to 1 when we should show char arrival times */
55int doprocess = 0;	/* set to 1 when we do processing analogous to driver */
56#ifdef CHULDISC
57int usechuldisc = 0;	/* set to 1 when CHU line discipline should be used */
58#endif
59#ifdef STREAM
60int usechuldisc = 0;	/* set to 1 when CHU line discipline should be used */
61#endif
62
63struct timeval lasttv;
64struct chucode chudata;
65
66extern u_long ustotslo[];
67extern u_long ustotsmid[];
68extern u_long ustotshi[];
69
70/*
71 * main - parse arguments and handle options
72 */
73int
74main(
75	int argc,
76	char *argv[]
77	)
78{
79	int c;
80	int errflg = 0;
81	extern int ntp_optind;
82	extern char *ntp_optarg;
83	void init_chu();
84
85	progname = argv[0];
86	while ((c = ntp_getopt(argc, argv, "cdfpt")) != EOF)
87	    switch (c) {
88		case 'c':
89#ifdef STREAM
90		    usechuldisc = 1;
91		    break;
92#endif
93#ifdef CHULDISC
94		    usechuldisc = 1;
95		    break;
96#endif
97#ifndef STREAM
98#ifndef CHULDISC
99		    (void) fprintf(stderr,
100				   "%s: CHU line discipline not available on this machine\n",
101				   progname);
102		    exit(2);
103#endif
104#endif
105		case 'd':
106		    ++debug;
107		    break;
108		case 'f':
109		    dofilter = 1;
110		    break;
111		case 'p':
112		    doprocess = 1;
113		case 't':
114		    showtimes = 1;
115		    break;
116		default:
117		    errflg++;
118		    break;
119	    }
120	if (errflg || ntp_optind+1 != argc) {
121#ifdef STREAM
122		(void) fprintf(stderr, "usage: %s [-dft] tty_device\n",
123			       progname);
124#endif
125#ifdef CHULDISC
126		(void) fprintf(stderr, "usage: %s [-dft] tty_device\n",
127			       progname);
128#endif
129#ifndef STREAM
130#ifndef CHULDISC
131		(void) fprintf(stderr, "usage: %s [-cdft] tty_device\n",
132			       progname);
133#endif
134#endif
135		exit(2);
136	}
137
138	(void) gettimeofday(&lasttv, (struct timezone *)0);
139	c = openterm(argv[ntp_optind]);
140	init_chu();
141#ifdef STREAM
142	if (usechuldisc)
143	    process_ldisc(c);
144	else
145#endif
146#ifdef CHULDISC
147	    if (usechuldisc)
148		process_ldisc(c);
149	    else
150#endif
151		process_raw(c);
152	/*NOTREACHED*/
153}
154
155
156/*
157 * openterm - open a port to the CHU clock
158 */
159int
160openterm(
161	char *dev
162	)
163{
164	int s;
165	struct sgttyb ttyb;
166
167	if (debug)
168	    (void) fprintf(stderr, "Doing open...");
169	if ((s = open(dev, O_RDONLY, 0777)) < 0)
170	    error("open(%s)", dev, "");
171	if (debug)
172	    (void) fprintf(stderr, "open okay\n");
173
174	if (debug)
175	    (void) fprintf(stderr, "Setting exclusive use...");
176	if (ioctl(s, TIOCEXCL, (char *)0) < 0)
177	    error("ioctl(TIOCEXCL)", "", "");
178	if (debug)
179	    (void) fprintf(stderr, "done\n");
180
181	ttyb.sg_ispeed = ttyb.sg_ospeed = B300;
182	ttyb.sg_erase = ttyb.sg_kill = 0;
183	ttyb.sg_flags = EVENP|ODDP|RAW;
184	if (debug)
185	    (void) fprintf(stderr, "Setting baud rate et al...");
186	if (ioctl(s, TIOCSETP, (char *)&ttyb) < 0)
187	    error("ioctl(TIOCSETP, raw)", "", "");
188	if (debug)
189	    (void) fprintf(stderr, "done\n");
190
191#ifdef CHULDISC
192	if (usechuldisc) {
193		int ldisc;
194
195		if (debug)
196		    (void) fprintf(stderr, "Switching to CHU ldisc...");
197		ldisc = CHULDISC;
198		if (ioctl(s, TIOCSETD, (char *)&ldisc) < 0)
199		    error("ioctl(TIOCSETD, CHULDISC)", "", "");
200		if (debug)
201		    (void) fprintf(stderr, "okay\n");
202	}
203#endif
204#ifdef STREAM
205	if (usechuldisc) {
206
207		if (debug)
208		    (void) fprintf(stderr, "Poping off streams...");
209		while (ioctl(s, I_POP, 0) >=0) ;
210		if (debug)
211		    (void) fprintf(stderr, "okay\n");
212		if (debug)
213		    (void) fprintf(stderr, "Pushing CHU stream...");
214		if (ioctl(s, I_PUSH, "chu") < 0)
215		    error("ioctl(I_PUSH, \"chu\")", "", "");
216		if (debug)
217		    (void) fprintf(stderr, "okay\n");
218	}
219#endif
220	return s;
221}
222
223
224/*
225 * process_raw - process characters in raw mode
226 */
227int
228process_raw(
229	int s
230	)
231{
232	u_char c;
233	int n;
234	struct timeval tv;
235	struct timeval difftv;
236
237	while ((n = read(s, &c, sizeof(char))) > 0) {
238		(void) gettimeofday(&tv, (struct timezone *)0);
239		if (dofilter)
240		    raw_filter((unsigned int)c, &tv);
241		else {
242			difftv.tv_sec = tv.tv_sec - lasttv.tv_sec;
243			difftv.tv_usec = tv.tv_usec - lasttv.tv_usec;
244			if (difftv.tv_usec < 0) {
245				difftv.tv_sec--;
246				difftv.tv_usec += 1000000;
247			}
248			(void) printf("%02x\t%lu.%06lu\t%lu.%06lu\n",
249				      c, tv.tv_sec, tv.tv_usec, difftv.tv_sec,
250				      difftv.tv_usec);
251			lasttv = tv;
252		}
253	}
254
255	if (n == 0) {
256		(void) fprintf(stderr, "%s: zero returned on read\n", progname);
257		exit(1);
258	} else
259	    error("read()", "", "");
260}
261
262
263/*
264 * raw_filter - run the line discipline filter over raw data
265 */
266int
267raw_filter(
268	unsigned int c,
269	struct timeval *tv
270	)
271{
272	static struct timeval diffs[10] = { 0 };
273	struct timeval diff;
274	l_fp ts;
275	void chufilter();
276
277	if ((c & 0xf) > 9 || ((c>>4)&0xf) > 9) {
278		if (debug)
279		    (void) fprintf(stderr,
280				   "character %02x failed BCD test\n");
281		chudata.ncodechars = 0;
282		return;
283	}
284
285	if (chudata.ncodechars > 0) {
286		diff.tv_sec = tv->tv_sec
287			- chudata.codetimes[chudata.ncodechars].tv_sec;
288		diff.tv_usec = tv->tv_usec
289			- chudata.codetimes[chudata.ncodechars].tv_usec;
290		if (diff.tv_usec < 0) {
291			diff.tv_sec--;
292			diff.tv_usec += 1000000;
293		} /*
294		    if (diff.tv_sec != 0 || diff.tv_usec > 900000) {
295		    if (debug)
296		    (void) fprintf(stderr,
297		    "character %02x failed time test\n");
298		    chudata.ncodechars = 0;
299		    return;
300		    } */
301	}
302
303	chudata.codechars[chudata.ncodechars] = c;
304	chudata.codetimes[chudata.ncodechars] = *tv;
305	if (chudata.ncodechars > 0)
306	    diffs[chudata.ncodechars] = diff;
307	if (++chudata.ncodechars == 10) {
308		if (doprocess) {
309			TVTOTS(&chudata.codetimes[NCHUCHARS-1], &ts);
310			ts.l_ui += JAN_1970;
311			chufilter(&chudata, &chudata.codetimes[NCHUCHARS-1]);
312		} else {
313			register int i;
314
315			for (i = 0; i < chudata.ncodechars; i++) {
316				(void) printf("%x%x\t%lu.%06lu\t%lu.%06lu\n",
317					      chudata.codechars[i] & 0xf,
318					      (chudata.codechars[i] >>4 ) & 0xf,
319					      chudata.codetimes[i].tv_sec,
320					      chudata.codetimes[i].tv_usec,
321					      diffs[i].tv_sec, diffs[i].tv_usec);
322			}
323		}
324		chudata.ncodechars = 0;
325	}
326}
327
328
329/* #ifdef CHULDISC*/
330/*
331 * process_ldisc - process line discipline
332 */
333int
334process_ldisc(
335	int s
336	)
337{
338	struct chucode chu;
339	int n;
340	register int i;
341	struct timeval diff;
342	l_fp ts;
343	void chufilter();
344
345	while ((n = read(s, (char *)&chu, sizeof chu)) > 0) {
346		if (n != sizeof chu) {
347			(void) fprintf(stderr, "Expected %d, got %d\n",
348				       sizeof chu, n);
349			continue;
350		}
351
352		if (doprocess) {
353			TVTOTS(&chu.codetimes[NCHUCHARS-1], &ts);
354			ts.l_ui += JAN_1970;
355			chufilter(&chu, &ts);
356		} else {
357			for (i = 0; i < NCHUCHARS; i++) {
358				if (i == 0)
359				    diff.tv_sec = diff.tv_usec = 0;
360				else {
361					diff.tv_sec = chu.codetimes[i].tv_sec
362						- chu.codetimes[i-1].tv_sec;
363					diff.tv_usec = chu.codetimes[i].tv_usec
364						- chu.codetimes[i-1].tv_usec;
365					if (diff.tv_usec < 0) {
366						diff.tv_sec--;
367						diff.tv_usec += 1000000;
368					}
369				}
370				(void) printf("%x%x\t%lu.%06lu\t%lu.%06lu\n",
371					      chu.codechars[i] & 0xf, (chu.codechars[i]>>4)&0xf,
372					      chu.codetimes[i].tv_sec, chu.codetimes[i].tv_usec,
373					      diff.tv_sec, diff.tv_usec);
374			}
375		}
376	}
377	if (n == 0) {
378		(void) fprintf(stderr, "%s: zero returned on read\n", progname);
379		exit(1);
380	} else
381	    error("read()", "", "");
382}
383/*#endif*/
384
385
386/*
387 * error - print an error message
388 */
389void
390error(
391	char *fmt,
392	char *s1,
393	char *s2
394	)
395{
396	(void) fprintf(stderr, "%s: ", progname);
397	(void) fprintf(stderr, fmt, s1, s2);
398	(void) fprintf(stderr, ": ");
399	perror("");
400	exit(1);
401}
402
403/*
404 * Definitions
405 */
406#define	MAXUNITS	4	/* maximum number of CHU units permitted */
407#define	CHUDEV	"/dev/chu%d"	/* device we open.  %d is unit number */
408#define	NCHUCODES	9	/* expect 9 CHU codes per minute */
409
410/*
411 * When CHU is operating optimally we want the primary clock distance
412 * to come out at 300 ms.  Thus, peer.distance in the CHU peer structure
413 * is set to 290 ms and we compute delays which are at least 10 ms long.
414 * The following are 290 ms and 10 ms expressed in u_fp format
415 */
416#define	CHUDISTANCE	0x00004a3d
417#define	CHUBASEDELAY	0x0000028f
418
419/*
420 * To compute a quality for the estimate (a pseudo delay) we add a
421 * fixed 10 ms for each missing code in the minute and add to this
422 * the sum of the differences between the remaining offsets and the
423 * estimated sample offset.
424 */
425#define	CHUDELAYPENALTY	0x0000028f
426
427/*
428 * Other constant stuff
429 */
430#define	CHUPRECISION	(-9)		/* what the heck */
431#define	CHUREFID	"CHU\0"
432
433/*
434 * Default fudge factors
435 */
436#define	DEFPROPDELAY	0x00624dd3	/* 0.0015 seconds, 1.5 ms */
437#define	DEFFILTFUDGE	0x000d1b71	/* 0.0002 seconds, 200 us */
438
439/*
440 * Hacks to avoid excercising the multiplier.  I have no pride.
441 */
442#define	MULBY10(x)	(((x)<<3) + ((x)<<1))
443#define	MULBY60(x)	(((x)<<6) - ((x)<<2))	/* watch overflow */
444#define	MULBY24(x)	(((x)<<4) + ((x)<<3))
445
446/*
447 * Constants for use when multiplying by 0.1.  ZEROPTONE is 0.1
448 * as an l_fp fraction, NZPOBITS is the number of significant bits
449 * in ZEROPTONE.
450 */
451#define	ZEROPTONE	0x1999999a
452#define	NZPOBITS	29
453
454/*
455 * The CHU table.  This gives the expected time of arrival of each
456 * character after the on-time second and is computed as follows:
457 * The CHU time code is sent at 300 bps.  Your average UART will
458 * synchronize at the edge of the start bit and will consider the
459 * character complete at the center of the first stop bit, i.e.
460 * 0.031667 ms later.  Thus the expected time of each interrupt
461 * is the start bit time plus 0.031667 seconds.  These times are
462 * in chutable[].  To this we add such things as propagation delay
463 * and delay fudge factor.
464 */
465#define	CHARDELAY	0x081b4e80
466
467static u_long chutable[NCHUCHARS] = {
468	0x2147ae14 + CHARDELAY,		/* 0.130 (exactly) */
469	0x2ac08312 + CHARDELAY,		/* 0.167 (exactly) */
470	0x34395810 + CHARDELAY,		/* 0.204 (exactly) */
471	0x3db22d0e + CHARDELAY,		/* 0.241 (exactly) */
472	0x472b020c + CHARDELAY,		/* 0.278 (exactly) */
473	0x50a3d70a + CHARDELAY,		/* 0.315 (exactly) */
474	0x5a1cac08 + CHARDELAY,		/* 0.352 (exactly) */
475	0x63958106 + CHARDELAY,		/* 0.389 (exactly) */
476	0x6d0e5604 + CHARDELAY,		/* 0.426 (exactly) */
477	0x76872b02 + CHARDELAY,		/* 0.463 (exactly) */
478};
479
480/*
481 * Keep the fudge factors separately so they can be set even
482 * when no clock is configured.
483 */
484static l_fp propagation_delay;
485static l_fp fudgefactor;
486static l_fp offset_fudge;
487
488/*
489 * We keep track of the start of the year, watching for changes.
490 * We also keep track of whether the year is a leap year or not.
491 * All because stupid CHU doesn't include the year in the time code.
492 */
493static u_long yearstart;
494
495/*
496 * Imported from the timer module
497 */
498extern u_long current_time;
499extern struct event timerqueue[];
500
501/*
502 * Time conversion tables imported from the library
503 */
504extern u_long ustotslo[];
505extern u_long ustotsmid[];
506extern u_long ustotshi[];
507
508
509/*
510 * init_chu - initialize internal chu driver data
511 */
512void
513init_chu(void)
514{
515
516	/*
517	 * Initialize fudge factors to default.
518	 */
519	propagation_delay.l_ui = 0;
520	propagation_delay.l_uf = DEFPROPDELAY;
521	fudgefactor.l_ui = 0;
522	fudgefactor.l_uf = DEFFILTFUDGE;
523	offset_fudge = propagation_delay;
524	L_ADD(&offset_fudge, &fudgefactor);
525
526	yearstart = 0;
527}
528
529
530void
531chufilter(
532	struct chucode *chuc,
533	l_fp *rtime
534	)
535{
536	register int i;
537	register u_long date_ui;
538	register u_long tmp;
539	register u_char *code;
540	int isneg;
541	int imin;
542	int imax;
543	u_long reftime;
544	l_fp off[NCHUCHARS];
545	l_fp ts;
546	int day, hour, minute, second;
547	static u_char lastcode[NCHUCHARS];
548	extern u_long calyearstart();
549	extern char *mfptoa();
550	void chu_process();
551	extern char *prettydate();
552
553	/*
554	 * We'll skip the checks made in the kernel, but assume they've
555	 * been done.  This means that all characters are BCD and
556	 * the intercharacter spacing isn't unreasonable.
557	 */
558
559	/*
560	 * print the code
561	 */
562	for (i = 0; i < NCHUCHARS; i++)
563	    printf("%c%c", (chuc->codechars[i] & 0xf) + '0',
564		   ((chuc->codechars[i]>>4) & 0xf) + '0');
565	printf("\n");
566
567	/*
568	 * Format check.  Make sure the two halves match.
569	 */
570	for (i = 0; i < NCHUCHARS/2; i++)
571	    if (chuc->codechars[i] != chuc->codechars[i+(NCHUCHARS/2)]) {
572		    (void) printf("Bad format, halves don't match\n");
573		    return;
574	    }
575
576	/*
577	 * Break out the code into the BCD nibbles.  Only need to fiddle
578	 * with the first half since both are identical.  Note the first
579	 * BCD character is the low order nibble, the second the high order.
580	 */
581	code = lastcode;
582	for (i = 0; i < NCHUCHARS/2; i++) {
583		*code++ = chuc->codechars[i] & 0xf;
584		*code++ = (chuc->codechars[i] >> 4) & 0xf;
585	}
586
587	/*
588	 * If the first nibble isn't a 6, we're up the creek
589	 */
590	code = lastcode;
591	if (*code++ != 6) {
592		(void) printf("Bad format, no 6 at start\n");
593		return;
594	}
595
596	/*
597	 * Collect the day, the hour, the minute and the second.
598	 */
599	day = *code++;
600	day = MULBY10(day) + *code++;
601	day = MULBY10(day) + *code++;
602	hour = *code++;
603	hour = MULBY10(hour) + *code++;
604	minute = *code++;
605	minute = MULBY10(minute) + *code++;
606	second = *code++;
607	second = MULBY10(second) + *code++;
608
609	/*
610	 * Sanity check the day and time.  Note that this
611	 * only occurs on the 31st through the 39th second
612	 * of the minute.
613	 */
614	if (day < 1 || day > 366
615	    || hour > 23 || minute > 59
616	    || second < 31 || second > 39) {
617		(void) printf("Failed date sanity check: %d %d %d %d\n",
618			      day, hour, minute, second);
619		return;
620	}
621
622	/*
623	 * Compute seconds into the year.
624	 */
625	tmp = (u_long)(MULBY24((day-1)) + hour);	/* hours */
626	tmp = MULBY60(tmp) + (u_long)minute;		/* minutes */
627	tmp = MULBY60(tmp) + (u_long)second;		/* seconds */
628
629	/*
630	 * Now the fun begins.  We demand that the received time code
631	 * be within CLOCK_WAYTOOBIG of the receive timestamp, but
632	 * there is uncertainty about the year the timestamp is in.
633	 * Use the current year start for the first check, this should
634	 * work most of the time.
635	 */
636	date_ui = tmp + yearstart;
637	if (date_ui < (rtime->l_ui + CLOCK_WAYTOOBIG)
638	    && date_ui > (rtime->l_ui - CLOCK_WAYTOOBIG))
639	    goto codeokay;	/* looks good */
640
641	/*
642	 * Trouble.  Next check is to see if the year rolled over and, if
643	 * so, try again with the new year's start.
644	 */
645	date_ui = calyearstart(rtime->l_ui);
646	if (date_ui != yearstart) {
647		yearstart = date_ui;
648		date_ui += tmp;
649		(void) printf("time %u, code %u, difference %d\n",
650			      date_ui, rtime->l_ui, (long)date_ui-(long)rtime->l_ui);
651		if (date_ui < (rtime->l_ui + CLOCK_WAYTOOBIG)
652		    && date_ui > (rtime->l_ui - CLOCK_WAYTOOBIG))
653		    goto codeokay;	/* okay this time */
654	}
655
656	ts.l_uf = 0;
657	ts.l_ui = yearstart;
658	printf("yearstart %s\n", prettydate(&ts));
659	printf("received %s\n", prettydate(rtime));
660	ts.l_ui = date_ui;
661	printf("date_ui %s\n", prettydate(&ts));
662
663	/*
664	 * Here we know the year start matches the current system
665	 * time.  One remaining possibility is that the time code
666	 * is in the year previous to that of the system time.  This
667	 * is only worth checking if the receive timestamp is less
668	 * than CLOCK_WAYTOOBIG seconds into the new year.
669	 */
670	if ((rtime->l_ui - yearstart) < CLOCK_WAYTOOBIG) {
671		date_ui = tmp + calyearstart(yearstart - CLOCK_WAYTOOBIG);
672		if ((rtime->l_ui - date_ui) < CLOCK_WAYTOOBIG)
673		    goto codeokay;
674	}
675
676	/*
677	 * One last possibility is that the time stamp is in the year
678	 * following the year the system is in.  Try this one before
679	 * giving up.
680	 */
681	date_ui = tmp + calyearstart(yearstart + (400*24*60*60)); /* 400 days */
682	if ((date_ui - rtime->l_ui) >= CLOCK_WAYTOOBIG) {
683		printf("Date hopelessly off\n");
684		return;		/* hopeless, let it sync to other peers */
685	}
686
687    codeokay:
688	reftime = date_ui;
689	/*
690	 * We've now got the integral seconds part of the time code (we hope).
691	 * The fractional part comes from the table.  We next compute
692	 * the offsets for each character.
693	 */
694	for (i = 0; i < NCHUCHARS; i++) {
695		register u_long tmp2;
696
697		off[i].l_ui = date_ui;
698		off[i].l_uf = chutable[i];
699		tmp = chuc->codetimes[i].tv_sec + JAN_1970;
700		TVUTOTSF(chuc->codetimes[i].tv_usec, tmp2);
701		M_SUB(off[i].l_ui, off[i].l_uf, tmp, tmp2);
702	}
703
704	/*
705	 * Here is a *big* problem.  What one would normally
706	 * do here on a machine with lots of clock bits (say
707	 * a Vax or the gizmo board) is pick the most positive
708	 * offset and the estimate, since this is the one that
709	 * is most likely suffered the smallest interrupt delay.
710	 * The trouble is that the low order clock bit on an IBM
711	 * RT, which is the machine I had in mind when doing this,
712	 * ticks at just under the millisecond mark.  This isn't
713	 * precise enough.  What we can do to improve this is to
714	 * average all 10 samples and rely on the second level
715	 * filtering to pick the least delayed estimate.  Trouble
716	 * is, this means we have to divide a 64 bit fixed point
717	 * number by 10, a procedure which really sucks.  Oh, well.
718	 * First compute the sum.
719	 */
720	date_ui = 0;
721	tmp = 0;
722	for (i = 0; i < NCHUCHARS; i++)
723	    M_ADD(date_ui, tmp, off[i].l_ui, off[i].l_uf);
724	if (M_ISNEG(date_ui, tmp))
725	    isneg = 1;
726	else
727	    isneg = 0;
728
729	/*
730	 * Here is a multiply-by-0.1 optimization that should apply
731	 * just about everywhere.  If the magnitude of the sum
732	 * is less than 9 we don't have to worry about overflow
733	 * out of a 64 bit product, even after rounding.
734	 */
735	if (date_ui < 9 || date_ui > 0xfffffff7) {
736		register u_long prod_ui;
737		register u_long prod_uf;
738
739		prod_ui = prod_uf = 0;
740		/*
741		 * This code knows the low order bit in 0.1 is zero
742		 */
743		for (i = 1; i < NZPOBITS; i++) {
744			M_LSHIFT(date_ui, tmp);
745			if (ZEROPTONE & (1<<i))
746			    M_ADD(prod_ui, prod_uf, date_ui, tmp);
747		}
748
749		/*
750		 * Done, round it correctly.  Prod_ui contains the
751		 * fraction.
752		 */
753		if (prod_uf & 0x80000000)
754		    prod_ui++;
755		if (isneg)
756		    date_ui = 0xffffffff;
757		else
758		    date_ui = 0;
759		tmp = prod_ui;
760		/*
761		 * date_ui is integral part, tmp is fraction.
762		 */
763	} else {
764		register u_long prod_ovr;
765		register u_long prod_ui;
766		register u_long prod_uf;
767		register u_long highbits;
768
769		prod_ovr = prod_ui = prod_uf = 0;
770		if (isneg)
771		    highbits = 0xffffffff;	/* sign extend */
772		else
773		    highbits = 0;
774		/*
775		 * This code knows the low order bit in 0.1 is zero
776		 */
777		for (i = 1; i < NZPOBITS; i++) {
778			M_LSHIFT3(highbits, date_ui, tmp);
779			if (ZEROPTONE & (1<<i))
780			    M_ADD3(prod_ovr, prod_uf, prod_ui,
781				   highbits, date_ui, tmp);
782		}
783
784		if (prod_uf & 0x80000000)
785		    M_ADDUF(prod_ovr, prod_ui, (u_long)1);
786		date_ui = prod_ovr;
787		tmp = prod_ui;
788	}
789
790	/*
791	 * At this point we have the mean offset, with the integral
792	 * part in date_ui and the fractional part in tmp.  Store
793	 * it in the structure.
794	 */
795	/*
796	 * Add in fudge factor.
797	 */
798	M_ADD(date_ui, tmp, offset_fudge.l_ui, offset_fudge.l_uf);
799
800	/*
801	 * Find the minimun and maximum offset
802	 */
803	imin = imax = 0;
804	for (i = 1; i < NCHUCHARS; i++) {
805		if (L_ISGEQ(&off[i], &off[imax])) {
806			imax = i;
807		} else if (L_ISGEQ(&off[imin], &off[i])) {
808			imin = i;
809		}
810	}
811
812	L_ADD(&off[imin], &offset_fudge);
813	if (imin != imax)
814	    L_ADD(&off[imax], &offset_fudge);
815	(void) printf("mean %s, min %s, max %s\n",
816		      mfptoa(date_ui, tmp, 8), lfptoa(&off[imin], 8),
817		      lfptoa(&off[imax], 8));
818}
819