Sarnold Acronym for Multiple Precision Arithmetics. Library in pure Tcl available at : http://sarnold.free.fr/
Released with a BSDish License. The fact is that the floating-point numbers are still experimental, even if it 'works'. (that means, the precision of fp. numbers is a subject i personnaly do not master) I guess it is not obvious to do such computations without losing precision at each step.
What's new : in tcllib there is a bignum implementation that is much faster than mine. For your information it uses lists as an internal representation of bignums.
Version 1.2 dated 12th November 2004 with many performance enhancements
AM (12 november 2004) Do you care to contribute this to Tcllib?
Sarnold Yes, but slowly, as I am not very active recently. I aim to do that, but doing all the stuff with selftests and manual pages, and more, treating with CVS on a 56K modem w/ Win ME is horrible :(.
AM Good - if you need help, just let me know.
Examples
with a 333MHz PII 100MHz PCI Bus PC with Windows Me.
(mpa-v0.11) 50 % package require mpa
0.11
(mpa-v0.11) 51 % mpa::int::add 111111111111111111111 98237827463784678346478
98348938574895789457589
(mpa-v0.11) 52 % mpa::int::add 1000000000000000000000 99999999999999999
1000099999999999999999
(mpa-v0.11) 53 % mpa::int::mul 1001 9009
9018009
(mpa-v0.11) 54 % mpa::int::mul 10010001 9009009
90180189099009
(mpa-v0.11) 55 % mpa::float::pi 20
# Pi constant with 20 decimals, the result is cached in a namespace variable
3.14159265358979323846
(mpa-v0.11) 56 % time {mpa::float::pi 20};# result is cached in memory
363 microseconds per iteration
(mpa-v0.11) 57 % time {mpa::float::pi 21};# result is now recomputed
90643 microseconds per iteration
(mpa-v0.11) 58 % time {puts [mpa::float::pi 21]};# result is cached
3.141592653589793238462
4962 microseconds per iteration
(mpa-v0.11) 59 % mpa::float::add 1.0 2.0
3.0~2 # the precision with 1.0 ans 2.0 is assumed to be 0.1
# and the precision of the sum is the sum of the precisions
(mpa-v0.11) 60 % mpa::float::add 1.0000 2.0000
3.0000~2
(mpa-v0.11) 62 % mpa::float::format [mpa::float::mul 1.0000000 2.0000000];
# the precision of the product is approximatly A*Pb+B*Pa = 3e-7
2.000000 # now the number is formatted to be put on screen
(mpa-v0.11) 63 % mpa::float::mul 1.0000000 2.0000000
2.0000000~4 # precision is computed to be always the smallest ,
# even if such computations are complex to handle and loss of CPUSarnold This package is especially good when you want it to represent arbitrary floating-point numbers. See also : bignum in pure tcl.
AM Yesterday I was reminded of a Fortran 90 library that deals with quadruple-precision variables written by Alan Miller [L1 ] and I decided to try how this would work in Tcl.
The motivation: when double precision is not enough, quadruple precision may be a useful alternative to a full MPA package. Such packages are inherently slow - they simply have to do a lot of work.
I am not sure that all is working correctly - there are a few tricky bits and the results are not completely clear ... that is: as there is no proper string conversion as yet, I can not quite see whether the last test case comes up with the proper results.
# quadprec.tcl --
# Experiments with quadruple-precision in Tcl
# The code below consists of a transcription of Fortran 90
# code by Alan Miller (<http://.../quad.html>)
#
# A quadruple-precision number is implemented as a list of two
# elements, the higher-order and the lower-order parts of the
# full number
#
# Note:
# The parentheses are essential - they should cause the
# finite-precision arithmetic to generate proper intermediate
# results.
#
namespace eval ::Quadprecision {
namespace export qadd qmult
#
# This constant may be machine-dependent: it reflects the
# number of binary digits that the intermediate results
# have!
# The term 11 reflects the extended precision on Intel
# machines ...
#
variable nbits 53
#variable constant [expr {pow(2.0,($nbits+11-$nbits/2)) + 1.0}]
#variable constant [expr {pow(2.0,($nbits-$nbits/2)) + 1.0}]
variable constant 1.0
for { set i 0 } { $i < $nbits-$nbits/2 } { incr i } {
set constant [expr {$constant * 2.0}]
}
set constant [expr {$constant+1.0}]
}
# qadd --
# Quadruple-precision procedure for adding two numbers
# Arguments:
# qp1 First quadruple-precision operand
# qp2 Second quadruple-precision operand
# Result:
# List of two numbers, representing the quadruple-precision sum
#
proc ::Quadprecision::qadd { qp1 qp2 } {
foreach {qp1_hi qp1_lo} $qp1 {break}
foreach {qp2_hi qp2_lo} $qp2 {break}
set z [expr {$qp1_hi+$qp2_hi}]
set q [expr {$qp1_hi-$z}]
set zz [expr { ((($q+$qp2_hi) + ($qp1_hi-($q+$z))) + $qp1_lo ) \
+ $qp2_lo}]
set r_hi [expr {$z + $zz}]
set r_lo [expr {($z - $r_hi) + $zz}]
return [list $r_hi $r_lo]
}
# exactmul2 --
# Quadruple-precision procedure for multiplying two
# double-precision numbers exactly
# Arguments:
# dp1 First double-precision operand
# dp2 Second double-precision operand
# Result:
# List of two numbers, representing the quadruple-precision product
#
proc ::Quadprecision::Exactmul2 {dp1 dp2 } {
variable constant
set t [expr {$constant*$dp1}]
set a1 [expr {($dp1-$t) + $t}]
set a2 [expr {$dp1-$a1}]
set t [expr {$constant*$dp2}]
set c1 [expr {($dp2-$t) + $t}]
set c2 [expr {$dp2-$c1}]
set r_hi [expr {$dp1 * $dp2}]
set r_lo [expr {((($a1 * $c1 - $r_hi) + $a1 * $c2) + $c1 * $a2) \
+ $c2 *$a2}]
return [list $r_hi $r_lo]
}
# qmult --
# Quadruple-precision procedure for multiplying two numbers
# Arguments:
# qp1 First quadruple-precision operand
# qp2 Second quadruple-precision operand
# Result:
# List of two numbers, representing the quadruple-precision product
#
proc ::Quadprecision::qmult { qp1 qp2 } {
foreach {qp1_hi qp1_lo} $qp1 {break}
foreach {qp2_hi qp2_lo} $qp2 {break}
foreach {z_hi z_lo} [Exactmul2 $qp1_hi $qp2_hi] {break}
set zz [expr {(($qp1_hi+$qp1_lo) * $qp2_lo + $qp1_lo * $qp2_hi) + $z_lo }]
set r_hi [expr {$z_hi + $zz}]
set r_lo [expr {($z_hi-$r_hi) + $zz}]
return [list $r_hi $r_lo]
}
#
# A few simple tests
#
set tcl_precision 17
set qp1 [list 1.0 1.0e-20]
set qp2 [list 1.0 2.0e-21]
set qp3 [list 1.00000000001 0.0]
set qp4 [list 1.00000000001 0.0]
puts "Sum of '$qp1' and '$qp2':"
puts [::Quadprecision::qadd $qp1 $qp2]
puts "Product of '$qp3' and '$qp4':"
puts [::Quadprecision::qmult $qp3 $qp4]
set qp1 [list 1.0 -1.0e-20]
set qp2 [list 1.0 1.0e-20]
puts "Sum of '$qp1' and '$qp2':"
puts [::Quadprecision::qadd $qp1 $qp2]
set qp1 [list 1.00000000001 -1.0e-20]
set qp2 [list 1.00000000001 1.0e-20]
puts "Sum of '$qp1' and '$qp2':"
puts [::Quadprecision::qadd $qp1 $qp2]
# 12345678901
set qp3 [list 10.0000000001 0.0e-20]
set qp4 [list 0.10000000001 0.0e-20]
puts "Product of '$qp3' and '$qp4':"
puts [set qp5 [::Quadprecision::qmult $qp3 $qp4]]
puts [::Quadprecision::qadd $qp5 {-1.0 0.0}]
puts "Just a crude comparison: [expr {10.0001*0.10001}]"Sarnold IMHO when 64 bit processors will generalize, they all will provide "long double" like GCC does for some platforms. Take a look at [L2 ] (french site).
AM What I gather from postings in comp.lang.fortran is that many 64-bits platforms or compilers for these platforms concentrate on the use of many gigabytes of memory, rather than on quadruple-precision numbers ...