Diophantine Approximation, Ostrowski Numeration and the Double-Base Number System

A partition of x > 0 of the form x = P i 2 a i 3 b i with distinct parts is called a double-base expansion of x . Such a representation can be obtained using a greedy approach, assuming one can efﬁciently compute the largest { 2 , 3 } -integer, i.e., a number of the form 2 a 3 b , less than or equal to x . In order to solve this problem, we propose an algorithm based on continued fractions in the vein of the Ostrowski number system, we prove its correctness and we analyse its complexity. In a second part, we present some experimental results on the length of double-base expansions when only a few iterations of our algorithm are performed.


Introduction
In this paper, we consider representations of integers in two coprime bases.More exactly, we analyze decompositions of positive integers in the form with a i , b i ≥ 0 and (a i , b i ) = (a j , b j ) if i = j.
Apart from its practical interest in real life applications such as digital signal processing [13] and cryptography [10,14], this so-called double-base number system (DBNS) has many intriguing properties and leads to interesting Diophantine approximation problems.We focus on bases 2 and 3 but most of the results and algorithms presented in this paper remain valid for arbitrary bases.
Clearly, such a decomposition always exists, the binary expansion is a special case of (1).In fact, this number system is even extremely redundant.For a given integer x > 0, the number f (x) of DBNS 1365-8050 c 2009 Discrete Mathematics and Theoretical Computer Science (DMTCS), Nancy, France representations is given by (see [11] for more details): otherwise.
For large numbers, such as those used in cryptographic applications, finding a representation of minimal length, in a reasonable amount of time, seems very hard.Fortunately, one can use a greedy approach (see Algorithm 1) to find a fairly sparse representation very quickly, based on the determination of the best default approximation of x (i.e., the largest integer ≤ x) of the form z = 2 a 3 b .

Algorithm 1 Greedy decomposition
Find the best default approximation of x of the form z = 2 a 3 b

3:
Print (a, b) 4: x := x − z Although it sometimes fails in finding a minimal representation (the smallest example is 41: the minimal representation is 32 + 9, whereas Algorithm 1 returns 41 = 36 + 4 + 1), it is very easy to implement.More importantly, it guarantees an expansion of sublinear length ℓ.Indeed, an important theorem from Dimitrov [12] states that ℓ is in O(log x/ log log x).The proof is based on a result by Tijdeman [26], which says that there exists an absolute constant C such that there is always a number of the form 2 a 3 b in the interval x − x/(log x) C , x .See [12] for a complete proof.
In this paper we investigate the problem of finding the best default approximation of x of the form z = 2 a 3 b .This operation is clearly of crucial importance for Algorithm 1.We present an algorithm as well as inhomogeneous approximation results in the vein of a number system from Ostrowski [20], which uses the series of convergents of the continued fraction expansion of an irrational number.
We introduce the problem in Diophantine terms in Section 2 and recall basic facts on Ostrowski's numeration in Section 3. In Section 4, we describe our algorithm and prove its correctness.Moreover, we show that for all x of bounded binary length, our algorithm terminates in O(log log x) iterations.The natural extension to signed digits and some implementation considerations are discussed in Section 5. We then conclude by presenting some numerical experiments in Section 6.The present paper is an extended version of [7].

Problem and notation
Let x ≥ 2 be a given positive integer.We want to find two integers a, b ≥ 0 such that 2 a 3 b ≤ x and among the solutions to this problem, 2 a 3 b is the largest possible value.In other words ( Equivalently, one can solve the following Diophantine inequality in non-negative integers In the following, {t} denotes the fractional part of t.Graphically, the solutions to cα + d ≤ log 3 x are the points (c, d) with non-negative integer coordinates located under the line of equation These are the integer points in the gray area of Figure 1.Among those points, the best left approximation of x we are looking for is the point (a, b), which has the smallest vertical distance to the line, that we denote by δ(u) = {−αu + log 3 x}, for u ≥ 0.
A naive approach consists in computing δ(u) for every integer 0 ≤ u ≤ log 2 x and to keep the point (u, v) which gives the smallest value.Since log 3 x < log 2 x, a little more efficient solution is to exchange the coordinate axis, i.e., to consider the line v ′ = u ′ /α − log 2 x, and to keep the closest point among all integers 0 ≤ u ′ ≤ log 3 x (we perform here the change of coordinates system Let us note that in some cases, it may be interesting to compute the best right approximation, i.e., the smallest integer of the form 2 a 3 b greater than x.Observe that it can be easily obtained by symmetry, by looking for the point (a ′ , b ′ ) located under and with the smallest vertical distance to the line of equation v ′ = −αu ′ + α⌈log 2 x⌉ + ⌈log 3 x⌉ − log 3 (x) (the dashed line in Figure 1), and to apply the change of The best left approximation is the point of coordinates (2,4).Similarly the closest point below the dashed line is the point (2,5), which by symmetry gives the best right approximation: (9 − 2, 6 − 5) = (7, 1).We easily verify that 2α + 4 ≈ 5.26 < log 3 x < 7α + 1 ≈ 5.41.
The complexity of all those graphical approaches is O(log x).In the next sections, we give an algorithm which solves the problem in O(log log x).
In the following, we define β = {log 3 x}.We assume 0 < β < 1.Indeed, if β = 0, then we take a = 0 and b = log 3 x ∈ N. Clearly, if (a, b) is a maximal solution of (3), then we have Let us briefly describe the strategy followed in the present paper.We set k = a and l = ⌊log 3 x⌋ − b.We have 0 ≤ kα − l ≤ β.
We thus are looking for (k, l) ∈ N 2 such that 0 ≤ kα − l ≤ β and We shall define two increasing sequences (k n ) n≥0 , (l n ) n≥0 such that, for some j > 0, k j = k, l j = l satisfy (4).We then set a = k and b = ⌊log 3 x⌋ − l to get the solution of (2).
We will consider this classical inhomogeneous approximation problem in Section 4, where the sequence of inhomogeneous best approximations of β by points of the form N α will be explicitly given.For that purpose, we first introduce Ostrowski's number system.

Ostrowski's number system
In this section we introduce a number system due to Ostrowski [20].Let us first recall some basic facts about continued fractions (see [17] for more details).
A simple continued fraction is an expression of the form , where a 0 = ⌊α⌋ and a 1 , a 2 , . . .are integers ≥ 1.The sequence (a n ) n∈N of partial quotients can either be finite or infinite.
Every rational number α can be expressed as a finite simple continued fraction, denoted by α = [a 0 , a 1 , . . ., a n ] in the more convenient compact notation.Similarly, every irrational number α can be expressed uniquely as an infinite simple continued fraction.We write α = [a 0 , a 1 , a 2 , . . .].For example, the continued fraction decomposition of log 2/ log 3 = 0.630929753571 . . . is The rational number obtained by restricting the continued fraction of α to its first n + 1 partial quotients is called the nth convergent of α.The integers p n , q n are easily computable; we have p −1 = 1, q −1 = 0, p 0 = a 0 , q 0 = 1, and It is well known that the sequences (p n ) n≥0 , (q n ) n≥0 satisfy lim n→∞ p n /q n = α and p n+1 q n −p n q n+1 = (−1) n .Moreover, simple continued fractions provide the sequence of best rational approximations of an irrational number.
Ostrowski's number system [20] (see [2] for an introduction) is associated with the numeration scale (q n ) n≥0 of denominators of the convergents of the continued fraction expansion of an irrational number 0 < α < 1.The following proposition holds.
Ostrowski's representation of integers can be extended to real numbers (see e.g. the survey [5]).The base is given by the sequence (θ n ) n≥0 , where θ n = q n α − p n .Note that the sign of θ n is equal to (−1) n .
Proposition 2 Let 0 < α < 1 be an irrational number.Every real number −α ≤ β < 1 − α can be written uniquely in the form , and b k = a k for infinitely many odd integers.
Proposition 2 can be used to approximate β modulo 1 by numbers of the form N α, with N ∈ N. Indeed, one verifies that the integers provides a series of arbitrarily good approximations to β on both sides, i.e., the fractional part of N n α can either be greater or smaller than β.
If we are only interested in the best left approximations to β, which is our case, we can express β in base (|θ n |) n≥0 (see Section 4 below for the exact definition of best left approximations).The following proposition holds.Proposition 3 Let 0 < α < 1 be an irrational number.Every real number 0 ≤ β < 1 can be written uniquely in the form for infinitely many even and odd integers.
In this case, the sequence of best left approximations is more difficult to define due to the alternating signs of θ n .Indeed, the corresponding numeration system on integers is defined with respect to the numeration scale ((−1) n q n ) n≥0 .Nevertheless, none of these expansions provides the sequence of best left approximations (see Remark 3 below for more details).In the next sections, we present an algorithm inspired by [22] and [6], we prove its correctness and analyse its complexity.

The sequence of inhomogeneous best approximations of β
From now on we assume α irrational, 0 < α < 1 and 0 < β ≤ 1. Inhomogeneous left approximations of β are numbers of the form kα + l ≤ β, with k, l integers.Clearly, there exist infinitely many such approximations.We want to define two increasing sequences of integers and, furthermore, ∀k < k n+1 , k = k n , and ∀l ∈ Z such that 0 ≤ kα − l ≤ β, we have For simplicity, we define f n = |θ n | for all n.We have One has f n > 0 for all n, since α is irrational.Hence, the sequence (f n−1 + f n ) n≥0 is strictly decreasing and tends to zero.Since 0 < β ≤ 1, there exists a unique integer n ≥ 0 such that Before we give the algorithm that defines the series of inhomogeneous best left approximations of β, we need to prove the following two lemmas.
Lemma 1 Let 0 < β ≤ 1 and (f n ) n≥0 be defined as above.There exists a unique integer n ≥ 0, a unique integer c ≥ 1, and a unique e ∈ R such that with 0 < e ≤ f n , and and n is uniquely determined since (f n + f n+1 ) n≥−1 is strictly decreasing.Now, let n ∈ N be the unique integer that satisfies Let 0 < β ≤ 1 and (f n ) n≥0 be defined as above, and let n, c, e be the unique numbers satisfying (7).We define k, l ∈ N by setting Then we have 0 < kα − l < β.
Proof: Let us prove that 0 < β − (kα − l) < β.Assume first that n is even.In this case, Algorithm 2 Inhomogeneous best left approximations to β Input : Two real numbers 0 < α < 1 and 0 < β ≤ 1, with α irrational Output : The infinite sequence (k i , l i ) i≥1 of best left approximations to β with 0 < k i α − l i < β, ∀i 1: (k 0 , l 0 ) := (0, 0) 2: while true do if n i is even then 5: else 7: This algorithm is inspired by [22]; similar ones can be found in [23,24,25].Note that β−(k i+1 α−l i+1 ) is equal to e i if n i is odd, and to (c i − 1)f ni + f ni+1 + e i , if n i is even.Hence, we may have n i+1 = n i .This happens if and only if n i is even and c i > 1; it therefore happens (c i − 1) times before the sequence (n i ) keeps growing towards +∞.We thus have The following proposition proves the correctness of Algorithm 2.
Proposition 4 Let 0 < α < 1 and 0 < β ≤ 1 be given.We assume α irrational.The increasing sequences of integers (k i ) i≥0 and (l i ) i≥0 provided by Algorithm 2 satisfy and furthermore, for all i, for all k such that k i < k < k i+1 , and for all l ∈ Z such that 0 Proof: We first prove (9).From Lemma 2, we know that, for all i, 0 < k i α − l i < β.We consider two cases: Let us now prove (10).Let us assume that k i < k < k i+1 , and What we prove in the next two cases depending on the parity of 1. Assume first that n i is even.We have Thus, what remains to be proved is that one has |k i+1 α − l i+1 − kα + l| > f ni (we use the best approximation property of continued fractions, see e.g.[8]).Next we show that (k i+1 α − l i+1 − kα + l) cannot be negative by considering two cases.Note first that, from (5) and line 3 of Algorithm 2, we have 2. If we now assume that n i is odd, we have Here, what remains to be proved is that Proposition 5 Let x ≥ 2 be a fixed positive integer.Let α = log 3 2; note that α is irrational and 0 < α < 1.Let β = {log 3 x}.We assume furthermore that β = 0. Let (k n , l n ) n∈N be the sequence of inhomogeneous best left approximations of β.Let v be the uniquely defined integer that satisfies k v ≤ ⌊log 2 x⌋ < k v+1 .Then Proof: From the proof of Proposition 4, we know that Algorithm 2 returns the infinite sequence (k i , l i ) i≥0 of best left approximations to β = {log 3 x}.In order to find (k, l) as above, we only need to perform a finite number of iterations.In fact, this number of iterations is given by the smallest integer v such that ⌊log 2 x⌋ < k v , or equivalently, the smallest integer v such that ⌊log 3 x⌋ < l v .We use this bound in our complexity analysis (see proof of Proposition 6).✷ Remark 1 As pointed out in the proof of Proposition 5, we only need to perform a finite number of iterations of Algorithm 2. Remark that this number of iterations depends on the sequence of partial quotients in the continued fraction expansion of α.For the same reason, we can also use a rational approximation of α = log 3 2, say p w /q w , given by the wth convergent.Note that the required precision for the rational approximation of α depends on the binary length of the input x.This is why we state our complexity result for integers x with bounded binary length.We analyse the required precision in Remark 2.
Proposition 6 Let m > 0. For all x having at most m bits, the finite version of Algorithm 2 (see Remark 4,above) terminates in O(log log x) iterations.
Proof: We use the notation of Algorithm 2. According to (8), we obtain and thus Let w(x) = max i {n i }.Since a ni+1 − c i ≥ 0, for all n i ≥ 1, and c i ≥ 1, for all n i , we have Moreover, we also know that the sequence (q i ) i≥0 grows at least as fast as the sequence of Fibonacci numbers given by the denominators of the continued fraction expansion of (1 + √ 5)/2 = [1, 1, 1, 1, . . .] (all the partial quotients are equal to 1).Therefore, from Binet's formula (see [28]), we deduce that there exists a constant C > 0 such that ∀x, w(x) < C log log x.
Consequently, there exists C ′ such that for all x having at most m bits, w(x) < C ′ log m.We now set A = max i=1... C ′ log m {a i }.For any m-bit integer x, the number of iterations s(x) satisfies s(x) ≤ Aw(x).This concludes the proof.✷

Corollary 1
The greedy algorithm combined with Algorithm 2 is optimal.
Proof: Since the greedy algorithm produces decompositions with O(log x/ log log x) terms according to [12], its overall asymptotic complexity is O(log x), i.e., it reaches the lower bound for a recoding algorithm (all the bits/digits of x have to be scanned).✷ Remark 2 Let us analyse the required precision for α in terms of the bitlength of x.From the proof of Proposition 6, we have and c i ≥ 1 for all n i .Therefore, if x is an m-bit integer, we only need the denominators q i such that q i ≤ m for all i ≥ 1.Let w = max{i : q i ≤ m}.Table 1 gives some useful numerical values for w.Now, we need to know how many bits of precision µ are given by α = pw qw , the wth rational approximation of α.It is well known (see [17]) that for all n From (11), we can e.g.deduce that Tab.1: Number w of partial quotients ai, and convergents pi/qi to be computed based on the size of input x.The last column gives the number of bits of precision provided by pw/qw, the wth rational approximation of α.

Remark 3
We have seen that Algorithm 2 provides an expansion of β of the form Let us stress the fact that sequences (n i ) i≥0 and (c i ) i≥0 cannot be directly derived from the sequences of digits provided by Proposition 2 and 3 by performing finitely many operations.Note that these latter sequences of digits can be generated in an effective way by maps defined on the unit square as skew products of the Gauss map, such as described in [15] and [16].

Signed expansions
It is natural to consider signed expansions of the form The following proposition shows that it is not possible to define a greedy algorithm based on a similar logarithm reduction approach, that is, by working with the sequence of two-sided inhomogeneous best approximations of β.More precisely, Proposition 7 below shows that it may give an erroneous solution infinitely often; this happens only when the logarithm reduction approach returns a solution greater than log x whereas the best approximation of x is possibly less than x.
Proposition 7 Let (K n ) n∈N be the increasing sequence of integers of the form 2 a 3 b , with a, b ∈ N.For any positive integer n, we define N n as the number of positive integers x such that The sequence (N n ) n∈N tends towards infinity.
For all n, we set δ n = K n+1 K n .One has According to [26], there exist C, C ′ with C > 0, C ′ > 0 such that for all n We thus deduce that lim n→+∞ N n = +∞.✷ Nevertheless, we can apply the following strategy in the signed case.Let x ≥ 2 be a given positive integer.We first find a, b ≥ 0 such that 2 a 3 b ≤ x and among the solutions to this problem, 2 a 3 b is the largest possible value: We similarly find a ′ , b ′ ≥ 0 such that 2 a ′ 3 b ′ ≥ x and among the solutions to this problem, 2 a ′ 3 b ′ is the smallest possible value: For this latter problem, we can either deal with an analogous algorithm providing the sequence of inhomogeneous best right approximations of β, or work with 1 −α and 1 −β, with the notation of the previous sections.It then remains to compare x − 2 a 3 b and 2 a ′ 3 b ′ − x, and to take the smallest of both values.

An approximate greedy algorithm
The greedy algorithm is not optimal, in the sense that it does not necessarily produce a DBNS representation of minimal length.However, within this algorithm, we proved that finding the largest number of the form 2 a 3 b less than equal to x is optimal in complexity.It seems then natural to investigate the potential advantages of an approximate greedy algorithm, where we only perform a few iterations of Algorithm 2 to define a "good" integer of the form 2 a 3 b , although not the largest, less than or equal to x.
The greedy decomposition described in Algorithm 1 can easily be adapted to compute a DBNS expansion of x, where only d iterations of Algorithm 2 are performed at each step to define the term z ≤ x.In the following, we shall refer to d as the depth.
In order to visualize DBNS expansions, we use a two-dimensional array where the columns represent the powers of 2 while the rows represent the powers of 3. A black square at position (i, j) indicates that the term 2 i 3 j is part of the DBNS decomposition.(The upper-left corner corresponds to the term 2 0 3 0 = 1.) Figure 2 shows the DBNS representations obtained for x = 23832098195 when the greedy decomposition algorithm is applied at depth 1,2, and 3. (Note that depth 3 is equivalent to the complete greedy decomposition.) As expected, the length of the DBNS decomposition decreases as the depth increases; we get 12 terms at depth 1; 8 terms at depth 2; and 7 terms at depth 3. We note also that the largest binary exponent (the number of columns) is smaller for small depths.These results were to be expected.However, the following questions seem difficult to answer in the general case.By using a greedy decomposition algorithm at a given depth d, how many summands are required to represent x compared to the solution provided by the complete greedy approach?What decrease (resp.increase) can we expect on the largest binary (resp.ternary) exponents?In order to answer those questions, we have implemented the algorithm at different depths, for a thousand randomly chosen numbers of sizes ranging from 128 to 512 bits.The difference in length between the full-greedy approach and the decomposition at fixed depths is shown in Figures 3 to 5. Our experiments are summarized in Table 4.
We first note that, for all our tests (up to 512-bit integers), d = 5 is equivalent to the complete greedy algorithm; which means that we perform 6 iterations of Algorithm 2 to get an approximation which corresponds to a negative ternary exponent (the stop condition, cf.Remark 4).With the greedy algorithm at fixed depth d, it is possible that the stop condition occurs before we reach the desired depth, or would occur at the (d + 1)th iteration.Only in those two cases, we get the optimal solution; i.e, the largest  Tab.4: Average distance to the full-greedy solution at various depths.Note that depth 5 is equivalent to the fullgreedy algorithm in all our experiments 2 a 3 b ≤ x.It is also interesting to note that at medium depth (d = 3 or 4), the length of the representations can be very close to the full greedy ones, and in some cases even shorter.
For some applications, we might not need the DBNS decomposition given by the complete greedy approach, but rather satisfy with an approximated one, with slightly more terms, especially if the time required for the conversion is reduced.Another consequence of using the greedy decomposition at fixed depth is that the binary exponent is likely to be smaller at small depths than with the complete greedy decomposition.It is clear that the ternary exponent increases at the same time, but clearly not as fast (the ratio is α ≃ 0.63).

Discussions
The question of the generalization to more than two bases, such as 2 a 3 b 5 c is natural.Nevertheless, there is no canonical generalization of Ostrowski numeration to higher dimensions.This is mainly due to the fact that there is no canonical notion of multidimensional continued fractions.To remedy to the lack of a satisfactory tool replacing continued fractions, several approaches are possible.One can consider e.g.lattice reduction algorithms: let us quote the computation of the n-th Hamming number in [21] based on results from [4].Let us recall that Hamming numbers are numbers of the form 2 a 3 b 5 c for a, b, c ∈ N, and that they are called after Hamming who asked for an efficient algorithm to generate in order the list of these numbers.The problem was popularized by E. W. Dijkstra in [9].
Algorithm 2 relies on the so-called three-gap theorem (see [22] and also the survey [5]).In the same vein, an algorithm is given in [18] which computes the points in Z 2 located at a distance smaller than a given ε from a segment.This latter algorithm is based on the three-distance theorem, which can be considered as a dual version of the three-gap theorem.As an application, an algorithm which produces worst cases when trying to round the elementary functions in floating-point arithmetic is given in [19].
Instead of working with integer bases, one can also consider numeration in rational base (2/3).The dynamical and arithmetical behaviour of such a numeration system differs drastically.Let us quote [1] which considers connections to number theory and applications to the problem of the distribution of the fractional part of the powers of rational numbers.
with the constraint that no other integers c, d ≥ 0 give a better left approximation to log x; meaning that if (a, b) is a maximal solution of (3), then for all integers c, d with (c, d) = (a, b), we have cα + d < aα + b ≤ log 3 x, where α = log 3 2.