[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 :(.
----

'''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 CPU

-----
[Sarnold] 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 [http://users.bigpond.net.au/amiller/quad.html] 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. I am investigating that point.


-----
[Category Mathematics] |
[Category Acronym] |
[Category Package]