Shanks' logarithm on the Stern–Brocot tree

Shanks' logarithm

Logarithms were rolled out in Europe in the early 17th century. In the form of tables and later slide rulers, they were primary tools for high-precision calculations in navigation and engineering until they moved on to the the electronic calculator in the 20th century.

Daniel Shanks was rightly surprised to discover an apparently new, very simple algorithm for calculating logarithms in 1954 — three centuries into their history.

It is a recursive algorithm, similar to Euclid's algorithm, that successively outputs terms of the continued fraction of the logarithm.

Let's use it to calculate the logarithm of to the base .

We first look which power of is the largest below . In this case it is the second power, , and we got out first term of the continued fraction, the exponent :

We now recur as if we'd want to calculate the logarithm of to the base . That means, we now need to find the largest power of below 2. This is , because so we can output the next term

And we recur again now looking for the logarithm of to the base , which would output 9.

And so on.

The interested reader can try out the small application below. The app is preset to calculate the logarithm of to the base as we did above. Pressing ⏯︎ starts or pauses the calculation. Pressing ⏮︎ will reset the algorithm. You can enter other numbers in the two input fields.

⚠ If the base is smaller than the argument (a<b), this browser tab will crash.^⚠

We can now succesively calculate the elements of the continued fraction of any logarithm of two integers .

For some logarithms, like the, the algorithm will return a 3 and stop. But in most cases logarithms are irrational numbers, as in the example . Then the continued fractions will never end.

The case and, more general, can be covered using the identities

However, it is not obvious how to use Shanks' algorithm on irrational numbers themselves.

Suppose we want to compute the first three elements of the continued fraction of . The continued fractions of and are known, but we don't not know in advance how much of them we need to get the desired precision. A trial-and-error method could be used to successively use more terms of and until the 3rd element of the resulting sequence no longer changes.

It would be more elegant, if there were a version of Shanks' algorithm that worked only on continued fractions, without having to switch back and forth between the continued and the reduced fractions.

Fortunately, two decades after Shanks' discovery, another very simple and beautiful algorithm was discovered — lazy continued fraction arithmetic.

In the following, I will show how an extension of Shanks' logarithm to the real numbers comes as a free lunch using lazy arithmetic.

Adaptation to lazy arithmetic

How to do arithmetic directly on continued fractions in a lazy way was figured out in 1972 by Bill Gosper — considered by some to be the first hacker.^g

Let's say we want to calculate the product of two, potentially infinite, continued fractions. Gosper's algorithm will output the product term by term, using only as many terms from the operands as needed for the next term of the result.

Explaining the algorithm in detail is beyond my scope here, but to get an idea of how it works, I'd like to point the reader to a web visualization. ^v

Continued fractions are isomorphic to paths in the Stern–Brocot tree

and thus the arithmetic can be adapted to them as well. This was done by Milad Niqui in 2004 [Niqui].

The generalization of Shanks' logarithm to the Stern–Brocot tree is remarkably simple, as can be seen in a stripped-down Clojure version of the algorithm.^c

  (defn shanks
  	"- `mem` holds the power of a
  	 - `dir` is either L or R, the value of the next output bit"
    [a b mem dir]
  	(lazy-seq
  	  (cond (< (* a mem) b)
  		      ;; output L or R
  		  		(cons dir (shanks a b (* a mem) dir))

  					(> (* a mem) b)
  					;; flip direction
  				  (shanks (/ b mem) a 0 (flip dir))))

We simply use lazy multiplication and comparison inside the algorithm.

To my surprise, this simple adjustment was sufficient to compute logarithms of irrational numbers. Here, for example, the first 23 bits of .

In summary, by using lazy arithmetic within Shanks' logarithm algorithm, the logarithm algorithm itself becomes lazy.

The following table shows 5 values of computed up to 32bit for 5 different values.

The log of 1 is catched beforehand since Shanks' algorithm does not terminate in that case. For and other powers of two, however, it does. The best five-bit approximation of has a dual logarithm very close to .

And here are some values for the natural logarithm. Full tables for numbers up to 5 bit are at the end of the article.

As you can see from the following plot, the logarithms are approximated very accurately. This is despite the fact that is only computed to 16 bits.

Exponentiation

The logarithm can now be used to define exponentiation of arbitrary irrationals in a greedy manner. A stripped down version of a Clojure program to calculate might look like this.

  (defn exp
    [a r b]
    (lazy-seq
      (cond (< (log a r) b)
            (cons R (exp a b (conj r R)))

            (> (log a r) b)
            (cons L (exp a b (conj r L))))))

We start with and calulate . If it's less than , we increase by appending an , if it's greater, we append an . Otherwise we stop, because we've found the solution.

Surprisingly, high precision calculations of exponentials seem to be much cheaper than logarithms, even though exponentiation uses logarithms.

This suggests that not many bits of the logarithms are actually needed to compute the exponential functions, but for now a detailed analysis of the complexity is an open topic, and I hope to come back to it in the future.

References

• [Shanks] Daniel Shanks, A Logarithm Algorithm, Mathematics of Computation 8, no. 46 (1954): 60–64

• [Gosper] Bill Gosper, HAKMEM, 1972, 37-44 and Appendix 2: Continued Fraction Arithmetic

• [Niqui] Milad Niqui, Exact Arithmetic on the Stern–Brocot Tree, Journal of Discrete Algorithms, 2004 Symposium on String Processing and Information Retrieval, 5, no. 2 (June 1, 2007): 356–79

Appendix: Tables

For a better bird's eye view, we replace with ● and with ○.