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 
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.
Programming Praxis 
[…] Praxis – Calculating Logarithms By Remco Niemeijer In today’s Programming Praxis problem we have to implement Euler’s method for calculating logarithms […]
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) / 2Mine, 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))[…] 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 […]
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))))))))
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