## Calculating Logarithms

### December 18, 2009

Our solution is based on a description by Ed Sandifer. We define `epsilon` as the desired accuracy of the final result:

`(define epsilon 1e-7)`

Newton’s method uses an initial estimate of 1.0, then calculates a new estimate using the formula for the derivative:

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

The implementation of Euler’s method merges two loops into one. The first pesudo-loop, at the `(< hi n)` condition, calculates the smallest power of `base` greater than `n`. The second pseudo-loop performs Euler’s algorithm; `next` is the geometric mean, `log-next` is the arithmetic mean, and the `if` chooses the sub-interval of the bisection:

```(define (euler base n)   (let loop ((lo 1.0) (log-lo 0.0) (hi 1.0) (log-hi 0.0))     (cond ((< hi n) (loop lo log-lo (* hi base) (+ log-hi 1)))           ((< (abs (- (/ log-lo log-hi) 1)) epsilon)             (/ (+ log-lo log-hi) 2))           (else (let* ((next (newton (* lo hi)))                        (log-next (/ (+ log-lo log-hi) 2)))                   (if (< next n)                       (loop next log-next hi log-hi)                       (loop lo log-lo next log-next)))))))```

Here are some sample calculations:

```> (newton 2) 1.4142135623746899 > (euler 10 612) 2.7867514193058014 > (expt 10 2.7867514193058014) 611.9999959982614```

You can run the code at http://programmingpraxis.codepad.org/oOdXhS7g.

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