1 /***************************************************************************
2 * genutils.c
3 *
4 * Generic utility routines
5 *
6 * Written by Chad Trabant
7 * ORFEUS/EC-Project MEREDIAN
8 * IRIS Data Management Center
9 *
10 * modified: 2017.053
11 ***************************************************************************/
12
13 #include <errno.h>
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include <time.h>
18
19 #include "libmseed.h"
20
21 static hptime_t ms_time2hptime_int (int year, int day, int hour,
22 int min, int sec, int usec);
23
24 static struct tm *ms_gmtime_r (int64_t *timep, struct tm *result);
25
26 /* A constant number of seconds between the NTP and Posix/Unix time epoch */
27 #define NTPPOSIXEPOCHDELTA 2208988800LL
28
29 /* Global variable to hold a leap second list */
30 LeapSecond *leapsecondlist = NULL;
31
32 /***************************************************************************
33 * ms_recsrcname:
34 *
35 * Generate a source name string for a specified raw data record in
36 * the format: 'NET_STA_LOC_CHAN' or, if the quality flag is true:
37 * 'NET_STA_LOC_CHAN_QUAL'. The passed srcname must have enough room
38 * for the resulting string.
39 *
40 * Returns a pointer to the resulting string or NULL on error.
41 ***************************************************************************/
42 char *
ms_recsrcname(char * record,char * srcname,flag quality)43 ms_recsrcname (char *record, char *srcname, flag quality)
44 {
45 struct fsdh_s *fsdh;
46 char network[6];
47 char station[6];
48 char location[6];
49 char channel[6];
50
51 if (!record)
52 return NULL;
53
54 fsdh = (struct fsdh_s *)record;
55
56 ms_strncpclean (network, fsdh->network, 2);
57 ms_strncpclean (station, fsdh->station, 5);
58 ms_strncpclean (location, fsdh->location, 2);
59 ms_strncpclean (channel, fsdh->channel, 3);
60
61 /* Build the source name string including the quality indicator*/
62 if (quality)
63 sprintf (srcname, "%s_%s_%s_%s_%c",
64 network, station, location, channel, fsdh->dataquality);
65
66 /* Build the source name string without the quality indicator*/
67 else
68 sprintf (srcname, "%s_%s_%s_%s", network, station, location, channel);
69
70 return srcname;
71 } /* End of ms_recsrcname() */
72
73 /***************************************************************************
74 * ms_splitsrcname:
75 *
76 * Split srcname into separate components: "NET_STA_LOC_CHAN[_QUAL]".
77 * Memory for each component must already be allocated. If a specific
78 * component is not desired set the appropriate argument to NULL.
79 *
80 * Returns 0 on success and -1 on error.
81 ***************************************************************************/
82 int
ms_splitsrcname(char * srcname,char * net,char * sta,char * loc,char * chan,char * qual)83 ms_splitsrcname (char *srcname, char *net, char *sta, char *loc, char *chan,
84 char *qual)
85 {
86 char *id;
87 char *ptr, *top, *next;
88 int sepcnt = 0;
89
90 if (!srcname)
91 return -1;
92
93 /* Verify number of separating underscore characters */
94 id = srcname;
95 while ((id = strchr (id, '_')))
96 {
97 id++;
98 sepcnt++;
99 }
100
101 /* Either 3 or 4 separating underscores are required */
102 if (sepcnt != 3 && sepcnt != 4)
103 {
104 return -1;
105 }
106
107 /* Duplicate srcname */
108 if (!(id = strdup (srcname)))
109 {
110 ms_log (2, "ms_splitsrcname(): Error duplicating srcname string");
111 return -1;
112 }
113
114 /* Network */
115 top = id;
116 if ((ptr = strchr (top, '_')))
117 {
118 next = ptr + 1;
119 *ptr = '\0';
120
121 if (net)
122 strcpy (net, top);
123
124 top = next;
125 }
126 /* Station */
127 if ((ptr = strchr (top, '_')))
128 {
129 next = ptr + 1;
130 *ptr = '\0';
131
132 if (sta)
133 strcpy (sta, top);
134
135 top = next;
136 }
137 /* Location */
138 if ((ptr = strchr (top, '_')))
139 {
140 next = ptr + 1;
141 *ptr = '\0';
142
143 if (loc)
144 strcpy (loc, top);
145
146 top = next;
147 }
148 /* Channel & optional Quality */
149 if ((ptr = strchr (top, '_')))
150 {
151 next = ptr + 1;
152 *ptr = '\0';
153
154 if (chan)
155 strcpy (chan, top);
156
157 top = next;
158
159 /* Quality */
160 if (*top && qual)
161 {
162 /* Quality is a single character */
163 *qual = *top;
164 }
165 }
166 /* Otherwise only Channel */
167 else if (*top && chan)
168 {
169 strcpy (chan, top);
170 }
171
172 /* Free duplicated stream ID */
173 if (id)
174 free (id);
175
176 return 0;
177 } /* End of ms_splitsrcname() */
178
179 /***************************************************************************
180 * ms_strncpclean:
181 *
182 * Copy up to 'length' characters from 'source' to 'dest' while
183 * removing all spaces. The result is left justified and always null
184 * terminated. The destination string must have enough room needed
185 * for the non-space characters within 'length' and the null
186 * terminator, a maximum of 'length + 1'.
187 *
188 * Returns the number of characters (not including the null terminator) in
189 * the destination string.
190 ***************************************************************************/
191 int
ms_strncpclean(char * dest,const char * source,int length)192 ms_strncpclean (char *dest, const char *source, int length)
193 {
194 int sidx, didx;
195
196 if (!dest)
197 return 0;
198
199 if (!source)
200 {
201 *dest = '\0';
202 return 0;
203 }
204
205 for (sidx = 0, didx = 0; sidx < length; sidx++)
206 {
207 if (*(source + sidx) == '\0')
208 {
209 break;
210 }
211
212 if (*(source + sidx) != ' ')
213 {
214 *(dest + didx) = *(source + sidx);
215 didx++;
216 }
217 }
218
219 *(dest + didx) = '\0';
220
221 return didx;
222 } /* End of ms_strncpclean() */
223
224 /***************************************************************************
225 * ms_strncpcleantail:
226 *
227 * Copy up to 'length' characters from 'source' to 'dest' without any
228 * trailing spaces. The result is left justified and always null
229 * terminated. The destination string must have enough room needed
230 * for the characters within 'length' and the null terminator, a
231 * maximum of 'length + 1'.
232 *
233 * Returns the number of characters (not including the null terminator) in
234 * the destination string.
235 ***************************************************************************/
236 int
ms_strncpcleantail(char * dest,const char * source,int length)237 ms_strncpcleantail (char *dest, const char *source, int length)
238 {
239 int idx, pretail;
240
241 if (!dest)
242 return 0;
243
244 if (!source)
245 {
246 *dest = '\0';
247 return 0;
248 }
249
250 *(dest + length) = '\0';
251
252 pretail = 0;
253 for (idx = length - 1; idx >= 0; idx--)
254 {
255 if (!pretail && *(source + idx) == ' ')
256 {
257 *(dest + idx) = '\0';
258 }
259 else
260 {
261 pretail++;
262 *(dest + idx) = *(source + idx);
263 }
264 }
265
266 return pretail;
267 } /* End of ms_strncpcleantail() */
268
269 /***************************************************************************
270 * ms_strncpopen:
271 *
272 * Copy 'length' characters from 'source' to 'dest', padding the right
273 * side with spaces and leave open-ended. The result is left
274 * justified and *never* null terminated (the open-ended part). The
275 * destination string must have enough room for 'length' characters.
276 *
277 * Returns the number of characters copied from the source string.
278 ***************************************************************************/
279 int
ms_strncpopen(char * dest,const char * source,int length)280 ms_strncpopen (char *dest, const char *source, int length)
281 {
282 int didx;
283 int dcnt = 0;
284 int term = 0;
285
286 if (!dest)
287 return 0;
288
289 if (!source)
290 {
291 for (didx = 0; didx < length; didx++)
292 {
293 *(dest + didx) = ' ';
294 }
295
296 return 0;
297 }
298
299 for (didx = 0; didx < length; didx++)
300 {
301 if (!term)
302 if (*(source + didx) == '\0')
303 term = 1;
304
305 if (!term)
306 {
307 *(dest + didx) = *(source + didx);
308 dcnt++;
309 }
310 else
311 {
312 *(dest + didx) = ' ';
313 }
314 }
315
316 return dcnt;
317 } /* End of ms_strncpopen() */
318
319 /***************************************************************************
320 * ms_doy2md:
321 *
322 * Compute the month and day-of-month from a year and day-of-year.
323 *
324 * Year is expected to be in the range 1800-5000, jday is expected to
325 * be in the range 1-366, month will be in the range 1-12 and mday
326 * will be in the range 1-31.
327 *
328 * Returns 0 on success and -1 on error.
329 ***************************************************************************/
330 int
ms_doy2md(int year,int jday,int * month,int * mday)331 ms_doy2md (int year, int jday, int *month, int *mday)
332 {
333 int idx;
334 int leap;
335 int days[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
336
337 /* Sanity check for the supplied year */
338 if (year < 1800 || year > 5000)
339 {
340 ms_log (2, "ms_doy2md(): year (%d) is out of range\n", year);
341 return -1;
342 }
343
344 /* Test for leap year */
345 leap = (((year % 4 == 0) && (year % 100 != 0)) || (year % 400 == 0)) ? 1 : 0;
346
347 /* Add a day to February if leap year */
348 if (leap)
349 days[1]++;
350
351 if (jday > 365 + leap || jday <= 0)
352 {
353 ms_log (2, "ms_doy2md(): day-of-year (%d) is out of range\n", jday);
354 return -1;
355 }
356
357 for (idx = 0; idx < 12; idx++)
358 {
359 jday -= days[idx];
360
361 if (jday <= 0)
362 {
363 *month = idx + 1;
364 *mday = days[idx] + jday;
365 break;
366 }
367 }
368
369 return 0;
370 } /* End of ms_doy2md() */
371
372 /***************************************************************************
373 * ms_md2doy:
374 *
375 * Compute the day-of-year from a year, month and day-of-month.
376 *
377 * Year is expected to be in the range 1800-5000, month is expected to
378 * be in the range 1-12, mday is expected to be in the range 1-31 and
379 * jday will be in the range 1-366.
380 *
381 * Returns 0 on success and -1 on error.
382 ***************************************************************************/
383 int
ms_md2doy(int year,int month,int mday,int * jday)384 ms_md2doy (int year, int month, int mday, int *jday)
385 {
386 int idx;
387 int leap;
388 int days[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
389
390 /* Sanity check for the supplied parameters */
391 if (year < 1800 || year > 5000)
392 {
393 ms_log (2, "ms_md2doy(): year (%d) is out of range\n", year);
394 return -1;
395 }
396 if (month < 1 || month > 12)
397 {
398 ms_log (2, "ms_md2doy(): month (%d) is out of range\n", month);
399 return -1;
400 }
401 if (mday < 1 || mday > 31)
402 {
403 ms_log (2, "ms_md2doy(): day-of-month (%d) is out of range\n", mday);
404 return -1;
405 }
406
407 /* Test for leap year */
408 leap = (((year % 4 == 0) && (year % 100 != 0)) || (year % 400 == 0)) ? 1 : 0;
409
410 /* Add a day to February if leap year */
411 if (leap)
412 days[1]++;
413
414 /* Check that the day-of-month jives with specified month */
415 if (mday > days[month - 1])
416 {
417 ms_log (2, "ms_md2doy(): day-of-month (%d) is out of range for month %d\n",
418 mday, month);
419 return -1;
420 }
421
422 *jday = 0;
423 month--;
424
425 for (idx = 0; idx < 12; idx++)
426 {
427 if (idx == month)
428 {
429 *jday += mday;
430 break;
431 }
432
433 *jday += days[idx];
434 }
435
436 return 0;
437 } /* End of ms_md2doy() */
438
439 /***************************************************************************
440 * ms_btime2hptime:
441 *
442 * Convert a binary SEED time structure to a high precision epoch time
443 * (1/HPTMODULUS second ticks from the epoch). The algorithm used is
444 * a specific version of a generalized function in GNU glibc.
445 *
446 * Returns a high precision epoch time on success and HPTERROR on
447 * error.
448 ***************************************************************************/
449 hptime_t
ms_btime2hptime(BTime * btime)450 ms_btime2hptime (BTime *btime)
451 {
452 hptime_t hptime;
453 int shortyear;
454 int a4, a100, a400;
455 int intervening_leap_days;
456 int days;
457
458 if (!btime)
459 return HPTERROR;
460
461 shortyear = btime->year - 1900;
462
463 a4 = (shortyear >> 2) + 475 - !(shortyear & 3);
464 a100 = a4 / 25 - (a4 % 25 < 0);
465 a400 = a100 >> 2;
466 intervening_leap_days = (a4 - 492) - (a100 - 19) + (a400 - 4);
467
468 days = (365 * (shortyear - 70) + intervening_leap_days + (btime->day - 1));
469
470 hptime = (hptime_t) (60 * (60 * ((hptime_t)24 * days + btime->hour) + btime->min) + btime->sec) * HPTMODULUS + (btime->fract * (HPTMODULUS / 10000));
471
472 return hptime;
473 } /* End of ms_btime2hptime() */
474
475 /***************************************************************************
476 * ms_btime2isotimestr:
477 *
478 * Build a time string in ISO recommended format from a BTime struct.
479 *
480 * The provided isostimestr must have enough room for the resulting time
481 * string of 25 characters, i.e. '2001-07-29T12:38:00.0000' + NULL.
482 *
483 * Returns a pointer to the resulting string or NULL on error.
484 ***************************************************************************/
485 char *
ms_btime2isotimestr(BTime * btime,char * isotimestr)486 ms_btime2isotimestr (BTime *btime, char *isotimestr)
487 {
488 int month = 0;
489 int mday = 0;
490 int ret;
491
492 if (!isotimestr)
493 return NULL;
494
495 if (ms_doy2md (btime->year, btime->day, &month, &mday))
496 {
497 ms_log (2, "ms_btime2isotimestr(): Error converting year %d day %d\n",
498 btime->year, btime->day);
499 return NULL;
500 }
501
502 ret = snprintf (isotimestr, 25, "%4d-%02d-%02dT%02d:%02d:%02d.%04d",
503 btime->year, month, mday,
504 btime->hour, btime->min, btime->sec, btime->fract);
505
506 if (ret != 24)
507 return NULL;
508 else
509 return isotimestr;
510 } /* End of ms_btime2isotimestr() */
511
512 /***************************************************************************
513 * ms_btime2mdtimestr:
514 *
515 * Build a time string in month-day format from a BTime struct.
516 *
517 * The provided isostimestr must have enough room for the resulting time
518 * string of 25 characters, i.e. '2001-07-29 12:38:00.0000' + NULL.
519 *
520 * Returns a pointer to the resulting string or NULL on error.
521 ***************************************************************************/
522 char *
ms_btime2mdtimestr(BTime * btime,char * mdtimestr)523 ms_btime2mdtimestr (BTime *btime, char *mdtimestr)
524 {
525 int month = 0;
526 int mday = 0;
527 int ret;
528
529 if (!mdtimestr)
530 return NULL;
531
532 if (ms_doy2md (btime->year, btime->day, &month, &mday))
533 {
534 ms_log (2, "ms_btime2mdtimestr(): Error converting year %d day %d\n",
535 btime->year, btime->day);
536 return NULL;
537 }
538
539 ret = snprintf (mdtimestr, 25, "%4d-%02d-%02d %02d:%02d:%02d.%04d",
540 btime->year, month, mday,
541 btime->hour, btime->min, btime->sec, btime->fract);
542
543 if (ret != 24)
544 return NULL;
545 else
546 return mdtimestr;
547 } /* End of ms_btime2mdtimestr() */
548
549 /***************************************************************************
550 * ms_btime2seedtimestr:
551 *
552 * Build a SEED time string from a BTime struct.
553 *
554 * The provided seedtimestr must have enough room for the resulting time
555 * string of 23 characters, i.e. '2001,195,12:38:00.0000' + NULL.
556 *
557 * Returns a pointer to the resulting string or NULL on error.
558 ***************************************************************************/
559 char *
ms_btime2seedtimestr(BTime * btime,char * seedtimestr)560 ms_btime2seedtimestr (BTime *btime, char *seedtimestr)
561 {
562 int ret;
563
564 if (!seedtimestr)
565 return NULL;
566
567 ret = snprintf (seedtimestr, 23, "%4d,%03d,%02d:%02d:%02d.%04d",
568 btime->year, btime->day,
569 btime->hour, btime->min, btime->sec, btime->fract);
570
571 if (ret != 22)
572 return NULL;
573 else
574 return seedtimestr;
575 } /* End of ms_btime2seedtimestr() */
576
577 /***************************************************************************
578 * ms_hptime2tomsusecoffset:
579 *
580 * Convert a high precision epoch time to a time value in tenths of
581 * milliseconds (aka toms) and a microsecond offset (aka usecoffset).
582 *
583 * The tenths of milliseconds value will be rounded to the nearest
584 * value having a microsecond offset value between -50 to +49.
585 *
586 * Returns 0 on success and -1 on error.
587 ***************************************************************************/
588 int
ms_hptime2tomsusecoffset(hptime_t hptime,hptime_t * toms,int8_t * usecoffset)589 ms_hptime2tomsusecoffset (hptime_t hptime, hptime_t *toms, int8_t *usecoffset)
590 {
591 if (toms == NULL || usecoffset == NULL)
592 return -1;
593
594 /* Split time into tenths of milliseconds and microseconds */
595 *toms = hptime / (HPTMODULUS / 10000);
596 *usecoffset = (int8_t) (hptime - (*toms * (HPTMODULUS / 10000)));
597
598 /* Round tenths and adjust microsecond offset to -50 to +49 range */
599 if (*usecoffset > 49 && *usecoffset < 100)
600 {
601 *toms += 1;
602 *usecoffset -= 100;
603 }
604 else if (*usecoffset < -50 && *usecoffset > -100)
605 {
606 *toms -= 1;
607 *usecoffset += 100;
608 }
609
610 /* Convert tenths of milliseconds to be in hptime_t (HPTMODULUS) units */
611 *toms *= (HPTMODULUS / 10000);
612
613 return 0;
614 } /* End of ms_hptime2tomsusecoffset() */
615
616 /***************************************************************************
617 * ms_hptime2btime:
618 *
619 * Convert a high precision epoch time to a SEED binary time
620 * structure. The microseconds beyond the 1/10000 second range are
621 * truncated and *not* rounded, this is intentional and necessary.
622 *
623 * Returns 0 on success and -1 on error.
624 ***************************************************************************/
625 int
ms_hptime2btime(hptime_t hptime,BTime * btime)626 ms_hptime2btime (hptime_t hptime, BTime *btime)
627 {
628 struct tm tms;
629 int64_t isec;
630 int ifract;
631 int bfract;
632
633 if (btime == NULL)
634 return -1;
635
636 /* Reduce to Unix/POSIX epoch time and fractional seconds */
637 isec = MS_HPTIME2EPOCH (hptime);
638 ifract = (int)(hptime - (isec * HPTMODULUS));
639
640 /* BTime only has 1/10000 second precision */
641 bfract = ifract / (HPTMODULUS / 10000);
642
643 /* Adjust for negative epoch times, round back when needed */
644 if (hptime < 0 && ifract != 0)
645 {
646 /* Isolate microseconds between 1e-4 and 1e-6 precision and adjust bfract if not zero */
647 if (ifract - bfract * (HPTMODULUS / 10000))
648 bfract -= 1;
649
650 isec -= 1;
651 bfract = 10000 - (-bfract);
652 }
653
654 if (!(ms_gmtime_r (&isec, &tms)))
655 return -1;
656
657 btime->year = tms.tm_year + 1900;
658 btime->day = tms.tm_yday + 1;
659 btime->hour = tms.tm_hour;
660 btime->min = tms.tm_min;
661 btime->sec = tms.tm_sec;
662 btime->unused = 0;
663 btime->fract = (uint16_t)bfract;
664
665 return 0;
666 } /* End of ms_hptime2btime() */
667
668 /***************************************************************************
669 * ms_hptime2isotimestr:
670 *
671 * Build a time string in ISO recommended format from a high precision
672 * epoch time.
673 *
674 * The provided isostimestr must have enough room for the resulting time
675 * string of 27 characters, i.e. '2001-07-29T12:38:00.000000' + NULL.
676 *
677 * The 'subseconds' flag controls whenther the sub second portion of the
678 * time is included or not.
679 *
680 * Returns a pointer to the resulting string or NULL on error.
681 ***************************************************************************/
682 char *
ms_hptime2isotimestr(hptime_t hptime,char * isotimestr,flag subseconds)683 ms_hptime2isotimestr (hptime_t hptime, char *isotimestr, flag subseconds)
684 {
685 struct tm tms;
686 int64_t isec;
687 int ifract;
688 int ret;
689
690 if (isotimestr == NULL)
691 return NULL;
692
693 /* Reduce to Unix/POSIX epoch time and fractional seconds */
694 isec = MS_HPTIME2EPOCH (hptime);
695 ifract = (int)(hptime - (isec * HPTMODULUS));
696
697 /* Adjust for negative epoch times */
698 if (hptime < 0 && ifract != 0)
699 {
700 isec -= 1;
701 ifract = HPTMODULUS - (-ifract);
702 }
703
704 if (!(ms_gmtime_r (&isec, &tms)))
705 return NULL;
706
707 if (subseconds)
708 /* Assuming ifract has at least microsecond precision */
709 ret = snprintf (isotimestr, 27, "%4d-%02d-%02dT%02d:%02d:%02d.%06d",
710 tms.tm_year + 1900, tms.tm_mon + 1, tms.tm_mday,
711 tms.tm_hour, tms.tm_min, tms.tm_sec, ifract);
712 else
713 ret = snprintf (isotimestr, 20, "%4d-%02d-%02dT%02d:%02d:%02d",
714 tms.tm_year + 1900, tms.tm_mon + 1, tms.tm_mday,
715 tms.tm_hour, tms.tm_min, tms.tm_sec);
716
717 if (ret != 26 && ret != 19)
718 return NULL;
719 else
720 return isotimestr;
721 } /* End of ms_hptime2isotimestr() */
722
723 /***************************************************************************
724 * ms_hptime2mdtimestr:
725 *
726 * Build a time string in month-day format from a high precision
727 * epoch time.
728 *
729 * The provided mdtimestr must have enough room for the resulting time
730 * string of 27 characters, i.e. '2001-07-29 12:38:00.000000' + NULL.
731 *
732 * The 'subseconds' flag controls whenther the sub second portion of the
733 * time is included or not.
734 *
735 * Returns a pointer to the resulting string or NULL on error.
736 ***************************************************************************/
737 char *
ms_hptime2mdtimestr(hptime_t hptime,char * mdtimestr,flag subseconds)738 ms_hptime2mdtimestr (hptime_t hptime, char *mdtimestr, flag subseconds)
739 {
740 struct tm tms;
741 int64_t isec;
742 int ifract;
743 int ret;
744
745 if (mdtimestr == NULL)
746 return NULL;
747
748 /* Reduce to Unix/POSIX epoch time and fractional seconds */
749 isec = MS_HPTIME2EPOCH (hptime);
750 ifract = (int)(hptime - (isec * HPTMODULUS));
751
752 /* Adjust for negative epoch times */
753 if (hptime < 0 && ifract != 0)
754 {
755 isec -= 1;
756 ifract = HPTMODULUS - (-ifract);
757 }
758
759 if (!(ms_gmtime_r (&isec, &tms)))
760 return NULL;
761
762 if (subseconds)
763 /* Assuming ifract has at least microsecond precision */
764 ret = snprintf (mdtimestr, 27, "%4d-%02d-%02d %02d:%02d:%02d.%06d",
765 tms.tm_year + 1900, tms.tm_mon + 1, tms.tm_mday,
766 tms.tm_hour, tms.tm_min, tms.tm_sec, ifract);
767 else
768 ret = snprintf (mdtimestr, 20, "%4d-%02d-%02d %02d:%02d:%02d",
769 tms.tm_year + 1900, tms.tm_mon + 1, tms.tm_mday,
770 tms.tm_hour, tms.tm_min, tms.tm_sec);
771
772 if (ret != 26 && ret != 19)
773 return NULL;
774 else
775 return mdtimestr;
776 } /* End of ms_hptime2mdtimestr() */
777
778 /***************************************************************************
779 * ms_hptime2seedtimestr:
780 *
781 * Build a SEED time string from a high precision epoch time.
782 *
783 * The provided seedtimestr must have enough room for the resulting time
784 * string of 25 characters, i.e. '2001,195,12:38:00.000000\n'.
785 *
786 * The 'subseconds' flag controls whenther the sub second portion of the
787 * time is included or not.
788 *
789 * Returns a pointer to the resulting string or NULL on error.
790 ***************************************************************************/
791 char *
ms_hptime2seedtimestr(hptime_t hptime,char * seedtimestr,flag subseconds)792 ms_hptime2seedtimestr (hptime_t hptime, char *seedtimestr, flag subseconds)
793 {
794 struct tm tms;
795 int64_t isec;
796 int ifract;
797 int ret;
798
799 if (seedtimestr == NULL)
800 return NULL;
801
802 /* Reduce to Unix/POSIX epoch time and fractional seconds */
803 isec = MS_HPTIME2EPOCH (hptime);
804 ifract = (int)(hptime - (isec * HPTMODULUS));
805
806 /* Adjust for negative epoch times */
807 if (hptime < 0 && ifract != 0)
808 {
809 isec -= 1;
810 ifract = HPTMODULUS - (-ifract);
811 }
812
813 if (!(ms_gmtime_r (&isec, &tms)))
814 return NULL;
815
816 if (subseconds)
817 /* Assuming ifract has at least microsecond precision */
818 ret = snprintf (seedtimestr, 25, "%4d,%03d,%02d:%02d:%02d.%06d",
819 tms.tm_year + 1900, tms.tm_yday + 1,
820 tms.tm_hour, tms.tm_min, tms.tm_sec, ifract);
821 else
822 ret = snprintf (seedtimestr, 18, "%4d,%03d,%02d:%02d:%02d",
823 tms.tm_year + 1900, tms.tm_yday + 1,
824 tms.tm_hour, tms.tm_min, tms.tm_sec);
825
826 if (ret != 24 && ret != 17)
827 return NULL;
828 else
829 return seedtimestr;
830 } /* End of ms_hptime2seedtimestr() */
831
832 /***************************************************************************
833 * ms_time2hptime_int:
834 *
835 * Convert specified time values to a high precision epoch time. This
836 * is an internal version which does no range checking, it is assumed
837 * that checking the range for each value has already been done.
838 *
839 * Returns epoch time on success and HPTERROR on error.
840 ***************************************************************************/
841 static hptime_t
ms_time2hptime_int(int year,int day,int hour,int min,int sec,int usec)842 ms_time2hptime_int (int year, int day, int hour, int min, int sec, int usec)
843 {
844 BTime btime;
845 hptime_t hptime;
846
847 memset (&btime, 0, sizeof (BTime));
848 btime.day = 1;
849
850 /* Convert integer seconds using ms_btime2hptime */
851 btime.year = (int16_t)year;
852 btime.day = (int16_t)day;
853 btime.hour = (uint8_t)hour;
854 btime.min = (uint8_t)min;
855 btime.sec = (uint8_t)sec;
856 btime.fract = 0;
857
858 hptime = ms_btime2hptime (&btime);
859
860 if (hptime == HPTERROR)
861 {
862 ms_log (2, "ms_time2hptime(): Error converting with ms_btime2hptime()\n");
863 return HPTERROR;
864 }
865
866 /* Add the microseconds */
867 hptime += (hptime_t)usec * (1000000 / HPTMODULUS);
868
869 return hptime;
870 } /* End of ms_time2hptime_int() */
871
872 /***************************************************************************
873 * ms_time2hptime:
874 *
875 * Convert specified time values to a high precision epoch time. This
876 * is essentially a frontend for ms_time2hptime that does range
877 * checking for each input value.
878 *
879 * Expected ranges:
880 * year : 1800 - 5000
881 * day : 1 - 366
882 * hour : 0 - 23
883 * min : 0 - 59
884 * sec : 0 - 60
885 * usec : 0 - 999999
886 *
887 * Returns epoch time on success and HPTERROR on error.
888 ***************************************************************************/
889 hptime_t
ms_time2hptime(int year,int day,int hour,int min,int sec,int usec)890 ms_time2hptime (int year, int day, int hour, int min, int sec, int usec)
891 {
892 if (year < 1800 || year > 5000)
893 {
894 ms_log (2, "ms_time2hptime(): Error with year value: %d\n", year);
895 return HPTERROR;
896 }
897
898 if (day < 1 || day > 366)
899 {
900 ms_log (2, "ms_time2hptime(): Error with day value: %d\n", day);
901 return HPTERROR;
902 }
903
904 if (hour < 0 || hour > 23)
905 {
906 ms_log (2, "ms_time2hptime(): Error with hour value: %d\n", hour);
907 return HPTERROR;
908 }
909
910 if (min < 0 || min > 59)
911 {
912 ms_log (2, "ms_time2hptime(): Error with minute value: %d\n", min);
913 return HPTERROR;
914 }
915
916 if (sec < 0 || sec > 60)
917 {
918 ms_log (2, "ms_time2hptime(): Error with second value: %d\n", sec);
919 return HPTERROR;
920 }
921
922 if (usec < 0 || usec > 999999)
923 {
924 ms_log (2, "ms_time2hptime(): Error with microsecond value: %d\n", usec);
925 return HPTERROR;
926 }
927
928 return ms_time2hptime_int (year, day, hour, min, sec, usec);
929 } /* End of ms_time2hptime() */
930
931 /***************************************************************************
932 * ms_seedtimestr2hptime:
933 *
934 * Convert a SEED time string (day-of-year style) to a high precision
935 * epoch time. The time format expected is
936 * "YYYY[,DDD,HH,MM,SS.FFFFFF]", the delimiter can be a dash [-],
937 * comma [,], colon [:] or period [.]. Additionally a [T] or space
938 * may be used to seprate the day and hour fields. The fractional
939 * seconds ("FFFFFF") must begin with a period [.] if present.
940 *
941 * The time string can be "short" in which case the omitted values are
942 * assumed to be zero (with the exception of DDD which is assumed to
943 * be 1): "YYYY,DDD,HH" assumes MM, SS and FFFF are 0. The year is
944 * required, otherwise there wouldn't be much for a date.
945 *
946 * Ranges are checked for each value.
947 *
948 * Returns epoch time on success and HPTERROR on error.
949 ***************************************************************************/
950 hptime_t
ms_seedtimestr2hptime(char * seedtimestr)951 ms_seedtimestr2hptime (char *seedtimestr)
952 {
953 int fields;
954 int year = 0;
955 int day = 1;
956 int hour = 0;
957 int min = 0;
958 int sec = 0;
959 float fusec = 0.0;
960 int usec = 0;
961
962 fields = sscanf (seedtimestr, "%d%*[-,:.]%d%*[-,:.Tt ]%d%*[-,:.]%d%*[-,:.]%d%f",
963 &year, &day, &hour, &min, &sec, &fusec);
964
965 /* Convert fractional seconds to microseconds */
966 if (fusec != 0.0)
967 {
968 usec = (int)(fusec * 1000000.0 + 0.5);
969 }
970
971 if (fields < 1)
972 {
973 ms_log (2, "ms_seedtimestr2hptime(): Error converting time string: %s\n", seedtimestr);
974 return HPTERROR;
975 }
976
977 if (year < 1800 || year > 5000)
978 {
979 ms_log (2, "ms_seedtimestr2hptime(): Error with year value: %d\n", year);
980 return HPTERROR;
981 }
982
983 if (day < 1 || day > 366)
984 {
985 ms_log (2, "ms_seedtimestr2hptime(): Error with day value: %d\n", day);
986 return HPTERROR;
987 }
988
989 if (hour < 0 || hour > 23)
990 {
991 ms_log (2, "ms_seedtimestr2hptime(): Error with hour value: %d\n", hour);
992 return HPTERROR;
993 }
994
995 if (min < 0 || min > 59)
996 {
997 ms_log (2, "ms_seedtimestr2hptime(): Error with minute value: %d\n", min);
998 return HPTERROR;
999 }
1000
1001 if (sec < 0 || sec > 60)
1002 {
1003 ms_log (2, "ms_seedtimestr2hptime(): Error with second value: %d\n", sec);
1004 return HPTERROR;
1005 }
1006
1007 if (usec < 0 || usec > 999999)
1008 {
1009 ms_log (2, "ms_seedtimestr2hptime(): Error with fractional second value: %d\n", usec);
1010 return HPTERROR;
1011 }
1012
1013 return ms_time2hptime_int (year, day, hour, min, sec, usec);
1014 } /* End of ms_seedtimestr2hptime() */
1015
1016 /***************************************************************************
1017 * ms_timestr2hptime:
1018 *
1019 * Convert a generic time string to a high precision epoch time. The
1020 * time format expected is "YYYY[/MM/DD HH:MM:SS.FFFF]", the delimiter
1021 * can be a dash [-], comma[,], slash [/], colon [:], or period [.].
1022 * Additionally a 'T' or space may be used between the date and time
1023 * fields. The fractional seconds ("FFFFFF") must begin with a period
1024 * [.] if present.
1025 *
1026 * The time string can be "short" in which case the omitted values are
1027 * assumed to be zero (with the exception of month and day which are
1028 * assumed to be 1): "YYYY/MM/DD" assumes HH, MM, SS and FFFF are 0.
1029 * The year is required, otherwise there wouldn't be much for a date.
1030 *
1031 * Ranges are checked for each value.
1032 *
1033 * Returns epoch time on success and HPTERROR on error.
1034 ***************************************************************************/
1035 hptime_t
ms_timestr2hptime(char * timestr)1036 ms_timestr2hptime (char *timestr)
1037 {
1038 int fields;
1039 int year = 0;
1040 int mon = 1;
1041 int mday = 1;
1042 int day = 1;
1043 int hour = 0;
1044 int min = 0;
1045 int sec = 0;
1046 float fusec = 0.0;
1047 int usec = 0;
1048
1049 fields = sscanf (timestr, "%d%*[-,/:.]%d%*[-,/:.]%d%*[-,/:.Tt ]%d%*[-,/:.]%d%*[-,/:.]%d%f",
1050 &year, &mon, &mday, &hour, &min, &sec, &fusec);
1051
1052 /* Convert fractional seconds to microseconds */
1053 if (fusec != 0.0)
1054 {
1055 usec = (int)(fusec * 1000000.0 + 0.5);
1056 }
1057
1058 if (fields < 1)
1059 {
1060 ms_log (2, "ms_timestr2hptime(): Error converting time string: %s\n", timestr);
1061 return HPTERROR;
1062 }
1063
1064 if (year < 1800 || year > 5000)
1065 {
1066 ms_log (2, "ms_timestr2hptime(): Error with year value: %d\n", year);
1067 return HPTERROR;
1068 }
1069
1070 if (mon < 1 || mon > 12)
1071 {
1072 ms_log (2, "ms_timestr2hptime(): Error with month value: %d\n", mon);
1073 return HPTERROR;
1074 }
1075
1076 if (mday < 1 || mday > 31)
1077 {
1078 ms_log (2, "ms_timestr2hptime(): Error with day value: %d\n", mday);
1079 return HPTERROR;
1080 }
1081
1082 /* Convert month and day-of-month to day-of-year */
1083 if (ms_md2doy (year, mon, mday, &day))
1084 {
1085 return HPTERROR;
1086 }
1087
1088 if (hour < 0 || hour > 23)
1089 {
1090 ms_log (2, "ms_timestr2hptime(): Error with hour value: %d\n", hour);
1091 return HPTERROR;
1092 }
1093
1094 if (min < 0 || min > 59)
1095 {
1096 ms_log (2, "ms_timestr2hptime(): Error with minute value: %d\n", min);
1097 return HPTERROR;
1098 }
1099
1100 if (sec < 0 || sec > 60)
1101 {
1102 ms_log (2, "ms_timestr2hptime(): Error with second value: %d\n", sec);
1103 return HPTERROR;
1104 }
1105
1106 if (usec < 0 || usec > 999999)
1107 {
1108 ms_log (2, "ms_timestr2hptime(): Error with fractional second value: %d\n", usec);
1109 return HPTERROR;
1110 }
1111
1112 return ms_time2hptime_int (year, day, hour, min, sec, usec);
1113 } /* End of ms_timestr2hptime() */
1114
1115 /***************************************************************************
1116 * ms_nomsamprate:
1117 *
1118 * Calculate a sample rate from SEED sample rate factor and multiplier
1119 * as stored in the fixed section header of data records.
1120 *
1121 * Returns the positive sample rate.
1122 ***************************************************************************/
1123 double
ms_nomsamprate(int factor,int multiplier)1124 ms_nomsamprate (int factor, int multiplier)
1125 {
1126 double samprate = 0.0;
1127
1128 if (factor > 0)
1129 samprate = (double)factor;
1130 else if (factor < 0)
1131 samprate = -1.0 / (double)factor;
1132 if (multiplier > 0)
1133 samprate = samprate * (double)multiplier;
1134 else if (multiplier < 0)
1135 samprate = -1.0 * (samprate / (double)multiplier);
1136
1137 return samprate;
1138 } /* End of ms_nomsamprate() */
1139
1140 /***************************************************************************
1141 * ms_readleapseconds:
1142 *
1143 * Read leap seconds from a file indicated by the specified
1144 * environment variable and populate the global leapsecondlist.
1145 *
1146 * Returns positive number of leap seconds read, -1 on file read
1147 * error, and -2 when the environment variable is not set.
1148 ***************************************************************************/
1149 int
ms_readleapseconds(char * envvarname)1150 ms_readleapseconds (char *envvarname)
1151 {
1152 char *filename;
1153
1154 if ((filename = getenv (envvarname)))
1155 {
1156 return ms_readleapsecondfile (filename);
1157 }
1158
1159 return -2;
1160 } /* End of ms_readleapseconds() */
1161
1162 /***************************************************************************
1163 * ms_readleapsecondfile:
1164 *
1165 * Read leap seconds from the specified file and populate the global
1166 * leapsecondlist. The file is expected to be standard IETF leap
1167 * second list format. The list is usually available from:
1168 * https://www.ietf.org/timezones/data/leap-seconds.list
1169 *
1170 * Returns positive number of leap seconds read on success and -1 on error.
1171 ***************************************************************************/
1172 int
ms_readleapsecondfile(char * filename)1173 ms_readleapsecondfile (char *filename)
1174 {
1175 FILE *fp = NULL;
1176 LeapSecond *ls = NULL;
1177 LeapSecond *lastls = NULL;
1178 int64_t expires;
1179 char readline[200];
1180 char *cp;
1181 int64_t leapsecond;
1182 int TAIdelta;
1183 int fields;
1184 int count = 0;
1185
1186 if (!filename)
1187 return -1;
1188
1189 if (!(fp = fopen (filename, "rb")))
1190 {
1191 ms_log (2, "Cannot open leap second file %s: %s\n", filename, strerror (errno));
1192 return -1;
1193 }
1194
1195 while (fgets (readline, sizeof (readline) - 1, fp))
1196 {
1197 /* Guarantee termination */
1198 readline[sizeof (readline) - 1] = '\0';
1199
1200 /* Terminate string at first newline character if any */
1201 if ((cp = strchr (readline, '\n')))
1202 *cp = '\0';
1203
1204 /* Skip empty lines */
1205 if (!strlen (readline))
1206 continue;
1207
1208 /* Check for and parse expiration date */
1209 if (!strncmp (readline, "#@", 2))
1210 {
1211 expires = 0;
1212 fields = sscanf (readline, "#@ %" SCNd64, &expires);
1213
1214 if (fields == 1)
1215 {
1216 /* Convert expires to Unix epoch */
1217 expires = expires - NTPPOSIXEPOCHDELTA;
1218
1219 /* Compare expire time to current time */
1220 if (time (NULL) > expires)
1221 {
1222 char timestr[100];
1223 ms_hptime2mdtimestr (MS_EPOCH2HPTIME (expires), timestr, 0);
1224 ms_log (1, "Warning: leap second file (%s) has expired as of %s\n",
1225 filename, timestr);
1226 }
1227 }
1228
1229 continue;
1230 }
1231
1232 /* Skip comment lines */
1233 if (*readline == '#')
1234 continue;
1235
1236 fields = sscanf (readline, "%" SCNd64 " %d ", &leapsecond, &TAIdelta);
1237
1238 if (fields == 2)
1239 {
1240 if ((ls = malloc (sizeof (LeapSecond))) == NULL)
1241 {
1242 ms_log (2, "Cannot allocate LeapSecond, out of memory?\n");
1243 return -1;
1244 }
1245
1246 /* Convert NTP epoch time to Unix epoch time and then to HPT */
1247 ls->leapsecond = MS_EPOCH2HPTIME ((leapsecond - NTPPOSIXEPOCHDELTA));
1248 ls->TAIdelta = TAIdelta;
1249 ls->next = NULL;
1250
1251 /* Add leap second to global list */
1252 if (!leapsecondlist)
1253 {
1254 leapsecondlist = ls;
1255 lastls = ls;
1256 }
1257 else
1258 {
1259 lastls->next = ls;
1260 lastls = ls;
1261 }
1262 }
1263 else
1264 {
1265 ms_log (1, "Unrecognized leap second file line: '%s'\n", readline);
1266 }
1267 }
1268
1269 if (ferror (fp))
1270 {
1271 ms_log (2, "Error reading leap second file (%s): %s\n", filename, strerror (errno));
1272 }
1273
1274 fclose (fp);
1275
1276 return count;
1277 } /* End of ms_readleapsecondfile() */
1278
1279 /***************************************************************************
1280 * ms_reduce_rate:
1281 *
1282 * Reduce the specified sample rate into two "factors" (in some cases
1283 * the second factor is actually a divisor).
1284 *
1285 * Integer rates between 1 and 32767 can be represented exactly.
1286 *
1287 * Integer rates higher than 32767 will be matched as closely as
1288 * possible with the deviation becoming larger as the integers reach
1289 * (32767 * 32767).
1290 *
1291 * Non-integer rates between 32767.0 and 1.0/32767.0 are represented
1292 * exactly when possible and approximated otherwise.
1293 *
1294 * Non-integer rates greater than 32767 or less than 1/32767 are not supported.
1295 *
1296 * Returns 0 on success and -1 on error.
1297 ***************************************************************************/
1298 int
ms_reduce_rate(double samprate,int16_t * factor1,int16_t * factor2)1299 ms_reduce_rate (double samprate, int16_t *factor1, int16_t *factor2)
1300 {
1301 int num;
1302 int den;
1303 int32_t intsamprate = (int32_t) (samprate + 0.5);
1304
1305 int32_t searchfactor1;
1306 int32_t searchfactor2;
1307 int32_t closestfactor;
1308 int32_t closestdiff;
1309 int32_t diff;
1310
1311 /* Handle case of integer sample values. */
1312 if (ms_dabs (samprate - intsamprate) < 0.0000001)
1313 {
1314 /* If integer sample rate is less than range of 16-bit int set it directly */
1315 if (intsamprate <= 32767)
1316 {
1317 *factor1 = intsamprate;
1318 *factor2 = 1;
1319 return 0;
1320 }
1321 /* If integer sample rate is within the maximum possible nominal rate */
1322 else if (intsamprate <= (32767 * 32767))
1323 {
1324 /* Determine the closest factors that represent the sample rate.
1325 * The approximation gets worse as the values increase. */
1326 searchfactor1 = (int)(1.0 / ms_rsqrt64 (samprate));
1327 closestdiff = searchfactor1;
1328 closestfactor = searchfactor1;
1329
1330 while ((intsamprate % searchfactor1) != 0)
1331 {
1332 searchfactor1 -= 1;
1333
1334 /* Track the factor that generates the closest match */
1335 searchfactor2 = intsamprate / searchfactor1;
1336 diff = intsamprate - (searchfactor1 * searchfactor2);
1337 if (diff < closestdiff)
1338 {
1339 closestdiff = diff;
1340 closestfactor = searchfactor1;
1341 }
1342
1343 /* If the next iteration would create a factor beyond the limit
1344 * we accept the closest factor */
1345 if ((intsamprate / (searchfactor1 - 1)) > 32767)
1346 {
1347 searchfactor1 = closestfactor;
1348 break;
1349 }
1350 }
1351
1352 searchfactor2 = intsamprate / searchfactor1;
1353
1354 if (searchfactor1 <= 32767 && searchfactor2 <= 32767)
1355 {
1356 *factor1 = searchfactor1;
1357 *factor2 = searchfactor2;
1358 return 0;
1359 }
1360 }
1361 }
1362 /* Handle case of non-integer less than 16-bit int range */
1363 else if (samprate <= 32767.0)
1364 {
1365 /* For samples/seconds, determine, potentially approximate, numerator and denomiator */
1366 ms_ratapprox (samprate, &num, &den, 32767, 1e-8);
1367
1368 /* Negate the factor2 to denote a division operation */
1369 *factor1 = (int16_t)num;
1370 *factor2 = (int16_t)-den;
1371 return 0;
1372 }
1373
1374 return -1;
1375 } /* End of ms_reduce_rate() */
1376
1377 /***************************************************************************
1378 * ms_genfactmult:
1379 *
1380 * Generate an appropriate SEED sample rate factor and multiplier from
1381 * a double precision sample rate.
1382 *
1383 * If the samplerate > 0.0 it is expected to be a rate in SAMPLES/SECOND.
1384 * If the samplerate < 0.0 it is expected to be a period in SECONDS/SAMPLE.
1385 *
1386 * Results use SAMPLES/SECOND notation when sample rate >= 1.0
1387 * Results use SECONDS/SAMPLE notation when samles rates < 1.0
1388 *
1389 * Returns 0 on success and -1 on error or calculation not possible.
1390 ***************************************************************************/
1391 int
ms_genfactmult(double samprate,int16_t * factor,int16_t * multiplier)1392 ms_genfactmult (double samprate, int16_t *factor, int16_t *multiplier)
1393 {
1394 int16_t factor1;
1395 int16_t factor2;
1396
1397 if (!factor || !multiplier)
1398 return -1;
1399
1400 /* Convert sample period to sample rate */
1401 if (samprate < 0.0)
1402 samprate = -1.0 / samprate;
1403
1404 /* Handle special case of zero */
1405 if (samprate == 0.0)
1406 {
1407 *factor = 0;
1408 *multiplier = 0;
1409 return 0;
1410 }
1411 /* Handle sample rates >= 1.0 with the SAMPLES/SECOND representation */
1412 else if (samprate >= 1.0)
1413 {
1414 if (ms_reduce_rate (samprate, &factor1, &factor2) == 0)
1415 {
1416 *factor = factor1;
1417 *multiplier = factor2;
1418 return 0;
1419 }
1420 }
1421 /* Handle sample rates < 1 with the SECONDS/SAMPLE representation */
1422 else
1423 {
1424 /* Reduce rate as a sample period and invert factor/multiplier */
1425 if (ms_reduce_rate (1.0 / samprate, &factor1, &factor2) == 0)
1426 {
1427 *factor = -factor1;
1428 *multiplier = -factor2;
1429 return 0;
1430 }
1431 }
1432
1433 return -1;
1434 } /* End of ms_genfactmult() */
1435
1436 /***************************************************************************
1437 * ms_ratapprox:
1438 *
1439 * Find an approximate rational number for a real through continued
1440 * fraction expansion. Given a double precsion 'real' find a
1441 * numerator (num) and denominator (den) whose absolute values are not
1442 * larger than 'maxval' while trying to reach a specified 'precision'.
1443 *
1444 * Returns the number of iterations performed.
1445 ***************************************************************************/
1446 int
ms_ratapprox(double real,int * num,int * den,int maxval,double precision)1447 ms_ratapprox (double real, int *num, int *den, int maxval, double precision)
1448 {
1449 double realj, preal;
1450 char pos;
1451 int pnum, pden;
1452 int iterations = 1;
1453 int Aj1, Aj2, Bj1, Bj2;
1454 int bj = 0;
1455 int Aj = 0;
1456 int Bj = 1;
1457
1458 if (real >= 0.0)
1459 {
1460 pos = 1;
1461 realj = real;
1462 }
1463 else
1464 {
1465 pos = 0;
1466 realj = -real;
1467 }
1468
1469 preal = realj;
1470
1471 bj = (int)(realj + precision);
1472 realj = 1 / (realj - bj);
1473 Aj = bj;
1474 Aj1 = 1;
1475 Bj = 1;
1476 Bj1 = 0;
1477 *num = pnum = Aj;
1478 *den = pden = Bj;
1479 if (!pos)
1480 *num = -*num;
1481
1482 while (ms_dabs (preal - (double)Aj / (double)Bj) > precision &&
1483 Aj < maxval && Bj < maxval)
1484 {
1485 Aj2 = Aj1;
1486 Aj1 = Aj;
1487 Bj2 = Bj1;
1488 Bj1 = Bj;
1489 bj = (int)(realj + precision);
1490 realj = 1 / (realj - bj);
1491 Aj = bj * Aj1 + Aj2;
1492 Bj = bj * Bj1 + Bj2;
1493 *num = pnum;
1494 *den = pden;
1495 if (!pos)
1496 *num = -*num;
1497 pnum = Aj;
1498 pden = Bj;
1499
1500 iterations++;
1501 }
1502
1503 if (pnum < maxval && pden < maxval)
1504 {
1505 *num = pnum;
1506 *den = pden;
1507 if (!pos)
1508 *num = -*num;
1509 }
1510
1511 return iterations;
1512 }
1513
1514 /***************************************************************************
1515 * ms_bigendianhost:
1516 *
1517 * Determine the byte order of the host machine. Due to the lack of
1518 * portable defines to determine host byte order this run-time test is
1519 * provided. The code below actually tests for little-endianess, the
1520 * only other alternative is assumed to be big endian.
1521 *
1522 * Returns 0 if the host is little endian, otherwise 1.
1523 ***************************************************************************/
1524 int
ms_bigendianhost(void)1525 ms_bigendianhost (void)
1526 {
1527 int16_t host = 1;
1528 return !(*((int8_t *)(&host)));
1529 } /* End of ms_bigendianhost() */
1530
1531 /***************************************************************************
1532 * ms_dabs:
1533 *
1534 * Determine the absolute value of an input double, actually just test
1535 * if the input double is positive multiplying by -1.0 if not and
1536 * return it.
1537 *
1538 * Returns the positive value of input double.
1539 ***************************************************************************/
1540 double
ms_dabs(double val)1541 ms_dabs (double val)
1542 {
1543 if (val < 0.0)
1544 val *= -1.0;
1545 return val;
1546 } /* End of ms_dabs() */
1547
1548 /***************************************************************************
1549 * ms_rsqrt64:
1550 *
1551 * An optimized reciprocal square root calculation from:
1552 * Matthew Robertson (2012). "A Brief History of InvSqrt"
1553 * https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf
1554 *
1555 * Further reference and description:
1556 * https://en.wikipedia.org/wiki/Fast_inverse_square_root
1557 *
1558 * Modifications:
1559 * Add 2 more iterations of Newton's method to increase accuracy,
1560 * specifically for large values.
1561 * Use memcpy instead of assignment through differing pointer types.
1562 *
1563 * Returns 0 if the host is little endian, otherwise 1.
1564 ***************************************************************************/
1565 double
ms_rsqrt64(double val)1566 ms_rsqrt64 (double val)
1567 {
1568 uint64_t i;
1569 double x2;
1570 double y;
1571
1572 x2 = val * 0.5;
1573 y = val;
1574 memcpy (&i, &y, sizeof(i));
1575 i = 0x5fe6eb50c7b537a9ULL - (i >> 1);
1576 memcpy (&y, &i, sizeof(y));
1577 y = y * (1.5 - (x2 * y * y));
1578 y = y * (1.5 - (x2 * y * y));
1579 y = y * (1.5 - (x2 * y * y));
1580
1581 return y;
1582 } /* End of ms_rsqrt64() */
1583
1584 /***************************************************************************
1585 * ms_gmtime_r:
1586 *
1587 * An internal version of gmtime_r() that is 64-bit compliant and
1588 * works with years beyond 2038.
1589 *
1590 * The original was called pivotal_gmtime_r() by Paul Sheer, all
1591 * required copyright and other hoohas are below. Modifications were
1592 * made to integrate the original to this code base, avoid name
1593 * collisions and formatting so I could read it.
1594 *
1595 * Returns a pointer to the populated tm struct on success and NULL on error.
1596 ***************************************************************************/
1597
1598 /* pivotal_gmtime_r - a replacement for gmtime/localtime/mktime
1599 that works around the 2038 bug on 32-bit
1600 systems. (Version 4)
1601
1602 Copyright (C) 2009 Paul Sheer
1603
1604 Redistribution and use in source form, with or without modification,
1605 is permitted provided that the above copyright notice, this list of
1606 conditions, the following disclaimer, and the following char array
1607 are retained.
1608
1609 Redistribution and use in binary form must reproduce an
1610 acknowledgment: 'With software provided by http://2038bug.com/' in
1611 the documentation and/or other materials provided with the
1612 distribution, and wherever such acknowledgments are usually
1613 accessible in Your program.
1614
1615 This software is provided "AS IS" and WITHOUT WARRANTY, either
1616 express or implied, including, without limitation, the warranties of
1617 NON-INFRINGEMENT, MERCHANTABILITY or FITNESS FOR A PARTICULAR
1618 PURPOSE. THE ENTIRE RISK AS TO THE QUALITY OF THIS SOFTWARE IS WITH
1619 YOU. Under no circumstances and under no legal theory, whether in
1620 tort (including negligence), contract, or otherwise, shall the
1621 copyright owners be liable for any direct, indirect, special,
1622 incidental, or consequential damages of any character arising as a
1623 result of the use of this software including, without limitation,
1624 damages for loss of goodwill, work stoppage, computer failure or
1625 malfunction, or any and all other commercial damages or losses. This
1626 limitation of liability shall not apply to liability for death or
1627 personal injury resulting from copyright owners' negligence to the
1628 extent applicable law prohibits such limitation. Some jurisdictions
1629 do not allow the exclusion or limitation of incidental or
1630 consequential damages, so this exclusion and limitation may not apply
1631 to You.
1632
1633 */
1634
1635 const char pivotal_gmtime_r_stamp_lm[] =
1636 "pivotal_gmtime_r. Copyright (C) 2009 Paul Sheer. Terms and "
1637 "conditions apply. Visit http://2038bug.com/ for more info.";
1638
1639 static const int tm_days[4][13] = {
1640 {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
1641 {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
1642 {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365},
1643 {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366},
1644 };
1645
1646 #define TM_LEAP_CHECK(n) ((!(((n) + 1900) % 400) || (!(((n) + 1900) % 4) && (((n) + 1900) % 100))) != 0)
1647 #define TM_WRAP(a, b, m) ((a) = ((a) < 0) ? ((b)--, (a) + (m)) : (a))
1648
1649 static struct tm *
ms_gmtime_r(int64_t * timep,struct tm * result)1650 ms_gmtime_r (int64_t *timep, struct tm *result)
1651 {
1652 int v_tm_sec, v_tm_min, v_tm_hour, v_tm_mon, v_tm_wday, v_tm_tday;
1653 int leap;
1654 long m;
1655 int64_t tv;
1656
1657 if (!timep || !result)
1658 return NULL;
1659
1660 tv = *timep;
1661
1662 v_tm_sec = ((int64_t)tv % (int64_t)60);
1663 tv /= 60;
1664 v_tm_min = ((int64_t)tv % (int64_t)60);
1665 tv /= 60;
1666 v_tm_hour = ((int64_t)tv % (int64_t)24);
1667 tv /= 24;
1668 v_tm_tday = (int)tv;
1669
1670 TM_WRAP (v_tm_sec, v_tm_min, 60);
1671 TM_WRAP (v_tm_min, v_tm_hour, 60);
1672 TM_WRAP (v_tm_hour, v_tm_tday, 24);
1673
1674 if ((v_tm_wday = (v_tm_tday + 4) % 7) < 0)
1675 v_tm_wday += 7;
1676
1677 m = (long)v_tm_tday;
1678
1679 if (m >= 0)
1680 {
1681 result->tm_year = 70;
1682 leap = TM_LEAP_CHECK (result->tm_year);
1683
1684 while (m >= (long)tm_days[leap + 2][12])
1685 {
1686 m -= (long)tm_days[leap + 2][12];
1687 result->tm_year++;
1688 leap = TM_LEAP_CHECK (result->tm_year);
1689 }
1690
1691 v_tm_mon = 0;
1692
1693 while (m >= (long)tm_days[leap][v_tm_mon])
1694 {
1695 m -= (long)tm_days[leap][v_tm_mon];
1696 v_tm_mon++;
1697 }
1698 }
1699 else
1700 {
1701 result->tm_year = 69;
1702 leap = TM_LEAP_CHECK (result->tm_year);
1703
1704 while (m < (long)-tm_days[leap + 2][12])
1705 {
1706 m += (long)tm_days[leap + 2][12];
1707 result->tm_year--;
1708 leap = TM_LEAP_CHECK (result->tm_year);
1709 }
1710
1711 v_tm_mon = 11;
1712
1713 while (m < (long)-tm_days[leap][v_tm_mon])
1714 {
1715 m += (long)tm_days[leap][v_tm_mon];
1716 v_tm_mon--;
1717 }
1718
1719 m += (long)tm_days[leap][v_tm_mon];
1720 }
1721
1722 result->tm_mday = (int)m + 1;
1723 result->tm_yday = tm_days[leap + 2][v_tm_mon] + m;
1724 result->tm_sec = v_tm_sec;
1725 result->tm_min = v_tm_min;
1726 result->tm_hour = v_tm_hour;
1727 result->tm_mon = v_tm_mon;
1728 result->tm_wday = v_tm_wday;
1729
1730 return result;
1731 } /* End of ms_gmtime_r() */
1732