source: trunk/third/xntp/clockstuff/chutest.c @ 10832

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