PageRenderTime 57ms CodeModel.GetById 15ms app.highlight 35ms RepoModel.GetById 1ms app.codeStats 0ms

/contrib/ntp/clockstuff/chutest.c

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