Big Floating-Point Numbers

Tcl, with its dependency on the C programming language, usually uses IEEE floating numbers, which have a very specific value range. IEEE double-precision values have an 11 bit exponent, allowing for exponents between -1024 and 1023, base 2. Therefore, the minimum positive value larger than zero is 2^-1024 or approximately 5*10^-309, and the maximum value is about 2*10^308.

If you feel constrained by that narrow value range, there are some options. First, there is bignum, which is a compiled extension. Second, there is BigFloat for Tcl, which is slow (e.g., multiplication is O(n^2), with n being the logarithm of the larger operand). Third, I give you Big Floating-Point Numbers, or bigfp.

This package is complete (includes log, exp, and pow), fast (all operations, but the factorial, are constant time), and small (see below). It's value range tops out at around 10^(2^63) or just shy of 10^(10^19).

Have you ever wondered what (10000!)^42 is? Trivial:

 % bigfp::pow [bigfp::fac 10000] 42
 1.2009636788737112e1497697

The catch, you ask? No catch. It's just that bigfp neither promises, nor delivers, arbitrary precision. While having a vastly increased value range, bigfp is no more precise than Tcl's own floating point numbers.

This package generalizes an old idea which I've used for years to estimate large factorial values.

The code below introduces the bigfp namespace, with the following operations:

  • bigfp::+ a b - Returns the sum of a and b.
  • bigfp::- a [b] - With one parameter, returns -a. With two parameters, substracts b from a.
  • bigfp::* a b - Multiplies the two numbers.
  • bigfp::/ a b - Divides a by b.
  • bigfp::log a - Returns the natural logarithm of a.
  • bigfp::exp a - Returns e to the power of a.
  • bigfp::pow a b - Computes the value of a raised to the power of b. a must be positive.
  • bigfp::fac n - Computes the factorial of n, which must be a positive integer.

All parameters (but for the parameter to the factorial operator) may be floating-point values that may well exceed Tcl's usual floating-point range, e.g., "1.23e9876".

Have fun! FPX


 #
 # ----------------------------------------------------------------------
 # Big Floating-Point Numbers. By Frank Pilhofer.
 # ----------------------------------------------------------------------
 #

 namespace eval bigfp {}

 #
 # ----------------------------------------------------------------------
 # Internal interface
 # ----------------------------------------------------------------------
 #

 #
 # Maximum exponent for which to print "normal" numbers, instead of
 # using exponential notation. E.g, when set to 3, "1234" will show
 # up as "1234", but "12345" will show up as "1.2345e4".
 #

 namespace eval bigfp::int {
    variable maxnormexp 3
 }

 #
 # Read a number from a string into its internal representation.
 #

 proc bigfp::int::in {on} {
    set n [string trim $on]

    if {$n == ""} {
        error "not a number: \"$on\""
    }

    if {[set ei [string first "e" $n]] == -1} {
        set ei [string first "E" $n]
    }

    if {$ei != -1} {
        set m [string range $n 0 [expr {$ei-1}]]
        set e [string range $n [expr {$ei+1}] end]
    } else {
        set m $n
        set e 0
    }

    if {[catch {set m [expr {double($m)}]}] || \
            [catch {set e [expr {wide($e)}]}]} {
        error "not a number: \"$on\""
    }

    return [norm [list $m $e]]
 }

 #
 # Print a number.
 #

 proc bigfp::int::out {n} {
    variable maxnormexp
    foreach {m e} $n {}
    if {$e >= -$maxnormexp && $e <= $maxnormexp} {
        set m [expr {$m*pow(10.,$e)}]
        set e 0
    }
    if {$e == 0} {
        set r $m
    } else {
        set r "${m}e$e"
    }
    return $r
 }

 #
 # Normalize a number.
 #

 proc bigfp::int::norm {n} {
    foreach {m e} $n {}
    if {$m < 0} {
        set s -1
        set m [expr {-$m}]
    } else {
        set s 1
    }
    if {[catch {
        set lm [expr {log10($m)}]
    }]} {
        return [list 0 0]
    }
    if {$lm < 0 || $lm >= 1} {
        set de [expr {wide(floor($lm))}]
        set lm [expr {$lm-$de}]
        incr e $de
        set m [expr {$s*pow(10.,$lm)}]
    } else {
        set m [expr {$s*$m}]
    }
    return [list $m $e]
 }

 #
 # Addition
 #

 proc bigfp::int::add {a b} {
    foreach {am ae} $a {}
    foreach {bm be} $b {}

    set ed [expr {$ae-$be}]

    if {$ed > 0} {
        set e $ae
        if {[catch {set bm [expr {$bm/pow(10.,$ed)}]}]} {
            set bm 0
        }
    } elseif {$ed < 0} {
        set e $be
        if {[catch {set am [expr {$am/pow(10.,-$ed)}]}]} {
            set am 0
        }
    } else {
        set e $ae
    }

    set m [expr {$am+$bm}]
    return [norm [list $m $e]]
 }

 #
 # Subtraction
 #

 proc bigfp::int::sub {a b} {
    foreach {bm be} $b {}
    return [add $a [list [expr {-$bm}] $be]]
 }

 #
 # Multiplication
 #

 proc bigfp::int::mul {a b} {
    foreach {am ae} $a {}
    foreach {bm be} $b {}
    set m [expr {$am*$bm}]
    set e [expr {$ae+$be}]
    return [norm [list $m $e]]
 }

 #
 # Division
 #

 proc bigfp::int::div {a b} {
    foreach {am ae} $a {}
    foreach {bm be} $b {}
    if {$bm == 0} {
        error "division by zero"
    }
    set m [expr {$am/$bm}]
    set e [expr {$ae-$be}]
    return [norm [list $m $e]]
 }

 #
 # Natural logarithm
 #

 proc bigfp::int::log {n} {
    foreach {m e} $n {}
    if {$m <= 0} {
        error "invalid argument"
    }
    set rm [expr {log($m)+2.3025850929940459*$e}]
    return [norm [list $rm 0]]
 }

 #
 # Exponentiation, base e
 #

 proc bigfp::int::exp {n} {
    foreach {m e} $n {}
    set lm [expr {$m*pow(10.,$e)*0.43429448190325182}]
    set re [expr {wide(floor($lm))}]
    set lm [expr {$lm-$re}]
    set rm [expr {pow(10.,$lm)}]
    return [list $rm $re]
 }

 #
 # Exponentiation, arbitrary base
 #

 proc bigfp::int::pow {a b} {
    foreach {am ae} $a {}
    foreach {bm be} $b {}
    if {$am < 0} {
        error "invalid argument"
    }
    if {[catch {set lm [expr {$bm*($ae+log10($am))*pow(10.,$be)}]}]} {
        return [list 0 0]
    }
    set re [expr {wide(floor($lm))}]
    set lm [expr {$lm-$re}]
    set rm [expr {pow(10.,$lm)}]
    return [list $rm $re]
 }

 #
 # ----------------------------------------------------------------------
 # Public interface
 # ----------------------------------------------------------------------
 #
 # These functions use regular floating-point numbers as input parameters
 # and result. The numbers may well exceed Tcl's numeric range, such as
 # 1e2345.
 #

 proc bigfp::+ {a b} {
    return [bigfp::int::out \
                [bigfp::int::add [bigfp::int::in $a] \
                     [bigfp::int::in $b]]]
 }

 proc bigfp::- {a {b ""}} {
    if {$b == ""} {
        return [bigfp::int::out \
                    [bigfp::int::sub [list 0 0] \
                         [bigfp::int::in $a]]]
    }
    return [bigfp::int::out \
                [bigfp::int::sub [bigfp::int::in $a] \
                     [bigfp::int::in $b]]]
 }

 proc bigfp::* {a b} {
    return [bigfp::int::out \
                [bigfp::int::mul [bigfp::int::in $a] \
                     [bigfp::int::in $b]]]
 }

 proc bigfp::/ {a b} {
    return [bigfp::int::out \
                [bigfp::int::div [bigfp::int::in $a] \
                     [bigfp::int::in $b]]]
 }

 proc bigfp::pow {a b} {
    return [bigfp::int::out \
                [bigfp::int::pow [bigfp::int::in $a] \
                     [bigfp::int::in $b]]]
 }

 proc bigfp::log {a} {
    return [bigfp::int::out [bigfp::int::log [bigfp::int::in $a]]]
 }

 proc bigfp::exp {a} {
    return [bigfp::int::out [bigfp::int::exp [bigfp::int::in $a]]]
 }

 #
 # Factorial
 #

 proc bigfp::fac {n} {
    if {[catch {set n [expr {int($n)}]}] || $n < 0} {
        error "invalid argument"
    }
    set r [list 1 0]
    for {set i 2} {$i<=$n} {incr i} {
        set r [bigfp::int::mul $r [list $i 0]]
    }
    return [bigfp::int::out $r]
 }