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.

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