#
 # ----------------------------------------------------------------------
 # 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]
 }