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.

About these ads

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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 645 other followers

%d bloggers like this: