ON SHANKS’ ALGORITHM FOR COMPUTING THE CONTINUED FRACTION OF log_{b} a

TERENCE JACKSON^{1} AND KEITH MATTHEWS^{2}

1 Department of Mathematics University of York, Heslington York YO105DD, England UK

                              thj1@york.ac.uk

2 Department of Mathematics University of Queensland Brisbane, Australia, 4072

                            krm@maths.uq.edu.au

Abstract. We give a more practical variant of Shanks’ 1954 algorithm for computing the continued fraction of log_{b} a, for integers a > b > 1, using the floor and ceiling functions and an integer parameter c > 1. The variant, when repeated for a few values of c = 10^{r}, enables one to guess if log_{b} a is rational and to find approximately r partial quotients.

1. Shanks’ algorithm

In his article [1], Shanks gave an algorithm for computing the partial quotients of log_{b} a, where a > b are positive integers greater than 1. Construct two sequences a_{0} = a, a_{1} = b, a_{2}, … and n_{0}, n_{1}, n_{2}, …, where the a_{i} are positive rationals and the n_{i} are positive integers, by the following rule: If i ≥ 1 and a_{i−1} > a_{i} > 1, then

a_{i}^{n_{i−1}} ≤ a_{i−1} < a_{i}^{n_{i−1}+1}  (1.1)
a_{i+1} = a_{i−1}/a_{i}^{n_{i−1}}.  (1.2)

Clearly (1.1) and (1.2) imply a_{i} > a_{i+1} ≥ 1. Also (1.1) implies a_{i} ≤ a_{i−1}^{1/n_{i−1}} for i ≥ 1 and hence by induction on i ≥ 0,

a_{i+1} ≤ a_{0}^{1/n_{0}···n_{i}}.  (1.3)

Also by induction on j ≥ 0, we get

a_{2j} = a_{0}^{r}/a_{1}^{s}, a_{2j+1} = a_{1}^{u}/a_{0}^{v},  (1.4)

where r and u are positive integers and s and v are non–negative integers.

Two possibilities arise:

2 TERENCE JACKSON^{1} AND KEITH MATTHEWS^{2}

(i) a_{r+1} = 1 for some r ≥ 1. Then equations (1.4) imply a relation a_{0}^{q} = a_{1}^{p} for positive integers p and q and so log_{a_{1}} a_{0} = p/q.

(ii) a_{i+1} > 1 for all i. In this case the decreasing sequence {a_{i}} tends to a ≥ 1. Also (1.3) implies a = 1, unless perhaps n_{i} = 1 for all sufficiently large i; but then (1.2) becomes a_{i+1} = a_{i−1}/a_{i} and hence a = a/a = 1.

If a_{i−1} > a_{i} > 1, then from (1.1) we have

n_{i−1} = ⌊\begin{matrix}log & a_{i−1} \\ log & a_{i} \end{matrix} ⌋.  (1.5)

Let x_{i} = log_{a_{i+1}} a_{i} if a_{i+1} > 1. Then we have

Lemma 1. If a_{i+2} > 1, then

x_{i} = n_{i} + 1/x_{i+1}.  (1.6)

Proof. From (1.2), we have

log a_{i+2} = log a_{i} − n_{i} log a_{i+1}  (1.7)
1 = \begin{matrix}log & a_{i} \\ log & a_{i+1} \end{matrix} · \begin{matrix}log & a_{i+1} \\ log & a_{i+2} \end{matrix} − n_{i} · \begin{matrix}log & a_{i+1} \\ log & a_{i+2} \end{matrix}  (1.8)
= x_{i}x_{i+1} − n_{i}x_{i+1},  (1.9)

from which (1.6) follows. □

From Lemma 1.1 and (1.5), we deduce

Lemma 2. (a) If log_{a_{1}} a_{0} is irrational, then

x_{i} = n_{i} + 1/x_{i+1} for all i ≥ 0.

(b) If log_{a_{1}} a_{0} is rational, with a_{r+1} = 1, then

x_{i} = \{ \begin{matrix}n_{i} + 1/x_{i+1}, & if 0 ≤ i < r − 1; \\ n_{r−1}, & if i = r − 1. \end{matrix}

In view of the equation log_{a_{1}} a_{0} = x_{0}, Lemma 2 leads immediately to

Corollary 1.

log_{a_{1}} a_{0} = \{ \begin{matrix}[n_{0}, n_{1}, …], & if log_{a_{1}} a_{0} is irrational; \\ [n_{0}, n_{1}, …, n_{r−1}] , & if log_{a_{1}} a_{0} is rational and a_{r+1} = 1. \end{matrix}  (1.10)

Remark. It is an easy exercise to show that for j ≥ 0,

a_{2j} = a_{0}^{q_{2j−2}}/a_{1}^{p_{2j−2}}, a_{2j+1} = a_{1}^{p_{2j−1}}a_{0}^{q_{2j−1}}  (1.11)

where p_{k}/q_{k} is the k–th convergent to log_{a_{1}} a_{0}.

Example 1. log_{2} 10: Here a_{0} = 10, a_{1} = 2. Then 2^{3} < 10 < 2^{4}, so n_{0} = 3 and a_{2} = 10/2^{3} = 1.25.

Further, 1.25^{3} < 2 < 1.25^{4}, so n_{1} = 3 and a_{3} = 2/1.25^{3} = 1.024.

Also, 1.024^{9} < 1.25 < 1.024^{10}, so n_{2} = 9 and

a_{4} = 1.25/1.024^{9}
= 1250000000000000000000000000/1237940039285380274899124224
= 1.0097419586 ⋯

Continuing in this fashion, we obtain Table 1 and log_{2} 10 = [3, 3, 9, 2, 2, 4, 6, 2, 1, 1, …].

\begin{matrix}i & n_{i} & a_{i} & p_{i}/q_{i} \\ 0 & 3 & 10 & 3/1 \\ 1 & 3 & 2 & 10/3 \\ 2 & 9 & 1.25 & 93/28 \\ 3 & 2 & 1.024 & 196/59 \\ 4 & 2 & 1.0097419586 ⋯ & 485/146 \\ 5 & 4 & 1.0043362776 ⋯ & 2136/643 \\ 6 & 6 & 1.0010415475 ⋯ & 13301/4004 \\ 7 & 2 & 1.0001628941 ⋯ & 28738/8651 \\ 8 & 1 & 1.0000637223 ⋯ & 42039/12655 \\ 9 & 1 & 1.0000354408 ⋯ & 70777/21306 \\ 10 & & 1.0000282805 ⋯ \\ 11 & & 1.0000071601 ⋯ \end{matrix} ; Table 1.

2. Some Pseudocode

In Table 2 we present pseudocode for the Shanks algorithm.

It soon becomes impractical to perform the calculations in multiprecision arithmetic, as the numerators and denominators a_{i} grow rapidly. If we truncate the decimal expansions of the a[i] to r places and represent a positive rational a as g(a)/10^{r}, where g(a) = ⌊10^{r}a⌋, the ratio aa/bb will be calculated as ⌊10^{r}g(aa)/g(bb)⌋. Working explicitly in integers, using the g(a), then results in algorithm 1, also depicted in Table 2, with c = 10^{r}, where int(x,y) equals ⌊x/y⌋, when x and y are integers.

As shown in the next section, the A[i] decrease strictly until they reach c. Also m[0]=n[0] and we can expect a number of the initial m[i] will be partial quotients. Naturally, the larger we take c, the more partial quotients will be produced.

Shanks’ algorithm  algorithm 1

input:  integers a>b>1  input:  integers a>b>1, c> 1
output:  n[0],n[1],. . . output:  m[0],m[1],. . .
s:= 0                    s:= 0
a[0]:= a; a[1]:= b       A[0]:= a*c; A[1]:= b*c
aa:= a[0]; bb:= a[1]     aa:= A[0]; bb:= A[1]
while(bb > 1){  while(bb > c){
i:=0                     i:=0
while(aa ≥ bb){  while(aa ≥ bb){
  aa:= aa/bb               aa:= int(aa*c,bb)
  i:= i+1                  i:= i+1

}  }

a[s+2]:= aa              A[s+2]:= aa
n[s]:= i                 m[s]:= i
t:= bb                   t:= bb
bb:= aa                  bb:= aa
aa:= t                   aa:= t
s:= s+1                  s:= s+1

}  }

Table 2.

3. Formal description of algorithm 1

We show in Theorem 2.1 below, that algorithm 1 will give the correct partial quotients when log_{a_{1}} a_{0} is rational and otherwise gives a parameterised sequence of integers which tend to the correct partial quotients when log_{a_{1}} a_{0} is irrational.

Algorithm 1 is now explicitly described. We define two integer sequences {A_{i,c}}, i = 0, …, l(c) and {m_{j,c}}, j = 0, …, l(c) − 2, as follows.

Let A_{0,c} = c · a_{0}, A_{1,c} = c · a_{1}. Then if i ≥ 1 and A_{i−1,c} > A_{i,c} > c, we define m_{i−1,c} and A_{i+1,c} by means of an intermediate sequence {B_{i,r,c}}, defined for r ≥ 0, by B_{i,0,c} = A_{i−1,c} and

B_{i,r+1,c} = ⌊\frac{cB_{i,r,c}}{A_{i,c}} ⌋ , r ≥ 0.  (3.1)

Then c ≤ B_{i,r+1,c} < B_{i,r,c}, if B_{i,r,c} ≥ A_{i,c} > c and hence there is a unique integer m = m_{i−1,c} ≥ 1 such that

B_{i,m,c} < A_{i,c} ≤ B_{i,m−1,c}.

Then we define A_{i+1,c} = B_{i,m,c}. Hence A_{i+1,c} ≥ c and the sequence {A_{i,c}} decreases strictly until A_{l(c),c} = c.

There are two possible outcomes, depending on whether or not log_{b}(a) is rational:

Theorem 2. (1) If log_{a_{1}} a_{0} is a rational number p/q with p > q ≥ 1 and gcd(p, q) = 1, then

(a) a_{0} = d^{p}, a_{1} = d^{q} for some positive integer d; (b) if p/q = [n_{0}, …, n_{r−1}], where n_{r−1} > 1 if r > 1, then

(i) A_{r+1,c} = c, a_{r+1} = 1;
(ii) A_{i,c} = c · a_{i} for 0 ≤ i ≤ r + 1;
(iii) m_{i,c} = n_{i} for 0 ≤ i ≤ r − 1.

(2) If log_{a_{1}} a_{0} is irrational, then (a) m_{0,c} = n_{0};

(b) l(c) → ∞ and for fixed i, A_{i,c}/c → a_{i} as c → ∞ and m_{i,c} = n_{i} for all large c.

Proof. 1(a) follows from the equation a_{1}^{p} = a_{0}^{q}.

1(b) is also straightforward on noticing that a_{i} is a power of d and that we are implicitly performing Euclid’s algorithm on the pair (p, q).

For 2(a), we have

a_{1}^{n_{0}} < a_{0} < a_{1}^{n_{0}+1}  (3.2)

and A_{0,c} = c · a_{0}, A_{1,c} = c · a_{1}. Also by induction on 0 ≤ r ≤ n_{0},

B_{1,r,c} ≥ ca_{1}^{n_{0}−r},  (3.3)
B_{1,r,c} ≤ \frac{ca_{0}}{a_{1}^{r}} .  (3.4)

Inequality (3.3) with r ≤ n_{0} − 1 gives B_{1,r,c} ≥ A_{1,c}, while inequality (3.4) with r = n_{0} gives

B_{1,n_{0},c} ≤ \frac{ca_{0}}{a_{1}^{n_{0}}} < ca_{1} = A_{1,c},

by inequality (3.2). Hence m_{0,c} = n_{0}.

For 2(b), we use induction on i ≥ 1 and assume l(c) ≥ i holds for all large c and that A_{i−1,c}/c → a_{i−1} and A_{i,c}/c → a_{i} as c → ∞. This is clearly true when i = 1.

By properties of the integer part symbol, equation (3.1) gives

\frac{c^{r}A_{i−1,c}}{A_{i,c}^{r}} − \frac{(1 − \frac{c^{r}}{A_{i,c}^{r}})}{1 − \frac{c}{A_{i,c}}} < B_{i,r,c} ≤ \frac{c^{r}A_{i−1,c}}{A_{i,c}^{r}} .  (3.5)

for r ≥ 0.

Hence for r < n_{i−1}, inequalities (3.5) give

B_{i,r,c}/c → a_{i−1}/a_{i}^{r} ≥ a_{i−1}/a_{i}^{n_{i−1}−1} > a_{i}.

Then, because A_{i,c}/c → a_{i}, it follows that B_{i,r,c} > A_{i,c} for all large c.

Also B_{i,n_{i−1},c}/c → a_{i−1}/a_{i}^{n_{i−1}} < a_{i}, so B_{i,n_{i−1},c} < A_{i,c} for all large c. Hence m_{i−1,c} = n_{i−1} for all large c. Also A_{i+1,c} = B_{i,n_{i−1},c} > c, so l(c) > i + 1 for all large c. Moreover A_{i+1,c}/c → a_{i−1}/a_{i}^{n_{i−1}} = a_{i+1} and the induction goes through. □

Example 3. Table 3 lists the sequences m_{0,c}, …, m_{l(c)−2,c} for c = 2^{u}, u = 1, …, 30, when a_{0} = 3, a_{1} = 2.

1,1,
1,1,1,
1,1,1,1,
1,1,1,2,
1,1,1,2,
1,1,1,2,3,
1,1,1,2,2,2,
1,1,1,2,2,2,1,
1,1,1,2,2,2,1,2,
1,1,1,2,2,3,2,3,
1,1,1,2,2,3,2,
1,1,1,2,2,3,1,2, 1, 1,1, 2,
1,1,1,2,2,3,1,3, 1, 1,3, 1,
1,1,1,2,2,3,1,4, 3, 1,
1,1,1,2,2,3,1,4, 1, 9,1,
1,1,1,2,2,3,1,5,24, 1,2,
1,1,1,2,2,3,1,5, 3, 1,1, 2,7,
1,1,1,2,2,3,1,5, 2, 1,1, 5,3, 1,
1,1,1,2,2,3,1,5, 2, 2,1, 3,1,16,
1,1,1,2,2,3,1,5, 2,15,1, 6,2
1,1,1,2,2,3,1,5, 2, 9,5, 1,2,
1,1,1,2,2,3,1,5, 2,13,1, 1,1, 6,  1, 2, 2,
1,1,1,2,2,3,1,5, 2,17,2, 7,8,
1,1,1,2,2,3,1,5, 2,19,1,49,2, 1,
1,1,1,2,2,3,1,5, 2,22,4, 8,3, 4,  1,
1,1,1,2,2,3,1,5, 2,22,2, 1,3, 1,  3, 8,
1,1,1,2,2,3,1,5, 2,22,1, 6,3, 1,  1, 3, 4,  2,
1,1,1,2,2,3,1,5, 2,23,2, 1,1, 2,  1,12,17,
1,1,1,2,2,3,1,5, 2,23,3, 2,2, 2,  2, 1, 3,  2,
1,1,1,2,2,3,1,5, 2,23,2, 1,7, 2,  2,14, 1,  1, 6,

Table 3.

In fact log_{2} 3 = [1, 1, 1, 2, 2, 3, 1, 5, 2, 23, 2, …].

4. A heuristic algorithm

We can replace the ⌊x⌋ function in equation (3.1) by ⌈x⌉, the least integer exceeding x. This produces an algorithm with similar properties to algorithm 1, with integer sequences

{A′_{i,c}}, i = 0, …, l′(c) and {m′_{j,c}}, j = 0, …, l′(c)−2. Here A_{0,c} = A′_{0,c} = a_{0} ·c, A_{1,c} = A′_{1,c} = a_{1} · c and m_{0,c} = m′_{0,c} = n_{0}. Then if i ≥ 1 and A′_{i−1,c} > A′_{i,c} > c, we define m′_{i−1,c} and A′_{i+1,c} by means of an intermediate sequence {B′_{i,r,c}}, defined for r ≥ 0, by B′_{i,0,c} = A′_{i−1,c} and

B′_{i,r+1,c} = ⌈\frac{cB′_{i,r,c}}{A′_{i,c}} ⌉ , r ≥ 0.  (4.1)

Then c ≤ B′_{i,r+1,c} < B′_{i,r,c}, if B′_{i,r,c} ≥ A′_{i,c} > c.

For

B′_{i,r+1,c} ≤ \frac{cB′_{i,r,c}}{A′_{i,c}} + 1

and

\frac{cB′_{i,r,c}}{A′_{i,c}} + 1 ≤ B′_{i,r,c} ⇔ cB′_{i,r,c} + A′_{i,c} ≤ A′_{i,c}B′_{i,r,c}

⇔ \frac{A′_{i,c}}{A′_{i,c} − c} ≤ B′_{i,r,c}.

The last inequality is certainly true if B′_{i,r,c} ≥ A′_{i,c} > c.

Hence there is a unique integer m′ = m′_{i−1,c} ≥ 1 such that

B′_{i,m′,c} < A′_{i,c} ≤ B′_{i,m′−1,c}.

Then we define A′_{i+1,c} = B′_{i,m′,c}. Hence A′_{i+1,c} ≥ c and the sequence {A′_{i,c}} decreases strictly until A′_{l′(c),c} = c.

If we perform the two computations simultaneously, the common initial elements of the sequences {m_{j,c}} and {m′_{k,c}} are likely to be partial quotients of log_{b}(a). With c = 10^{r} we expect roughly r partial quotients to be produced.

If l(c) = l′(c) and A_{j,c} = A′_{j,c} and m_{j,c} = m′_{j,c} for j = 0, …, l(c) − 2, then log_{b} a is likely to be rational.

In practice, to get a feeling of certainty regarding the output when c = 10^{r}, we also run the algorithm for c = 10^{t}, r − 5 ≤ t ≤ r + 5.

Example 4. Table 4 lists the common values of m_{i,c} and m′_{i,c}, when a = 3, b = 2 and c = 2^{r}, 1 ≤ r ≤ 31. It seems likely that only partial quotients are produced for all r ≥ 1.

 1: 1
 2: 1
 3: 1,1,1
 4: 1,1,1
 5: 1,1,1,2
 6: 1,1,1,2
 7: 1,1,1,2,2
 8: 1,1,1,2,2
 9: 1,1,1,2,2
10: 1,1,1,2,2
11: 1,1,1,2,2
12: 1,1,1,2,2
13: 1,1,1,2,2,3,1
14: 1,1,1,2,2,3,1
15: 1,1,1,2,2,3,1
16: 1,1,1,2,2,3,1,5
17: 1,1,1,2,2,3,1,5
18: 1,1,1,2,2,3,1,5
19: 1,1,1,2,2,3,1,5,2
20: 1,1,1,2,2,3,1,5
21: 1,1,1,2,2,3,1,5,2
22: 1,1,1,2,2,3,1,5,2
23: 1,1,1,2,2,3,1,5,2
24: 1,1,1,2,2,3,1,5,2
25: 1,1,1,2,2,3,1,5,2
26: 1,1,1,2,2,3,1,5,2
27: 1,1,1,2,2,3,1,5,2
28: 1,1,1,2,2,3,1,5,2,23
29: 1,1,1,2,2,3,1,5,2,23
30: 1,1,1,2,2,3,1,5,2,23,2
31: 1,1,1,2,2,3,1,5,2,23,2

Table 4. a = 3, b = 2, c = 2^{r}, 1 ≤ r ≤ 31.

Example 5. Table 5 lists the common values of m_{i,c} and m′_{i,c}, when a = 34, b = 2 and c = 10^{r}, 1 ≤ r ≤ 20. Partial quotients are not always produced, as is seen from lines 9,14 and 17.

1: 1,2,2
2: 1,2,2,1,1
3: 1,2,2,1,1,2
4: 1,2,2,1,1,2
5: 1,2,2,1,1,2,3,1
6: 1,2,2,1,1,2,3,1,8,1
7: 1,2,2,1,1,2,3,1,8,1,1
8: 1,2,2,1,1,2,3,1,8,1,1,2
9: 1,2,2,1,1,2,3,1,8,1,1,2,2,1,13,3,2,32,7
10:1,2,2,1,1,2,3,1,8,1,1,2,2,1
11:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1
12:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1
13:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13
14:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,3
15:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2
16:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2,2
17:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2,2,18,1,1,1,1,1
18:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2,2,17,1
19:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2,2,17,1
20:1,2,2,1,1,2,3,1,8,1,1,2,2,1,12,1,13,3,2,2,17,1

Table 5. a = 34, b = 12, c = 10^{r}, r = 1, …, 20.

5. Acknowledgement

The second author is grateful for the hospitality provided by the School of Mathematical Sciences, ANU, where research for part of this paper was carried out.

References

1. D. Shanks, A logarithm algorithm, Math. Tables and Other Aids to Computation 8 (1954), 60–64.

2000 Mathematics Subject Classification: 11D09. Keywords: Shanks’ algorithm, continued fraction, log, heuristic algorithm

(Concerned with sequence \underline{A028507}.)

Received November 19, 2002; revised version received December 6, 2002. Published in Journal of Integer Sequences December 10, 2002.

Return to Journal of Integer Sequences home page.