Reliable Computing, scripting language perspective

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.