Euclidean Affine Functions and Applications to Calendar Algorithms
aa r X i v : . [ c s . D S ] F e b Euclidean Affine Functions and Applications to CalendarAlgorithms ∗ Cassio Neri † Lorenz Schneider ‡ February 16, 2021
Abstract
We study properties of Euclidean affine functions (EAFs), namely those of the form f ( r ) = ( α · r + β ) /δ ,and their closely related expression ˚ f ( r ) = ( α · r + β )% δ , where r , α , β and δ are integers, and where / and % respectively denote the quotient and remainder of Euclidean division. We derive algebraic relations andnumerical approximations that are important for the efficient evaluation of these expressions in modernCPUs. Since simple division and remainder are particular cases of EAFs (when α = 1 and β = 0 ), theoptimisations proposed in this paper can also be appplied to them. Such expressions appear in some ofthe most common tasks in any computer system, such as printing numbers, times and dates. We usecalendar calculations as the main application example because it is richer with respect to the number ofEAFs employed. Specifically, the main application presented in this article relates to Gregorian calendaralgorithms. We will show how they can be implemented substantially more efficiently than is currentlythe case in widely used C, C++, C Keywords:
Euclidean Affine Functions, Integer Division, Calendar Algorithms
Mathematics Subject Classification (2020):
Divisions by constants appear in some of the most common tasks performed by software systems, includingprinting decimal numbers (division by powers of ) and working with times (division by , and ).Since division is the slowest of the four basic arithmetical operations, various authors [1, 5, 10, 16, 22] haveproposed strength reduction optimisations ( i.e., the replacement of instructions with alternatives that aremathematically equivalent but faster) for integer divisions when divisors are constants known by the compiler.The algorithms proposed by Granlund and Montgomery [10] have been implemented by major compilers.A typical example that will appear later in this article is the division n/ · n/ , ∀ n ∈ [0 ,
28 825 529[ . ∗ Submitted to the editors February 16, 2021 † JPMorgan Chase & Co, [email protected] (Opinions expressed in this paper are those of the authors, and do notnecessarily reflect the view of JP Morgan.) ‡ EMLYON Business School, [email protected] . , the multiplier 2 939 745, the exponent , and theboundary
28 825 529 of the interval in which equality holds. Clearly, there are many possible solutions.Remainder calculation is a closely related problem, also arising in the tasks mentioned above, although thecited works do not consider optimisations for this operation. Compilers are content to apply strength reductionto obtain the quotient r/δ and the remainder evaluating r − δ · ( r/δ ) . For this reason, [13, 23] consideredthe problem of directly obtaining the remainder without first calculating the quotient.This paper expands previous works in two ways. Firstly, it considers the more general setting of Euclideanaffine functions (EAFs) f ( r ) = ( α · r + β ) /δ where α , β and δ are integer constants and r is an integer variable.Secondly, it suggests optimisations that are even more effective in applications where both the quotient andthe remainder need to be evaluated.We derive EAF-related equalities that provide alternative ways of evaluating expressions commonly used inapplications. (For instance, r − δ · ( r/δ ) , which is used to obtain the remainder as explained above.) Theseequalities underpin optimisations, other than strength reduction, that take into account aspects of modernCPUs. Specifically, they foster instruction-level parallelism implemented by superscalar processors and theyprofit from the backward compatibility features that drove the design of the x86 64 instruction set.Calendar calculations are a rich application field for our EAF results. Performed by various software systems,they tackle questions such as: What is the date today? For how many days do interest rates accrue over theperiod of a loan? Your mobile phone no doubt performed a similar calculation while you read this paragraph.Our algorithms are substantially faster than those in widely used open source implementations, as shown bybenchmarks against counterparts in glibc [7, 8] (Linux contains a similar implementation [14]), Boost [4],libc++ [15], .NET [17] and OpenJDK [19] (Android contains the same implementations [2]). Our algorithmsare also faster than others found in the academic literature [3, 6, 11, 12, 20, 21].The critical point setting apart implementations of calendar algorithms is how they deal with non-linearitiescaused by irregular month and year lengths. Some authors [14, 8, 17] resort to look-up tables, which canbe costly when the L1-level cache is cold. They conduct linear searches on these tables, entailing branchingthat can cause stalls in the processor’s execution pipeline. Others [3, 6, 11, 12, 4, 15, 20, 21] tackle the issueentirely through EAFs. Nevertheless, they do not go far enough and do not use the mathematical propertiesderived in this paper. To the best of our knowledge, the most successful attempt at achieving this efficiencyis Baum [3]. This author provides some explanations, but appears to resort to trial and error when it comesto pre-calculating certain magic numbers used in his algorithm. We go further by providing a systematic andgeneral framework for such calculations.Our paper and its main contributions are organised as follows.Section 2 introduces EAFs and derives some of their properties. Although they are common in applications,we were unable to find any systematic coverage of them. Hatcher [11] and Richards [21] appear to be awareof some of the results of Theorem 2.5, but they do not point to any proof.Section 3 concerns efficient evaluations of EAFs. Theorems 3.2 and 3.4 generalise prior-known results on di-vision to EAFs. Section 3.1 revisits those results and brings geometric insights to the problem. Theorem 3.11concerns the efficient evaluation of residuals: they are to EAFs what remainders are to division. The opti-misation proposed by Theorem 3.11 is fundamentally different from other optimisations, both in the presentpaper and in prior works, since it does not involve strength reduction. Indeed, this theorem shows how tobreak a data dependency present in the instructions currently emitted by compilers, enabling instruction-levelparallelism, as we shall see in Example 3.12.Section 4 provides a generic mathematical framework for calendars, which we use to study the Gregoriancalendar. Efficient algorithms for this calendar are derived in Sections 5 and 6. Section 7 presents performanceanalysis. 2e finish this introduction by setting the notation used throughout and recalling some well-known results.On the totally ordered Euclidean domain of integer numbers ( Z , + , · , ≤ ) , forward slash / and percent % respectively denote the quotient and remainder of Euclidean division. More precisely, given n, δ ∈ Z with δ = 0 , there exist unique q, r ∈ Z , with ≤ r < | δ | , such that n = q · δ + r . Then n/δ = q and n % δ = r .(These concepts, for non-negative operands, match the usual operators / and % of many programminglanguages in the C-family.)We follow the usual algebraic order of operation rules and, in accordance with C-style languages, we give % the same precedence as · and / . Hence, a · b % c = ( a · b )% c , a % b · c = ( a % b ) · c and a/b % c = ( a/b )% c , a % b/c = ( a % b ) /c . Moreover, a + b % c = a + ( b % c ) , and a − b % c = a − ( b % c ) . In general, a · b/c = a · ( b/c ) and ( a + b ) /c = a/c + b/c . This indicates that these precedence rules are even more important than whenworking in algebraic fields. Nevertheless, we shall drop unnecessary parentheses henceforth.The set of non-negative integer numbers is denoted by Z + = { x ∈ Z ; x ≥ } .On the totally ordered field of rational numbers ( Q , + , · , ≤ ) , the reciprocal of δ ∈ Q \ { } is denoted by δ − .For n , δ ∈ Z with δ = 0 , we shall write n · δ − ∈ Q to distinguish it from n/δ ∈ Z . If the former belongs to Z , then n · δ − = n/δ . This distinction is particularly important when dealing with inequalities. On the onehand, for x, y, δ ∈ Z with δ > , we have x · δ − ≥ y · δ − if, and only if, x ≥ y ; and on the another hand x/δ ≥ y/δ if x ≥ y , but the converse does not hold.For any x ∈ Q , there exist unique ⌊ x ⌋ ∈ Z and ⌈ x ⌉ ∈ Z such that ⌊ x ⌋ ≤ x < ⌊ x ⌋ + 1 and ⌈ x ⌉ − < x ≤ ⌈ x ⌉ .For any finite set X , the number of elements of X is denoted by X . Multiplication by converts hours to minutes. Conversely, division by converts minutes to hours. Ifthe amount of minutes in an hour were variable, these problems would be more complex. Since monthsand years have variable numbers of days, such complexities arise in calendar calculations. Different waysof tackling this problem appear across implementations with various degrees of performance. Nevertheless,similarities between calendar calculations and the previous linear problems are stronger than they might appear,as revealed by our study of EAFs. Definition 2.1.
A function f : Z → Z is a Euclidean affine function , EAF for short, if it has the form f ( r ) = ( α · r + β ) /δ for all r ∈ Z and fixed α , β , δ ∈ Z with δ = 0 . This terminology is ours. The analogues of EAFs in higher dimensions appear in Discrete Geometry andare called quasi-affine transformations. That area focuses on periodicity, tiling and other geometric aspects,whereas we are concerned with efficient calculations. We therefore use the term EAF to distinguish betweenthese two approaches.
Definition 2.2.
Let f be the EAF f ( r ) = ( α · r + β ) /δ . The residual function of f is ˚ f ( r ) = ( α · r + β )% δ or, equivalently, ˚ f ( r ) = α · r + β − δ · f ( r ) . An important function related to f ( r ) = r/δ is its right inverse ˆ f ( q ) = δ · q . If δ > and r ∈ [ ˆ f ( q ) , ˆ f ( q +1)[ =[ δ · q, δ · ( q +1)[ , then f ( r ) = q and ˚ f ( r ) = r − δ · ( r/δ ) = r − ˆ f ( f ( r )) = r − ˆ f ( q ) . As we shall see, Theorem 2.5generalises these results and supports the following terminology. Definition 2.3.
Let f be the EAF f ( r ) = ( α · r + β ) /δ , with δ ≥ α > . The minimal right inverse of f is the EAF ˆ f ( q ) = ( δ · q + α − β − /α . emma 2.4. Let f be the EAF f ( r ) = ( α · r + β ) /δ , with δ > , and g : Z → Z . For any q ∈ Z such that f ( g ( q ) − < f ( g ( q )) , we have: r ∈ Z and f ( r ) = f ( g ( q )) = ⇒ ˚ f ( r ) /α = r − g ( q ) . Proof.
Since ˚ f ( r ) = α · r + β − δ · f ( r ) , we have ˚ f ( g ( q )) = α · g ( q ) + β − δ · f ( g ( q )) . Hence, ˚ f ( r ) = α · r + β − δ · f ( r ) = α · ( r − g ( q )) + α · g ( q ) + β − δ · f ( g ( q ))= α · ( r − g ( q )) + ˚ f ( g ( q )) . We will show that ≤ ˚ f ( g ( q )) < α and it will follow that ˚ f ( r ) /α = r − g ( q ) . By definition of ˚ f and % weobtain ˚ f ( g ( q )) ≥ . Suppose, by contradiction, that ˚ f ( g ( q )) ≥ α . Then, α · ( g ( q ) −
1) + β = α · g ( q ) + β − α = δ · f ( g ( q )) + ˚ f ( g ( q )) − α ≥ δ · f ( g ( q )) . Dividing the above by δ gives f ( g ( q ) − ≥ f ( g ( q )) , which contradicts the assumption on q . Theorem 2.5.
Let f be the EAF f ( r ) = ( α · r + β ) /δ with δ ≥ α > . Then, for any r and q ∈ Z we have:(i). f ( ˆ f ( q )) = q and f ( ˆ f ( q ) −
1) = q − ;(ii). f ( r ) = q if, and only if, r ∈ [ ˆ f ( q ) , ˆ f ( q + 1)[ ;(iii). r ∈ [ ˆ f ( f ( r )) , ˆ f ( f ( r ) + 1)[ and ˚ f ( r ) /α = r − ˆ f ( f ( r )) .Proof. Let q ∈ Z and set n = δ · q + α − β − so that ˆ f ( q ) = ( δ · q + α − β − /α = n/α .(i): Since ≤ n % α ≤ α − and α ≤ δ , we have: ≤ α − − n % α ≤ δ − − n % α < δ. (1)From ˆ f ( q ) = n/α , we obtain α · ˆ f ( q ) = n − n % α and α · ˆ f ( q ) + β = n + β − n % α = δ · q + α − − n % α. This and Equation (1) yield ( α · ˆ f ( q ) + β ) /δ = q , that is, f ( ˆ f ( q )) = q . Furthermore, the equality above gives: α · ( ˆ f ( q ) −
1) + β = δ · q − − n % α = δ · ( q −
1) + δ − − n % α. This and Equation (1) yield ( α · ( ˆ f ( q ) −
1) + β ) /δ = q − , that is, f ( ˆ f ( q ) −
1) = q − .(ii): Since δ ≥ α > , f is non-decreasing. If r ≤ ˆ f ( q ) − , then f ( r ) ≤ f ( ˆ f ( q ) −
1) = q − . If ˆ f ( q ) ≤ r ≤ ˆ f ( q + 1) − , then q = f ( ˆ f ( q )) ≤ f ( r ) ≤ f ( ˆ f ( q + 1) −
1) = q (the last equality follows from (i)applied to q + 1 ) and thus, f ( r ) = q . Finally, if r ≥ ˆ f ( q + 1) , then f ( r ) ≥ f ( ˆ f ( q + 1)) = q + 1 (again from(i) applied to q + 1 ).(iii): Since δ ≥ α > , ˆ f is strictly increasing and there thus exists p ∈ Z such that r ∈ [ ˆ f ( p ) , ˆ f ( p + 1)[ . Using(i) and (ii) (for q = p ) gives f ( ˆ f ( p )) = p , f ( ˆ f ( p ) −
1) = p − and f ( r ) = p . In particular, f ( ˆ f ( p ) − Theorem 2.5 for f ( r ) = (4 · r + 3) / and ˆ f ( q ) = 146097 · q/ gives: r ∈ [ ˆ f ( q ) , ˆ f ( q + 1)[ = ⇒ (4 · r + 3) / q and (4 · r + 3)%146097 / r − ˆ f ( q ) . Each side of last equality provides a way to evaluate the same quantity. To analyse performance trade-offs,we introduce auxiliary variables referring to partial results and compare both methods side-by-side:Calculation of (4 · r + 3)%146097 / n = 4 · r + 3 ,q = n/ ,u = n − · q, ( i.e., u = n %146097) r = u/ , Calculation of r − ˆ f ( q ) : n = 4 · r + 3 ,q = n/ ,v = 146097 · q/ i.e., v = ˆ f ( q )) ,r = r − v. For fairness of comparison, u is set to n − · q rather than n %146097 , since the former better representsthe CPU instructions typically emitted by compilers to evaluate the latter. Both listings show exactly thesame operations, the only difference being the order of the last two operations. Our preference for the leftlisting will be explained in the next steps of the Gregorian calendar algorithm. Indeed, the expressions aboveappear in these calculations and are followed by the evaluation of · r + 3 . If r is calculated as shown on theleft, then · r + 3 = 4 · ( u/ 4) + 3 and smart compilers reduce this expression to u | , where | denotes bitwise or . Example 2.7. Theorem 2.5 for f ( r ) = (4 · r + 3) / and ˆ f ( q ) = 1461 · q/ gives: r ∈ [ ˆ f ( q ) , ˆ f ( q + 1)[ = ⇒ (4 · r + 3) / q and (4 · r + 3)%1461 / r − ˆ f ( q ) . Example 2.8. Theorem 2.5 for f ( r ) = (5 · r + 461) / and ˆ f ( m ) = (153 · m − / gives: r ∈ [ ˆ f ( m ) , ˆ f ( m + 1)[ = ⇒ (5 · r + 461) / 153 = m and (5 · r + 461)%153 / r − ˆ f ( m ) . This section covers optimisations for ( α · r + β ) /δ . The particular case α = 1 and β = 0 has been consideredby many authors [1, 5, 10, 16, 22], who have derived strength reductions whereby r/δ is reduced to amultiplication and cheaper operations. Major compilers implement the algorithms of [10]. Faster algorithmsexist but compilers are not able to use them. For instance, r might be restricted to a small interval but,unaware of this fact, the compiler must assume that r can take any allowed value for its type, usually in aninterval of form [0 , w [ . We change focus and instead of looking for an algorithm on a given interval, westart with the best algorithm we know and search the largest interval on which it can be applied. If such aninterval is satisfactory for our application, then we use this algorithm.When δ is a power of two, the evaluation of r/δ reduces to a bitwise shift, which is very cheap in binaryCPUs. The Fundamental Theorem of Arithmetic implies that δ is a power of two if, and only if, k % δ = 0 for some k ∈ Z + . For ( α · r + β ) /δ , a trivial reduction is also available when k · α % δ = 0 . Indeed, we have k · α = δ · (2 k · α/δ ) , and setting α ′ = 2 k · α/δ and β ′ = 2 k · β/δ gives k · α = δ · α ′ and ( α · r + β ) /δ = (cid:0) k · α · r + 2 k · β (cid:1) / (cid:0) k · δ (cid:1) = (cid:0) δ · α ′ · r + 2 k · β (cid:1) / (cid:0) k · δ (cid:1) = ( α ′ · r + β ′ ) / k . ( α · r + β ) /δ we multiply r by α · δ − and add the result to β · δ − . Assume, forthe time being, that β = 0 . We take an approximation α ′ ∈ Z of k · α · δ − , where k ∈ Z + is carefullychosen, and evaluate α ′ · r/ k . Two natural choices for α ′ are (cid:6) k · α · δ − (cid:7) and (cid:4) k · α · δ − (cid:5) , covered byTheorems 3.2 and 3.4, respectively. Lemma 3.1. Let f be the EAF f ( r ) = ( α · r + β ) /δ . Then, α · ( r/δ ) = f ( r ) − f ( r % δ ) , ∀ r ∈ Z . Proof. f ( r ) = ( α · ( δ · ( r/δ ) + r % δ ) + β ) /δ = α · ( r/δ ) + ( α · ( r % δ ) + β ) /δ = α · ( r/δ ) + f ( r % δ ) .The next theorem does not assume k · α % δ = 0 , but when equality holds a simpler reduction can bemade as seen above. The result becomes more interesting when we do have k · α % δ = 0 , in which case (cid:6) k · α · δ − (cid:7) = 2 k · α/δ + 1 . Theorem 3.2. Let k ∈ Z + and f be the EAF f ( r ) = ( α · r + β ) /δ with δ > . Set α ′ = 2 k · α/δ + 1 , ε = δ − k · α % δ and β ′ = − min { α ′ · r − k · f ( r ) ; r ∈ [0 , δ [ } . For r ∈ [0 , δ [ define: q ( r ) = min { p ∈ Z + ; ε · p + α ′ · r + β ′ − k · f ( r ) ≥ k } and M ( r ) = δ · q ( r ) + r. Let N = min { M ( r ) ; r ∈ [0 , δ [ } . Then, ( α · r + β ) /δ = ( α ′ · r + β ′ ) / k , ∀ r ∈ [0 , N [ . Proof. First note that q ( r ) is well-defined since ε > . So are M ( r ) and N . Also, δ > implies M ( r ) ≥ and, consequently, N ≥ . Nevertheless, we do not exclude the possibility that N = 0 , in which case thistheorem’s conclusion is vacuously true. The sequel assumes that N > and r ∈ Z . Note that: ε = δ − (cid:0) k · α − δ · (cid:0) k · α/δ (cid:1)(cid:1) = δ · (cid:0) k · α/δ (cid:1) − k · α = δ · α ′ − k · α. (2)It follows that δ · α ′ · ( r/δ ) − k · α · ( r/δ ) = ε · ( r/δ ) . Hence, α ′ · r − k · α · ( r/δ ) = α ′ · ( δ · ( r/δ ) + r % δ ) − k · α · ( r/δ ) = ε · ( r/δ ) + α ′ · ( r % δ ) . From the above and Lemma 3.1 we conclude that: α ′ · r − k · f ( r ) = ε · ( r/δ ) + α ′ · ( r % δ ) − k · f ( r % δ ) . (3)Suppose further that r ∈ [0 , N [ . By definition, N ≤ M ( r % δ ) and hence r < M ( r % δ ) , i.e. , δ · ( r/δ ) + r % δ <δ · q ( r % δ ) + r % δ . It follows that r/δ < q ( r % δ ) . From the definition of q ( r % δ ) and the fact that r/δ ≥ ,we obtain: ε · ( r/δ ) + α ′ · ( r % δ ) + β ′ − k · f ( r % δ ) < k . (4)The definition of β ′ yields α ′ · ( r % δ ) + β ′ − k · f ( r % δ ) ≥ . This and ε · ( r/δ ) ≥ mean that the leftside of Equation (4) is non-negative. Therefore, Equation (3) yields ≤ α ′ · r + β ′ − k · f ( r ) < k and theconclusion follows from Euclidean division by k . Example 3.3. Theorem 3.2 for k = 5 and f ( m ) = (153 · m − / yields α ′ = 980 , β ′ = − and N = 12 , that is, (153 · m − / · m − / , ∀ m ∈ [0 , . The next theorem does assume that k · α % δ = 0 and thus, k · α/δ = (cid:4) k · α · δ − (cid:5) .6 heorem 3.4. Let k ∈ Z + and f be the EAF f ( r ) = ( α · r + β ) /δ with δ > and k · α % δ > . Set α ′ = 2 k · α/δ , ε = 2 k · α % δ and β ′ = min { k − − ( α ′ · r − k · f ( r )) ; r ∈ [0 , δ [ } . For r ∈ [0 , δ [ define: q ( r ) = min { p ∈ Z + ; − ε · p + α ′ · r + β ′ − k · f ( r ) < } and M ( r ) = δ · q ( r ) + r. Let N = min { M ( r ) ; r ∈ [0 , δ [ } . Then, ( α · r + β ) /δ = ( α ′ · r + β ′ ) / k , ∀ r ∈ [0 , N [ . Proof. The assumption k · α % δ > means that ε > and ensures that q ( r ) is well-defined, as are M ( r ) and N . Also, δ > implies M ( r ) ≥ and, consequently, N ≥ . Nevertheless, we do not exclude the possibilitythat N = 0 , in which case this theorem’s conclusion is vacuously true. The sequel assumes that N > and r ∈ Z .Note that: ε = 2 k · α % δ = 2 k · α − δ · (cid:0) k · α/δ (cid:1) = 2 k · α − δ · α ′ . (5)It follows that k · α · ( r/δ ) − δ · α ′ · ( r/δ ) = ε · ( r/δ ) . Hence, α ′ · r − k · α · ( r/δ ) = α ′ · ( δ · ( r/δ ) + r % δ ) − k · α · ( r/δ ) = − ε · ( r/δ ) + α ′ · ( r % δ ) . From the above and Lemma 3.1 we conclude that: α ′ · r − k · f ( r ) = − ε · ( r/δ ) + α ′ · ( r % δ ) − k · f ( r % δ ) . (6)Suppose further that r ∈ [0 , N [ . By definition N ≤ M ( r % δ ) and hence, r < M ( n % δ ) , i.e. , δ · ( r/δ ) + r % δ <δ · q ( r % δ ) + r % δ . Therefore, r/δ < q ( r % δ ) . From the definition of q ( r % δ ) and the fact that r/δ ≥ , weobtain: − ε · ( r/δ ) + α ′ · ( r % δ ) + β ′ − k · f ( r % δ ) ≥ . (7)The definition of β ′ yields α ′ · ( r % δ ) + β ′ − k · f ( r % δ ) ≤ k − < k . This and ε · ( r/δ ) ≥ mean thatthe left side of Equation (7) is less than k . Therefore, Equation (6) yields ≤ α ′ · r + β ′ − k · f ( r ) < k and the conclusion follows from Euclidean division by k . Example 3.5. Theorem 3.4 for k = 5 and f ( m ) = (153 · m − / yields α ′ = 979 , β ′ = − and N = 34 , namely: (153 · m − / · m − / , ∀ m ∈ [0 , . (8)Examples 3.3 and 3.5 give two fast alternatives for the same EAF. In this case, the latter has the advantageof being valid on a larger range, [0 , as opposed to [0 , , but in other cases the opposite is true. ForEAFs known prior to compilation, one can find both alternatives and select the one that is valid in the largerrange. Compilers can do the same for EAFs found in source code. However, if the costs of finding the twoalternatives need to be avoided, then a simple and fast-selecting criterion is based on the following heuristics.To maximise N we seek to maximise q ( r ) , which is the minimum value of p for which ε · p reaches a certainthreshold. The smaller ε is, the larger p must be for this to happen. Hence, we choose the alternative withthe smallest ε amongst its two possible values, namely, ε = δ − k · α % δ and ε = 2 k · α % δ . (If δ is odd, then ε = ε .) Another interpretation of this criterion is that it sets α ′ to the best approximation of k · α · δ − given by (cid:6) k · α · δ − (cid:7) or (cid:4) k · α · δ − (cid:5) . Indeed, if k · α % δ = 0 , then Equations (2) and (5) give: ε = δ · (cid:0) k · α/δ + 1 (cid:1) − k · α = δ · (cid:6) k · α · δ − (cid:7) − k · α,ε = 2 k · α − δ · (cid:0) k · α/δ (cid:1) = 2 k · α − δ · (cid:4) k · α · δ − (cid:5) . It follows that ε < ε if, and only if, (cid:6) k · α · δ − (cid:7) − k · α · δ − < k · α · δ − − (cid:4) k · α · δ − (cid:5) .7 xample 3.6. Theorem 3.4 for k = 16 and f ( r ) = (5 · r + 461) / yields α ′ = 2141 , β ′ = 197913 and N = 734 , namely: (5 · r + 461) / 153 = (2141 · r + 197913) / , ∀ r ∈ [0 , . (9)In Theorems 3.2 and 3.4, α ′ and ε are obtained in O (1) operations and β ′ is found by an O ( δ ) search. Foreach r ∈ [0 , δ [ , q ( r ) and N ( r ) are calculated in O (1) operations. N and the overall search for a fast EAF thenhave O ( δ ) time complexity. Since many EAFs featuring in calendar calculations have small divisors, findingmore efficient alternatives for them is reasonably fast. α = 1 , β = 0 This section examines r/δ with δ > , i.e. , the particular case of the EAF f ( r ) = ( α · r + β ) /δ where α = 1 , β = 0 and δ > .For α = 1 , Theorems 3.2 and 3.4 set α ′ to (cid:6) k · δ − (cid:7) and (cid:4) k · δ − (cid:5) , respectively. The former is the choicetaken by [1, 5, 10] and the latter by [16]. Finally, [22] and an appendix to [5], consider both and, in thisparticular case, rigourously justify the heuristics we have suggested for choosing between the two approaches.This section follows a more direct path but most of its results can be obtained from the above for α = 1 and β = 0 . For instance, for α ′ > , let β ′ = − min { α ′ · r − k · f ( r ) ; r ∈ [0 , δ [ } as in Theorem 3.2. Since f ( r ) = r/δ = 0 for r ∈ [0 , δ [ , we have β ′ = − min { α ′ · r ; r ∈ [0 , δ [ } = 0 . Hence, in contrast to the generalcase, there is no need for an O ( δ ) search to obtain β ′ . Similarly, β ′ = min { k − − ( α ′ · r − k · f ( r )) ; r ∈ [0 , δ [ } , as defined by Theorem 3.4, can be proven to be α ′ + (2 k % δ − , the same value found in [16].Theorem 3.4 assumes that k % δ ≥ and, in the definition of β ′ , had we subtracted k % δ instead of theproof would still work but β ′ would have a smaller value, namely, β ′ = α ′ . In this case, the final reductionwould be r/δ = α ′ · ( r + 1) / k as found in [5, 22]. In addition, by making β ′ smaller, q ( r ) , M ( r ) and N alsodecrease. Hence, Theorem 3.4 obtains a range of validity that is no smaller than the one obtained in [5, 22].The remainder of this section focuses on the round-up approach but the round-down alternative could besimilarly considered. Our reasons for this are that the round-up approach is mostly used by compilers andthat the appearance of the term r + 1 can potentially lead to an overflow. Lemma 3.7. Let δ, k ∈ Z + with δ > . Set α ′ = 2 k /δ + 1 and ε = δ − k % δ . Then, α ′ · δ = 2 k + ε and forany n ∈ Z we have: α ′ · n/ k = n/δ if, and only if, ≤ ε · ( n/δ ) + α ′ · ( n % δ ) < k . (10) Proof. Taking α = 1 in Theorem 3.2 gives the same α ′ and ε as here. Equation (2) then reads ε = α ′ · δ − k .Similarly, from Equation (3) with f ( n ) = n/δ and f ( n % δ ) = n % δ/δ = 0 we obtain: α ′ · n − k · ( n/δ ) = ε · ( n/δ ) + α ′ · ( n % δ ) . The result follows from Euclidean division by k .Since α ′ , ε ∈ Z + , we have ≤ ε · ( n/δ ) + α ′ · ( n % δ ) for all n ∈ Z + . There is an illuminating geometricinterpretation for the second part of the double inequality in Equation (10) that follows from mapping each n ∈ Z + to P n = ( n/δ, n % δ ) in the xy plane. Figure 1 shows, for δ = 5 , some points of particular interestlabelled by their circled value. Notably, the origin P = (0 / , , P = (1 / , and P = (5 / , .Inequality ε · ( n/δ ) + α ′ · ( n % δ ) < k states that P n is below the line ε · x + α ′ · y = 2 k , which is representedby the dotted line in Figure 1, for k = 8 , ε = 4 and α ′ = 52 . Hence, Equation (10) means that a point below8 %d n/d Figure 1: The geometry of replacing n/δ with α ′ · n/ k , where δ = 5 and k = 8 , α ′ = 2 k /δ + 1 = 52 .the line corresponds to n ∈ Z + such that n/δ = α ′ · n/ k and a point above or on the line is related to n such that n/δ = α ′ · n/ k .If the slope of the dotted line is not too steep (more precisely, − ε · α ′− ≥ − , then the graph makes itobvious that the smallest N for which P N is above or on the dotted line must lie on the line y = δ − .In our case N = 64 . For Theorem 3.2, with α = 1 and β = 0 , this means that N = M ( δ − . In particular,we are not interested in either M ( r ) or q ( r ) when r = δ − . This and β ′ = 0 turn finding a fast EAF andits interval of applicability [0 , N [ into an O (1) calculation rather than an O ( δ ) search as in the general case.These geometric ideas are present although algebraically disguised in the proof of Theorem 3.8.The result of Theorem 3.8 also appears in [5]. In addition to providing the aforementioned geometric insightsand an arguably simpler proof, we favour faster algorithms over range of applicability. In other words, wereduce division to multiplication and bitwise shift and find the largest N for which this optimisation yieldscorrect results for all dividends in [0 , N [ . (In [5] N is called critical value and is denoted by N cr .) Theorem 3.8. Let δ, k ∈ Z + with δ > . Set α ′ = 2 k /δ + 1 , ε = δ − k % δ and N = (cid:6) α ′ · ε − (cid:7) · δ − . If ε ≤ α ′ , then: n/δ = α ′ · n/ k , ∀ n ∈ [0 , N [ . Proof. Let N ′ = N − δ = (cid:0)(cid:6) α ′ · ε − (cid:7) − (cid:1) · δ + δ − . We have N ′ /δ = (cid:6) α ′ · ε − (cid:7) − and N ′ % δ = δ − .Recall from Lemma 3.7 that α ′ · δ − ε = 2 k . Therefore, using ε · (cid:6) α ′ · ε − (cid:7) < α ′ + ε we obtain: ε · ( N ′ /δ ) + α ′ · ( N ′ % δ ) = ε · (cid:6) α ′ · ε − (cid:7) − · ε + α ′ · δ − α ′ < k . (11)Let n ∈ Z + with n < N . We shall show that ≤ ε · ( n/δ ) + α ′ · ( n % δ ) < k and the result will follow fromLemma 3.7. Since ε , α ′ and r are non-negative, it is sufficient to show that ε · ( n/δ ) + α ′ · ( n % δ ) < k .Assume first that n ≤ N ′ . As a result, n/δ ≤ N ′ /δ . Then use Equation (11) to obtain: ε · ( n/δ ) + α ′ · ( n % δ ) ≤ ε · ( N ′ /δ ) + α ′ · ( δ − 1) = ε · ( N ′ /δ ) + α ′ · ( N ′ % δ ) < k . Now assume that N ′ + 1 ≤ n < N . Then (cid:0)(cid:6) α ′ · ε − (cid:7) − (cid:1) · δ ≤ n < (cid:0)(cid:6) α ′ · ε − (cid:7) − (cid:1) · δ + δ − . Hence, n/δ = (cid:6) α ′ · ε − (cid:7) − N ′ /δ + 1 = N/δ . Since n < N , we must have n % δ ≤ N % δ − δ − . Therefore, ε · ( n/δ ) + α ′ · ( n % δ ) ≤ ε · ( N ′ /δ + 1) + α ′ · ( δ − 2) = ε · ( N ′ /δ ) + α ′ · ( δ − 1) + ε − α ′ ≤ ε · ( N ′ /δ ) + α ′ · ( N ′ % δ ) . The result follows from Equation (11). Remark 3.9. If δ is not a power of two and k ≥ δ · ( δ − , then ε ≤ α ′ . Indeed, under these hypotheseswe have k /δ ≥ δ − and k % δ ≥ . It follows that δ − k % δ ≤ δ − ≤ k /δ + 1 , i.e. , ε ≤ α ′ . xample 3.10. Consider the division n / . For k = 39 the constants given by Theorem 3.8 are α ′ =376 287 347 , ε = 79 and N = 6 958 934 390 . Since ε ≤ α ′ , this theorem gives: n / · n / , ∀ n ∈ [0 , . (12) Following [10], major compilers replace the expression on the left with the one on the right when n is a32-bit unsigned integer. (It is worth mentioning that k = 39 is the smallest value for which N ≥ .)For k = 32 , the constants given by Theorem 3.8 are α ′ = 2 939 745 , ε = 149 and N = 28 825 529 . Hence, n / · n / , ∀ n ∈ [0 , 28 825 529[ . (13) Although the interval is smaller than the one for k = 39 , it is large enough for the Gregorian calendar algorithmthat we present later on. Example 3.12 will unveil the advantage of Equation (13) over Equation (12) . Theorems 3.2 and 3.4 suggest optimisations for EAFs, generalising the current practice for divisions revisitedby Theorem 3.8. Our next result extends the optimisation to residual functions. Theorem 3.11. Let f and f ′ be EAFs given by f ( r ) = ( α · r + β ) /δ and f ′ ( r ) = ( α ′ · r + β ′ ) /δ ′ , with δ ≥ α > , and assume that f ≡ f ′ on [ a, b [ . If a = ˆ f ( f ( a )) and f ′ ( a − < f ′ ( a ) , then: ˚ f ( r ) /α = ˚ f ′ ( r ) /α ′ ∀ r ∈ [ a, b [ . Proof. Let r ∈ [ a, b [ and set q = f ( r ) . From Theorem 2.5-(iii), we obtain ˚ f ( r ) /α = r − ˆ f ( q ) . We will finishthe proof by using Lemma 2.4 to show that ˚ f ′ ( r ) /α ′ = r − ˆ f ( q ) .Since δ ≥ α > , both f and ˆ f are non-decreasing. Hence, from a ≤ r ≤ b − we obtain ˆ f ( f ( a )) ≤ ˆ f ( f ( r )) ≤ ˆ f ( f ( b − . By assumption, a = ˆ f ( f ( a )) and Theorem 2.5-(iii) (for r = b − ) gives b − ∈ [ ˆ f ( f ( b − , ˆ f ( f ( b − 1) + 1)[ , in particular ˆ f ( f ( b − ≤ b − . Therefore, a ≤ ˆ f ( f ( r )) ≤ b − , in otherwords, ˆ f ( q ) ∈ [ a, b [ .Since f ≡ f ′ in [ a, b [ and ˆ f ( q ) is in this interval, we have f ′ ( ˆ f ( q )) = f ( ˆ f ( q )) . Theorem 2.5-(i) yields f ( ˆ f ( q )) = q and thus f ′ ( ˆ f ( q )) = q = f ( r ) .Similarly, it can be shown that if ˆ f ( q ) − ∈ [ a, b [ , then f ′ ( ˆ f ( q ) − 1) = q − , specifically, f ′ ( ˆ f ( q ) − Consider the EAFs f ( n ) = n / , ˆ f ( q ) = q · and f ′ ( n ) = 2 939 745 · n / .Equation (13) shows that f ( n ) = f ′ ( n ) for all n ∈ [ a, b [ = [0 , 28 825 529[ . Simple calculations give ˆ f ( f (0)) = 0 and f ′ ( − 1) = − < f ′ (0) . Hence, it follows from Equation (13) and Theorem 3.11 that: n / · n / and n %1461 = 2 939 745 · n %2 / , ∀ n ∈ [0 , 28 825 529[ . Revisiting Example 2.7 we show three alternatives that can be used to obtain q = (4 · r + 3) / and r = r − · q / side-by-side. The first evaluates r as expressed. The second uses the expression · r + 3)%1461 / seen in Example 2.7. The third is similar but replaces n / and n %1461 with theaforementioned alternatives. n = 4 · r + 3 ,q = n / ,r = r − · q / , n = 4 · r + 3 ,q = n / ,r = ( n − · q ) / , ( i.e., r = n %1461 / n = 4 · r + 3 ,u = 2 939 745 · n ,q = u / ,r = u %2 / / . We are interested in instructions emitted by compilers for the three listings above but, for clarity, we did notperform typical compiler optimisations ( e.g. , substituting division by and with bitwise shifts). However,in the middle listing we did substitute n %1461 with n − · q for greater resemblance to emittedinstructions and fairer comparison with the first listing. Compilers reduce division by to multiplicationand bitwise shift, the same operations used to get u and q in the third listing. Therefore, up to and includingthe calculation of q , the instructions for the three columns are the same and only the constants differ.The calculation of r in the first two listings requires three instructions: one multiplication, one subtractionand one bitwise shift ( / ). For the third listing, there appear to be four: one bitwise and ( %2 ), onemultiplication and one bitwise shift (compilers’ reduction of / ), and another bitwise shift ( / ). Infact, these two bitwise shifts collapse into one. Hence, there are also three instructions for the third listing.We conclude that the number of instructions is exactly the same for the three listings, but this is only partof the story. The breakthrough of the third listing is removing the dependency of r on q , seen in the firsttwo listings, which forces the CPU to wait for the calculation of q to finish before starting that of r . Inlisting three, once u is obtained, the evaluations of q and r can start concurrently. There is a small butworthwhile price to pay for this parallelisation and even that can be avoided in some platforms as we shall seenow.The calculations of q and r need u to be stored in two different registers, which implies a mov from theregister on which u was originally obtained. This does not necessarily increase the number of instructionsbecause instead of leaving the compiler to reduce division by to multiplication and a bitwise shift of bits (see Equation (12) ), we perform our own reduction (see Equation (13) ) using bits. For backwardcompatibility, x86-64 CPUs provide mov instructions that reset the upper bits of registers, effectivelyperforming a mov and a bitwise and at the same time. Hence, the operation %2 in the calculation of r might come for free. Example 3.13. Consider the EAFs f ( r ) = (5 · r + 461) / , ˆ f ( m ) = (153 · m − / , and f ′ ( r ) =(2141 · r + 197913) / . Equation (9) states that f ( r ) = f ′ ( r ) for all r ∈ [ a, b [ = [0 , . Simplecalculations give ˆ f ( f (0)) = 0 and f ′ ( − 1) = 2 < f ′ (0) . It follows from Equation (9) and Theorem 3.11that: (5 · r + 461) / 153 = (2141 · r + 197913) / and (5 · r + 461)%153 / · r + 197913)%2 / , ∀ r ∈ [0 , . Example 2.8 shows that the left side of the last equation matches d = r − (153 · m − / . Hence, as inExample 3.12, we have three alternatives for evaluating this quantity. The analysis of Example 3.12 also holdshere and once n = 2141 · r + 197913 is obtained the calculations of m = n / and d = n %2 / can start concurrently. Furthermore, in some platforms, the operation %2 might come for free. Example 3.14. In line with the previous examples, we have n/ · n/ , n %3600 = 1 193 047 · n %2 / , ∀ n ∈ [0 , n/ 60 = 71 582 789 · n/ , n %60 = 71 582 789 · n %2 / 71 582 789 , ∀ n ∈ [0 , 97 612 919[; n/ 10 = 429 496 730 · n/ , n %10 = 429 496 730 · n %2 / 429 496 730 , ∀ n ∈ [0 , . ompilers create a data dependency when evaluating the left side of the equalities above but not when theexpressions on the right are used. The first two lines can be used in conversions of the seconds elapsed sincethe start of the day, a quantity in [0 , , to hours, minutes and seconds. The third line can be used inconversions of non-negative integers up to digits into their decimal representations. α = 1 , β = 0 Again for this particular case, the EAF f ( n ) = ( α · n + β ) /δ simplifies to f ( n ) = n/δ and its residual functionsimplifies to ˚ f ( n ) = n % δ . In this section we use Theorem 3.8 and Theorem 3.11 to derive an efficient wayto calculate remainders.Formally, Theorem 3.8 states how to replace n/δ with α ′ · n/ k , where α ′ ≈ k · δ − and k ∈ Z + . Theorem 3.11then gives the equality n % δ = ( α ′ · n %2 k ) /α ′ and Theorem 3.8, again, suggests replacing division by α ′ withmultiplication by an approximation of k · α ′− and division by k . It turns out that δ ≈ k · α ′− is theapproximation we need and we obtain n % δ = δ · ( α ′ · n %2 k ) / k . This is the idea behind our next theorem. Theorem 3.15. Let δ, k ∈ Z + with δ > . Set α ′ = 2 k /δ + 1 , ε = δ − k % δ and M = (cid:6) k · ε − (cid:7) . If ε ≤ α ′ ,then: n % δ = δ · ( α ′ · n %2 k ) / k , ∀ n ∈ [0 , M [ . (14) Proof. Set N = (cid:6) α ′ · ε − (cid:7) · δ − ≥ α ′ · ε − · δ − . Recall from Lemma 3.7 that α ′ · δ = 2 k + ε and thus N ≥ (2 k + ε ) · ε − − k · ε − . It follows that N ≥ M . Theorem 3.8 gives: n/δ = α ′ · n/ k , ∀ n ∈ [0 , N [ . (15)Hence, for the EAFs f ( n ) = n/δ and f ′ ( n ) = α ′ · n/ k , Equation (15) states that f ≡ f ′ on [0 , N [ . Simplecalculations give ˆ f ( f (0)) = δ · f (0) = 0 and f ′ ( − < f ′ (0) . Theorem 3.11 (for α = 1 and [ a, b [ = [0 , N [ )therefore yields: n % δ = ( α ′ · n %2 k ) /α ′ , ∀ n ∈ [0 , N [ . (16)Let n ∈ [0 , M [ and set m = ε · ( n/δ )+ α ′ · ( n % δ ) . Since n < N , we have n/δ ≤ N/δ = (cid:6) α ′ · ε − (cid:7) − < α ′ · ε − .It follows that ≤ ε · ( n/δ ) < α ′ . Hence, m/α ′ = n % δ and m % α ′ = ε · ( n/δ ) . Moreover, Equation (15)and Lemma 3.7 give ≤ m < k . Since α ′ · n = α ′ · ( δ · ( n/δ ) + n % δ ) = 2 k · ( n/δ ) + ε · ( n/δ ) + α ′ · ( n % δ ) =2 k · ( n/δ ) + m , we conclude that α ′ · n %2 k = m .Since n < M , we have ε · n < k , otherwise n ≥ k · ε − and then, n ≥ M . Therefore: ≤ ε · ( m/α ′ ) + δ · ( m % α ′ ) = ε · ( n % δ ) + δ · ε · ( n/δ ) = ε · n < k . (17)Using α ′ · δ = 2 k + ε again gives k = δ · α ′ − ε = ( δ − · α ′ + α ′ − ε and, since ≤ α ′ − ε < α ′ , we obtain k /α ′ = δ − and k % α ′ = α ′ − ε , that is, δ = 2 k /α ′ + 1 and ε = α ′ − k % α ′ . Therefore, Equation (17) andLemma 3.7 (with α ′ and δ interchanged) give m/α ′ = δ · m/ k , namely, ( α ′ · n %2 k ) /α ′ = δ · ( α ′ · n %2 k ) / k .Finally, Equation (16) yields n % δ = ( α ′ · n %2 k ) /α ′ = δ · ( α ′ · n %2 k ) / k . Example 3.16. Revisiting Example 3.14, Theorem 3.15 gives n %3600 = 3600 · (1 193 047 · n %2 ) / , ∀ n ∈ [0 , n %60 = 60 · (71 582 789 · n %2 ) / , ∀ n ∈ [0 , 97 612 894[; n %10 = 10 · (429 496 730 · n %2 ) / , ∀ n ∈ [0 , . he expressions on the right side of the equals sign provide efficient ways of evaluating remainders. However,greater benefits are achieved when they are used in conjunction with the quotient expressions presented inExample 3.14. The equality in Equation (14) also appears in [13]. As in other works, it focuses on obtaining the value k ∈ Z + for which the equality holds on an interval of the form [0 , w [ or, in other words, w ≤ M as the nextresult shows. Corollary 3.17. Let δ, l, w ∈ Z + with < δ < w and δ − w + l % δ ≤ l . Set α ′ = 2 w + l /δ + 1 . Then, n % δ = δ · ( α ′ · n %2 w + l ) / w + l , ∀ n ∈ [0 , w [ . Proof. Set k = w + l , ε = δ − k % δ and M = (cid:6) k · ε − (cid:7) . From Theorem 3.15, it is sufficient to show that ε ≤ α ′ and k ≤ M . Since δ < w and ε ≤ l , we have δ · ε < w + l = 2 k and thus, ε ≤ k /δ < α ′ . We alsohave M ≥ k · ε − ≥ k · − l = 2 w . We begin this section with a mathematical framework that can be applied to all calendars, before specificallyturning our attention to the Gregorian calendar. Definition 4.1. A date is an ordered pair ( y, x ) ∈ Z × X and a calendar is a non-empty set of dates C ⊂ Z × X such that for all y ∈ Z the set { x ∈ X ; ( y, x ) ∈ C } is finite. Example 4.2. In the most common interpretation, y is called year and x is the day of the year . Moreover, X = Z and x = ( m, d ) is broken down into sub-coordinates, namely month m and day (of the month) d . Example 4.3. In another interpretation the coordinates of a date ( c, x ) ∈ Z × Z are called century c and day of the century x and the sub-coordinates of x = ( y, m, d ) are the year of the century y , month m and day d . We mostly use the interpretation of Example 4.2, although that of Example 4.3 appears, in passing, inPropositions 5.4 and 5.5. To ease notation, a date ( y, ( m, d )) ∈ Z × Z is identified with ( y, m, d ) ∈ Z .Under the lexicographical order of Z (considered throughout), any calendar C ⊂ Z is totally ordered andany bounded interval in C is finite. This observation supports the following definition. Definition 4.4. Let C ⊂ Z be a calendar and e ∈ C . The function ρ : C → Z given by ρ ( x ) = e, x [ , if x ≥ e , and ρ ( x ) = − x, e [ , if x < e , is called the rata die function with epoch e . The term rata die is usually [20] applied to the particular case where the epoch is 31 December 0000 in theproleptic Gregorian calendar but we shall use it regardless of the epoch.Rata die functions are strictly increasing and thus invertible on their images. We are interested in developingalgorithms to evaluate rata die functions and their inverses. This calendar is used by most countries and is familiar to most readers. Nevertheless, we recall some of itsproperties. Although it is meaningless to refer to dates in the Gregorian calendar prior to its introduction13n , we can extrapolate the calendar backwards indefinitely, yielding the so called proleptic Gregoriancalendar . (See Richards [21].) For the sake of brevity, we drop the adjective proleptic and refer to theGregorian calendar even when dealing with dates prior to .A year y ∈ Z is said to be a leap year if y %4 = 0 and y %100 = 0 or y %400 = 0 . The Gregorian calendar isdefined by G = { ( y, m, d ) ∈ Z ; m ∈ [1 , and d ∈ [1 , L ( y, m )] } , where L ( y, m ) is the length of month m of year y as given by Table 1.Table 1: Months of the Gregorian calendar. m Name L ( y, m ) ( a ) m Name L ( y, m ) m Name L ( y, m ) ( a ) L ( y, 2) = 29 if y is a leap year, or otherwise. From the definition of length of month we obtain: X m =1 L ( y, m ) = ( , if y is a leap year , , otherwise . (18) Borrowing Hatcher’s [11, 12] terminology, we introduce a computational calendar that allows efficient eval-uations of rata die functions and their inverses. This calendar is derived from G by rotating the months toplace February last and setting its minimum to e = (0 , , .Let P : Z → Z be the map defined by: P ( y , m , d ) = ( y − { m ≤ } , m + 12 · { m ≤ } , d − , ∀ ( y , m , d ) ∈ Z . Remark 4.5. It is easy to see that P is a strictly increasing bijection with P − ( y , m , d ) = ( y + { m ≥ } , m − · { m ≥ } , d + 1) , ∀ ( y , m , d ) ∈ G . (19)The computational calendar is G = { P ( x ) ; x ∈ G and P ( x ) ≥ e } . By setting e = (0 , , ∈ G and noting that P ( e ) = e , the monotonicity of P gives G = P ( G ) , where G = { x ∈ G ; x ≥ e } .Figure 2 illustrates the relationships between years and months in the Gregorian and computational calendars.Month m = 14 of year (generally, y ) of the computational calendar corresponds to m = 2 (February)of year (generally, y + 1 ) of the Gregorian calendar. Hence, if y = y + 1 is a leap year in the Gregoriancalendar, then from Equation (18), year y of the computational calendar will have days; otherwise it willhave days.If ( y , m , d ) ∈ G , then m ∈ [3 , and d ∈ [0 , L ( y , m )[ where ( y , m , d ) = P − ( y , m , d ) . (SeeTable 2.) Remark 4.6. Table 2 shows that, except for m = 14 , m L ( y , m ) is a periodic function and there are 153 = 31 + 30 + 31 + 30 + 31 days in each -month period. y GregorianComputational ...... 1 2 3 4 5 6 7 8 9 10 111213 4 5 6 7 8 9 10 1112 13141 ...... ......01 2 3 4 5 6 7 8 9 10 1112m y Figure 2: Years and in the Gregorian and computational calendars.Table 2: Months of the computational calendar and their lengths. m Name L ( y , m ) ( a ) m Name L ( y , m ) ( a ) m Name L ( y , m ) ( b ) 13 January 3114 February 28 or 29 ( c ) ( a ) y = y , m = m . ( b ) y = y + 1 , m = m − . ( c ) L ( y , m ) = 29 if y = y + 1 is a leap year, or L ( y , m ) = 28 otherwise. e , the minimum date in G , is undoubtedly the most natural epoch for a rata die function ρ on G , ρ ( x ) = e , x [ for all x ∈ G .For any x = ( y , m , d ) ∈ G , integers y , m , d and ρ ( x ) are non-negative. Hence, implementationsof G can work exclusively on unsigned integer types, which are usually faster than signed ones.Let x = ( y , m , d ) ∈ G and split [ e , x [ into three disjoint intervals: [ e , x [ = [(0 , , , ( y , , ∪ [( y , , , ( y , m , ∪ [( y , m , , ( y , m , d )[ . Hence, ρ ( x ) is the sum of the number of elements of these three intervals, which we consider separately.The year count , defined by y c ( y ) = , , , ( y , , is the number of dates prior to year y : y c ( y ) = 365 · y + { y ′ ∈ [0 , y [ ; y ′ has 366 days } Recalling that year y ′ of the computational calendar has 366 days if, and only if, y ′ + 1 is a leap year of theGregorian calendar gives: y c ( y ) = 365 · y + { y ′ ∈ [0 , y [ ; y ′ + 1 is a leap year } = 365 · y + { y ∈ [1 , y + 1[ ; y is a leap year } = 365 · y + { y ∈ [1 , y ] ; y %4 = 0 and y %100 = 0 or y %400 = 0 } = 365 · y + y / − y / 100 + y / . (20) Remark 5.1. For all y , z ∈ Z with y + 400 · z ≥ , we have y c ( y + 400 · z ) = 365 · ( y + 400 · z ) + ( y + 400 · z ) / − ( y + 400 · z ) / 100 + ( y + 400 · z ) / · y + y / − y / 100 + y / 400 + 146000 · z + 100 · z − · z + z = y c ( y ) + 146097 · z. Hence, by adding · z years to a date, its rata die increases by · z . There are dates in any -year interval. Such intervals are called leap cycles . month count defined by m c ( y , m ) = y , , , ( y , m , is the number of days between the 1stday of year y and the 1st day of month m of year y . It can be obtained by adding the lengths of previousmonths in year y . Since m = 14 is the last month of the year, its length is not included in the summationand becomes irrelevant to m c . Therefore, the year y also becomes irrelevant to m c , because the last monthis the only one whose length depends on the year. In other words, m c ( y , m ) = m c ( m ) as per Table 3.Table 3: Month count. m L ( y , m ) m c ( m ) m L ( y , m ) m c ( m ) m L ( y , m ) m c ( m ) 13 31 30614 (irrelevant) 337Some implementations [14, 8, 17] use look-up arrays to recover m c ( m ) , but it can be expressed by an EAF: m c ( m ) = (153 · m − / , ∀ m ∈ [3 , . (21)Variations of this formula have been found by many authors [3, 11, 12, 4, 15, 17, 20, 21], most often througha trial-and-error line fitting to the set of points (cid:8) ( m , m c ( m )) ∈ Z ; m ∈ [3 , (cid:9) . Remark 4.6 clarifies thematter. However, a simple textbook linear regression on the set (cid:8) ( m , m c ( m ) + 0 . ∈ Q ; m ∈ [3 , (cid:9) finds m c ( m ) = (26256 · m − / . Using this EAF in Theorem 3.4 produces the same fasteralternative seen in Example 3.5.The day count defined by d c ( y , m , d ) = y , m , , ( y , m , d )[ is the number of dates since the 1stof month m of year y . Trivially, d c ( y , m , d ) = d .Putting together the above results yields the next proposition. Proposition 5.2. Let ( y , m , d ) ∈ G . Set: y c = 365 · y + y / − y / 100 + y / , m c = (153 · m − / , d c = d . Then, ρ ( y , m , d ) = y c + m c + d c . The calculations of Proposition 5.2 can be performed more efficiently by first computing the century. Corollary 5.3. Let ( y , m , d ) ∈ G . Set: q = y / , y c = 1461 · y / − q + q / , m c = (979 · m − / , d c = d . Then, ρ ( y , m , d ) = y c + m c + d c .Proof. It suffices to show that the following alternatives for the expressions seen in Proposition 5.2 hold.(i). · y / · y + y / : This follows from · y = 4 · (365 · y + y / y %4 and ≤ y %4 < .(ii). (153 · m − / · m − / : This follows from m ∈ [3 , and Equation (8).(iii). y / 400 = q / : This follows from y / / y / .16 .1 Inverting the rata die function on the computational calendar The computational calendar G is unbounded from above and, thus, for any r ∈ Z + a unique x ∈ G existssuch that ρ ( x ) = r . The objective of this section is to find x given that r .For x = ( y , m , d ) ∈ G , the quotient q = y / is its century . Since · q ≤ y < · ( q + 1) , wehave (100 · q, , ≤ x < (100 · ( q + 1) , , . The quantity · q, , , x [ is the day of the century of x . The next proposition retrieves the century and the day of the century from a rata die. Proposition 5.4. Let x ∈ G and r be its rata die. Set: n = 4 · r + 3 , q = n / , r = n %146097 / . Then, q is the century of x and r is its day of the century.Proof. Write x = ( y , m , d ) and set q = y / and r = · q, , , x [ . We must prove that q = q and r = r .For p ∈ Z + , let g ( p ) = e , (100 · p, , . The definition of y c gives g ( p ) = y c (100 · p ) and Equation (20)yields g ( p ) = 365 · · p + 100 · p/ − · p/ 100 + 100 · p/ 400 = 36524 · p + p/ 4= 146097 · p/ . (The last equality is obtained as Item (i) in the proof of Corollary 5.3, using instead of .)Since r = e , x [ = e , (100 · q, , · q, , , x [ , we have r = g ( q ) + r ≥ g ( q ) . Now, x < (100 · ( q + 1) , , and r = e , x [ < e , (100 · ( q + 1) , , g ( q + 1) . Hence, r < g ( q + 1) ,in other words, r ∈ [ g ( q ) , g ( q + 1)[ . As seen in Example 2.6 (where ˆ f = g ): (4 · r + 3) / q and (4 · r + 3)%146097 / r − g ( q ) . These equalities mean that q = q and r = r .For x = ( y , m , d ) ∈ G , the remainder y %100 is its year of the century . We have ( y , , ≤ x < ( y + 1 , , and the quantity y , , , x [ is called the day of the year of x . The next propositionstarts where Proposition 5.4 left off, i.e. , at the day of the century, to obtain the year of the century and theday of the year. Proposition 5.5. Let x ∈ G and r be its day of the century. Set: n = 4 · r + 3 , q = n / , r = n %1461 / . Then, q is the year of the century of x and r is its day of the year.Proof. Write x = ( y , m , d ) and set q = y / , q = y %100 and r = y , , , x [ . We must provethat q = q and that r = r .For p ∈ [0 , , it is easy to show that (100 · q + p ) / 100 = 100 · q / and (100 · q + p ) / 400 = q / . Let g ( p ) = · q , , , (100 · q + p, , . The definition of y c gives g ( p ) = y c (100 · q + p ) − y c (100 · q ) and Equation (20) yields: g ( p ) = 365 · (100 · q + p ) + (100 · q + p ) / − (100 · q + p ) / 100 + (100 · q + p ) / − y c (100 · q )= 365 · p + p/ · p/ . y = 100 · q + q and r = · q , , , x [ , we have r = · q , , , ( y , , y , , , x [= · q , , , (100 · q + q, , y , , , x [ = g ( q ) + r ≥ g ( q ) . Since x < ( y +1 , , 0) = (100 · q + q +1 , , , we have r < · q , , , (100 · q + q +1 , , g ( q +1) and it follows that r ∈ [ g ( q ) , g ( q + 1)[ . As seen in Example 2.7 (where ˆ f = g ): (4 · r + 3) / q and (4 · r + 3)%1461 / r − g ( q ) . These equalities mean that q = q and that r = r .Like Proposition 5.5, the following result provides calculations for year of century and day of year, but in acomputationally more efficient way. Corollary 5.6. Let x ∈ G and r be its day of the century. Set: n = 4 · r + 3 , u = 2 939 745 · n , q = u / , r = u %2 / / . Then, q is the year of the century of x and r is its day of the year.Proof. Example 3.12 shows that the q and r set above match those of Proposition 5.5, provided that n ∈ [0 , 28 825 529[ , which we shall prove now. Proposition 5.4 provides that r = n %146097 / for some n ∈ Z . From ≤ n %146097 < , it easily follows that ≤ · r + 3 < , in other words, n ∈ [3 , .The following proposition extracts the month and day from the day of the year. Proposition 5.7. Let x ∈ G and r be its day of the year. Set: n = 5 · r + 461 , q = n / , r = n %153 / . Then, q is the month of x and r is its day.Proof. Let x = ( y , m , d ) . We must prove that q = m and that r = d .For p ∈ [3 , , let g ( p ) = (153 · p − / . Equation (21) provides that if p ∈ [3 , , then g ( p ) = m c ( p ) = y , , , ( y , p, . By definition, r = y , , , x [ can be written as: r = y , , , ( y , m , y , m , , ( y , m , d )[ = g ( m ) + d ≥ g ( m ) . We shall now prove that r < g ( m + 1) . If m < , then x < ( y, m + 1 , and we obtain r = y , , , x [ < y , , , ( y , m + 1 , m c ( m + 1) = g ( m + 1) . Now, if m = 14 , then x ≤ ( y , , L ( y , and it follows that r ≤ y , , , ( y , , L ( y , ≤ < 367 = g (15) = g ( m + 1) .Hence, r < g ( m + 1) , in other words, r ∈ [ g ( m ) , g ( m + 1)[ . As seen in Example 2.8 (where ˆ f = g ): (5 · r + 461) / 153 = m and (5 · r + 461)%153 / r − g ( m ) . These equalities mean that q = m and that r = d . Remark 5.8. In the proof of Proposition 5.7, we have shown that if r is the day of the year of ( y , m , d ) ∈ G and m < , then r ∈ [ m c ( m ) , m c ( m + 1)[ . It follows that m ≥ is equivalent to r ≥ m c (13) =306 . Corollary 5.9. Let x ∈ G and r be its day of the year. Set: n = 2141 · r + 197913 , q = n / , r = n %2 / . Then, q is the month of x and r is its day.Proof. Example 3.13 shows that the q and r set above match those of Proposition 5.7, provided that r ∈ [0 , , which we shall prove now. Proposition 5.5 provides that r = n %1461 / for some n ∈ Z .From ≤ n %1461 < , it easily follows that ≤ r ≤ . Section 5 covered evaluation of a rata die function and its inverse on the computational calendar G . Thissection adapts them to the Gregorian calendar. G allows implementations to work exclusively on unsigned integer types, which might be faster than signedones. To avoid compromising this advantage, we only consider subsets of G that have a minimum. Thesecalendars, their respective rata die functions and epochs are summarised in Table 4. For ease of reference,the first row is a reminder for G . Table 4: Gregorian-related calendars.Calendar Rata die Epoch G = P ( G ) , where G is defined below. ρ : G → Z e = (0 , , 0) = P ( e ) G = { x ∈ G ; x ≥ e } ρ : G → Z e = (0 , , G = { x ∈ G ; x ≥ e } ρ : G → Z e = ( z , , , where z is a multiple of . G (as above) ρ : G → Z e ∈ G .The rows of Table 4 show increasing degrees of freedom of how epochs and minima can be set. For G and ρ there are no choices: epoch and the minimum date are set to e = (0 , , . A match between the epochand the minimum implies that ρ is non-negative. Since e ’s year is zero, dates in G have non-negativeyears. The same holds for G , ρ and e = (0 , , . For G and ρ , epoch and minimum also match andare set to e = ( z , , , but z can be any multiple of . Therefore, ρ is also non-negative, althoughnegative years are possible when z < . Finally, ρ provides a more flexible setting where any epoch andminimum can be chosen.As we shall see, ρ = ρ ◦ P . In general, maps from one calendar to another provide a way to express ratadie functions. Conversely, rata die functions provide a way to map one calendar into another [12, 20, 21]. Proposition 6.1. Let C and C ′ be calendars and let ρ and ρ ′ be their respective rata die functions withepochs e ∈ C and e ′ ∈ C ′ . If P : C → C ′ is a strictly increasing bijection and P ( e ) = e ′ , then ρ = ρ ′ ◦ P and ρ − = P − ◦ ρ ′− on Im ( ρ ′ ) = Im ( ρ ) . Conversely, if Im ( ρ ) = Im ( ρ ′ ) , then ρ ′− ◦ ρ is a strictly increasingbijection from C to C ′ .Proof. Let x ∈ C and assume that x ≥ e . Since P is a strictly increasing bijection, we have ρ ( x ) = e, x [ = P ([ e, x [) = P ( e ) , P ( x )[) = e ′ , P ( x )[ = ρ ′ ( P ( x )) . The case x < e is treated analogously and we19onclude that ρ = ρ ′ ◦ P . From this and well-known results on functions and their inverses, we obtain ρ − = P − ◦ ρ ′− on Im ( ρ ′ ) = Im ( ρ ) .By definition, ρ ( x ) ∈ Im ( ρ ) = Im ( ρ ′ ) and ρ ′− ( ρ ( x )) ∈ C ′ , in other words, ρ ′− ◦ ρ is a map from C to C ′ .Since ρ : C → Z and ρ ′− : Im ( ρ ′ ) → C ′ are strictly increasing surjective functions, hence bijections, theircomposition is again a strictly increasing bijection. This section concerns the rata die function ρ on G = { x ∈ G ; x ≥ e } with epoch e = (0 , , .From Remark 4.5 and Proposition 6.1 applied to C = G , C ′ = G , ρ = ρ , ρ ′ = ρ , e = e , e ′ = e and P = P we obtain: ( y , m , d ) ∈ G = ⇒ P ( y , m , ρ ) ∈ G and ρ ( y , m , d ) = ρ ( P ( y , m , d )) (22)and, conversely, ( y , m , d ) ∈ G = ⇒ P − ( y , m , ρ ) ∈ G and ρ ( y , m , d ) = ρ (cid:0) P − ( y , m , d ) (cid:1) . (23) Let z ∈ Z be a multiple of and e = ( z , , . This section concerns two rata die functions: ρ on G = { x ∈ G ; x ≥ e } with epoch e , and ρ also on G but with epoch e ∈ G arbitrarily chosen.Let P : Z → Z be the map P ( y , m , d ) = ( y − z , m , d ) , which is a strictly increasing bijection with P − ( y , m , d ) = ( y + z , m , d ) and P ( e ) = (0 , , 1) = e .Let ( y , m , d ) ∈ G and set ( y , m , d ) = P ( y , m , d ) . We have m = m ∈ [1 , , d = d ∈ [1 , L ( y , m )] and, since z %400 = 0 , y − z is a leap year if, and only if, y is a leap year. Hence, L ( y , m ) = L ( y − z , m ) = L ( y , m ) and it follows that d ∈ [1 , L ( y , m )] . Therefore, ( y , m , d ) ∈ G .From ( y , m , d ) ≥ e , we obtain ( y , m , d ) ≥ e and conclude that ( y , m , d ) ∈ G . Similarly, we canshow that P − ( y , m , d ) ∈ G for all ( y , m , d ) ∈ G and thus, P is a bijection from G to G .Applying Proposition 6.1 to C = G , C ′ = G , ρ = ρ , ρ ′ = ρ , e = e , e ′ = e and P = P gives: ( y , m , d ) ∈ G = ⇒ ( y − z , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( y − z , m , d ) (24)and, reciprocally, ( y , m , d ) ∈ G = ⇒ ( y + z , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( y + z , m , d ) . (25)Fix e ∈ G and let ρ : G → Z be the rata die function with epoch e . Let x ∈ G . On the onehand, if x ≥ e , then e , x [ = e , e [ + e , x [ , in other words, ρ ( x ) = ρ ( e ) + ρ ( x ) . On theother hand, if x < e , then e , e [ = e , x [ + x , e [ , in other words, ρ ( e ) = ρ ( x ) − ρ ( x ) .Therefore, x ∈ G = ⇒ ρ ( x ) = ρ ( x ) − ρ ( e ) . (26)The next two propositions prove the correctness of our Gregorian calendar algorithms. Proposition 6.2 con-siders the calculation of the rata die function for any chosen epoch and Proposition 6.3 considers its inverse.20 roposition 6.2. Let ρ , ρ be rata die functions on G with epochs e = ( z , , and e , respectively,where z ∈ Z is a multiple of . Given ( y , m , d ) ∈ G , set: y = y − z ,m = m ,d = d , y = y − { m ≤ } ,m = m + 12 · { m ≤ } ,d = d − , q = y / ,y c = 1461 · y / − q + q / ,m c = (979 · m − / ,d c = d . Then, ρ ( y , m , d ) = y c + m c + d c − ρ ( e ) .Proof. The first column and Equation (24) give ( y , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( y , m , d ) .The second column gives ( y , m , d ) = P ( y , m , d ) and Equation (22) yields ( y , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( y , m , d ) . Hence, ρ ( y , m , d ) = ρ ( y , m , d ) , and from Corollary 5.3 we obtain ρ ( y , m , d ) = y c + m c + d c . The result follows from Equation (26). Proposition 6.3. Let ρ , ρ be rata die functions on G with epochs e = ( z , , and e ≥ e , respectively,where z ∈ Z is a multiple of . Given r ∈ Z with r ≥ − ρ ( e ) , set r = r + ρ ( e ) and: n = 4 · r + 3 ,q = n / ,r = n %146097 / , n = 4 · r + 3 ,u = 2 939 745 · n ,q = u / ,r = u %2 / / , n = 2141 · r + 197913 ,q = n / ,r = n %2 / . and: y = 100 · q + q ,m = q ,d = r , y = y + { r ≥ } ,m = m − · { r ≥ } ,d = d + 1 . y = y + z ,m = m ,d = d . (27) Then, ρ ( y , m , d ) = r .Proof. Since r ≥ − ρ ( e ) , we have r = r + ρ ( e ) ≥ and thus x ∈ G exists such that ρ ( x ) = r .Proposition 5.4 provides that q is the century of x and that r is its day of the century. Corollary 5.6 thenyields that q is the year of the century of x and r is its day of the year, while Corollary 5.9 states that q is the month of x and r is its day.Hence, the first column of Equation (27) reads x = ( y , m , d ) . From Remark 5.8 we can substitute { r ≥ } with { m ≥ } and then, by Equation (19) the second column of Equation (27) becomes ( y , m , d ) = P − ( x ) and Equation (23) gives ( y , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( x ) . Finally, Equation (25) andthe third column of Equation (27) yield ( y , m , d ) ∈ G and ρ ( y , m , d ) = ρ ( y , m , d ) = ρ ( x ) = r = r + ρ ( e ) .We have shown that r = ρ ( y , m , d ) − ρ ( e ) . Therefore, Equation (26) gives ρ ( y , m , d ) = r . We benchmarked our algorithms from Propositions 6.2 and 6.3 against counterparts in five of the most widelyused C, C++, C libc The GNU C Library [7, 8]. (The Linux Kernel contains a similar implementation [14].) Boost The Boost C++ libraries [4]. libc++ LLVM’s implementation of the C++ Standard Library [15]. .NET Microsoft .NET framework [17]. OpenJDK Oracle’s open source implementation of the Java Platform SE [19]. (Android uses the samecode [2].)We used source files as publicly available on 2 May 2020. Non-C++ implementations have been ported to thislanguage and have all been slightly modified to achieve consistent (a) function signatures; (b) storage types(for years, months, days and rata dies); and (c) epoch (Unix epoch, i.e. , 1 January 1970). Some originals dealwith date and time but our variants work on dates only. (Given the uniform durations of days, hours, minutesand seconds, it would be trivial to incorporate the time component to any dates-only algorithm. Moreover,Example 3.14 suggests improvements that are not used by major implementations.)We did not include Microsoft’s C++ Standard Library because in May 2020 it did not yet implement thesefunctionalities. For the same reason, libstdc++, the GNU implementation of the C++ Standard Library, isalso absent. Furthermore, a forthcoming release of libstdc++ is expected to implement our algorithms.We also considered our own implementations of algorithms described in the academic literature, namely, Baum[3], Fliegel and Flandern [6], Hatcher [11, 12, 21] and Reingold and Dershovitz [20].Timings were obtained with the help of the Google Benchmark library [9], to which we delegated the task ofproducing statistically relevant results. The code was run on an Intel i7-10510U CPU at . GHz. Sourcecode, available at [18], was compiled by GCC 10.2.0 at optimisation level -O3 . Results vary across platformsbut maintain qualitative consistency between them.The table in Figure 3 shows the time taken by each algorithm to evaluate ρ at 16 384 pseudo-random dates,uniformly distributed in [(1570 , , , (2370 , , ( i.e., Unix epoch ± years). They encompass the timespent scanning the array of dates (also shown). Subtracting the scanning time from that of each algorithmgives a fairer account of the time spent by the algorithm itself. The chart plots these adjusted timings relativeto ours (Proposition 6.2). ReingoldDershowitz 3.1OpenJDK 2.9glibc 2.4.NET 2.3FliegelFlandern 2.3Hatcher 1.8Boost 1.7libc++ 1.5Baum 1.0NeriSchneider Algorithm Time (ns)scanning only 3 430.3Reingold Dershowitz 76 523.4OpenJDK 69 906.5glibc 64 891.2.NET 55 179.5Hatcher 52 494.6Fliegel Flandern 53 893.4Boost 41 281.2libc++ 39 309.4Baum 35 263.8Neri Schneider 24 963.3 Figure 3: Relative and absolute timings of rata die evaluations.Similarly, the table in Figure 4 shows the time taken by each algorithm to evaluate ρ − at pseudo-random integer numbers uniformly distributed in [ − , (again, Unix epoch ± years as perRemark 5.1). The chart displays adjusted times relative to ours (Proposition 6.3).22 .8ReingoldDershowitz 7.3glibc 4.3.NET 2.4Hatcher 2.4FliegelFlandern 2.2libc++ 2.2OpenJDK 1.5Baum 1.3Boost 1.0NeriSchneider Algorithm Time (ns)scanning only 3 429.1Reingold Dershowitz 372 984.0glibc 350 648.0.NET 206 279.0Hatcher 119 239.0Fliegel Flandern 117 635.0libc++ 107 537.0OpenJDK 106 850.0Baum 76 516.2Boost 65 564.4Neri Schneider 50 770.3 Figure 4: Relative and absolute timings of rata die inverse evaluations. References [1] R. Alverson. Integer division using reciprocals. In [1991] Proceedings 10th IEEE Symposium on ComputerArithmetic , pages 186–190, 1991. doi.org/10.1109/ARITH.1991.145558.[2] Android. ojluni/src/main/java/java/time/LocalDate.java. tinyurl.com/ycsltdnq, April 2017.[3] Peter Baum. Date algorithms. tinyurl.com/y44rgx2j, 1998.[4] Boost C++ Libraries. include/boost/date time/gregorian calendar.ipp. tinyurl.com/y4buxmmf, March2020.[5] D. Cavagnino and A. E. Werbrouck. Efficient algorithms for integer division by constants using multipli-cation. The Computer Journal , 51(4):470–480, 11 2007. doi.org/10.1093/comjnl/bxm082.[6] Henry F. Fliegel and Thomas C. van Flandern. Letters to the editor: A machine algo-rithm for processing calendar dates. Communications of the ACM , 11(10):657–658, 10 1968.doi.org/10.1145/364096.364097.[7] GNU C Library. time/mktime.c. tinyurl.com/y42pkbhp, January 2020.[8] GNU C Library. time/offtime.c. tinyurl.com/yyr7uazb, January 2020.[9] Google. Benchmark v1.5.2. tinyurl.com/y29mh7q5, September 2020.[10] Torbj¨orn Granlund and Peter L. Montgomery. Division by invariant integers using multiplication. In Proceedings of the ACM SIGPLAN 1994 Conference on Programming Language Design and Imple-mentation , PLDI ’94, pages 61–72, New York, NY, USA, 1994. Association for Computing Machinery.doi.org/10.1145/178243.178249.[11] D. A. Hatcher. Simple formulae for julian day numbers and calendar dates. Quarterly Journal of theRoyal Astronomical Society , 25(1):55–53, 03 1984. tinyurl.com/y2orpwfr.[12] D. A. Hatcher. Generalized equations for julian day numbers and calendar dates. Quarterly Journal ofthe Royal Astronomical Society , 26(2):151–155, 06 1985. tinyurl.com/y6ec7t3h.[13] Daniel Lemire, Owen Kaser, and Nathan Kurz. Faster remainder by direct computation: Applica-tions to compilers and software libraries. Software: Practice and Experience , 49(6):953–970, 2019.doi.org/10.1002/spe.2689. 2314] Linux Kernel. kernel/time/timeconv.c. tinyurl.com/y6x8so24, October 2018.[15] LLVM Project. libcxx/include/chrono. tinyurl.com/yytw67zb, April 2019.[16] D. J. Magenheimer, L. Peters, K. W. Pettis, and D. Zuras. Integer multiplication and division on the hpprecision architecture. IEEE Transactions on Computers , 37(8):980–990, 1988. doi.org/10.1109/12.2248.[17] Microsoft .NET. src/libraries/System.Private.CoreLib/src/System/DateTime.cs. tinyurl.com/y4kej3mm,April 2020.[18] Cassio Neri and Lorenz Schneider. Gregorian calendar algorithms. github.com/cassioneri/calendar, 2020.[19] OpenJDK. jdk/src/java.base/share/classes/java/time/localdate.java. tinyurl.com/y92svzxw, August2019.[20] Edward M. Reingold and Nachum Dershowitz. Calendrical Calculations: The Ultimate Edition . CambridgeUniversity Press, USA, 4th edition, 2018. doi.org/10.1017/9781107415058.[21] Edward Graham Richards. Mapping Time: The CALENDAR and its HISTORY . Oxford University Press,1998. tinyurl.com/y3uegvsu.[22] A. D. Robison. N-bit unsigned division via n-bit multiply-add. In , pages 131–139. IEEE, 2005. doi.org/10.1109/ARITH.2005.31.[23] Henry S. Warren. Hacker’s Delight . Addison-Wesley Professional, 2nd edition, 2013.tinyurl.com/y23o57kr.