Formulas and algorithms for the length of a Farey sequence

This paper proves several novel formulas for the length of a Farey sequence of order n. The formulas use different trade-offs between iteration and recurrence and they range from simple to more complex. The paper also describes several iterative algorithms for computing the length of a Farey sequence based on these formulas. The algorithms are presented from the slowest to the fastest in order to explain the improvements in computational techniques from one version to another. The last algorithm in this progression runs in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(n^{2/3})$$\end{document}O(n2/3) time and uses only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\sqrt{n})$$\end{document}O(n) memory, which makes it the most efficient algorithm for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|F_n|$$\end{document}|Fn| described to date. With this algorithm we were able to compute the length of the Farey sequence of order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{18}$$\end{document}1018.


S1. ALGORITHM FOR ENUMERATING AND COUNTING THE ELEMENTS OF A FAREY SEQUENCE
For the Farey sequence F n , the mediant property implies that the numerator p and the denominator q of the third fraction in a triple of adjacent fractions ( a /b, c /d, p /q) can be expressed in terms of the values of a, b, c, and d as follows: This pair of formulas enables an elegant algorithm 11 for enumerating the elements of F n . The main idea is to use a window of adjacent fractions ( a /b, c /d) that 'slides' over the elements of F n until it reaches 1 /1. This window is initialized with the pair of fractions ( 0 /1, 1 /n) that are the first two elements of the sequence F n . Algorithm S1 uses this approach to compute the length of the Farey sequence of order n by explicitly enumerating its elements and counting them. Because the length of F n grows quadratically with n, this algorithm runs in O(n 2 ) time, which makes it too slow for large values of n. Nevertheless, for small n, it can be used to cross-check the results of the faster algorithms described in the main paper. This algorithm uses a tiny amount of memory because it only prints the sequence elements and does not need to store them in memory. In other words, its space complexity is O(1).
Algorithm S1. Enumerate all elements of the Farey sequence F n and return |F n |. Runs in O(n 2 ) time and uses O(1) memory. The main paper describes five algorithms for computing the length of a Farey sequence. To distinguish between them, the algorithms are denoted with the first five letters of the alphabet. These letters are also used as suffixes in the function name for each algorithm, e.g., FAREYLENGTHA or FAREYLENGTHE.
The algorithms are presented from the slowest to the fastest. This is done in order to explain the computational techniques and optimizations that were needed to derive the most efficient algorithm, i.e., algorithm E, which runs in O(n 2/3 ) time and uses O( √ n ) memory. Figure Figure S1: Visual summary of the time and space complexities of the five algorithms for computing the length of the Farey sequence of order n that are described in the main paper.

S3. COMPUTING EULER'S TOTIENT FUNCTION
This section defines Euler's totient function ϕ(n) and proves some of its properties. It also gives the pseudo-code for two algorithms that can be used to compute it. An example that illustrates these algorithms is also provided.
Theorem 1. For each n > 1, the length of the Farey sequence of order n can be expressed as the sum of ϕ(n) and the length of the Farey sequence of order n − 1. That is, |F n | = |F n−1 | + ϕ(n), for each n > 1.
Proof. Let D n be a set that contains the elements of the sequence F n that are not members of the sequence F n−1 , i.e., D n = p /q ∈ F n s.t. p /q ∈ F n−1 . Then, It remains to show that |D n | = ϕ(n). For each irreducible fraction p /q ∈ D n the denominator q must be equal to n, otherwise p /q would be in F n−1 . Thus, the cardinality of D n is equal to the number of positive integers p between 1 and n for which the fraction p /n is irreducible. This fraction is irreducible if and only if p is coprime with n. Therefore, |D n | = ϕ(n), which completes the proof.
The following theorem gives a formula for computing the value of ϕ(p · k), where p is a prime number and k is any positive integer. This formula is used as an update rule by several algorithms described in this manuscript.
Theorem 3. Let p be a prime number and let k be a positive integer. Then, the value of Euler's totient ϕ(p · k) can be expressed as follows: Proof. Euler's product formula allows us to express the values of ϕ(k) and ϕ(p · k) as shown below: where the products run over the distinct prime numbers q that divide k and k · p, respectively. If p is a factor of k, then ϕ(p · k) can be expressed as follows: which proves the first branch in formula (22). If p is not a factor of k, then formulas (23) and (24) imply that ϕ(p · k) has the following value: which proves the second branch in formula (22).
One way to compute Euler's totients for all integers between 1 and N is to use formula (22) with the array L p that the linear sieve algorithm 18 generates. Algorithm S2 gives the pseudo-code for the linear sieve. Its main goal is to compute a list P that consists of all prime numbers in the interval [1, N ]. The algorithm solves this problem in O(N ) time and uses O(N ) memory. In addition to P , the linear sieve algorithm also fills an array L p that contains the smallest prime factor for each integer in the interval [1, N ] (the algorithm sets the value of L p [1] to 1 as a special case). The array L p is sufficient for computing the Euler's totients ϕ(1), ϕ(2), . . . , ϕ(N ) using a special case of formula (22) as described below.
Algorithm S3 fills the array ϕ such that its elements ϕ[1], ϕ[2], . . . , ϕ[N ] are equal to Euler's totients ϕ(1), ϕ(2), . . . , ϕ(N ). The algorithm starts by setting ϕ[1] to ϕ(1), which is equal to 1. Then, the algorithm iterates over m from 2 to N and sets ϕ[m] to Euler's totient ϕ(m). Each iteration of this for-loop uses the following version of formula (22): where p = L p [m] is the smallest prime number that divides m and k = m/p = m/p (we use the floor function to indicate that integer division should be used here). In other words, m = k · p. Because p is the smallest prime number that divides m, the prime factorization of m can be derived from the prime factorization of k by adding one more factor p to it. This also implies that the smallest prime factor of k must be greater than or equal to p, i.e., return ϕ; 14: end function Figure S2 visualizes a sample run of the linear sieve (i.e., Algorithm S2) with N = 5. It shows the state of the array of smallest prime factors L p and the list of prime numbers P for k ∈ {1, 2, 3, 4, 5}. The elements of L p and P that change after each iteration are colored in red. All elements of L p are initially set to null, which is indicated with ∅ in the figure. The list P starts empty.  Figure S2: Illustration of the linear sieve algorithm for N = 5.
The algorithm starts by setting L p [1] to 1. Then, it proceeds to the for-loop that starts at k = 2. Because L p [k] is null at the beginning of this iteration, the algorithm marks k = 2 as a prime number by setting L p [2] to 2 and appending 2 to the list P . Subsequently, the algorithm iterates over each p in P = (2), sets the variable q to the product p · k, and then sets L p [q] to p until q exceeds N . In our example, this leads to setting L p [4] to 2. Next, the algorithm proceeds to k = 3. Because L p [k] is null here, this value of k is also a prime number, which leads to setting L p [3] to 3 and appending 3 to P . No other elements of L p change because 2 · 3 > N = 5, which leads to breaking out of the inner for-loop on its first iteration. Subsequently, the algorithm continues to k = 4. Because L p [k] is not null, this value of k is not a prime number and the list P remains unchanged here. The inner for-loop also terminates early and does not modify any elements of L p . Finally, the algorithm advances to k = 5. The value of L p [k] is null here, which leads to setting L p [5] to 5 and appending 5 to the list P (i.e., 5 is a prime number). The inner for loop, again, terminates early and the function returns the final state of P and L p . Figure S3 shows an example for Algorithm S3 with N = 5. This algorithm computes Euler's totients ϕ(1), ϕ(2), ϕ(3), ϕ(4), and ϕ(5) using the array L p computed by the linear sieve algorithm in the previous example, i.e., L p = [1, 2, 3, 2, 5]. Algorithm S3 stores the computed totients in the array ϕ, which starts uninitialized. The uninitialized elements of ϕ are denoted with the symbol ∅ in the figure.  Figure S3: Illustration of the algorithm for computing Euler's totients with N = 5.
The algorithm starts by setting ϕ[1] to 1. Then, it continues to the for-loop, which starts at m = 2. It sets the value of the variable p to L p [m], which in this case is equal to 2. The variable k is set to m/p = 2/2 = 1. In this iteration, L p [k] = p and therefore the value of ϕ[m] is set to (p − 1) · ϕ[k] = (2 − 1) · ϕ[1] = 1. Next, the algorithm proceeds to m = 3. It sets the variable p to L p [m], which is equal to 3 and sets the variable k to m/p , which is equal to 1. The if-statement, again, sets the value of ϕ[m] to (p − 1) · ϕ[k] = (3 − 1) · 1 = 2. Subsequently, the algorithm proceeds to m = 4. In this iteration, p becomes 2 because L p [4] = 2. The variable k is also set to 2, because m/p = 4/2 = 2. In this case, L p [k] = p and therefore the first branch of the if-statement is executed, which sets ϕ[m] to p · ϕ[k] = 2 · ϕ[2] = 2. Finally, the algorithm continues to m = 5. Here, p becomes equal to 5 = L p [m] and k becomes equal to 1 = 5/5 . Because L p [k] = p, the value of ϕ[m] is set to (p − 1) · ϕ[k] = (5 − 1) · 1 = 4. After the last iteration, the function returns the final state of the array ϕ.

S4. LEMMAS USED IN THE PROOFS
The following lemma is used in the proof of Theorem 13.
Lemma 4. Let k and n be two positive integers and let S(k, n) be a set of all integers m such that n /m = k, i.e., S(k, n) = m ∈ N s.t. n /m = k . (28) Then, the set S(k, n) consists of all integers in the interval n k+1 , n k , i.e., S(k, n) = n k + 1 + 1, n k + 1 + 2, . . . , n k .
Moreover, if k 1 = k 2 , then the sets S(k 1 , n) and S(k 2 , n) are disjoint, i.e., Finally, the union of the sets S(k, n) over all k ∈ {1, 2, . . . , √ n } is equal to the set of integers that fall in the interval Proof. First, we will prove formula (29). Let m ∈ S(k, n). Then, n /m = k, i.e., k ≤ n m < k + 1.
Reciprocating the three terms leads to: Multiplying this inequality by n leads to the following interval for the value of m: Because m is an integer, the floor function can be applied to the lower and upper bounds for the value of m in the previous inequality, i.e., Therefore, m ∈ n k+1 , n k , which implies that S(k, n) ⊆ n k+1 , n k . Conversely, let m be an integer that lies in the interval n k+1 , n k . That is, Dividing by n and reciprocating leads to: n Because n /k ≤ n /k, the lower bound can be stated as: On the other hand, n k+1 + 1 > n k+1 . Therefore, the upper bound in (37) can be expressed as: Inequalities (38) and (39) imply that k ≤ n /m < k + 1 and, thus, n /m = k. Therefore, N ∩ n k+1 , n k ⊆ S(k, n), which proves formula (29).
Formula (30) is proven by contradiction. Suppose that there is a pair of integers k 1 and k 2 such that k 1 = k 2 and that the intersection of S(k 1 , n) and S(k 2 , n) is not empty. Then, there is an integer m that is an element of both sets, which implies that n /m = k 1 and n /m = k 2 . Thus, k 1 = k 2 , which is a contradiction that proves formula (30). Formula (31) is proven by showing that the union of the sets in the left-hand side forms a contiguous range of integers that matches the set in the right-hand side. In this case, larger values of k correspond to smaller elements of S(k, n). From formula (29) we can derive that the maximum element in the set S(k + 1, n) can be obtained by decrementing the minimum element in the set S(k, n). More formally, The minimum element in the set S( √ n , n) is equal to n √ n +1 + 1 = u(n) + 1. The maximum element in the set S(1, n) is equal to n. Thus, all integers between u(n) + 1 and n are elements of the union.
Many formulas in this manuscript use the function u(n) = n/( √ n + 1) to simplify the notation. Depending on the value of the integer n, this function evaluates to either √ n − 1 or √ n . For example, if n = 4, then it evaluates to 1 = √ 4 − 1. For n = 6, however, it evaluates to 2 = √ 6 . The next lemma formally proves this property.
Lemma 5. For any positive integer n, the value of the function u(n) Proof. For the fraction in the definition of u(n), the value of the denominator lies between √ n and √ n + 1, i.e., Reciprocating all three values changes the direction of the inequality and leads to: Multiplying by n leads to the following bounds for the fraction in u(n) that is used as the argument of the floor function: The value of the term n/( √ n + 1) in the left-hand side can be expressed as follows: Therefore, Combining inequality (43) with inequality (45) leads to the following formula: Applying the floor function to the previous inequality makes the middle term equal to u(n). It also changes the inequality such that it is no longer strict when √ n is an integer. This completes the proof, i.e., where √ n − 1 = √ n − 1 follows from the definition of the floor function.
The following two lemmas are used in the proof of Theorem 12. They express the value of n+1 k in two ways depending on whether k divides n + 1 or not.
Lemma 6. Let n be a positive integer and let k be an integer between 1 and n. If k does not divide n + 1, then the value of n+1 k is equal to n k . That is, Proof. The properties of integer division imply that there is exactly one pair of integers x and y such that the value of n + 1 can be expressed in the following form: where 0 ≤ y < k. In this case, x = n+1 k is the quotient and y = (n+1) mod k is the remainder. The lemma states that k does not divide n + 1, which implies that y cannot be zero, i.e., 1 ≤ y < k.
Subtracting 1 from both sides of equation (49) leads to the following expression for the value of n: The value of the term y − 1 is an integer that lies in the interval [0, k − 1). Once again, the uniqueness of integer division implies that y − 1 = n mod k and that x = n k , which completes the proof.
Lemma 7. Let n be a positive integer. Also, let k be an integer between 1 and n. If k divides n + 1, then the value of n+1 k is a whole number equal to n k + 1, i.e., n + 1 Proof. Similarly to the previous proof, the properties of integer division imply that: where x = n+1 k and y = (n+1) mod k = 0. In this case, y = 0 because n + 1 is evenly divisible by k. Once again, we can subtract 1 from both sides of (52), which leads to the following formula for the value of n: Furthermore, we can express x k as (x − 1)k + k, which leads to Here the uniqueness of integer division implies that n k = x − 1.
Recalling that x = n+1 k completes the proof.
The next lemma proves an interesting property of the floor function, which is used in the correctness proof for Algorithm C.
Lemma 8. Let k, j, and n be three positive integers. Then, Proof. Similarly to the previous two lemmas, we start by using integer division to express the value of n/k as follows: where x = n/k j and y = n/k mod j. The integer y lies in the closed interval [0, j − 1]. The value of n/k is equal to the sum of its integer and fractional parts, i.e., where f = n/k − n/k . Because n and k are integers, the value of f is an element of the closed interval 0, k−1 k , i.e., it never reaches 1.
Combining (56) and (57) allows us to derive the following formula for the value of n/k j : This formula, in turn, implies that the following upper and lower bounds hold for n/k j : In other words, x ≤ n/k j < x + 1. Applying the floor function to all three terms leads to x ≤ n/k j < x + 1, since x is an integer. Because the term in the middle must be an integer and because the second inequality remains strict, it follows that n/k j must be equal to x = n/k j . This completes the proof.
Lemma 9. Let k, j, and n be three positive integers such that n ≥ kj. Then, the lengths of the following two Farey sequences are equal: Proof. The proof follows from Lemma 8, which proved that the values of the two subindices in (60) are equal. Also, when n ≥ kj, the value of n/(kj) is an integer that is greater than or equal to 1, i.e., it identifies an order of a Farey sequence.
The next two lemmas are used to prove the computational complexities of FAREYLENGTHL and FAREYLENGTHC (i.e., Algorithm S10 and Algorithm 5), respectively. Lemma 10. Let N be a positive integer. Then, Proof. Let f (x) = 1 x . This function is monotonically decreasing in the interval (0, ∞). Therefore, for each k ∈ N = {1, 2, . . . }. This implies that we can find lower and upper bounds for the sum in formula (61) as follows: The left-most integral is equal to: The right-most term in this inequality evaluates to: Therefore, the bounds for the sum are: which completes the proof.
Lemma 11. Let N be a positive integer. Then, Proof. Let f (x) = 1 √ x . This function is monotonically decreasing in the interval (0, ∞). Thus, for each k ∈ N = {1, 2, . . . } the following inequality holds: From this inequality we can derive lower and upper bounds for the value of the sum 1/ The left-most integral evaluates to: Similarly, the right-most integral has the following value: Combining (69) with (70) and (71) leads to: from which (67) follows.
S5. PROOF OF FORMULA (2) This section proves formula (2), which expresses the length of the Farey sequence F n recursively in terms of the lengths of the Farey sequences of lower orders. This formula is well-known, but its proof is hard to find. Our proof uses only basic algebra and mathematical induction. The proof also uses Euler's totient function ϕ(n), which is defined in Section S3.
Theorem 12. The length of the Farey sequence of order n can be computed recursively as follows: where x denotes the largest integer that does not exceed x.
Proof. The proof is by mathematical induction. The base case of the induction is formed by the length of F n when n = 1, i.e., To prove the inductive step, we will introduce the term T n that is defined as follows: Our goal is to prove that We have already established that this is true for n = 1. Assuming that this formula holds for some n ≥ 1, our next goal is to prove that it also holds for T n+1 . In other words, we would like to show that the value of T n+1 can be expressed as follows: Using Theorem 1 and Definition 1 (see Section S3), we can express the term T n+1 as shown below: Lemmas 6 and 7 (see Section S4) imply that in (78) the value of each term in the sum from 2 to n can be expressed as follows: Moreover, for the second case in (79), Theorem 1 and Lemma 7 imply that: Let D be the set of all positive integer divisors of n+1 that lie between 2 and n. That is, Then, the sum in equation (78) can be expressed as follows: If an integer k is an element of the set D, then the value of (n+1) /k is also an element of D and vice versa. Therefore, Thus, equation (82) can be restated as: Plugging this result into (78) leads to the following formula for the value of T n+1 : Because ϕ(1) = 1, it follows that ϕ(1) + 1 = 2. This allows us to rearrange the terms in the last equation as follows: In this formula, k|n+1 ϕ(k) denotes the sum of the values of Euler's totient function ϕ(k) for all k that are positive integer divisors of n + 1. In 1798, Gauss proved 14 that this sum is equal to n+1. Thus, T n+1 = T n + (n + 1) + 1 = T n + n + 2, which proves equation (77) as required.
S6. PROOF OF FORMULA (5) This section proves another recursive formula for the value of |F n |. This formula is derived from formula (2) by splitting the sum into two segments and grouping the repeated terms in the second segment. The resulting formula is longer, but it needs to add fewer terms. It leads to Algorithm 5 (i.e., FAREYLENGTHC), which is described in the main paper.
Theorem 13. For each n > 1, the length of the Farey sequence F n can be expressed using the following recursive formula: where u(n) = n √ n +1 . Proof. This formula is derived from formula (2), which is replicated below: The first step is to split the summation over k at u(n) into two separate sums: And then change the index variable from k to m in the second sum: The next step is to express the last sum in formula (91) so that it has only √ n terms. Lemma 4 proves that the set of the indices m used in that sum is equal to the following union: where each set S(k, n) consists of all positive integers m for which n /m = k, i.e., S(k, n) = n k + 1 + 1, n k + 1 + 2, . . . , n k .
Lemma 4 also proves that all sets S(k, n) in this union are disjoint. Each of these sets has n k − n k+1 elements. Therefore, for each n > 1 the last sum in (91) can be expressed as follows: Substituting this expression into (91) and that result into (89) proves the theorem.
Formula (88) cannot express F 1 because the last sum leads to an infinite recursion when n = 1. Also, in this special case the sum in formula (89) is degenerate. In contrast to the cases when n ≥ 2, the sum cannot be split in a way that leads to formula (90) when n = 1. Therefore, any algorithm that implements formula (88) should handle n = 1 as a special case. For example, F 1 can be hard-coded to be equal to 2.

√
n AND 3 √ n 2 EXACTLY Algorithms S4 and S5 give the pseudo-code for two helper functions that are used by several of the algorithms. The first function computes √ n exactly using Newton's method and integer arithmetic. 28 The second function computes the value of 3 √ n exactly using a similar approach. All algorithms, however, call ICBRT with an argument that is equal to n 2 instead of n. Thus, it is used to compute x ← n; 3: y ← n /2 ; 4: while y < x do 5: x ← y; x ← n; 3: y ← n /3 ; 4: while y < x do 5: x ← y; This appendix describes an implementation for a lookup table that satisfies these time and space complexity constraints. The implementation is specific to our algorithms, i.e., it works for the special case when all keys (or indices) that are used are in the set {1, 2, . . . , d} ∪ n n/d , . . . , n 2 , n 1 . Alternatively, it is possible to use a generic hash table as a lookup table (e.g., the hash table implemented in the standard dictionary class in the Python language). In that case the computational complexity of setting or getting entries is in O(1) on average. Even though this is a minor technical point, this appendix shows that the desired performance can be achieved in all cases, i.e., not just on average.
Algorithm S6 gives the pseudo-code for initializing the lookup table. The function has two parameters: n and d. The value of n determines the size of the problem, i.e., the order of the Farey sequence F n for which we want to compute the length. The parameter d controls the split between the two sets of keys that the table supports. The table is organized as two arrays. The first array, X, stores the entries for the keys in the set {1, 2, . . . , d}. The second array, Y , stores the entries for the keys in the set n n/d , . . . , n 2 , n 1 . In most cases, d = √ n + 1, except for Algorithm S13, where it is set to 3 √ n 2 + 1.
Algorithm S7 gives the pseudo-code for getting the value of an element from the lookup table, given a key k. This function is called when the square bracket notation is used with F in our algorithms, i.e., F [k] translates to GETITEM(F, k); Algorithm S8 gives the pseudo-code for setting an element of the lookup table. Once again, it is assumed that an assignment F [k] ← x in our algorithms translates to a call to SETITEM(F, k, x).
Finally, Algorithm S9 gives the pseudo-code for a function that checks if the value for the key k is set in the lookup table F . This check is used in some of our algorithms for diagnostic purposes. That is, it is assumed that a check k in F translates to a function call CONTAINS(F, k).  (2) and (5) from the main paper, which are reproduced below: Algorithm S10 gives the pseudo-code for the first algorithm, which is given a suffix letter L. The helper function in this case implements formula (95). The main algorithm, however, structures the calls to this function in the same way as Algorithm 5, which is based on formula (96). Theorem 14 (see below) proves that the computational complexity of this algorithm is in O(n log n). The space complexity is in O( √ n ). This is easy to show because FAREYLENGTHL is structured similarly to FAREYLENGTHC, except for using a different helper function (see also the proofs in Section S10).
Algorithm S10. Compute the length of the Farey sequence F n . Runs in O(n log n) time and uses O( √ n ) memory.

22: end function
Algorithm S11 is a recursive version of Algorithm S10. This algorithm also runs in O(n log n) time and uses O( √ n ) memory. This can be established by unpacking the recursion, which leads to the iterative version presented above. This is the shortest algorithm described in this paper. Unfortunately, it is also one of the slowest. For small values of n, however, it could be very useful. For example, it can be used to verify the output of the faster algorithms.
Algorithm S11. Recursive algorithm for computing the length of F n . Runs in O(n log n) time and uses O( √ n ) memory. The main algorithm calls the helper function from within two for-loops, but they are not nested. Therefore, the computational complexity of the algorithm can be expressed as follows: The first sum in the right-hand side of this formula is in O(n). That is, The second sum in the right-hand side of (97) is in O(n log n). This can be shown as follows: In this formula, the left-most inequality follows from Lemma 5, which implies that u(n) ≤ √ n . The derivation can be continued using Lemma 10, which implies that: Therefore, Combining formulas (98) and (101) with formula (97) proves that Algorithm S10 runs in O(n log n) time.

S10. COMPUTATIONAL COMPLEXITY ANALYSIS OF FAREYLENGTHC
This section analyzes the computational complexity of FAREYLENGTHC (i.e., Algorithm 5) and proves that it runs in O(n 3/4 ) time and uses O( √ n ) memory. It also proves that the elements of the lookup table are computed in such a way that they are always available when the helper function needs them.
The following theorem gives the time complexity of the helper function UPDATELOOKUPTABLE (i.e., Algorithm 6).
Proof. This algorithm implements a helper function that uses formula (5) to compute the value of F m . That is, This function has two for-loops, but they are not nested. These loops implement the first and the second summation in the formula, respectively. One call to this function computes exactly one new entry in the lookup table F . The main algorithm (i.e., Algorithm 5) ensures that all calls to this function are performed in the correct order such that the values from F that it uses are already set because they were computed earlier (see the proofs on the next page). The computational complexity of each iteration in either of the two loops is in O(1). Thus, the computational complexity of the helper function can be determined as follows: This result follows from Lemma 5, which proves that u(m) ≤ √ m . This, in turn, implies that u(m) − 1 + √ m < 2 √ m . Therefore, this function runs in O( √ m) time, where m is its argument, i.e., not the n for which the main algorithm is called.
The next theorem establishes the time and space complexity of FAREYLENGTHC (i.e., Algorithm 5). Proof. The algorithm calls the helper function from inside two for-loops, which are not nested. The first loop iterates over m from 2 to √ n . The second loop iterates over j from u(n) down to 1. Thus, the overall computational complexity can be expressed as follows: The first sum in the right-hand side of (104) can be bounded from above as follows: This result implies that The second sum in the right-hand side of (104) can also be bounded from above: which follows from Lemmas 5 and 11. To summarize, we showed that both terms in the right-hand side of (104) are in O(n 3/4 ). Therefore, the computational complexity of FAREYLENGTHC is also in O(n 3/4 ). The algorithm fills √ n entries of the lookup table F in the first for-loop and u(n) entries in the second for-loop. Because u(n) ≤ √ n , the space complexity of the algorithm is in O( √ n ).
Algorithm 5 uses the helper function to compute the value of |F m | for each integer m in two lists. The first list has √ n − 1 elements and consists of the numbers 2, 3, . . . , √ n . The second list has u(n) elements and consists of the numbers n u(n) , n u(n)−1 , ..., n 1 . The algorithm processes the first list sequentially starting from 2. After that, it processes the second list sequentially starting from n/u(n) . In total, the algorithm processes √ n − 1 + u(n) ≤ 2 √ n values of m. Let Q(m) be the set of keys in the lookup table F that are used by the first for-loop of the helper function (i.e., Algorithm 6). For a given value of m this set is equal to: Let K(m) be the set of keys in the lookup table F that are used by the second for-loop of the helper function. That is, The following theorem helps us to prove that all entries in the lookup table F that the helper function uses to compute the value of F [m] are already set when the function is called in the first for-loop of Algorithm 5.
Theorem 17. Let m be an integer that is greater than 1. Then, both Q(m) and K(m) are subsets of the set of all integers in the interval [1, m − 1]. More formally, Let Y (j) be the set of keys in the lookup table F that are already computed by Algorithm 5 before it calls the helper function in the second for-loop with the parameter m equal to n/j . More formally, , n u(n) − 1 , . . . , n j + 1 .
Then, the following theorem helps us to prove that the second for-loop in Algorithm 5 also uses only those entries of F that have already been set.
This implies that k ∈ Y (j), which shows that Q( n/j ) ⊆ Y (j) and completes the proof of the theorem.
Corollary 20. For each iteration of the second for-loop in Algorithm 5, all entries in the lookup table F that are necessary for computing the value of F n/j are already available when the helper function is called to compute it.
Proof. This corollary follows from Theorem 19. After the end of the first for-loop of Algorithm 5 the values of |F m | for m ∈ {1, 2, . . . , √ n } are already computed. The second for-loop iterates over values of j from the list (u(n), u(n) − 1, . . . , 1) in decreasing order starting from u(n). Thus, the entries for keys from the set Y (j) are already set when the control flow enters the helper function with m = n/j . The function uses keys from the sets Q( n/j ) and K( n/j ) to compute F n/j . Because both Q( n/j ) and K( n/j ) are subsets of Y (j), all required entries are available and the helper function can compute F n/j successfully.

S11. RECURSIVE VERSION OF FAREYLENGTHC
Algorithm S12 gives the pseudo-code for a recursive version of Algorithm 5 (i.e., FAREYLENGTHC). This version is a direct implementation of formula (5), which is replicated below: return F [n]; 20: end function S12. COMPUTATIONAL COMPLEXITY ANALYSIS OF FAREYLENGTHD FAREYLENGTHD can be viewed as a mixture between FAREYLENGTHA and FAREYLENGTHC. It picks an optimal point to divide the computation between these two algorithms as described below.
In order to compute |F n |, FAREYLENGTHC (i.e., Algorithm 5) computes |F m | for values of m that can be split into two sets M 1 and M 2 that correspond to the two for-loops in the algorithm. That is, the first for-loop computes the length of F m for m ∈ M 1 and the second for-loop computes the length of F m for m ∈ M 2 . The set M 1 includes all integers between 1 and √ n , i.e., The set M 2 is defined by the following formula: Because u(n) = n √ n +1 , it follows that n u(n) ≥ √ n + 1. Thus, each element of M 2 falls between √ n + 1 and n. This set is sparse because it has u(n) elements, but the total number of integers between √ n + 1 and n is equal to n − √ n . In other words, u(n) ≤ √ n (see Lemma 5), but n − √ n > √ n for each n > 2. n . Nevertheless, it is possible to change the split between the sets M 1 and M 2 so that the sieve-based approach processes more values of m, including some values that are greater than √ n . Because that approach computes the values of |F m | for each integer m in the designated interval, it would still compute |F m | for each m ∈ M 1 . In other words, the modified algorithm can run faster despite computing more values of |F m | than required for computing |F n | using the recursive formula.
Let M * 1 (x) be a function that maps a real number x to the following set: Also, letû(x) be the following function:û Furthermore, let M * 2 (x) be another function that maps a real number to a set, which is defined as follows: is the following function: The value of g(x) can be bounded from above as follows: Lemma 11 implies that the sum in the right-hand side of this inequality is in O( n/x ). In other words, The time complexity of computing |F m | for each m ∈ M * 1 (x) ∪ M * 2 (x) by an algorithm that uses the sieving approach to process the set M * 1 (x) and the recursive approach to process M * 2 (x) is in O(f (x)), where f (x) is the following function: What is the optimal value of x for which f (x) is minimized? Differentiating f (x) leads to the following formula for its first derivative: The equation f (x) = 0 has only one root at x = ( n /2) 2 /3 . Multiplying this solution by a constant factor that is independent of n does not affect the time complexity class of the overall computation because any constant factor is subsumed by the big-O notation. For simplicity, the optimal split x * can be set to n 2/3 . FAREYLENGTHD is the result of the modifications to FAREYLENGTHC described above. It replaces the first for-loop with FAREYLENGTHA in order to process the set M * 1 ( n 2/3 ). It also modifies the second for-loop to start from v = n / 3 √ n 2 +1 in order to compute |F m | for each m in the set M * 2 ( n 2/3 ). The resulting algorithm runs in O(n 2/3 ) time because each of these two sub-parts requires O(n 2/3 ) time. The space complexity of this algorithm is determined by the sieve-based approach that processes the set M * 1 ( n 2/3 ). Because this step requires holding arrays of length n 2/3 in memory, FAREYLENGTHD uses O(n 2/3 ) memory.

S13. ALTERNATIVE VERSION OF FAREYLENGTHD
Algorithm S13 gives the pseudo-code for an alternative version of Algorithm 7 (i.e., FAREYLENGTHD). This version has the same run-time complexity and the same memory complexity as the one described in the main paper. It is shorter than the other algorithm, but updates more entries of the lookup table than necessary to compute |F n |. It also obscures the link to Algorithm 8 (i.e., FAREYLENGTHE), which improves the space complexity. Additional details are provided in the main paper. ϕ ← COMPUTETOTIENTS(c, L p ); The following theorem proves the time and space complexity of FAREYLENGTHE (i.e., Algorithm 8).
Proof. The time complexity of Algorithm 8 is in the same complexity class as Algorithm 7, i.e., O(n 2/3 ). More specifically, FAREYLENGTHE starts by processing the integers in the interval [1, √ n ] using the linear sieve and Algorithm S3. Next, the algorithm enumerates all √ n -smooth numbers in the interval [ √ n + 1, n 2/3 ] using Algorithm 3 and adds their totients to the corresponding elements of the array B. This step requires O(n 2/3 ) time and uses O( √ n ) memory. Next, the algorithm enumerates all numbers that are not √ n -smooth in the interval [ √ n + 1, n 2/3 ] using Algorithm 4 and also adds their totients to the corresponding elements of B. This step requires O(n 2/3 ) time and O( √ n ) memory because it uses the sieve of Atkin 19 to find the prime numbers between √ n + 1 and n 2/3 . Finally, the algorithm computes the value of |F m | for each m in the set The following three theorems prove that the computation in Algorithm 8 is correct. They show that the set of integers in the interval [α, β] that the algorithm enumerates is the correct set, where α = √ n + 1 and β = n/(v(n) + 1) . The theorems also show that the formula for mapping k to the index i of the array element B[i] is also correct for each integer k in [α, β].
The first theorem shows that the set S(u(n), n), which consists of all integers j such that n j = u(n), includes the value of √ n + 1. This implies that α = √ n + 1 is a suitable split point for switching between the linear sieve approach and enumeration of smooth and non-smooth integers. For this value of α the algorithm computes the value of B[0] correctly because the split point is an element of S(u(n), n) so that all integers in S(u(n), n) are accounted for either by the computation of |F 1 |, |F 2 |, . . . , |F √ n | using the linear sieving approach or by enumerating smooth and non-smooth numbers. Theorem 22. Let S(k, n) be a set of all positive integers m such that the value of n m is equal to k. That is, Then, √ n + 1 ∈ S(u(n), n). Moreover, the set S(u(n), n) is a contiguous range of integers such that min S(u(n), n) = n u(n)+1 + 1 ≤ √ n + 1 ≤ max S(u(n), n) = n u(n) .
Proof. Formula (129) follows from Lemma 4 after plugging u(n) as the value of k in formula (29). Formula (130) follows from (129) and the definitions for the set S(u(n), n) and the function u(n). That is, √ n+1 ∈ S(u(n), n) because S(u(n), n) is a contiguous range of integers. Thus, √ n + 1 lies between its minimum and maximum elements.
The second theorem shows that the endpoint β = n v(n)+1 that Algorithm 8 uses as the upper limit for enumerating smooth and non-smooth numbers is the correct endpoint. That is, it shows that the enumeration process covers all elements of the sets S(u(n) − 1, n), S(u(n) − 2, n), . . . , S(v(n) + 1, n), which is required for computing the elements of the array B correctly.
Theorem 23. The union of the sets S(u(n) − 1, n), S(u(n) − 2, n), . . . , S(v(n) + 1, n) is equal to the contiguous range of integers between n u(n) + 1 and Proof. This theorem follows from equation (29) in the statement of Lemma 4, which implies that the set S(u(n) − 1, n) is a contiguous range of integers that starts from n u(n) + 1, the set S(u(n) − 2, n) is a contiguous range of integers that follows S(u(n) − 1, n) without any gaps in-between, and so forth until S(v(n) + 1, n). The lemma also shows that the maximum element of S(v(n) + 1, n) is equal to n v(n)+1 . The last theorem proves that the formula i = u(n) − n m that Algorithm 8 uses to compute the value of the index for the array B in the visitor function correctly identifies the zero-based indices of these elements. This theorem implies that i = 0 corresponds to S(u(n), n), i = 1 corresponds to S(u(n) − 1, n), etc., until i = w(n), which maps to S(v(n) + 1, n), where w(n) = u(n) − v(n) − 1.
Theorem 24. Let m and n be two positive integers and let m be an element of the set S(k, n). Then, k = u(n) − i, where i = u(n) − n m . Proof. The proof follows from the definition of the set S(k, n) in Lemma 4. That is, m ∈ S(k, n) implies that n m = k. Therefore, i = u(n) − n m = u(n) − k. Thus, k = u(n) − i.

S15. RESULTS FOR THE NUMBER OF CPU INSTRUCTIONS
In addition to the run time and memory usage results described in the main paper, we also measured the number of CPU instructions that the code executed (see Methods). Figure S4a visualizes these additional results for the C implementations of algorithms C, D, and E. Figure S4b shows the results obtained with the Python versions of these three algorithms. In both plots the value of n was varied from 10 0 to 10 14 in increments of 0.5 on the decimal logarithm scale. Figure S4c lists the slopes and intercepts for six lines that were fitted to the six curves in (a) and (b) using least squares. The fitting procedure used the region between n = 10 10 and n = 10 14 .
These results show that the number of CPU instructions agrees with the theoretical time complexities. Each slope is within 0.01 of the corresponding power of n in the time complexity class, i.e., 3 /4 for algorithm C and 2 /3 for algorithms D and E. More specifically, the slopes for the two programming language implementations of algorithm C are equal to 0.7494 and 0.7520, respectively. For algorithm D they are equal to 0.6672 and 0.6754. For algorithm E they are equal to 0.6655 and 0.6744.
This precise agreement with the theory suggests that the run time results reported in Figure 8 may be slightly higher than the theoretical predictions due to the practical aspects of running the code on modern processors. One of them could be increasingly inaccurate branch prediction as n increases. Another could be a greater number of cache misses for larger n. Because the time spent waiting for the cache to be updated does not affect the CPU instruction counter for the current process, this metric agrees with the theory very well.