On precision in Math::Pari

tlm has asked for the wisdom of the Perl Monks concerning the following question:

Dear monks,

This node is somewhere in between a SoPW question and a meditation-like cautionary tale. The questions are at the end.

I just got burned pretty good by the following (keep an eye out for those decimal points):

use strict;
use warnings;
use Math::Pari;

my %zero;
my %zero_dot;
my %delta;

@zero    { qw( perl pari ) } = ( 0  ) x 2;
@zero_dot{ qw( perl pari ) } = ( 0. ) x 2;

$delta{ perl } =       1./7E300  ;
$delta{ pari } = PARI( 1./7E300 );

$zero    { perl } += $delta{ perl };
$zero_dot{ perl } += $delta{ perl };

$zero    { pari } += $delta{ pari };
$zero_dot{ pari } += $delta{ pari };

print '     $delta{ perl }: ', $delta   { perl }, $/;
print '0  + $delta{ perl }: ', $zero    { perl }, $/;
print '0. + $delta{ perl }: ', $zero_dot{ perl }, $/;
print $/;
print '     $delta{ pari }: ', $delta   { pari }, $/;
print '0  + $delta{ pari }: ', $zero    { pari }, $/;
print '0. + $delta{ pari }: ', $zero_dot{ pari }, $/;
__END__
     $delta{ perl }: 1.42857142857143e-301
0  + $delta{ perl }: 1.42857142857143e-301
0. + $delta{ perl }: 1.42857142857143e-301

     $delta{ pari }: 1.428571428571428547E-301
0  + $delta{ pari }: 1.428571428571428547E-301
0. + $delta{ pari }: 0.E-92
[download]

Notice how by adding $delta{ pari } to 0., all the precision is lost. The vague instinct that goes something like this: "by adding a decimal point we somehow indicate that we want floating point operations, therefore we avoid any stray integer division anywhere, therefore we preserve precision as much as possible..." is dead wrong in this case. Unfortunately, my code is rife with numbers written like 0., 1., etc., because, as shown above, Perl handles them in a reasonably predictable way. Switching to Math::Pari has brought on a major headache, due to the need to identify all these innocuous-looking but catastrophic precision leaks.

Where is the behavior above documented? What other dangers await the unsuspecting user of Math::Pari? I looked in perldoc Math::Pari, but no luck.

And more to the...er, point: is there some other way to get the maximum precision one gets by using 0 instead of 0., without having scan through the source for stray decimal points? I tried Math::Pari::setprecision with various arguments, but I got the same results as above every time.

Thanks in advance!

the lowliest monk

Comment on On precision in Math::Pari Download Code

Replies are listed 'Best First'.
Re: On precision in Math::Pari by kvale (Monsignor) on May 05, 2005 at 02:55 UTC
Your vague instinct is correct. From the top of the perdoc for Math::Pari: By default the package exports functions PARI(), PARIcol(), PARIvar(), + PARImat() and PARImat_tr() which convert their argument(s) to a PARI + object. (In fact PARI() is just an alias for new Math::Pari). The fu +nction PARI() accepts following data as its arguments One integer Is converted to a PARI integer. One float Is converted to a PARI float. One string Is executed as a PARI expresion (so should not contain whitespace) +. PARI object Is passed unchanged. Reference to a Perl array Each element is converted using the same rules, PARI vector-row wi +th these elements is returned. ... Conflicts of rules in PARI() In deciding what rule of the above to apply the preference is give +n to the uppermost choice of those available now. If none matches, th +en the string rule is used. So PARI(1) returns integer, PARI(1.) retu +rns float, PARI("1") evaluates 1 as a PARI expression (well, the resu +lt is the same as PARI(1), only slower). [download] So numbers with a decimal point are always interpreted as floats. An alternative method for insuring integer objects is to `int` Perl variables before converting to Pari objects. -Mark	[reply] [d/l] [select]
Re: On precision in Math::Pari by Zaxo (Archbishop) on May 05, 2005 at 02:58 UTC
As you say, integers are treated as exact by pari, while floating-point versions of the values have finite precision. The default floating-point η is larger than the number you're adding on. Libpari and the Math::Pari package don't document this that I know of. Their documentation is generally specific to the functions available and sketchy for elementary properties of the native pari types. In particular, the properties of floating-point representations are assumed to be familiar to the user. They are notorious for their inaccuracies. Floats aren't associative for addition, so they don't even form a mathematical group, much less a (finite) representation of the field of real numbers. After Compline, Zaxo	[reply]