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.
Sarnold 20041202 (december 2004) Work in progress :
- abandon the integer package because tcllib provides bignum which is far faster
- recenter on big floats and
- move all code into ::math::bigfloat namespace (it is planned that it would be introduced into tcllib, and therefore, in only one namespace
- need some tutorial about doctools (also for moving it in tcllib) (Lars H: How about the doctools page?) Thanks, that is everything I needed. Thank you, and the wiki is great
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 enhancementsAM (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.
Exampleswith 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 [1] 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 [2] (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 ...Lars H: Moreover, this kind of "quad precision floats" does have some advantages over "long doubles". E.g. in a formula such as
1 - sqrt(1-$x)where $x is close to 0, any kind of straightforward float (be that single, double, quadruple, or whatever given precision) will suffer from cancellation which grows worse as the gap in magnitude between $x and 1 gets larger. Not so with pairs of floats--in the extreme case, the hi component of 1-$x simply gets the 1 and the lo component the -$x, so there is always at least simple float precision no matter how large the magnitude gap grows.Of course, at the bottom everything has to be done using floats, so there is nothing magical about it. Implementing sqrt for pairs of floats is pretty much equivalent to rewriting 1 - sqrt(1-$x) on a form that prevents cancellation, and it is prefectly possible to do that with ordinary floats to begin with. The difference is that with pairs of floats arithmetic, you don't have to rewrite your expressions to avoid cancellation, as this will be done at runtime for you.
Sarnold I think in most cases, quadruple-precision numbers should be enough. And I would like to say that the first aim I had when I came to develop MPA was to do decimal arithmetic, and followly I aim to the point that each result should have all his digits exact. But there is a problem again when working with large floats in MPA : You surely may have to write 3.00000000000 when you just mean 3. I realize this is a serious problem, but I do not see the solution. IMHO I think it is the responsability of the user to choose the appropriate library, knowing each library's goods and bads. Please correct me if I am wrong.
