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