Algorithms for Complementary Sequences
Full text
#A95 INTEGERS 25 (2025) ALGORITHMS FOR COMPLEMENTARY SEQUENCES Chai Wah Wu Mathematics of Computation, IBM Research, T. J. Watson Research Center, Yorktown Heights, New York [email protected] Received: 9/9/24, Accepted: 9/30/25, Published: 11/5/25 Abstract Finding the n-th positive square number is easy, as it is simply n2. But how do we find the complementary sequence, i.e., the n-th positive non-square number? For this case there is an explicit formula. However, for general constraints on numbers, a formula is harder to find. In this paper, we study how to compute the n-th integer that does (or does not) satisfy a certain condition. In particular, we consider it as a fixed point problem, relate it to the iterative method of Lambek and Moser, study a bisection approach to this problem, and provide novel formulas for various complementary sequences including the non-k-gonal numbers, non-kgonal-pyramidal numbers, non-k-simplex numbers, non-sum-of-k-th-powers, non-kth-powers, and non-Jacobsthal numbers. 1. Introduction For a positive integer n∈N+, the n-th positive square number is simply n2. Can we also easily find the complementary sequence? In other words, what is the n-th positive non-square number? It is quite remarkable that there exists an explicit formula for the n-th positive non-square number: n+⌊1 2+√n⌋=n+⌊pn+⌊√n⌋⌋ [9,13,15,16]. This can also be computed as n+⌊√n⌋+ 1 if n > ⌊√n⌋(⌊√n⌋+ 1) and as n+⌊√n⌋otherwise. These formulas are well-suited for implementation in a computer algorithm since many computer languages and number theory software packages include functions to compute ⌊√n⌋. For instance, the isqrt function in Python, Julia, and Maple all perform this calculation. These formulas have been extended to higher powers as well. In particular, the n-th non-k-th-power number is given by n+⌊k pn+⌊k √n⌋⌋ [13,17,22].1 DOI: 10.5281/zenodo.17535229 1The computation of ⌊k √n⌋for arbitrary integers kand n≥0 is readily available in symbolic computer algebra systems and software for number theory. For instance, ⌊k √n⌋can be computed with the integer_nthroot function in the sympy Python module which in turn uses the mpz_root function in the multiple precision library gmp. Although these computer operations assume that
INTEGERS: 25 (2025) 2 For Pa logical statement on the natural numbers, let us define fP(n) as the n-th positive natural number msuch that P(m) is true. For the case where P(m) denotes the logical statement “mis square”, fP(n) is easily determined, since the list of integers msuch that P(m) is true are easily enumerated. As noted in the example above, for this particular P, the function f¬P(n) can also be computed by an explicit formula. However, in general, the simplicity of fPdoes not imply the simplicity of f¬P. Furthermore, for more general statements P, the formula for fP(n) or f¬P(n) may not be readily available. Even if such explicit formulas are available, some of them require the use of floating point arithmetic and it can be difficult to use computationally to find fP(n) or f¬P(n), especially for large n. See for example the formulas for the n-th non-Fibonacci number in [4,5] which require logϕat high precision for large n. While there have been many studies of explicit formulas for such complementary sequences [4,5,9,13,15–17,22], there have not been much study in computer algorithms to calculate such sequences. The purpose of this paper is to discuss algorithms to compute fP(n) or f¬P(n). As a consequence, we show that the n-th non-k-gonal number is given by n+$r2n−2+⌊k+1 4⌋ k−2+1 2% and that the n-th non-second-hexagonal number is n+pn 2−1. 2. Finding fP(n) as the Solution to a Fixed Point Problem For an integer a, define the counting function CP(a) = |{b∈N: (1 ≤b≤a)∧P(b)}| as the number of positive integers less than or equal to asuch that P(a) is true. It is clear that CP(a) is increasing, 0 ≤CP(a)≤a, and 0 ≤C¬P(a) = a−CP(a)≤a. Furthermore, fP(n) is the smallest integer msuch that CP(m) = n. Also note that fPis strictly increasing and fP(n)≥n. Define gn(x) = n+C¬P(x) = n+x−CP(x). A fixed point xof gnsatisfies x=n+x−CP(x), i.e., CP(x) = n. Thus, the smallest fixed point of gnis equal to fP(n). Furthermore, a fixed point of gnthat is in the range of fPis equal to fP(n). In particular, if gnhas a unique fixed point, then it must necessarily be equal to fP(n). Finding a fixed point of f(x) is equivalent to finding a root of f(x)−x. Equivalently, we could define ˜gn(x) = n−CP(x) and find the roots of ˜gn. However, in the sequel we will consider the fixed point formulation as gnis defined with C¬P and has a more natural interpretation for complementary sequences. Furthermore, we show below that the function iteration method solving this fixed point problem is equivalent to the well-known Lambek-Moser method for defining complementary sequences. The function gn, viewed as a function on the real numbers, is a piecewise-linear nis an integer, they can be used to compute ⌊k √n⌋for all real n≥0 since ⌊k √n⌋=⌊k p⌊n⌋⌋ for n≥0 (see [6, Equation 3.9]).
INTEGERS: 25 (2025) 3 function. For each value of n, the function gncoincides with the identity function on the segment spanned by {m∈N:CP(m) = n}. It is clear that gn(m)> m if m < fP(n) and gn(m)≤mfor m≥fP(n). As an example, we show in Figure 1 the function gnfor the case where P(m) denotes the statement “mis prime” and nis equal to 4. Notice that the minimal fixed point is at x= 7 which corresponds to fP(4), i.e., the fourth prime number. Next, let us consider methods to find the smallest fixed point of such an increasing piecewise-linear function gnon the integers. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 x 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 g 4( x ) Figure 1: gn(x) when n= 4 and P(m) denote the logical statement “mis prime”. The minimal fixed point is at x= 7 which corresponds to fP(4), the fourth prime number. 2.1. Function Iteration Method The function iteration method to find a fixed point of a function is a classical method that dates back to at least Heron’s method for finding an approximation to the square root [7] and is used in general root finding algorithms. We first pick an initial condition less than or equal to fP(n). Since n≤fP(n), we can start
INTEGERS: 25 (2025) 4 with the initial condition x=nand apply the iteration x→gn(x) repeatedly until convergence (Algorithm 1). This is for example implemented in the FixedPoint function in Mathematica. Note that since gn(x)≥xinitially, at each step of the algorithm the value of xincreases, until it reaches a point where gnintersects with the identity function, which is the smallest fixed point, i.e., fP(n). For the initial condition x=n, it is easy to see that this method is equivalent to the Lambek-Moser method and [13] showed that it indeed converges to the smallest fixed point. Algorithm 1 Function iteration method on gn(x) to compute fP(n). Require: gn(x)▷computes the minimal fixed point of gn. m←n while gn(m)=mdo m←gn(m) end while return m While the Lambek-Moser method assumes the initial condition x=n, depending on Pwe may choose a more suitable initial condition. For instance, if P(m) denotes the statement “mis the product of kdistinct primes”, then the initial condition can be chosen as max(n, pk#) since fP(n)≥pk# where pk# is the k-th primorial. The number of steps needed for convergence is less than fP(n)−nand thus this algorithm is efficient when fP(n)−nis small with respect to n, i.e., when the numbers satisfying Pare dense. In particular, in [13] it is shown that if the difference function of f¬Phas at least a linear growth rate (implying that f¬Pgrows at least quadratically) then 2 steps suffice. More precisely, the following result is shown. Theorem 1. If f¬P(m+ 1) −f¬P(m)≥mfor all m, then fP(n) = gn(gn(n)) = n+C¬P(n+C¬P(n)). Sequences satisfying these conditions include non-k-th-powers or the non-powers of k. Thus, in these cases the computation of fP(n) requires at most 2 evaluations of the counting function C¬P. We next show that 1 evaluation suffices. Theorem 2. If f¬P(m+ 1) −f¬P(m)≥mfor all m, then fP(n) = (n+C¬P(n)+1 if n+C¬P(n)≥f¬P(C¬P(n) + 1) n+C¬P(n)otherwise. Proof. First note that f¬P(C¬P(n) + 1) > n. By hypothesis, f¬P(C¬P(n) + 2) ≥C¬P(n) + 1 + f¬P(C¬P(n) + 1) > C¬P(n) + n=gn(n).
INTEGERS: 25 (2025) 5 This means that C¬P(n)+2> C¬P(C¬P(n) + n) = gn(gn(n)) −n, i.e., gn(gn(n)) < n +C¬P(n)+2. Since C¬Pis increasing, gn(gn(n)) ≥n+C¬P(n). Finally, it is easy to see that the threshold where gn(gn(n)) changes from n+C¬P(n) to n+C¬P(n) + 1 is precisely given by n+C¬P(n)≥f¬P(C¬P(n) + 1). This can be more compactly expressed using the Iverson bracket [10,12], which we denote using JK: fP(n) = n+C¬P(n) + Jn+C¬P(n)≥f¬P(C¬P(n) + 1)K. Next we show several applications of Theorem 2. 2.1.1. Non-k-th-powers As an example of applying Theorem 2, we give the following formula for the nth non-k-th power for k > 1, which simplifies the formula given in Section 1by requiring only one evaluation of the integer k-th root function ⌊k √n⌋: a(n) = (n+⌊k √n⌋+ 1 if n+⌊k √n⌋ ≥ (⌊k √n⌋+ 1)k n+⌊k √n⌋otherwise. 2.1.2. Non-Mersenne Numbers Similarly, for the non-Mersenne numbers (i.e., numbers not of the form 2p−1 for p prime), we have for n > 1 the formula a(n) = n+s+ 1 if n+1+s≥2ps+1 n+sotherwise, where s=⌊log2π(n)⌋,pkis the k-th prime, and π(n) is the prime counting function that returns the number of prime numbers less than or equal to n. 2.1.3. Non-Fermat Numbers The Fermat numbers are defined as 22n−1+ 1 for n≥1, i.e., 3,5,17,257,65537, . . .. The first two non-Fermat numbers are 1 and 2 and the n-th non-Fermat numbers for n > 2 are a(n) = (n+⌊log2(log2(n−1))⌋+ 2 if n+⌊log2(log2(n−1))⌋ ≥ 42⌊log2(log2(n−1))⌋ n+⌊log2(log2(n−1))⌋+ 1 otherwise.
INTEGERS: 25 (2025) 6 2.1.4. Non-powers of k Similarly, for the non-powers of kwe have the formula a(n) = (n+1+⌊logkn⌋+ 1 if n+1+⌊logkn⌋ ≥ k⌊logkn⌋+1 n+1+⌊logkn⌋otherwise. 2.1.5. Non-Jacobsthal Numbers The Jacobsthal numbers 0,1,1,3,5,11,21, . . . are defined recursively as the sequence a(n) = nfor n≤1 and a(n) = a(n−1)+2a(n−2) otherwise. It can also be defined as the nearest integer to 2n 3. This sequence is sequence A001045 of the on-line encyclopedia of integer sequences (OEIS) [19]. Looking only at the positive integers and ignoring the duplicate 1, we obtain the sequence 1,3,5,11,21, . . . whose definition can be written as f¬P(n) = j2n+1+1 3k. The number of such numbers less than or equal to nis 1,1,2,2,3, ... which can be written as C¬P(n) = ⌊log2(3n+ 1)⌋ − 1. Applying this to Theorem 2results in the following formula for the non-Jacobsthal numbers (OEIS A147613): a(n) = (n+⌊log2(3n+ 1)⌋if n+⌊log2(3n+ 1)⌋>j2⌊log2(3n+1)⌋+1+1 3k n+⌊log2(3n+ 1)⌋−1 otherwise. (1) Note that ⌊log2(x)⌋+1 is equal to the number of bits in the binary representation of the integer xand this quantity is easily available in many programming languages. Contrast Equation (1) with the first formula listed in the entry of sequence A147613 in OEIS which requires floating point arithmetic due to the computation of loge(2) and the Lambert W function. 2.2. Interleaving Functions Typically, the computation of C¬Prequires the inversion of f¬P, which can be difficult to do. For instance, if f¬P(n) is a polynomial in nof degree 5 or more then by the Abel-Ruffini theorem it is in general not solvable in radicals. However, if the sequence {f¬P(n)}is interleaved with another sequence {α(n)}that is more easily invertible (either analytically or computationally), then we can leverage this to more efficiently find the complementary sequence fP. More precisely, assume that we are given a real-valued increasing function αsuch that α(1) = 1 and f¬P(m−1) ≤α(m)≤f¬P(m) for all m. Note that we do not require that αis integer-valued. We define the increasing and integer-valued function h(n) = max{m∈N+:α(m)≤n}. The idea
INTEGERS: 25 (2025) 7 is to choose αsuch that his easier to compute than inverting f¬P. We can then compute fP(n) with one evaluation of h(n). Theorem 3. If f¬Pis such that f¬P(m+ 1) −f¬P(m)≥mfor all m, and α,h are as defined above, then fP(n) = n+h(n)+1 if n+h(n)≥f¬P(h(n) + 1) n+h(n)−1if n+h(n)≤f¬P(h(n)) n+h(n)otherwise. (2) Proof. The conditions on f¬Pimply that we can apply Theorem 2. Next note that n+h(n)≥f¬P(h(n) + 1) ≥f¬P(h(n)) + h(n) implies f¬P(h(n)) ≤nand that n+h(n)≤f¬P(h(n)) implies n<f¬P(h(n)) since h(n)≥1. The result then follows from the fact that C¬P(n) = h(n) if f¬P(h(n)) ≤nand C¬P(n) = h(n)−1 otherwise. This result can be further simplified depending on how close f¬Pand αare. Corollary 1. Given the hypothesis of Theorem 3, if f¬P(m)−α(m)< m, then fP(n) = (n+h(n)+1 if n+h(n)≥f¬P(h(n) + 1) n+h(n)otherwise. Proof. Note that by definiton of h,α(h(n)) ≤n. Out of the three conditions in Equation (2), the second condition is never satisfied by hypothesis since otherwise it reaches the contradiction n+h(n)−1< f¬P(h(n)) ≤h(n) + α(h(n)) −1≤n+h(n)−1. Corollary 2. Given the hypothesis of Theorem 3, if f¬P(m)−α(m)≥mfor m > 1, then fP(n) = (n+h(n)if n+h(n)> f¬P(h(n)) n+h(n)−1otherwise. Proof. Note that by definiton of h,α(h(n) + 1) > n. Out of the three conditions in Equation (2), the first condition is never satisfied by hypothesis since otherwise it reaches the contradiction n+h(n)≥f¬P(h(n) + 1) ≥h(n) + α(h(n) + 1) > n +h(n).
INTEGERS: 25 (2025) 8 To illustrate, we will use these results to obtain formulas for the non-k-gonal numbers, the non-k-gonal-pyramidal numbers, the non-k-simplex numbers, the nonsum-of-k-th-powers and the non-centered-k-gonal numbers. Furthermore, we choose hand αsuch that these formulas can be implemented algorithmically using integer arithmetic. 2.2.1. Non-k-gonal Numbers The n-th k-gonal numbers are defined as T(k, n) = (k−2)n(n−1) 2+nfor k≥2 and n≥1. For k= 2, the 2-gonal numbers are simply the natural numbers and thus there are no non-2-gonal numbers. In this section we give a general formula for the n-th non-k-gonal number for k≥3. Theorem 4. The n-th non-k-gonal number (k≥3) is given by a(n) = n+jq2n k−2k+ 1 if 2n > (k−2) jq2n k−2kjq2n k−2k+ 1 n+jq2n k−2kotherwise, (3) i.e., a(n) = n+$r2n k−2%+t2n > (k−2) $r2n k−2% $r2n k−2%+ 1!|. For 3≤k≤10, this can be written as a(n) = n+jq2n k−2+1 2k. Proof. Note that T(k, n + 1) −T(k, n)=(k−2)n+ 1 ≥nfor k≥3. Furthermore, T(k, n) is interleaved with α(n) = k−2 2(n−1)2with corresponding h(n) = $r2n k−2%+ 1. Since T(k, n)−α(n) = k 2(n−1) + 1 ≥n, we can apply Corollary 2and obtain Equation (3). Assume that k≤10. The inequality q2n k−2−jq2n k−2k≥1 2is equivalent to 2n≥(k−2) $r2n k−2% $r2n k−2%+ 1!+k−2 4. Since 2nand (k−2) jq2n k−2kjq2n k−2k+ 1are both even and k−2 4≤2 this is equivalent to 2n≥(k−2) jq2n k−2kjq2n k−2k+ 1+ 2 which in turn is equivalent to 2n > (k−2) jq2n k−2kjq2n k−2k+ 1which is exactly the condition in Equation (3).
INTEGERS: 25 (2025) 9 Recall that jq2n k−2kcan be computed as rj2n k−2k. For instance, in Python 3.x this can be implemented as isqrt(2*n//(k-2)). By setting k= 3 or k= 4 in Theorem 4, we get the following result. Corollary 3. The n-th non-triangular number is given by n+⌊√2n+1 2⌋. The n-th non-square number is given by n+⌊√n+1 2⌋. In discussing [13, Example 4], it was shown that the n-th non-square number is n+⌊√n+1 2⌋, while discussing [13, Example 6] it was shown that n+⌊pn+⌊√n⌋⌋ is a new formula for the n-th non-square number. Theorem 1along with Corollary 3 show the equivalence of these two formulas. In [13, Example 5] it was reported that the n-th non-triangular number is n+⌊√2n+1 2⌋which also follows from Corollary 3. Theorem 5. For k≥3, let rbe a real number such that −1≤r≤k−4. Then the n-th non-k-gonal number is given by a(n) = n+jq2n+r k−2k+ 1 if 2n > (k−2) jq2n+r k−2kjq2n+r k−2k+ 1 n+jq2n+r k−2kotherwise. (4) By choosing r=k+1 4−2, this can be simplified as a(n) = n+$r2n−2+⌊k+1 4⌋ k−2+1 2%. Proof. Let α(n) = k−2 2(n−1)2−r 2with corresponding h(n) = jq2n+r k−2k+ 1. Note that since r≥ −1, T(k, n)−α(n) = k 2(n−1)+1+r 2≥k 2(n−1) + 1 2≥nfor n > 1. Furthermore, α(n+ 1) −T(k, n) = (k−4)n 2−r 2≥0 so we can apply Corollary 2and obtain Equation (4). Next, q2n+r k−2−jq2n+r k−2k≥1 2is equivalent to 2n≥(k−2) $r2n+r k−2% $r2n+r k−2%+ 1!+k−2 4−r. By picking r=k−10 4=k+1 4−2, this ensures that k−2 4−r≤2. Furthermore, r≥ −1 for k≥3. Since 2nand (k−2) jq2n+r k−2kjq2n+r k−2k+ 1are both even this is equivalent to 2n≥(k−2) jq2n+r k−2kjq2n+r k−2k+ 1+ 2 which in turn is equivalent to 2n > (k−2) jq2n+r k−2kjq2n+r k−2k+ 1, i.e., the condition in Equation (4).
INTEGERS: 25 (2025) 16 Lemma 2. For each real number bsuch that ak−1 ak−1<b<ak−1 ak, ak(n+b)k≤a(n)≤ak(n+1+b)k for all sufficiently large n. If in addition k≥3, then ak(n+b)k+n≤a(n)≤ak(n+1+b)k for all sufficiently large n. Proof. Let ˜a(n) = Pk i=0 aini, then |a(n)−˜a(n)| ≤ 1. Next note that ak(n+b)k−˜a(n) = (akb−ak−1)nk−1+r1(n) where r1(n) is a polynomial of degree k−2 or less. Since akb−ak−1<0, it follows that ak(n+b)k−˜a(n)<−1 and thus ak(n+b)k−a(n)<0 for sufficiently large n. Similarly, ak(n+b+ 1)k−˜a(n)=(ak(b+ 1) −ak−1)nk−1+r2(n) for a polynomial r2(n) of degree k−2 or less. Since ak(b+1)−ak−1>0, this means that ak(n+b+ 1)k−˜a(n)>1 and thus ak(n+b+ 1)k−a(n)>0 for large enough n. Finally, if k≥3, the facts that ak(n+b)k+n−˜a(n)=(akb−ak−1)nk−1+n+r1(n), (akb−ak−1)<0, and n+r1(n) is a polynomial of degree k−2 or less imply that ak(n+b)k+n−˜a(n)<−1 and thus ak(n+b)k+n−a(n)<0 for large enough n. Note that a(n+ 1) −a(n)≥nfor large enough nsince ak>0. By choosing α(m) = ak(m+b)kwith corresponding h(n) = jk qn ak−bk, for large enough nthe n-th term of the complementary sequence to a(n) can be found using one evaluation of the function h. In particular, Theorem 3implies the following result. Theorem 13. Let {c(n)}be the complementary sequence to the sequence {a(n)}. If ak−1 ak−1<b<ak−1 akand h(n) = jk qn ak−bk, then there exists n0>0such that for all n≥n0, c(n) = n+h(n)+1 if n+h(n)≥a(h(n) + 1) n+h(n)−1if n+h(n)≤a(h(n)) n+h(n)otherwise. This can be implemented for the following special case using the integer k-th root function n→ ⌊ k √n⌋discussed in Section 1by choosing b=jak−1 akk.
INTEGERS: 25 (2025) 17 Corollary 6. Let {c(n)}be the complementary sequence to the sequence {a(n)}. If ak−1 akis not an integer, then there exists n0>0such that for all n≥n0, c(n) = n+jk qn akk−jak−1 akk+ 1 if n+jk qn akk−jak−1 akk≥ ajk qn akk−jak−1 akk+ 1 n+jk qn akk−jak−1 akk−1if n+jk qn akk−jak−1 akk≤ ajk qn akk−jak−1 akk n+jk qn akk−jak−1 akkotherwise. Similarly, Corollary 2implies the following result. Corollary 7. Let k≥3and let {c(n)}be the complementary sequence to the sequence {a(n)}. If ak−1 ak−1< b < ak−1 akand h(n) = jk qn ak−bk, then there exists n0>0such that for all n≥n0, c(n) = (n+h(n)if n+h(n)> a(h(n)) n+h(n)−1otherwise. Corollary 8. Let k≥3and let {c(n)}be the complementary sequence to the sequence {a(n)}. If ak−1 akis not an integer, then there exists n0>0such that for all n≥n0, c(n) = n+jk qn akk−jak−1 akkif n+jk qn akk−jak−1 akk> ajk qn akk−jak−1 akk n+jk qn akk−jak−1 akk−1otherwise. 2.2.8. Computing Characteristic Functions If his easily computable, then this can lead to an efficient algorithm to compute the characteristic function χ¬Pof f¬P. It is clear that if f¬P(m) = n, then h(n) = m. This implies that χ¬P(n) = 1 if and only if f¬P(h(n)) = n, i.e., χ¬P(n) = Jf¬P(h(n)) = nK.
INTEGERS: 25 (2025) 18 As an example, consider the characteristic function χ(n) of 4-simplex numbers, i.e., numbers of the form m 4for some m(OEIS A256436). Using h(n) = ⌊4 √24n⌋−1, we see that χ(n) = 1 if and only if n=⌊4 √24n⌋+2 4. This formula is simpler than the first formula in the entry for OEIS A256436 given by χ(n) = $p4√24n+ 1 + 5 2−3 2%− q4p24(n−1) + 1 + 5 2−3 2 . 2.3. Bisection Search Theorem 1shows that if f¬P(n) grows quadratically (or faster), then the number of steps needed to find fPusing the function iteration method is no more than two. Numerical experiments suggest a similar relationship for other grow rates. In particular, these experiments allow us to conjecture the following result. Conjecture 2. For k≥1, if f¬P(n) grows faster than n1+ 1 k, then the number of steps needed for the function iteration method to determine fP(n) is no more than k+ 1. For the complementary sequences where Theorem 1is applicable (and this include the sequences discussed in the above sections), the number of steps is bounded by a constant independent of n. For other types of sequences, this may not be the case. In these cases, a different method for finding fixed points is needed that is more efficient than the function iteration method. The function iteration method can take a large number of steps to converge when the set of integers that satisfy Pis sparse. In this case, it might be more optimal to use a bisection search method to find the fixed point of gn. To this end, we first find an interval [kmin, kmax] that bounds fP(n). Since we know that fP(n)≥n, we initially set kmin =n. If a better lower bound for fP(n) is known, this can be assigned to the initial kmin. Initially we can also set kmax =nunless we know a better lower or upper bound or an approximation of bfor fP(n) in which case we set kmax =b. We next double this initial value kmax repeatedly until gn(kmax)≤kmax. Then a bisection search is applied until the smallest fixed point is obtained. The pseudo code for this algorithm is shown in Algorithm 2. The number of steps needed to converge is on the order of log2(fP(n)) and is in general more efficient than the function iteration method in Section 2.1, especially when the numbers satisfying P are sparse. To illustrate this, we show in Figure 2the number of steps needed for these two methods to obtain fP(n) when P(m) denotes the logical statement “mis a product of exactly 6 distinct primes” (OEIS A067885). We see that the number of steps for the bisection method is much smaller than for the function iteration method.
INTEGERS: 25 (2025) 19 Algorithm 2 Bisection search on gn(x) to compute fP(n). Require: r∈N, gn(x)▷computes the minimal fixed point of gn. kmin ←n,kmax ←n. while gn(kmax)> kmax do kmax ←2kmax end while kmin ←max(kmin, kmax/2) while kmax −kmin >1do kmid =⌊(kmax +kmin)/2⌋ if g(kmid)≤kmid then kmax ←kmid else kmin ←kmid end if end while return kmax 100101102103104105106107108 n 102 103 104 number of steps P ( m ) denotes " m is a product of exactly 6 distinct primes" function iteration method bisection search Figure 2: Number of steps to find fP(n) when P(m) denotes the statement “mis a product of exactly 6 distinct primes”.
INTEGERS: 25 (2025) 20 2.4. Hybrid Method Since in several well known cases the number of steps that the function iteration method terminates in is small, we can take advantage of this by setting the initial kmin and kmax to g(r) n(n) for some small r, say r= 2. Here f(r)denotes the r-th iterate of the function f. The total number of steps is then rplus the number of steps of the bisection method. This is illustrated in Algorithm 3. Algorithm 3 Hybrid method. First riterations of the function iteration method is used to initialize the bisection search. Require: r∈N, gn(x)▷computes the minimal fixed point of gn. kmin ←g(r) n(n), kmax ←g(r) n(n). while gn(kmax)> kmax do kmax ←2kmax end while kmin ←max(kmin, kmax/2) while kmax −kmin >1do kmid =⌊(kmax +kmin)/2⌋ if g(kmid)≤kmid then kmax ←kmid else kmin ←kmid end if end while return kmax 3. The n-th Term of the Union or Difference of Two Sequences Consider the following scenario where {a(i)}and {b(i)}are disjoint sequences of integers with corresponding counting functions Caand Cb. For instance, the sequence acould be the set of square-free numbers and bthe perfect powers. The goal is to find the n-th element in the sorted list when aand bare sorted together. In this specific example of aand b, the joint sequence (denoted as c) is OEIS A304449. Other examples of such joint sequences are for instance, OEIS A000430,A006899, A089237,A126684,A168363,A174090, and A384419. Because the sequences aand bare disjoint, in this scenario, the counting function CP(n) for the joint sequence c is simply the sum of the counting functions of aand bgiven by Ca(n) + Cb(n) and the above algorithms can be used to find the n-th element of c. When the two sequences are not disjoint, CP(n) = Ca(n) + Cb(n)−Ca∩b(n) by the inclusion-exclusion principle and in some cases the intersection of the sequences can easily be determined. For instance, let pbe prime and consider the sequence
INTEGERS: 25 (2025) 21 of numbers ksuch that kkis a p-th power (OEIS A176693,A376379). If the prime factorization of kis k=Qipei i, then kk=Qipkei i. Thus, kkis a p-power if and only if kei≡0 (mod p). Since pis prime, the residue classes form a field, and this condition corresponds to when kis a multiple of por eiis a multiple of pfor all i, i.e., kis a p-th power. Thus, this sequence is the union of the multiples of pand the p-th powers with corresponding counting functions ⌊n/p⌋and ⌊p √n⌋. The counting function of their intersection is ⌊p √n/p⌋and thus the counting function of the union is ⌊n/p⌋+⌊p √n⌋−⌊p √n/p⌋. Similarly the counting function of the union of squares and powers of 2 (OEIS A188915) is given by ⌊√n⌋+⌊log2(n)⌋−⌊log2(n)/2⌋=⌊√n⌋+⌈log2(n)/2⌉. Similarly, the counting function of the sequence resulting from combining aand band removing their intersection is CP(n) = Ca(n) + Cb(n)−2Ca∩b(n). For an example see OEIS A377025. A complementary approach can be used to compute the difference of two sequences by noting that if bis a subsequence of a, then Ca\b(n) = Ca(n)−Cb(n). For instance, consider OEIS A388304 which is equal to OEIS A303606 without the terms of OEIS A177492. Similarly OEIS A388427 is equal to OEIS A303554 without the non-composites. In computer science, there is sometimes the need to find the n-th smallest element of an (large) unordered list of length lwithout having to sort the entire list. Typically, this is performed using a partial sort (see, e.g., the Quickselect algorithm [8]) which has an O(l) average performance. The above scenario can be considered a special case of this problem where the sequence can be decomposed into kdisjoint subsequences, the counting functions of the subsequences are efficiently computed, and the length of the list of elements is a priori unknown. Using the above algorithm leads to a running time depending on nunlike the partial sort algorithm which has a running time depending on the length lof the entire list. 4. Sequences of Repeated Terms Since fP(n), which is the complementary sequence of the sequence f¬P(n), can be viewed as the sequence of integers skipping the values of f¬P(n), we can consider the function fP(n)−nwhich lists consecutive integers, each one of which is repeated. For example, consider the sequence of non-square numbers (OEIS A000037) given by fP(n) = (2,3,5,6,7,8,10,11, . . .). The sequence fP(n)−nis (1,1,2,2,2,2,3,3, . . .), i.e., each integer mappears 2mtimes in the sequence (OEIS A000194). Since fP(n) = n+⌊√n+1 2⌋, this implies that the sequence (1,1,2,2,2,2,3,3, . . .) can be written as ⌊√n+1 2⌋. Similarly, the sequence fP(n) of the non-triangular numbers corresponds to a sequence fP(n)−nwhere each integer mappears mtimes and thus can be written explicitly as ⌊√2n+1 2⌋(OEIS A002024). This sequence has
INTEGERS: 25 (2025) 22 been studied in [11]. More generally, given a sequence of real numbers a1, a2, . . ., consider a sequence b(n) for n≥1 where each number am≥1 appears β(m) times consecutively: a1, a1,...a1 | {z } β(1) times , a2, a2,...a2 | {z } β(2) times , . . . The goal is to determine b(n) given n. The case of β(m) = md for a fixed dwas studied in [17]. We will consider the special case where ai=ias the approach to the general case is the same. Then fP(n) = b(n) + n−1 skips an integer after every additional β(i) numbers, meaning the i-th value skipped is β(i) + 1 plus the last number skipped. In other words f¬P(n) = Pn i=1(β(i) + 1) = n+Pn i=1 β(i). We can then apply the results and algorithms in the previous sections to find fP(n) and thus also b(n). Since f¬P(n)−f¬P(n−1) = β(n) + 1, if in addition β(m)≥m−2, then we can apply Theorems 1and 2to find b(n). To illustrate this approach, consider the sequence b(n) where each integer mappears m2times (OEIS A074279). Then f¬P(n) = n+Pn i=1 i2=n+n(n+ 1)(2n+ 1)/6 (OEIS A145066). Theorem 11 in Section 2.2.5 shows that fP(n) = (n+⌊3 √3n⌋if 6n > ⌊3 √3n⌋⌊3 √3n⌋+ 12⌊3 √3n⌋+ 1 n+⌊3 √3n⌋−1 otherwise. This implies that b(n) = fP(n)−n+1 = (⌊3 √3n⌋+ 1 if 6n > ⌊3 √3n⌋⌊3 √3n⌋+ 12⌊3 √3n⌋+ 1 ⌊3 √3n⌋otherwise. This formula is simpler than the formula for this sequence given in [21]. More generally, a sequence b(n) where each integer mrepeats mk−1times can be computed using the formula b(n) = (⌊k √kn⌋+ 1 if kn > Pk−1 j=0 k jB+ j⌊k √kn⌋k−j ⌊k √kn⌋otherwise. Again, these formulas require only one evaluation of ⌊k √kn⌋. This approach is used to find novel formulas for sequences of repeated integers such as OEIS A056556,A056557,A056558,A108581,A108582,A127321,A180447, A194847,A194848,A235463, and A360010. For instance, the sequence b(n) where each integer mis repeated m+3 3times (OEIS A127321) can be expressed as b(n) = ⌊4 p24(n+ 2)⌋−2 if n < ⌊4 √24(n+2)⌋+2 4 ⌊4 p24(n+ 2)⌋−1 otherwise.
INTEGERS: 25 (2025) 23 As mentioned above, the existence of a simple formula for fPdoes not imply a simple formula for f¬P. However, the existence of a simple formula for CPimplies the existence of a simple formula C¬Pwhich leads to efficient algorithms for both fPand f¬P. In Section 5, we look at other examples of logical statements P, mainly related to the prime factorizations of numbers, for which there are relatively efficient algorithms for computing CP(n) (and thus also C¬P(n)). 5. Some Other Explicitly Computable Counting Functions CP(n) The algorithm for computing the complementary sequence f¬P(n) requires the computation of the counting function CP(n). Note that CP(n) is the largest integer msuch that fP(m)≤n. If fP(n) is efficiently computed, then the fact that fPis strictly increasing means that CPcan be computed via a bisection search (similar to Algorithm 2for gn). An example of such functions fPis the case where P(m) denotes the statement “q(x) = mfor some positive integer x”, with qa strictly increasing polynomial with integer coefficients. In other cases, the function fP(n) is more easily enumerated sequentially in which case a naive approach would be to enumerate fP(n) in order to compute CP(n) and then solve for the fixed point of gn. However, this approach to find the complementary sequence f¬P(n) is less efficient than simply enumerating fP(n) and assigning the gaps between successive terms to f¬P(n). On the other hand, for several number theoretical statements P, computing CP(n) can be more efficient than enumerating fP(n). For instance, there exists algorithms for computing the prime counting function π(x) that are more efficient than enumerating all π(x) prime numbers less than or equal to x[2]. The formulas for the counting function of square-free numbers2given by ⌊√n⌋ X i=1 µ(i)⌊n/i2⌋ and the counting function of perfect powers3given by 1−⌊log2(x)⌋ X i=2 µ(i)(⌊i √x⌋−1) require the computation of the M¨obius function µ(n) which can be done efficiently using a sieve [3]. The counting functions of semiprimes and k-almost primes can be expressed using π(x) and are found in [1,24]. The counting function of square-free 2See [20] for a more efficient formula. 3See [18] for an alternative formula.
INTEGERS: 25 (2025) 24 k-almost primes (i.e., numbers that are products of kdistinct primes) has a similar form in terms of π(x) (see for instance, OEIS A067885). Another example where CP(n) is easily obtained is when P(m) denotes the statement “mis a repdigit in base b” (OEIS A139819). In this case, CP(n) = (b−1)⌊logb(n)⌋+(b−1)n b⌊logb(n)⌋+1 −1. 6. Conclusions We study algorithms to find the n-th integer that satisfies a certain condition Pvia a fixed point approach. We show that the function iteration method to solve this problem is equivalent to the Lambek-Moser method. Furthermore, we show that the two-step iteration of the Lambek-Moser algorithm requiring 2 evaluations of the counting function can be reduced to a single evaluation of the counting function. We show that the use of suitable simple interleaving functions can simplify the inversion of fPto obtain CP. We also present a bisection algorithm that is more efficient when the numbers satisfying Pare sparse. For a particular condition P, this approach is useful not only in finding the n-th term of the complementary sequences f¬Pbut also the n-th term of the original sequence fPitself when computing CP(n) is more efficient than enumerating fPup to n. For instance, we have used this approach in various Python programs to find the n-th term of OEIS sequences without enumerating all nterms. These sequences include perfect powers (OEIS A001597,A075090,A075109), prime powers (OEIS A000961,A246655,A246547,A025475), powerful numbers (OEIS A001694), k-full numbers (OEIS A036966,A036967,A069492,A069493), square-free numbers (OEIS A005117), semiprimes (OEIS A001358), square-free semiprimes (OEIS A006881), k-almost primes (OEIS A014612,A014613,A014614), square-free k-almost primes (OEIS A007304,A046386,A046387,A067885,A281222), Achilles numbers (OEIS A052486), orders of proper semifields and twisted fields (OEIS A088247,A088248), p-smooth numbers (OEIS A003586,A051037,A002473,A051038,A080197), numbers with at least one digit b−1 in base b(OEIS A074940,A337239,A337250), primes starting with digit b(OEIS A045707-A045715), sums of 3 squares (OEIS A000378), numbers with exactly kdivisors (OEIS A030513,A030515,A030626, A030627,A030632,A030633,A137484,A137485,A137488), selected sifting sequences [23] (OEIS A003159,A007417,A382744,A382745,A382746), characteristic functions (OEIS A256436,A387646), and complementary sequences (OEIS A004215,A024619,A052485,A029742,A007916,A002808,A013929,A100959, A139819,A374812,A090946,A057854,A185543,A138836,A138890,A325112, A059485,A279622,A145397,A302058,A376573,A183300,A387644,A388077). Finally, interested readers can obtain Python programs for the computation of
INTEGERS: 25 (2025) 25 the OEIS sequences discussed in this paper by accessing the entries of the corresponding OEIS sequences or the corresponding programs in the GitHub repository oeis-sequences [25]. References [1] D. Cri¸san and R. Erban, On the counting function of semiprimes, Integers 21 (2021), #A122. [2] M. Del´eglise and J. Rivat, Computing π(x): the Meissel, Lehmer, Lagarias, Miller, Odlyzko method, Math. Comp. 65 (213) (1996), 235–245. [3] M. Del´eglise and J. Rivat, Computing the summation of the M¨obius function, Exp. Math. 5 (4) (1996), 291–295. [4] B. Farhi, An explicit formula generating the non-Fibonacci numbers, preprint, arXiv: 1105.1127. [5] H. W. Gould, Non-Fibonacci numbers, Fibonacci Quart.,3(1965), 177–183. [6] R. L. Graham, D. E. Knuth, and O. Patashnik, Concrete Mathematics: A Foundation for Computer Science, Addison-Wesley, 2nd ed., 1994. [7] T. Heath, A History of Greek Mathematics, vol. 2, Clarendon Press, 1921. [8] C. A. R. Hoare, Algorithm 65: find, Commun. ACM 4(7) (1961), 321–322. [9] R. Honsberger, Essay 12, in Ingenuity in Mathematics, The Mathematical Association of America, 1970, 93–110. [10] K. E. Iverson, A Programming Language, Wiley, 1962. [11] D. E. Knuth, The Art of Computer Programming, Addison-Wesley, 1968. [12] D. E. Knuth, Two notes on notation, Amer. Math. Monthly 99 (5) (1992), 403–422. [13] J. Lambek and L. Moser, Inverse and complementary sequences of natural numbers, Amer. Math. Monthly 61 (7) (1954), 454–458. [14] K. MacMillan and J. Sondow, Proofs of power sum and binomial coefficient congruences via Pascal’s identity, Amer. Math. Monthly 118 (6) (2011), 549–551. [15] C. Mortici, Remarks on complementary sequences, Fibonacci Quart. 48 (4) (2010), 343–347. [16] R. D. Nelson, Sequences which omit powers, Math. Gaz. 72 (461) (1988), 208–211. [17] M. A. Nyblom, Some curious sequences involving floor and ceiling functions, Amer. Math. Monthly,109 (6) (2002), 559–564. [18] M. A. Nyblom, A counting function for the sequence of perfect powers, Austral. Math. Soc. Gaz. 33 (2006), 338–343. [19] The On-Line Encyclopedia of Integer Sequences, Available online at https://oeis.org. [20] J. Pawlewicz, Counting square-free numbers, preprint, arXiv:1107.4890. [21] B. Putievskiy, Integer sequences: Irregular arrays and intra-block permutations, preprint, arXiv:2310.18466.