Calculating Logarithms « Programming Praxis

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

Posted by programmingpraxis

Filed in Exercises

4 Comments »

4 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 [...]

Programming Praxis