Calculating Logarithms

December 18, 2009

Before computers, arithmetic was done by hand, and one of the primary tools of the practicing arithmetician was a table of logarithms, which was printed in a large book and used to perform multiplication, division and exponentiation. But how to calculate the logarithms?

Leonhard Euler, the Swiss mathematician, explained the algorithm in his 1748 book Introductio in analysin infinitorum. The basic idea is that, given two numbers a and b, a < n < b, the geometric mean of a and b is equal to the arithmetic mean of their logarithms. Thus, to calculate the logarithm of a given number n, find two powers of the logarithmic base that bound n, calculate the geometric mean (square root of the product) of the two numbers and the arithmetic mean (half of the sum) of the two numbers, and recur on the the smaller interval between one of the two numbers and the geometric mean that bounds n, continuing until the desired accuracy is reached.

But how did Euler calculate the square roots? Sir Isaac Newton, the English physicist, described a method based on derivatives. Given an initial estimate x of the square root of n, a better estimate is given by x - \frac{x^2 - n}{2x}. Again, this calculation is iterated until the desired accuracy is reached.

Your task is to write a function to compute square roots using Newton’s method, then use that function to compute logarithms using Euler’s method. When you are finished, you are welcome to read or run a suggested solution, or to post your solution or discuss the exercise in the comments below.


Pages: 1 2

6 Responses to “Calculating Logarithms”

  1. […] Praxis – Calculating Logarithms By Remco Niemeijer In today’s Programming Praxis problem we have to implement Euler’s method for calculating logarithms […]

  2. Remco Niemeijer said

    My Haskell solution (see for a version with comments):

    eps :: Double
    eps = 1e-7
    sqrt' :: Double -> Double
    sqrt' n = head . filter (\x -> abs (x^2 - n) < eps) $
        iterate (\x -> x - (x^2 - n) / (2*x)) 1
    log' :: Double -> Double -> Double
    log' b n = f 0 0 where
        f lo hi | b ** hi < n             = f lo (hi + 1)
                | abs (lo / hi - 1) < eps = arit
                | geom < n                = f arit hi
                | otherwise               = f lo arit
                where geom = sqrt' (b ** lo * b ** hi)
                      arit = (lo + hi) / 2
  3. Jaime said

    Mine, in Python

    import math
    #Maximum error
    DELTA = 0.0001
    def cmp_estimation(estimation, true):
        if (true - DELTA) <= estimation <= (true + DELTA):
            return True
            return False
    def newton(number):
        Returns the square root using Newton formula
        old_estimation = 0
        estimation = number / 2
        #check until next iteration gives us less precission than DELTA
        while not cmp_estimation(old_estimation, estimation):
            old_estimation = estimation
            estimation = estimation \
                        - (float(estimation ** 2 - number) / (2 * estimation))
        return estimation
    def euler(number, base):
        Calculate logarithms using Euler formula
        #initialize the numbers getting the closest powers
        upper = 1
        while base ** upper <= number:
            lower = upper
            upper += 1
        prev_arith_mean = -1
        arith_mean = 0
        geo_mean = upper
        while not cmp_estimation(arith_mean, prev_arith_mean):
            prev_arith_mean = arith_mean
            #Calculate geometric and arithmetic mean
            arith_mean = float(lower + upper) / 2
            geo_mean = newton((base ** lower) * (base ** upper))
            #Narrow the interval
            if geo_mean < number:
                lower = arith_mean
                upper = arith_mean
        return arith_mean
    #Some test to ensure everything is fine
    if __name__ == '__main__':
        for i in range(2, 100):
            assert cmp_estimation(newton(i), math.sqrt(i))
        for base in range(2, 11):
            for i in range(base, 100):
                assert cmp_estimation(euler(i, base), math.log(i, base))
  4. […] calculated the value of pi, and logarithms to seven digits, in two previous exercises. We continue that thread in the current exercise with a function that calculates […]

  5. Pat said

    I note that the approved solution calculates (- (* x x) x2)) twice. Is there some reason you don’t use a temporary?

    (define (newton x2)
    (let loop ((x 1.0))
    (let ((t (- (* x x) x2)))
    (if (< (abs t) epsilon) x
    (loop (- x (/ t (+ x x))))))))

  6. David said

    In Erlang.

    -export([log/2, sqrt/1]).
    -define(EPS, 1.0e-12).
    sqrt(A) -> sqrt(A, 1).
    sqrt(A, X) when abs(X*X - A) < ?EPS -> X;
    sqrt(A, X) -> sqrt(A, (X + A/X) / 2).
    log(B, A) when A > B -> 1 + log(B, A / B);
    log(B, A) when A == B -> 1;
    log(B, A) -> bisect(A, 0, 1, 1, B).
    bisect(N, LogLo, LogHi, Lo, Hi) ->
        Arith = (LogLo + LogHi) / 2,
        Geo = sqrt(Lo * Hi),
            abs(Lo / Hi - 1) < ?EPS -> Arith;
            N > Geo -> bisect(N, Arith, LogHi, Geo, Hi);
            true    -> bisect(N, LogLo, Arith, Lo, Geo)

    Interactive Erlang session

    11> calc:sqrt(2).
    12> calc:log(10, 612).
    13> math:pow(10, 2.786751422145585).

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: