I am happy that I attended the
Sun HPC Consortium Singapore 2006.
This was held in conjunction with GridAsia 2006.

One of the presentations (Reliable Computing and Petascale systems, by Ruud van der Pas, Sun Microsystems)
interested me to explore further, from the viewpoint of scripting languages. Although his
presentation material is not available, I managed to find a similar paper that he presented last year in
Australia.

The Future Of Floating-Point Computing Is In Reliable Computing, in Australian Partnership for Advanced Computing Conference and Exhibition on Advanced Computing, Grid Applications and eResearch, 26-30 September 2005.

The equation in question is:

F(a,b) = (333.75-a^2)*b^6 + a^2*(11*a^2*b^2-121*b^4-2) + 5.5*b^8 + a/(2*b)

a = 77617.0; b = 33096.0

The answer (using 32-bit precision): 1.172604. Question: How Good Is This Result?

In thr APAC 2005 paper, it said:

`
"If we bump up the precision and the
numbers match, we're okay"
Not Necessarily .......
`

% ./rump.exe

Computed Results

32-bit precision: 1.172604

64-bit precision: 1.1726039400531786

128-bit precision: 1.1726039400531786949244406059733592

% ./rump.exe

Computed Results

32-bit precision: 1.172604

64-bit precision: 1.1726039400531786

128-bit precision: 1.1726039400531786949244406059733592

Correct Results

32-bit precision: -0.82739603

64-bit precision: -0.8273960599468213

128-bit precision: -0.8273960599468213050755593940266408

Okay, let's look at how scripting languages handle the above equation on my
notebook running Cygwin:

**Perl** (5.8.7)

[ans=1.17260394005318] NOT CORRECT

$ cat rump.pl
#! /usr/bin/perl
$a=77617.0;
$b=33096.0;
$a2=$a*$a;
$b2=$b*$b;
$b4=$b2*$b2;
$b6=$b4*$b2;
$b8=$b6*$b2;
$ans=(333.75-$a2)*$b6 + $a2*(11*$a2*$b2-121*$b4-2) + 5.5*$b8 + $a/(2*$b);
print $ans;
$ perl rump.pl
1.17260394005318

**Tcl** (8.4.1)

[ans=1.17260394005] NOT CORRECT

$ cat rump.tcl
#! /usr/bin/tclsh
proc f { a b } {
set a2 [expr {$a*$a}]
set b2 [expr {$b*$b}]
set b4 [expr {$b2*$b2}]
set b6 [expr {$b4*$b2}]
set b8 [expr {$b6*$b2}]
set ans [expr {(333.75-$a2)*$b6 + $a2*(11*$a2*$b2-121*$b4-2) + 5.5*$b8 + $a/(2*$b)}]
return $ans
}
puts [f 77617.0 33096.0]
$ ./rump.tcl
1.17260394005

Apparently all these calculations are based on IEEE floating point standard.
If we were to calculate it long hand, it looks like this:

(333.75 - a^2)*b^6=
-7917110903377385049079188237280149504.00
a^2*(11a^2*b^2-121*b^4-2)=
-437291576312021946464244793346.00
5.5*b^8=
7917111340668961361101134701524942848.00
a/(2*b)=
1.17260394005317863185
-------------------------------------------------------------------------------
F(a,b)=
-0.8273960599468213681411650954798162919991
-7917110903377385049079188237280149504.00 +
-437291576312021946464244793346.00 +
7917111340668961361101134701524942848.00 = 2

Can we do it correctly with scripting languages (if we consider "bc" also) ? Let see

**bc** (1.06) [ANS=-.8273960599468213681411650954798162919991] CORRECT

$ cat rump.bc
scale=40
a=77617.0
b=33096.0
(333.75-a^2)*b^6+a^2*(11*a^2*b^2-121*b^4-2)+5.5*b^8+a/(2*b)
$ bc -l < rump.bc
-.8273960599468213681411650954798162919991

**Perl** (5.8.7) (same script as above, but run with Bignum module)

[ans=-0.827396059946821368141165095479816291999]

CORRECT
$ cat rump.pl
#! /usr/bin/perl
$a=77617.0;
$b=33096.0;
$a2=$a*$a;
$b2=$b*$b;
$b4=$b2*$b2;
$b6=$b4*$b2;
$b8=$b6*$b2;
$ans=(333.75-$a2)*$b6 + $a2*(11*$a2*$b2-121*$b4-2) + 5.5*$b8 + $a/(2*$b);
print $ans;
$ perl -Mbignum rump.pl
-0.827396059946821368141165095479816291999

**Tcl (math::bigfloat)**
[ans=-0.827396059946821368141165095479816291999]

CORRECT
The tricky part is how much precision one has to set to get the correct answer. Also, the implementation is very clumsy.

package require math::bigfloat
namespace import math::bigfloat::*
set precision 77
set a [fromdouble 77617.0 $precision]
set b [fromdouble 33096.0 $precision]
set c333_75 [fromdouble 333.75 $precision]
set c11 [fromdouble 11.0 $precision]
set c121 [fromdouble 121.0 $precision]
set c2 [fromdouble 2.0 $precision]
set c5_5 [fromdouble 5.5 $precision]
set a2 [mul $a $a]
set b2 [mul $b $b]
set b4 [mul $b2 $b2]
set b6 [mul $b4 $b2]
set b8 [mul $b6 $b2]
set term1 [mul [sub $c333_75 $a2] $b6]
set term2 [mul $a2 [sub [sub [mul [mul $c11 $a2] $b2] [mul $c121 $b4]] $c2]]
set term3 [mul $c5_5 $b8]
set term4 [div $a [mul $c2 $b]]
set ans [add [add [add $term1 $term2] $term3] $term4]
puts [tostr $ans]

**Tcl (Mpexpr)**
[ans=-0.827396059946821368141165095479816291999]

CORRECT
Tcl module (

Mpexr )
is able to do multiple precision arithmatic.
Below is Tcl 8.1.1 with Mpexr on Solaris, because I have problem compiling that with the latest Tcl8.4 on Cygwin

$ cat rump.tcl
package require Mpexpr
proc f { a b } {
set a2 [mpexpr {$a*$a}]
set b2 [mpexpr {$b*$b}]
set b4 [mpexpr {$b2*$b2}]
set b6 [mpexpr {$b4*$b2}]
set b8 [mpexpr {$b6*$b2}]
set ans [mpexpr {(333.75-$a2)*$b6 + $a2*(11*$a2*$b2-121*$b4-2) + 5.5*$b8 + $a/(2*$b)}]
return $ans
}
set mp_precision 40
puts [f 77617.0 33096.0]
$ tclsh rump.tcl
-0.827396059946821368141165095479816291999

## Conclusion

Scripting languages can be deployed in the area of reliable computing, but one has to be mindful in choosing the right modules for the job. BTW, I am very impressed by Perl's Bignum implementation which hides all the technicality from the user. The same perl script can run with or without Bignum module.

Labels: awk, bc, Perl, Tcl, unix