## 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

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
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))

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.

-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