logarithm

RS: Tcl's expr has logarithm to base e built-in:

% expr exp(1) ;# Euler's number, or e for short
2.71828182846
8 % expr log(exp(1))
1.0

Logarithms to most any other base can be had by division -- e.g for log10 of 1000:

% expr log(1000)/log(10)
3.0

Feynman apparently is in the domain literature for numerical analysis of an algorithm which computes the logarithm of a number between 1 and 2. Any such number is a unique product of factors of the form (1 + 2 ** -k) for integral k. Logarithms of all such factors fit in a table of modest size. Determination of the sequence of k-s can be done with shifts and subtracts.

Here's a Tcl coding that illustrates the algorithm:

...

See math::bigfloat::log for another algorithm. (with arbitrary precision) It uses this development:

log(1+x)=x + x^2/2 + x^3/3 ...

and does not require sqrt() to be present. It only needs arithmetics.

proc logarithm x {
    set eps 1e-[set ::tcl_precision]
    if {$x<=-1} {error {logarithm requires a positive argument}}
    set log $x
    set pow [expr {-$x*$x}]
    set i 2
    set add [expr {$pow/$i}]
    while {abs($add)>$eps*abs($log)} {
        set log [expr {$log+$add}]
        incr i
        set pow [expr {-$pow*$x}]
        set add [expr {$pow/$i}]
    }
    incr i -2
    #puts "$i iterations"
    return $log
}


proc log x {
    set doubles 0
    if {$x<=0} {error {logarithm requires a positive argument}}
    if {$x<1} {
        while {$x<=0.5} {
            # iterate : x=x*2 until 1/2<x<1
            incr doubles -1
            set x [expr {$x*2}]
        }
    } else  {
        while {$x>=2} {
            # iterate : x=x/2 until 1<x<2
            incr doubles
            set x [expr {$x/2}]
        }
    }
    set x [expr {double($x-1)}]
    global log2
    #puts x=$x
    return [expr {[logarithm $x]+$doubles*$log2}]
}
# compute log(2)
set log2 [expr {-[logarithm -0.2]*3-[logarithm [expr {-3./128.}]]}]


# for debugging purpose
proc close x {
    set eps 1e-[set ::tcl_precision]
    set my [log $x]
    set system [expr {log($x)}]
    if {abs($my-$system)>$eps*abs($system)} {
        puts "my log $x : $my"
        puts "system log($x) : $system"
        error "$x runs odd"
    }
    return $my
}

-- Sarnold

GS 2016-01-14: Until we have the solution for above, here is an algorithm which computes the logarithm of a number between 0 and 1.

# Description: Log evaluation between 0 and 1 with a recursion formula
#
#                                 /    2*u(i)     \
#              u(i+1) = u(i)*sqrt|-----------------|
#                                 \ u(i) + u(i-1) /
# Arguments: 
# . x : the value to be computed
# Results: 
# . y: the value of log(x)
# . n: then number of iterations

proc Log01 x {
    if {$x <= 0 || $x > 1} return 
    set eps .000000001
    set n 0
    set y 0
    set du 10000
    set x2 [expr {$x*$x}]
    set u0 [expr {($x2 - 1.0/$x2)/4}] 
    set u1 [expr {($x-1.0/$x)/2}]
   
    while {$du > $eps} {
         set u2 [expr {$u1*sqrt(2*$u1/($u1+$u0))}]
         set y $u2
         incr n
         set du [expr {abs($u2-$u1)}]
         set u0 $u1
         set u1 $u2 
    }
    return [list $y $n]
}

# Test
for {set i 2} {$i < 10} {incr i} {
   set x [expr {1.0/$i}]
   set r [Log01 $x]
   puts "Log($x) = [format "%.8f" [lindex $r 0]] in [lindex $r 1] iterations"
}

Test results:

Log(0.5) = -0.69314718 in 14 iterations
Log(0.333333333333) = -1.09861229 in 15 iterations
Log(0.25) = -1.38629436 in 16 iterations
Log(0.2) = -1.60943791 in 16 iterations
Log(0.166666666667) = -1.79175947 in 16 iterations
Log(0.142857142857) = -1.94591015 in 16 iterations
Log(0.125) = -2.07944154 in 17 iterations
Log(0.111111111111) = -2.19722458 in 17 iterations

Category Mathematics