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