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 (060114) 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