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.

Posted by programmingpraxis

Filed in Exercises

6 Comments »

6 Responses to “Calculating Logarithms”

Programming Praxis – Calculating Logarithms « Bonsai Code said
December 18, 2009 at 12:15 PM
[…] Praxis – Calculating Logarithms By Remco Niemeijer In today’s Programming Praxis problem we have to implement Euler’s method for calculating logarithms […]

Remco Niemeijer said

December 18, 2009 at 12:16 PM

My Haskell solution (see http://bonsaicode.wordpress.com/2009/12/18/programming-praxis-calculating-logarithms/ 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

Jaime said

December 18, 2009 at 12:31 PM

Mine, in Python

import math

#Maximum error
DELTA = 0.0001

def cmp_estimation(estimation, true):
    if (true - DELTA) <= estimation <= (true + DELTA):
        return True
    else:
        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
        else:
            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))

Calculating Sines « Programming Praxis said
January 12, 2010 at 9:01 AM
[…] 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 […]
Pat said
May 9, 2011 at 11:51 PM
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))))))))

David said

September 30, 2014 at 2:24 AM

In Erlang.

-module(calc).
-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),
    if
        abs(Lo / Hi - 1) < ?EPS -> Arith;
        N > Geo -> bisect(N, Arith, LogHi, Geo, Hi);
        true    -> bisect(N, LogLo, Arith, Lo, Geo)
    end.

Interactive Erlang session

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

Programming Praxis