summaryrefslogtreecommitdiffstats
path: root/contrib/ntp/libntp/ntp_calendar.c
blob: ff6ead364e5ae13a8fa0f57cd33f279914f8860c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
/*
 * ntp_calendar.c - calendar and helper functions
 *
 * Written by Juergen Perlinger (perlinger@ntp.org) for the NTP project.
 * The contents of 'html/copyright.html' apply.
 *
 * --------------------------------------------------------------------
 * Some notes on the implementation:
 *
 * Calendar algorithms thrive on the division operation, which is one of
 * the slowest numerical operations in any CPU. What saves us here from
 * abysmal performance is the fact that all divisions are divisions by
 * constant numbers, and most compilers can do this by a multiplication
 * operation.  But this might not work when using the div/ldiv/lldiv
 * function family, because many compilers are not able to do inline
 * expansion of the code with following optimisation for the
 * constant-divider case.
 *
 * Also div/ldiv/lldiv are defined in terms of int/long/longlong, which
 * are inherently target dependent. Nothing that could not be cured with
 * autoconf, but still a mess...
 *
 * Furthermore, we need floor division in many places. C either leaves
 * the division behaviour undefined (< C99) or demands truncation to
 * zero (>= C99), so additional steps are required to make sure the
 * algorithms work. The {l,ll}div function family is requested to
 * truncate towards zero, which is also the wrong direction for our
 * purpose.
 *
 * For all this, all divisions by constant are coded manually, even when
 * there is a joined div/mod operation: The optimiser should sort that
 * out, if possible. Most of the calculations are done with unsigned
 * types, explicitely using two's complement arithmetics where
 * necessary. This minimises the dependecies to compiler and target,
 * while still giving reasonable to good performance.
 *
 * The implementation uses a few tricks that exploit properties of the
 * two's complement: Floor division on negative dividents can be
 * executed by using the one's complement of the divident. One's
 * complement can be easily created using XOR and a mask.
 *
 * Finally, check for overflow conditions is minimal. There are only two
 * calculation steps in the whole calendar that suffer from an internal
 * overflow, and these conditions are checked: errno is set to EDOM and
 * the results are clamped/saturated in this case.  All other functions
 * do not suffer from internal overflow and simply return the result
 * truncated to 32 bits.
 *
 * This is a sacrifice made for execution speed.  Since a 32-bit day
 * counter covers +/- 5,879,610 years and the clamp limits the effective
 * range to +/-2.9 million years, this should not pose a problem here.
 *
 */

#include <config.h>
#include <sys/types.h>

#include "ntp_types.h"
#include "ntp_calendar.h"
#include "ntp_stdlib.h"
#include "ntp_fp.h"
#include "ntp_unixtime.h"

/* For now, let's take the conservative approach: if the target property
 * macros are not defined, check a few well-known compiler/architecture
 * settings. Default is to assume that the representation of signed
 * integers is unknown and shift-arithmetic-right is not available.
 */
#ifndef TARGET_HAS_2CPL
# if defined(__GNUC__)
#  if defined(__i386__) || defined(__x86_64__) || defined(__arm__)
#   define TARGET_HAS_2CPL 1
#  else
#   define TARGET_HAS_2CPL 0
#  endif
# elif defined(_MSC_VER)
#  if defined(_M_IX86) || defined(_M_X64) || defined(_M_ARM)
#   define TARGET_HAS_2CPL 1
#  else
#   define TARGET_HAS_2CPL 0
#  endif
# else
#  define TARGET_HAS_2CPL 0
# endif
#endif

#ifndef TARGET_HAS_SAR
# define TARGET_HAS_SAR 0
#endif

/*
 *---------------------------------------------------------------------
 * replacing the 'time()' function
 * --------------------------------------------------------------------
 */

static systime_func_ptr systime_func = &time;
static inline time_t now(void);


systime_func_ptr
ntpcal_set_timefunc(
	systime_func_ptr nfunc
	)
{
	systime_func_ptr res;

	res = systime_func;
	if (NULL == nfunc)
		nfunc = &time;
	systime_func = nfunc;

	return res;
}


static inline time_t
now(void)
{
	return (*systime_func)(NULL);
}

/*
 *---------------------------------------------------------------------
 * Get sign extension mask and unsigned 2cpl rep for a signed integer
 *---------------------------------------------------------------------
 */

static inline uint32_t
int32_sflag(
	const int32_t v)
{
#   if TARGET_HAS_2CPL && TARGET_HAS_SAR && SIZEOF_INT >= 4

	/* Let's assume that shift is the fastest way to get the sign
	 * extension of of a signed integer. This might not always be
	 * true, though -- On 8bit CPUs or machines without barrel
	 * shifter this will kill the performance. So we make sure
	 * we do this only if 'int' has at least 4 bytes.
	 */
	return (uint32_t)(v >> 31);
	
#   else

	/* This should be a rather generic approach for getting a sign
	 * extension mask...
	 */
	return UINT32_C(0) - (uint32_t)(v < 0);
	
#   endif
}

static inline uint32_t
int32_to_uint32_2cpl(
	const int32_t v)
{
	uint32_t vu;
	
#   if TARGET_HAS_2CPL

	/* Just copy through the 32 bits from the signed value if we're
	 * on a two's complement target.
	 */
	vu = (uint32_t)v;
	
#   else

	/* Convert from signed int to unsigned int two's complement. Do
	 * not make any assumptions about the representation of signed
	 * integers, but make sure signed integer overflow cannot happen
	 * here. A compiler on a two's complement target *might* find
	 * out that this is just a complicated cast (as above), but your
	 * mileage might vary.
	 */
	if (v < 0)
		vu = ~(uint32_t)(-(v + 1));
	else
		vu = (uint32_t)v;
	
#   endif
	
	return vu;
}

static inline int32_t
uint32_2cpl_to_int32(
	const uint32_t vu)
{
	int32_t v;
	
#   if TARGET_HAS_2CPL

	/* Just copy through the 32 bits from the unsigned value if
	 * we're on a two's complement target.
	 */
	v = (int32_t)vu;

#   else

	/* Convert to signed integer, making sure signed integer
	 * overflow cannot happen. Again, the optimiser might or might
	 * not find out that this is just a copy of 32 bits on a target
	 * with two's complement representation for signed integers.
	 */
	if (vu > INT32_MAX)
		v = -(int32_t)(~vu) - 1;
	else
		v = (int32_t)vu;
	
#   endif
	
	return v;
}

/* Some of the calculations need to multiply the input by 4 before doing
 * a division. This can cause overflow and strange results. Therefore we
 * clamp / saturate the input operand. And since we do the calculations
 * in unsigned int with an extra sign flag/mask, we only loose one bit
 * of the input value range.
 */
static inline uint32_t
uint32_saturate(
	uint32_t vu,
	uint32_t mu)
{
	static const uint32_t limit = UINT32_MAX/4u;
	if ((mu ^ vu) > limit) {
		vu    = mu ^ limit;
		errno = EDOM;
	}
	return vu;
}

/*
 *---------------------------------------------------------------------
 * Convert between 'time_t' and 'vint64'
 *---------------------------------------------------------------------
 */
vint64
time_to_vint64(
	const time_t * ptt
	)
{
	vint64 res;
	time_t tt;

	tt = *ptt;

#   if SIZEOF_TIME_T <= 4

	res.D_s.hi = 0;
	if (tt < 0) {
		res.D_s.lo = (uint32_t)-tt;
		M_NEG(res.D_s.hi, res.D_s.lo);
	} else {
		res.D_s.lo = (uint32_t)tt;
	}

#   elif defined(HAVE_INT64)

	res.q_s = tt;

#   else
	/*
	 * shifting negative signed quantities is compiler-dependent, so
	 * we better avoid it and do it all manually. And shifting more
	 * than the width of a quantity is undefined. Also a don't do!
	 */
	if (tt < 0) {
		tt = -tt;
		res.D_s.lo = (uint32_t)tt;
		res.D_s.hi = (uint32_t)(tt >> 32);
		M_NEG(res.D_s.hi, res.D_s.lo);
	} else {
		res.D_s.lo = (uint32_t)tt;
		res.D_s.hi = (uint32_t)(tt >> 32);
	}

#   endif

	return res;
}


time_t
vint64_to_time(
	const vint64 *tv
	)
{
	time_t res;

#   if SIZEOF_TIME_T <= 4

	res = (time_t)tv->D_s.lo;

#   elif defined(HAVE_INT64)

	res = (time_t)tv->q_s;

#   else

	res = ((time_t)tv->d_s.hi << 32) | tv->D_s.lo;

#   endif

	return res;
}

/*
 *---------------------------------------------------------------------
 * Get the build date & time
 *---------------------------------------------------------------------
 */
int
ntpcal_get_build_date(
	struct calendar * jd
	)
{
	/* The C standard tells us the format of '__DATE__':
	 *
	 * __DATE__ The date of translation of the preprocessing
	 * translation unit: a character string literal of the form "Mmm
	 * dd yyyy", where the names of the months are the same as those
	 * generated by the asctime function, and the first character of
	 * dd is a space character if the value is less than 10. If the
	 * date of translation is not available, an
	 * implementation-defined valid date shall be supplied.
	 *
	 * __TIME__ The time of translation of the preprocessing
	 * translation unit: a character string literal of the form
	 * "hh:mm:ss" as in the time generated by the asctime
	 * function. If the time of translation is not available, an
	 * implementation-defined valid time shall be supplied.
	 *
	 * Note that MSVC declares DATE and TIME to be in the local time
	 * zone, while neither the C standard nor the GCC docs make any
	 * statement about this. As a result, we may be +/-12hrs off
	 * UTC.  But for practical purposes, this should not be a
	 * problem.
	 *
	 */
#   ifdef MKREPRO_DATE
	static const char build[] = MKREPRO_TIME "/" MKREPRO_DATE;
#   else
	static const char build[] = __TIME__ "/" __DATE__;
#   endif
	static const char mlist[] = "JanFebMarAprMayJunJulAugSepOctNovDec";

	char		  monstr[4];
	const char *	  cp;
	unsigned short	  hour, minute, second, day, year;
 	/* Note: The above quantities are used for sscanf 'hu' format,
	 * so using 'uint16_t' is contra-indicated!
	 */

#   ifdef DEBUG
	static int        ignore  = 0;
#   endif

	ZERO(*jd);
	jd->year     = 1970;
	jd->month    = 1;
	jd->monthday = 1;

#   ifdef DEBUG
	/* check environment if build date should be ignored */
	if (0 == ignore) {
	    const char * envstr;
	    envstr = getenv("NTPD_IGNORE_BUILD_DATE");
	    ignore = 1 + (envstr && (!*envstr || !strcasecmp(envstr, "yes")));
	}
	if (ignore > 1)
	    return FALSE;
#   endif

	if (6 == sscanf(build, "%hu:%hu:%hu/%3s %hu %hu",
			&hour, &minute, &second, monstr, &day, &year)) {
		cp = strstr(mlist, monstr);
		if (NULL != cp) {
			jd->year     = year;
			jd->month    = (uint8_t)((cp - mlist) / 3 + 1);
			jd->monthday = (uint8_t)day;
			jd->hour     = (uint8_t)hour;
			jd->minute   = (uint8_t)minute;
			jd->second   = (uint8_t)second;

			return TRUE;
		}
	}

	return FALSE;
}


/*
 *---------------------------------------------------------------------
 * basic calendar stuff
 * --------------------------------------------------------------------
 */

/* month table for a year starting with March,1st */
static const uint16_t shift_month_table[13] = {
	0, 31, 61, 92, 122, 153, 184, 214, 245, 275, 306, 337, 366
};

/* month tables for years starting with January,1st; regular & leap */
static const uint16_t real_month_table[2][13] = {
	/* -*- table for regular years -*- */
	{ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 },
	/* -*- table for leap years -*- */
	{ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 }
};

/*
 * Some notes on the terminology:
 *
 * We use the proleptic Gregorian calendar, which is the Gregorian
 * calendar extended in both directions ad infinitum. This totally
 * disregards the fact that this calendar was invented in 1582, and
 * was adopted at various dates over the world; sometimes even after
 * the start of the NTP epoch.
 *
 * Normally date parts are given as current cycles, while time parts
 * are given as elapsed cycles:
 *
 * 1970-01-01/03:04:05 means 'IN the 1970st. year, IN the first month,
 * ON the first day, with 3hrs, 4minutes and 5 seconds elapsed.
 *
 * The basic calculations for this calendar implementation deal with
 * ELAPSED date units, which is the number of full years, full months
 * and full days before a date: 1970-01-01 would be (1969, 0, 0) in
 * that notation.
 *
 * To ease the numeric computations, month and day values outside the
 * normal range are acceptable: 2001-03-00 will be treated as the day
 * before 2001-03-01, 2000-13-32 will give the same result as
 * 2001-02-01 and so on.
 *
 * 'rd' or 'RD' is used as an abbreviation for the latin 'rata die'
 * (day number).  This is the number of days elapsed since 0000-12-31
 * in the proleptic Gregorian calendar. The begin of the Christian Era
 * (0001-01-01) is RD(1).
 */

/*
 * ==================================================================
 *
 * General algorithmic stuff
 *
 * ==================================================================
 */

/*
 *---------------------------------------------------------------------
 * Do a periodic extension of 'value' around 'pivot' with a period of
 * 'cycle'.
 *
 * The result 'res' is a number that holds to the following properties:
 *
 *   1)	 res MOD cycle == value MOD cycle
 *   2)	 pivot <= res < pivot + cycle
 *	 (replace </<= with >/>= for negative cycles)
 *
 * where 'MOD' denotes the modulo operator for FLOOR DIVISION, which
 * is not the same as the '%' operator in C: C requires division to be
 * a truncated division, where remainder and dividend have the same
 * sign if the remainder is not zero, whereas floor division requires
 * divider and modulus to have the same sign for a non-zero modulus.
 *
 * This function has some useful applications:
 *
 * + let Y be a calendar year and V a truncated 2-digit year: then
 *	periodic_extend(Y-50, V, 100)
 *   is the closest expansion of the truncated year with respect to
 *   the full year, that is a 4-digit year with a difference of less
 *   than 50 years to the year Y. ("century unfolding")
 *
 * + let T be a UN*X time stamp and V be seconds-of-day: then
 *	perodic_extend(T-43200, V, 86400)
 *   is a time stamp that has the same seconds-of-day as the input
 *   value, with an absolute difference to T of <= 12hrs.  ("day
 *   unfolding")
 *
 * + Wherever you have a truncated periodic value and a non-truncated
 *   base value and you want to match them somehow...
 *
 * Basically, the function delivers 'pivot + (value - pivot) % cycle',
 * but the implementation takes some pains to avoid internal signed
 * integer overflows in the '(value - pivot) % cycle' part and adheres
 * to the floor division convention.
 *
 * If 64bit scalars where available on all intended platforms, writing a
 * version that uses 64 bit ops would be easy; writing a general
 * division routine for 64bit ops on a platform that can only do
 * 32/16bit divisions and is still performant is a bit more
 * difficult. Since most usecases can be coded in a way that does only
 * require the 32-bit version a 64bit version is NOT provided here.
 * ---------------------------------------------------------------------
 */
int32_t
ntpcal_periodic_extend(
	int32_t pivot,
	int32_t value,
	int32_t cycle
	)
{
	uint32_t diff;
	char	 cpl = 0; /* modulo complement flag */
	char	 neg = 0; /* sign change flag	    */

	/* make the cycle positive and adjust the flags */
	if (cycle < 0) {
		cycle = - cycle;
		neg ^= 1;
		cpl ^= 1;
	}
	/* guard against div by zero or one */
	if (cycle > 1) {
		/*
		 * Get absolute difference as unsigned quantity and
		 * the complement flag. This is done by always
		 * subtracting the smaller value from the bigger
		 * one.
		 */
		if (value >= pivot) {
			diff = int32_to_uint32_2cpl(value)
			     - int32_to_uint32_2cpl(pivot);
		} else {
			diff = int32_to_uint32_2cpl(pivot)
			     - int32_to_uint32_2cpl(value);
			cpl ^= 1;
		}
		diff %= (uint32_t)cycle;
		if (diff) {
			if (cpl)
				diff = (uint32_t)cycle - diff;
			if (neg)
				diff = ~diff + 1;
			pivot += uint32_2cpl_to_int32(diff);
		}
	}
	return pivot;
}

/*
 *-------------------------------------------------------------------
 * Convert a timestamp in NTP scale to a 64bit seconds value in the UN*X
 * scale with proper epoch unfolding around a given pivot or the current
 * system time. This function happily accepts negative pivot values as
 * timestamps befor 1970-01-01, so be aware of possible trouble on
 * platforms with 32bit 'time_t'!
 *
 * This is also a periodic extension, but since the cycle is 2^32 and
 * the shift is 2^31, we can do some *very* fast math without explicit
 * divisions.
 *-------------------------------------------------------------------
 */
vint64
ntpcal_ntp_to_time(
	uint32_t	ntp,
	const time_t *	pivot
	)
{
	vint64 res;

#   if defined(HAVE_INT64)

	res.q_s = (pivot != NULL)
		      ? *pivot
		      : now();
	res.Q_s -= 0x80000000;		/* unshift of half range */
	ntp	-= (uint32_t)JAN_1970;	/* warp into UN*X domain */
	ntp	-= res.D_s.lo;		/* cycle difference	 */
	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */

#   else /* no 64bit scalars */

	time_t tmp;

	tmp = (pivot != NULL)
		  ? *pivot
		  : now();
	res = time_to_vint64(&tmp);
	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000);
	ntp -= (uint32_t)JAN_1970;	/* warp into UN*X domain */
	ntp -= res.D_s.lo;		/* cycle difference	 */
	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);

#   endif /* no 64bit scalars */

	return res;
}

/*
 *-------------------------------------------------------------------
 * Convert a timestamp in NTP scale to a 64bit seconds value in the NTP
 * scale with proper epoch unfolding around a given pivot or the current
 * system time.
 *
 * Note: The pivot must be given in the UN*X time domain!
 *
 * This is also a periodic extension, but since the cycle is 2^32 and
 * the shift is 2^31, we can do some *very* fast math without explicit
 * divisions.
 *-------------------------------------------------------------------
 */
vint64
ntpcal_ntp_to_ntp(
	uint32_t      ntp,
	const time_t *pivot
	)
{
	vint64 res;

#   if defined(HAVE_INT64)

	res.q_s = (pivot)
		      ? *pivot
		      : now();
	res.Q_s -= 0x80000000;		/* unshift of half range */
	res.Q_s += (uint32_t)JAN_1970;	/* warp into NTP domain	 */
	ntp	-= res.D_s.lo;		/* cycle difference	 */
	res.Q_s += (uint64_t)ntp;	/* get expanded time	 */

#   else /* no 64bit scalars */

	time_t tmp;

	tmp = (pivot)
		  ? *pivot
		  : now();
	res = time_to_vint64(&tmp);
	M_SUB(res.D_s.hi, res.D_s.lo, 0, 0x80000000u);
	M_ADD(res.D_s.hi, res.D_s.lo, 0, (uint32_t)JAN_1970);/*into NTP */
	ntp -= res.D_s.lo;		/* cycle difference	 */
	M_ADD(res.D_s.hi, res.D_s.lo, 0, ntp);

#   endif /* no 64bit scalars */

	return res;
}


/*
 * ==================================================================
 *
 * Splitting values to composite entities
 *
 * ==================================================================
 */

/*
 *-------------------------------------------------------------------
 * Split a 64bit seconds value into elapsed days in 'res.hi' and
 * elapsed seconds since midnight in 'res.lo' using explicit floor
 * division. This function happily accepts negative time values as
 * timestamps before the respective epoch start.
 * -------------------------------------------------------------------
 */
ntpcal_split
ntpcal_daysplit(
	const vint64 *ts
	)
{
	ntpcal_split res;
	uint32_t Q;

#   if defined(HAVE_INT64)
	
	/* Manual floor division by SECSPERDAY. This uses the one's
	 * complement trick, too, but without an extra flag value: The
	 * flag would be 64bit, and that's a bit of overkill on a 32bit
	 * target that has to use a register pair for a 64bit number.
	 */
	if (ts->q_s < 0)
		Q = ~(uint32_t)(~ts->Q_s / SECSPERDAY);
	else
		Q = (uint32_t)(ts->Q_s / SECSPERDAY);

#   else

	uint32_t ah, al, sflag, A;

	/* get operand into ah/al (either ts or ts' one's complement,
	 * for later floor division)
	 */
	sflag = int32_sflag(ts->d_s.hi);
	ah = sflag ^ ts->D_s.hi;
	al = sflag ^ ts->D_s.lo;

	/* Since 86400 == 128*675 we can drop the least 7 bits and
	 * divide by 675 instead of 86400. Then the maximum remainder
	 * after each devision step is 674, and we need 10 bits for
	 * that. So in the next step we can shift in 22 bits from the
	 * numerator.
	 *
	 * Therefore we load the accu with the top 13 bits (51..63) in
	 * the first shot. We don't have to remember the quotient -- it
	 * would be shifted out anyway.
	 */
	A = ah >> 19;
	if (A >= 675)
		A = (A % 675u);

	/* Now assemble the remainder with bits 29..50 from the
	 * numerator and divide. This creates the upper ten bits of the
	 * quotient. (Well, the top 22 bits of a 44bit result. But that
	 * will be truncated to 32 bits anyway.)
	 */
	A = (A << 19) | (ah & 0x0007FFFFu);
	A = (A <<  3) | (al >> 29);
	Q = A / 675u;
	A = A % 675u;

	/* Now assemble the remainder with bits 7..28 from the numerator
	 * and do a final division step.
	 */
	A = (A << 22) | ((al >> 7) & 0x003FFFFFu);
	Q = (Q << 22) | (A / 675u);

	/* The last 7 bits get simply dropped, as they have no affect on
	 * the quotient when dividing by 86400.
	 */

	/* apply sign correction and calculate the true floor
	 * remainder.
	 */
	Q ^= sflag;
	
#   endif
	
	res.hi = uint32_2cpl_to_int32(Q);
	res.lo = ts->D_s.lo - Q * SECSPERDAY;

	return res;
}

/*
 *-------------------------------------------------------------------
 * Split a 32bit seconds value into h/m/s and excessive days.  This
 * function happily accepts negative time values as timestamps before
 * midnight.
 * -------------------------------------------------------------------
 */
static int32_t
priv_timesplit(
	int32_t split[3],
	int32_t ts
	)
{
	/* Do 3 chained floor divisions by positive constants, using the
	 * one's complement trick and factoring out the intermediate XOR
	 * ops to reduce the number of operations.
	 */
	uint32_t us, um, uh, ud, sflag;

	sflag = int32_sflag(ts);
	us    = int32_to_uint32_2cpl(ts);

	um = (sflag ^ us) / SECSPERMIN;
	uh = um / MINSPERHR;
	ud = uh / HRSPERDAY;

	um ^= sflag;
	uh ^= sflag;
	ud ^= sflag;

	split[0] = (int32_t)(uh - ud * HRSPERDAY );
	split[1] = (int32_t)(um - uh * MINSPERHR );
	split[2] = (int32_t)(us - um * SECSPERMIN);
	
	return uint32_2cpl_to_int32(ud);
}

/*
 * ---------------------------------------------------------------------
 * Given the number of elapsed days in the calendar era, split this
 * number into the number of elapsed years in 'res.hi' and the number
 * of elapsed days of that year in 'res.lo'.
 *
 * if 'isleapyear' is not NULL, it will receive an integer that is 0 for
 * regular years and a non-zero value for leap years.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_split_eradays(
	int32_t days,
	int  *isleapyear
	)
{
	/* Use the fast cyclesplit algorithm here, to calculate the
	 * centuries and years in a century with one division each. This
	 * reduces the number of division operations to two, but is
	 * susceptible to internal range overflow. We make sure the
	 * input operands are in the safe range; this still gives us
	 * approx +/-2.9 million years.
	 */
	ntpcal_split res;
	int32_t	 n100, n001; /* calendar year cycles */
	uint32_t uday, Q, sflag;

	/* split off centuries first */
	sflag = int32_sflag(days);
	uday  = uint32_saturate(int32_to_uint32_2cpl(days), sflag);
	uday  = (4u * uday) | 3u;
	Q    = sflag ^ ((sflag ^ uday) / GREGORIAN_CYCLE_DAYS);
	uday = uday - Q * GREGORIAN_CYCLE_DAYS;
	n100 = uint32_2cpl_to_int32(Q);
	
	/* Split off years in century -- days >= 0 here, and we're far
	 * away from integer overflow trouble now. */
	uday |= 3;
	n001 = uday / GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;
	uday = uday % GREGORIAN_NORMAL_LEAP_CYCLE_DAYS;

	/* Assemble the year and day in year */
	res.hi = n100 * 100 + n001;
	res.lo = uday / 4u;

	/* Eventually set the leap year flag. Note: 0 <= n001 <= 99 and
	 * Q is still the two's complement representation of the
	 * centuries: The modulo 4 ops can be done with masking here.
	 * We also shift the year and the century by one, so the tests
	 * can be done against zero instead of 3.
	 */
	if (isleapyear)
		*isleapyear = !((n001+1) & 3)
		    && ((n001 != 99) || !((Q+1) & 3));
	
	return res;
}

/*
 *---------------------------------------------------------------------
 * Given a number of elapsed days in a year and a leap year indicator,
 * split the number of elapsed days into the number of elapsed months in
 * 'res.hi' and the number of elapsed days of that month in 'res.lo'.
 *
 * This function will fail and return {-1,-1} if the number of elapsed
 * days is not in the valid range!
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_split_yeardays(
	int32_t eyd,
	int     isleapyear
	)
{
	ntpcal_split    res;
	const uint16_t *lt;	/* month length table	*/

	/* check leap year flag and select proper table */
	lt = real_month_table[(isleapyear != 0)];
	if (0 <= eyd && eyd < lt[12]) {
		/* get zero-based month by approximation & correction step */
		res.hi = eyd >> 5;	   /* approx month; might be 1 too low */
		if (lt[res.hi + 1] <= eyd) /* fixup approximative month value  */
			res.hi += 1;
		res.lo = eyd - lt[res.hi];
	} else {
		res.lo = res.hi = -1;
	}

	return res;
}

/*
 *---------------------------------------------------------------------
 * Convert a RD into the date part of a 'struct calendar'.
 *---------------------------------------------------------------------
 */
int
ntpcal_rd_to_date(
	struct calendar *jd,
	int32_t		 rd
	)
{
	ntpcal_split split;
	int	     leapy;
	u_int	     ymask;

	/* Get day-of-week first. Since rd is signed, the remainder can
	 * be in the range [-6..+6], but the assignment to an unsigned
	 * variable maps the negative values to positive values >=7.
	 * This makes the sign correction look strange, but adding 7
	 * causes the needed wrap-around into the desired value range of
	 * zero to six, both inclusive.
	 */
	jd->weekday = rd % DAYSPERWEEK;
	if (jd->weekday >= DAYSPERWEEK)	/* weekday is unsigned! */
		jd->weekday += DAYSPERWEEK;

	split = ntpcal_split_eradays(rd - 1, &leapy);
	/* Get year and day-of-year, with overflow check. If any of the
	 * upper 16 bits is set after shifting to unity-based years, we
	 * will have an overflow when converting to an unsigned 16bit
	 * year. Shifting to the right is OK here, since it does not
	 * matter if the shift is logic or arithmetic.
	 */
	split.hi += 1;
	ymask = 0u - ((split.hi >> 16) == 0);
	jd->year = (uint16_t)(split.hi & ymask);
	jd->yearday = (uint16_t)split.lo + 1;

	/* convert to month and mday */
	split = ntpcal_split_yeardays(split.lo, leapy);
	jd->month    = (uint8_t)split.hi + 1;
	jd->monthday = (uint8_t)split.lo + 1;

	return ymask ? leapy : -1;
}

/*
 *---------------------------------------------------------------------
 * Convert a RD into the date part of a 'struct tm'.
 *---------------------------------------------------------------------
 */
int
ntpcal_rd_to_tm(
	struct tm  *utm,
	int32_t	    rd
	)
{
	ntpcal_split split;
	int	     leapy;

	/* get day-of-week first */
	utm->tm_wday = rd % DAYSPERWEEK;
	if (utm->tm_wday < 0)
		utm->tm_wday += DAYSPERWEEK;

	/* get year and day-of-year */
	split = ntpcal_split_eradays(rd - 1, &leapy);
	utm->tm_year = split.hi - 1899;
	utm->tm_yday = split.lo;	/* 0-based */

	/* convert to month and mday */
	split = ntpcal_split_yeardays(split.lo, leapy);
	utm->tm_mon  = split.hi;	/* 0-based */
	utm->tm_mday = split.lo + 1;	/* 1-based */

	return leapy;
}

/*
 *---------------------------------------------------------------------
 * Take a value of seconds since midnight and split it into hhmmss in a
 * 'struct calendar'.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_daysec_to_date(
	struct calendar *jd,
	int32_t		sec
	)
{
	int32_t days;
	int   ts[3];

	days = priv_timesplit(ts, sec);
	jd->hour   = (uint8_t)ts[0];
	jd->minute = (uint8_t)ts[1];
	jd->second = (uint8_t)ts[2];

	return days;
}

/*
 *---------------------------------------------------------------------
 * Take a value of seconds since midnight and split it into hhmmss in a
 * 'struct tm'.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_daysec_to_tm(
	struct tm *utm,
	int32_t	   sec
	)
{
	int32_t days;
	int32_t ts[3];

	days = priv_timesplit(ts, sec);
	utm->tm_hour = ts[0];
	utm->tm_min  = ts[1];
	utm->tm_sec  = ts[2];

	return days;
}

/*
 *---------------------------------------------------------------------
 * take a split representation for day/second-of-day and day offset
 * and convert it to a 'struct calendar'. The seconds will be normalised
 * into the range of a day, and the day will be adjusted accordingly.
 *
 * returns >0 if the result is in a leap year, 0 if in a regular
 * year and <0 if the result did not fit into the calendar struct.
 *---------------------------------------------------------------------
 */
int
ntpcal_daysplit_to_date(
	struct calendar	   *jd,
	const ntpcal_split *ds,
	int32_t		    dof
	)
{
	dof += ntpcal_daysec_to_date(jd, ds->lo);
	return ntpcal_rd_to_date(jd, ds->hi + dof);
}

/*
 *---------------------------------------------------------------------
 * take a split representation for day/second-of-day and day offset
 * and convert it to a 'struct tm'. The seconds will be normalised
 * into the range of a day, and the day will be adjusted accordingly.
 *
 * returns 1 if the result is in a leap year and zero if in a regular
 * year.
 *---------------------------------------------------------------------
 */
int
ntpcal_daysplit_to_tm(
	struct tm	   *utm,
	const ntpcal_split *ds ,
	int32_t		    dof
	)
{
	dof += ntpcal_daysec_to_tm(utm, ds->lo);

	return ntpcal_rd_to_tm(utm, ds->hi + dof);
}

/*
 *---------------------------------------------------------------------
 * Take a UN*X time and convert to a calendar structure.
 *---------------------------------------------------------------------
 */
int
ntpcal_time_to_date(
	struct calendar	*jd,
	const vint64	*ts
	)
{
	ntpcal_split ds;

	ds = ntpcal_daysplit(ts);
	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);
	ds.hi += DAY_UNIX_STARTS;

	return ntpcal_rd_to_date(jd, ds.hi);
}


/*
 * ==================================================================
 *
 * merging composite entities
 *
 * ==================================================================
 */

/*
 *---------------------------------------------------------------------
 * Merge a number of days and a number of seconds into seconds,
 * expressed in 64 bits to avoid overflow.
 *---------------------------------------------------------------------
 */
vint64
ntpcal_dayjoin(
	int32_t days,
	int32_t secs
	)
{
	vint64 res;

#   if defined(HAVE_INT64)

	res.q_s	 = days;
	res.q_s *= SECSPERDAY;
	res.q_s += secs;

#   else

	uint32_t p1, p2;
	int	 isneg;

	/*
	 * res = days *86400 + secs, using manual 16/32 bit
	 * multiplications and shifts.
	 */
	isneg = (days < 0);
	if (isneg)
		days = -days;

	/* assemble days * 675 */
	res.D_s.lo = (days & 0xFFFF) * 675u;
	res.D_s.hi = 0;
	p1 = (days >> 16) * 675u;
	p2 = p1 >> 16;
	p1 = p1 << 16;
	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);

	/* mul by 128, using shift */
	res.D_s.hi = (res.D_s.hi << 7) | (res.D_s.lo >> 25);
	res.D_s.lo = (res.D_s.lo << 7);

	/* fix sign */
	if (isneg)
		M_NEG(res.D_s.hi, res.D_s.lo);

	/* properly add seconds */
	p2 = 0;
	if (secs < 0) {
		p1 = (uint32_t)-secs;
		M_NEG(p2, p1);
	} else {
		p1 = (uint32_t)secs;
	}
	M_ADD(res.D_s.hi, res.D_s.lo, p2, p1);

#   endif

	return res;
}

/*
 *---------------------------------------------------------------------
 * get leap years since epoch in elapsed years
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_leapyears_in_years(
	int32_t years
	)
{
	/* We use the in-out-in algorithm here, using the one's
	 * complement division trick for negative numbers. The chained
	 * division sequence by 4/25/4 gives the compiler the chance to
	 * get away with only one true division and doing shifts otherwise.
	 */

	uint32_t sflag, sum, uyear;

	sflag = int32_sflag(years);
	uyear = int32_to_uint32_2cpl(years);
	uyear ^= sflag;

	sum  = (uyear /=  4u);	/*   4yr rule --> IN  */
	sum -= (uyear /= 25u);	/* 100yr rule --> OUT */
	sum += (uyear /=  4u);	/* 400yr rule --> IN  */

	/* Thanks to the alternation of IN/OUT/IN we can do the sum
	 * directly and have a single one's complement operation
	 * here. (Only if the years are negative, of course.) Otherwise
	 * the one's complement would have to be done when
	 * adding/subtracting the terms.
	 */
	return uint32_2cpl_to_int32(sflag ^ sum);
}

/*
 *---------------------------------------------------------------------
 * Convert elapsed years in Era into elapsed days in Era.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_days_in_years(
	int32_t years
	)
{
	return years * DAYSPERYEAR + ntpcal_leapyears_in_years(years);
}

/*
 *---------------------------------------------------------------------
 * Convert a number of elapsed month in a year into elapsed days in year.
 *
 * The month will be normalized, and 'res.hi' will contain the
 * excessive years that must be considered when converting the years,
 * while 'res.lo' will contain the number of elapsed days since start
 * of the year.
 *
 * This code uses the shifted-month-approach to convert month to days,
 * because then there is no need to have explicit leap year
 * information.	 The slight disadvantage is that for most month values
 * the result is a negative value, and the year excess is one; the
 * conversion is then simply based on the start of the following year.
 *---------------------------------------------------------------------
 */
ntpcal_split
ntpcal_days_in_months(
	int32_t m
	)
{
	ntpcal_split res;

	/* Add ten months and correct if needed. (It likely is...) */
	res.lo  = m + 10;
	res.hi  = (res.lo >= 12);
	if (res.hi)
		res.lo -= 12;

	/* if still out of range, normalise by floor division ... */
	if (res.lo < 0 || res.lo >= 12) {
		uint32_t mu, Q, sflag;
		sflag = int32_sflag(res.lo);
		mu    = int32_to_uint32_2cpl(res.lo);
		Q     = sflag ^ ((sflag ^ mu) / 12u);
		res.hi += uint32_2cpl_to_int32(Q);
		res.lo  = mu - Q * 12u;
	}
	
	/* get cummulated days in year with unshift */
	res.lo = shift_month_table[res.lo] - 306;

	return res;
}

/*
 *---------------------------------------------------------------------
 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
 * days in Gregorian epoch.
 *
 * If you want to convert years and days-of-year, just give a month of
 * zero.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_edate_to_eradays(
	int32_t years,
	int32_t mons,
	int32_t mdays
	)
{
	ntpcal_split tmp;
	int32_t	     res;

	if (mons) {
		tmp = ntpcal_days_in_months(mons);
		res = ntpcal_days_in_years(years + tmp.hi) + tmp.lo;
	} else
		res = ntpcal_days_in_years(years);
	res += mdays;

	return res;
}

/*
 *---------------------------------------------------------------------
 * Convert ELAPSED years/months/days of gregorian calendar to elapsed
 * days in year.
 *
 * Note: This will give the true difference to the start of the given year,
 * even if months & days are off-scale.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_edate_to_yeardays(
	int32_t years,
	int32_t mons,
	int32_t mdays
	)
{
	ntpcal_split tmp;

	if (0 <= mons && mons < 12) {
		years += 1;
		mdays += real_month_table[is_leapyear(years)][mons];
	} else {
		tmp = ntpcal_days_in_months(mons);
		mdays += tmp.lo
		       + ntpcal_days_in_years(years + tmp.hi)
		       - ntpcal_days_in_years(years);
	}

	return mdays;
}

/*
 *---------------------------------------------------------------------
 * Convert elapsed days and the hour/minute/second information into
 * total seconds.
 *
 * If 'isvalid' is not NULL, do a range check on the time specification
 * and tell if the time input is in the normal range, permitting for a
 * single leapsecond.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_etime_to_seconds(
	int32_t hours,
	int32_t minutes,
	int32_t seconds
	)
{
	int32_t res;

	res = (hours * MINSPERHR + minutes) * SECSPERMIN + seconds;

	return res;
}

/*
 *---------------------------------------------------------------------
 * Convert the date part of a 'struct tm' (that is, year, month,
 * day-of-month) into the RD of that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_tm_to_rd(
	const struct tm *utm
	)
{
	return ntpcal_edate_to_eradays(utm->tm_year + 1899,
				       utm->tm_mon,
				       utm->tm_mday - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * Convert the date part of a 'struct calendar' (that is, year, month,
 * day-of-month) into the RD of that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_date_to_rd(
	const struct calendar *jd
	)
{
	return ntpcal_edate_to_eradays((int32_t)jd->year - 1,
				       (int32_t)jd->month - 1,
				       (int32_t)jd->monthday - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * convert a year number to rata die of year start
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_year_to_ystart(
	int32_t year
	)
{
	return ntpcal_days_in_years(year - 1) + 1;
}

/*
 *---------------------------------------------------------------------
 * For a given RD, get the RD of the associated year start,
 * that is, the RD of the last January,1st on or before that day.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_rd_to_ystart(
	int32_t rd
	)
{
	/*
	 * Rather simple exercise: split the day number into elapsed
	 * years and elapsed days, then remove the elapsed days from the
	 * input value. Nice'n sweet...
	 */
	return rd - ntpcal_split_eradays(rd - 1, NULL).lo;
}

/*
 *---------------------------------------------------------------------
 * For a given RD, get the RD of the associated month start.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_rd_to_mstart(
	int32_t rd
	)
{
	ntpcal_split split;
	int	     leaps;

	split = ntpcal_split_eradays(rd - 1, &leaps);
	split = ntpcal_split_yeardays(split.lo, leaps);

	return rd - split.lo;
}

/*
 *---------------------------------------------------------------------
 * take a 'struct calendar' and get the seconds-of-day from it.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_date_to_daysec(
	const struct calendar *jd
	)
{
	return ntpcal_etime_to_seconds(jd->hour, jd->minute,
				       jd->second);
}

/*
 *---------------------------------------------------------------------
 * take a 'struct tm' and get the seconds-of-day from it.
 *---------------------------------------------------------------------
 */
int32_t
ntpcal_tm_to_daysec(
	const struct tm *utm
	)
{
	return ntpcal_etime_to_seconds(utm->tm_hour, utm->tm_min,
				       utm->tm_sec);
}

/*
 *---------------------------------------------------------------------
 * take a 'struct calendar' and convert it to a 'time_t'
 *---------------------------------------------------------------------
 */
time_t
ntpcal_date_to_time(
	const struct calendar *jd
	)
{
	vint64  join;
	int32_t days, secs;

	days = ntpcal_date_to_rd(jd) - DAY_UNIX_STARTS;
	secs = ntpcal_date_to_daysec(jd);
	join = ntpcal_dayjoin(days, secs);

	return vint64_to_time(&join);
}


/*
 * ==================================================================
 *
 * extended and unchecked variants of caljulian/caltontp
 *
 * ==================================================================
 */
int
ntpcal_ntp64_to_date(
	struct calendar *jd,
	const vint64    *ntp
	)
{
	ntpcal_split ds;

	ds = ntpcal_daysplit(ntp);
	ds.hi += ntpcal_daysec_to_date(jd, ds.lo);

	return ntpcal_rd_to_date(jd, ds.hi + DAY_NTP_STARTS);
}

int
ntpcal_ntp_to_date(
	struct calendar *jd,
	uint32_t	 ntp,
	const time_t	*piv
	)
{
	vint64	ntp64;

	/*
	 * Unfold ntp time around current time into NTP domain. Split
	 * into days and seconds, shift days into CE domain and
	 * process the parts.
	 */
	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
	return ntpcal_ntp64_to_date(jd, &ntp64);
}


vint64
ntpcal_date_to_ntp64(
	const struct calendar *jd
	)
{
	/*
	 * Convert date to NTP. Ignore yearday, use d/m/y only.
	 */
	return ntpcal_dayjoin(ntpcal_date_to_rd(jd) - DAY_NTP_STARTS,
			      ntpcal_date_to_daysec(jd));
}


uint32_t
ntpcal_date_to_ntp(
	const struct calendar *jd
	)
{
	/*
	 * Get lower half of 64-bit NTP timestamp from date/time.
	 */
	return ntpcal_date_to_ntp64(jd).d_s.lo;
}



/*
 * ==================================================================
 *
 * day-of-week calculations
 *
 * ==================================================================
 */
/*
 * Given a RataDie and a day-of-week, calculate a RDN that is reater-than,
 * greater-or equal, closest, less-or-equal or less-than the given RDN
 * and denotes the given day-of-week
 */
int32_t
ntpcal_weekday_gt(
	int32_t rdn,
	int32_t dow
	)
{
	return ntpcal_periodic_extend(rdn+1, dow, 7);
}

int32_t
ntpcal_weekday_ge(
	int32_t rdn,
	int32_t dow
	)
{
	return ntpcal_periodic_extend(rdn, dow, 7);
}

int32_t
ntpcal_weekday_close(
	int32_t rdn,
	int32_t dow
	)
{
	return ntpcal_periodic_extend(rdn-3, dow, 7);
}

int32_t
ntpcal_weekday_le(
	int32_t rdn,
	int32_t dow
	)
{
	return ntpcal_periodic_extend(rdn, dow, -7);
}

int32_t
ntpcal_weekday_lt(
	int32_t rdn,
	int32_t dow
	)
{
	return ntpcal_periodic_extend(rdn-1, dow, -7);
}

/*
 * ==================================================================
 *
 * ISO week-calendar conversions
 *
 * The ISO8601 calendar defines a calendar of years, weeks and weekdays.
 * It is related to the Gregorian calendar, and a ISO year starts at the
 * Monday closest to Jan,1st of the corresponding Gregorian year.  A ISO
 * calendar year has always 52 or 53 weeks, and like the Grogrian
 * calendar the ISO8601 calendar repeats itself every 400 years, or
 * 146097 days, or 20871 weeks.
 *
 * While it is possible to write ISO calendar functions based on the
 * Gregorian calendar functions, the following implementation takes a
 * different approach, based directly on years and weeks.
 *
 * Analysis of the tabulated data shows that it is not possible to
 * interpolate from years to weeks over a full 400 year range; cyclic
 * shifts over 400 years do not provide a solution here. But it *is*
 * possible to interpolate over every single century of the 400-year
 * cycle. (The centennial leap year rule seems to be the culprit here.)
 *
 * It can be shown that a conversion from years to weeks can be done
 * using a linear transformation of the form
 *
 *   w = floor( y * a + b )
 *
 * where the slope a must hold to
 *
 *  52.1780821918 <= a < 52.1791044776
 *
 * and b must be chosen according to the selected slope and the number
 * of the century in a 400-year period.
 *
 * The inverse calculation can also be done in this way. Careful scaling
 * provides an unlimited set of integer coefficients a,k,b that enable
 * us to write the calulation in the form
 *
 *   w = (y * a	 + b ) / k
 *   y = (w * a' + b') / k'
 *
 * In this implementation the values of k and k' are chosen to be
 * smallest possible powers of two, so the division can be implemented
 * as shifts if the optimiser chooses to do so.
 *
 * ==================================================================
 */

/*
 * Given a number of elapsed (ISO-)years since the begin of the
 * christian era, return the number of elapsed weeks corresponding to
 * the number of years.
 */
int32_t
isocal_weeks_in_years(
	int32_t years
	)
{	
	/*
	 * use: w = (y * 53431 + b[c]) / 1024 as interpolation
	 */
	static const uint16_t bctab[4] = { 157, 449, 597, 889 };

	int32_t  cs, cw;
	uint32_t cc, ci, yu, sflag;

	sflag = int32_sflag(years);
	yu    = int32_to_uint32_2cpl(years);
	
	/* split off centuries, using floor division */
	cc  = sflag ^ ((sflag ^ yu) / 100u);
	yu -= cc * 100u;

	/* calculate century cycles shift and cycle index:
	 * Assuming a century is 5217 weeks, we have to add a cycle
	 * shift that is 3 for every 4 centuries, because 3 of the four
	 * centuries have 5218 weeks. So '(cc*3 + 1) / 4' is the actual
	 * correction, and the second century is the defective one.
	 *
	 * Needs floor division by 4, which is done with masking and
	 * shifting.
	 */
	ci = cc * 3u + 1;
	cs = uint32_2cpl_to_int32(sflag ^ ((sflag ^ ci) / 4u));
	ci = ci % 4u;
	
	/* Get weeks in century. Can use plain division here as all ops
	 * are >= 0,  and let the compiler sort out the possible
	 * optimisations.
	 */
	cw = (yu * 53431u + bctab[ci]) / 1024u;

	return uint32_2cpl_to_int32(cc) * 5217 + cs + cw;
}

/*
 * Given a number of elapsed weeks since the begin of the christian
 * era, split this number into the number of elapsed years in res.hi
 * and the excessive number of weeks in res.lo. (That is, res.lo is
 * the number of elapsed weeks in the remaining partial year.)
 */
ntpcal_split
isocal_split_eraweeks(
	int32_t weeks
	)
{
	/*
	 * use: y = (w * 157 + b[c]) / 8192 as interpolation
	 */

	static const uint16_t bctab[4] = { 85, 130, 17, 62 };

	ntpcal_split res;
	int32_t  cc, ci;
	uint32_t sw, cy, Q, sflag;

	/* Use two fast cycle-split divisions here. This is again
	 * susceptible to internal overflow, so we check the range. This
	 * still permits more than +/-20 million years, so this is
	 * likely a pure academical problem.
	 *
	 * We want to execute '(weeks * 4 + 2) /% 20871' under floor
	 * division rules in the first step.
	 */
	sflag = int32_sflag(weeks);
	sw  = uint32_saturate(int32_to_uint32_2cpl(weeks), sflag);
	sw  = 4u * sw + 2;
	Q   = sflag ^ ((sflag ^ sw) / GREGORIAN_CYCLE_WEEKS);
	sw -= Q * GREGORIAN_CYCLE_WEEKS;
	ci  = Q % 4u;
	cc  = uint32_2cpl_to_int32(Q);

	/* Split off years; sw >= 0 here! The scaled weeks in the years
	 * are scaled up by 157 afterwards.
	 */ 
	sw  = (sw / 4u) * 157u + bctab[ci];
	cy  = sw / 8192u;	/* ws >> 13 , let the compiler sort it out */
	sw  = sw % 8192u;	/* ws & 8191, let the compiler sort it out */

	/* assemble elapsed years and downscale the elapsed weeks in
	 * the year.
	 */
	res.hi = 100*cc + cy;
	res.lo = sw / 157u;

	return res;
}

/*
 * Given a second in the NTP time scale and a pivot, expand the NTP
 * time stamp around the pivot and convert into an ISO calendar time
 * stamp.
 */
int
isocal_ntp64_to_date(
	struct isodate *id,
	const vint64   *ntp
	)
{
	ntpcal_split ds;
	int32_t      ts[3];
	uint32_t     uw, ud, sflag;

	/*
	 * Split NTP time into days and seconds, shift days into CE
	 * domain and process the parts.
	 */
	ds = ntpcal_daysplit(ntp);

	/* split time part */
	ds.hi += priv_timesplit(ts, ds.lo);
	id->hour   = (uint8_t)ts[0];
	id->minute = (uint8_t)ts[1];
	id->second = (uint8_t)ts[2];

	/* split days into days and weeks, using floor division in unsigned */
	ds.hi += DAY_NTP_STARTS - 1; /* shift from NTP to RDN */
	sflag = int32_sflag(ds.hi);
	ud  = int32_to_uint32_2cpl(ds.hi);
	uw  = sflag ^ ((sflag ^ ud) / DAYSPERWEEK);
	ud -= uw * DAYSPERWEEK;
	ds.hi = uint32_2cpl_to_int32(uw);
	ds.lo = ud;

	id->weekday = (uint8_t)ds.lo + 1;	/* weekday result    */

	/* get year and week in year */
	ds = isocal_split_eraweeks(ds.hi);	/* elapsed years&week*/
	id->year = (uint16_t)ds.hi + 1;		/* shift to current  */
	id->week = (uint8_t )ds.lo + 1;

	return (ds.hi >= 0 && ds.hi < 0x0000FFFF);
}

int
isocal_ntp_to_date(
	struct isodate *id,
	uint32_t	ntp,
	const time_t   *piv
	)
{
	vint64	ntp64;

	/*
	 * Unfold ntp time around current time into NTP domain, then
	 * convert the full time stamp.
	 */
	ntp64 = ntpcal_ntp_to_ntp(ntp, piv);
	return isocal_ntp64_to_date(id, &ntp64);
}

/*
 * Convert a ISO date spec into a second in the NTP time scale,
 * properly truncated to 32 bit.
 */
vint64
isocal_date_to_ntp64(
	const struct isodate *id
	)
{
	int32_t weeks, days, secs;

	weeks = isocal_weeks_in_years((int32_t)id->year - 1)
	      + (int32_t)id->week - 1;
	days = weeks * 7 + (int32_t)id->weekday;
	/* days is RDN of ISO date now */
	secs = ntpcal_etime_to_seconds(id->hour, id->minute, id->second);

	return ntpcal_dayjoin(days - DAY_NTP_STARTS, secs);
}

uint32_t
isocal_date_to_ntp(
	const struct isodate *id
	)
{
	/*
	 * Get lower half of 64-bit NTP timestamp from date/time.
	 */
	return isocal_date_to_ntp64(id).d_s.lo;
}

/* -*-EOF-*- */
OpenPOWER on IntegriCloud