## #166 rint, ceil, floor sometimes return the wrong value?

closed-works-for-me
None
5
2007-03-19
2007-03-19
John Harker
No

A very subtle bug can be introduced into a program if the writer is unaware that rint, ceil, and floor (and possibly others) do NOT return their result, but ONLY work inplace. Contrary to intuition, rint etc. actually return the INPUT value.

Here's the essence of the problem. This program:

===========================================
use PDL;
use PDL::Math;

\$Numerator = 0.25;

\$Min = 1 / 800;
\$Max = 1 / 400;
\$Step = ( \$Max - \$Min ) / ( 19 );

\$Denominator = rint( \$Numerator / \$Step );
print "Denominator = ".\$Denominator."\n";

\$Ratio1 = \$Numerator / \$Denominator;
\$Ratio2 = \$Numerator / 3800;

print "Ratio1 = ".\$Ratio1."\n";
print "Ratio2 = ".\$Ratio2."\n";

\$Divisor1 = \$Numerator / \$Ratio1;
\$Divisor2 = \$Numerator / \$Ratio2;

print "Divisor1 = ".\$Divisor1."\n";
print "Divisor2 = ".\$Divisor2."\n";

# die;

...will return this output...

===========================================
Denominator = 3800
Ratio1 = 6.57894706819206e-05
Ratio2 = 6.57894736842105e-05
Divisor1 = 3800.00024414062
Divisor2 = 3800
Died at test-02.pl line 25.
===========================================

You can see that the values for \$Ratio1 and \$Ratio2 are subtly different, and while \$Divisor1 might be reasonably expected to be an integer, it is not. Further, why does \$Denominator print as an integer when \$Divisor1 (which should be the same value) does not?

It is not obvious to me why rint etc. should return the input value. They should return the output value, or nothing. This oddity provides huge possibilities for "silent" errors that could potentially corrupt the results of any number of calculations.

Thanks!

Here's the full output of "perldl -V":

===========================================
perlDL shell v1.33
PDL comes with ABSOLUTELY NO WARRANTY. For details, see the file
'COPYING' in the PDL distribution. This is free software and you
are welcome to redistribute it under certain conditions, see
the same file for details.

Summary of my PDL configuration

VERSION: PDL v2.4.3

\$%PDL::Config = {
'WITH_PROJ' => undef,
'FFTW_TYPE' => 'double',
'FFTW_LIBS' => [
'/lib',
'/usr/lib',
'/usr/local/lib'
],
'WITH_FFTW' => undef,
'GSL_LIBS' => undef,
'WITH_IO_BROWSER' => '0',
'PROJ_INC' => undef,
'WHERE_PLPLOT_INCLUDE' => undef,
'WITH_KARMA' => '0',
'WHERE_KARMA' => undef,
'HTML_DOCS' => '1',
'WHERE_PLPLOT_LIBS' => undef,
'WITH_3D' => '0',
'FFTW_INC' => [
'/usr/include/',
'/usr/local/include'
],
'WITH_HDF' => undef,
'HDF_INC' => undef,
'OPENGL_LIBS' => undef,
'WITH_GD' => undef,
'FITS_LEGACY' => '1',
'WITH_SLATEC' => '1',
'TEMPDIR' => '/tmp',
'PROJ_LIBS' => undef,
'GD_LIBS' => undef,
'GSL_INC' => undef,
'GD_INC' => undef,
'OPTIMIZE' => undef,
'WITH_GSL' => undef,
'HDF_LIBS' => undef,
'MALLOCDBG' => {},
'WITH_PLPLOT' => '0'
};
Summary of my perl5 (revision 5 version 8 subversion 8) configuration:
Platform:
uname='linux rothera 2.6.15.7 #1 smp tue jun 27 18:34:43 utc 2006 i686 gnulinux '
config_args='-Dusethreads -Duselargefiles -Dccflags=-DDEBIAN -Dcccdlflags=-fPIC -Darchname=i486-linux-gnu -Dprefix=/usr -Dprivlib=/usr/share/perl/5.8 -Darchlib=/usr/lib/perl/5.8 -Dvendorprefix=/usr -Dvendorlib=/usr/share/perl5 -Dvendorarch=/usr/lib/perl5 -Dsiteprefix=/usr/local -Dsitelib=/usr/local/share/perl/5.8.8 -Dsitearch=/usr/local/lib/perl/5.8.8 -Dman1dir=/usr/share/man/man1 -Dman3dir=/usr/share/man/man3 -Dsiteman1dir=/usr/local/man/man1 -Dsiteman3dir=/usr/local/man/man3 -Dman1ext=1 -Dman3ext=3perl -Dpager=/usr/bin/sensible-pager -Uafs -Ud_csh -Uusesfio -Uusenm -Duseshrplib -Dlibperl=libperl.so.5.8.8 -Dd_dosuid -des'
hint=recommended, useposix=true, d_sigaction=define
useperlio=define d_sfio=undef uselargefiles=define usesocks=undef
use64bitint=undef use64bitall=undef uselongdouble=undef
usemymalloc=n, bincompat5005=undef
Compiler:
cc='cc', ccflags ='-D_REENTRANT -D_GNU_SOURCE -DTHREADS_HAVE_PIDS -DDEBIAN -fno-strict-aliasing -pipe -I/usr/local/include -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64',
optimize='-O2',
cppflags='-D_REENTRANT -D_GNU_SOURCE -DTHREADS_HAVE_PIDS -DDEBIAN -fno-strict-aliasing -pipe -I/usr/local/include'
ccversion='', gccversion='4.1.2 20060613 (prerelease) (Ubuntu 4.1.1-2ubuntu5)', gccosandvers=''
intsize=4, longsize=4, ptrsize=4, doublesize=8, byteorder=1234
d_longlong=define, longlongsize=8, d_longdbl=define, longdblsize=12
ivtype='long', ivsize=4, nvtype='double', nvsize=8, Off_t='off_t', lseeksize=8
alignbytes=4, prototype=define
ld='cc', ldflags =' -L/usr/local/lib'
libpth=/usr/local/lib /lib /usr/lib
libs=-lgdbm -lgdbm_compat -ldb -ldl -lm -lpthread -lc -lcrypt
libc=/lib/libc-2.4.so, so=so, useshrplib=true, libperl=libperl.so.5.8.8
gnulibc_version='2.4'
dlsrc=dl_dlopen.xs, dlext=so, d_dlsymun=undef, ccdlflags='-Wl,-E'
cccdlflags='-fPIC', lddlflags='-shared -L/usr/local/lib'
===========================================

## Discussion

• John Harker - 2007-03-19
• summary: rint, ceil, floor should return their expected result --> rint, ceil, floor sometimes return the wrong value?

• John Harker - 2007-03-19

Logged In: YES
user_id=1747641
Originator: YES

Slight update to my original bug post. I was confused- after experimenting I found that in MOST cases, rint etc. actually DO return their output value. In fact, the problem does not appear to be related to inplace vs. returned values at all.

# Instead, for certain special cases, rint etc. somehow seem to produce the wrong value. Here's an example:

use PDL;
use PDL::Math;

\$Numerator = 0.25;

#\$Step = (1/400 - 1/900)/19;
\$Step = (1/400 - 1/800)/19;

\$Denominator = rint( \$Numerator / \$Step );
rint( \$Denominator );

print "Denominator = ".\$Denominator."\n";

\$Ratio = \$Numerator / \$Denominator;

print "Ratio = ".\$Ratio."\n";

\$Divisor = \$Numerator / \$Ratio;

print "Divisor = ".\$Divisor."\n";

# die;

This program would produce this output:

========================================
Denominator = 3800
Ratio = 6.57894706819206e-05
Divisor = 3800.00024414062
Died at test-03.pl line 22.
========================================

which has an error. Switching the commented-out line would produce this output:

========================================
Denominator = 3420
Ratio = 7.30994142941199e-05
Divisor = 3420
Died at test-03.pl line 22.
========================================

which is correct.

• Craig DeForest - 2007-03-19

Logged In: YES
user_id=20200
Originator: NO

This is not a problem with rint, ceil, floor -- it is a problem with mixing perl scalars with PDL variables. Instead try the following:

use PDL;
use PDL::Math;
\$Num = pdl(0.25);
\$Stp = pdl(1/400-1/800)/19;
\$Den = rint(\$Num/\$Stp);
print "Den=\$Den\n";
\$Rat = \$Num/\$Den;
print "Rat=\$Rat\n";
\$Div=\$Num/\$Rat;
p "Div=\$Div\n";
die "All done!\n";

That outputs:

Den=3800
Div=3800
All done!

on my machine.

In the original demo, the arithmetic was all being performed by the perl engine -- perhaps in single precision?

Cheers,
Craig

• Craig DeForest - 2007-03-19
• assigned_to: nobody --> zowie
• status: open --> closed-works-for-me

• Chris Marshall - 2007-03-19

Logged In: YES
user_id=44920
Originator: NO

Using Craig's code but setting the PDL data types to float
I get the same "incorrect" values as originally reported.

• John Harker - 2007-03-19

Logged In: YES
user_id=1747641
Originator: YES

Ah, yes, I see what you mean. Still... Compare these two:

use PDL;
use PDL::Math;
\$Num = 0.25;
\$Stp = (1/400 - 1/800)/19;

\$Den = \$Num / \$Stp;
rint(\$Den);

\$Ratio = \$Num / \$Den;
\$Div = \$Num / \$Ratio;
print "Den = \$Den\n";
print "Ratio = \$Ratio\n";
print "Div = \$Div\n";
die "All done!\n";

That outputs:

Den = 3800
Ratio = 6.57894736842105e-05
Div = 3800
All done!

...which is correct, even though the first division is apparently being done with perl scalars. However, one small change produces the error:

use PDL;
use PDL::Math;
\$Num = 0.25;
\$Stp = (1/400 - 1/800)/19;

\$Den = \$Num / \$Stp;
\$Den = rint(\$Den); # Changed

\$Ratio = \$Num / \$Den;
\$Div = \$Num / \$Ratio;
print "Den = \$Den\n";
print "Ratio = \$Ratio\n";
print "Div = \$Div\n";
die "All done!\n";

produces

Den = 3800
Ratio = 6.57894706819206e-05
Div = 3800.00024414062
All done!

Why should the results of the second be different from the first?

• John Harker - 2007-03-19

Logged In: YES
user_id=1747641
Originator: YES

Ah, sorry, I composed my last reply before Marshall's was posted, so I missed it. It appears that it's definitely related to single vs. double precision. This'll teach me to make sure everything is declared as a piddle to start with!

Many thanks!
John

• Craig DeForest - 2007-03-19

Logged In: YES
user_id=20200
Originator: NO

Ah, interesting. The "problem" here is that rint, etc. default to *float* when pdlifying Perl scalars. Check this out:

use PDL;
printf("%s is a %s\n", "rint(5.5)", type rint(5.5));
printf("%s is a %s\n", "rint(pdl(5.5))", type rint(pdl(5.5)));

outputs:
rint(5.5) is a float
rint(pdl(5.5)) is a double

So this could conceivably be considered a bug -- after all it violates the principle of least surprise. Accordingly I've re-opened this but with low priority. The floor is now open for discussion on what rint Should do... :-)

• Craig DeForest - 2007-03-19
• labels: --> 101696
• priority: 5 --> 2
• status: closed-works-for-me --> open

• Chris Marshall - 2007-03-19
• priority: 2 --> 5
• labels: 101696 -->
• status: open --> closed-works-for-me

• Chris Marshall - 2007-03-19

Logged In: YES
user_id=44920
Originator: NO

Well, since I thought perl floating point values were all
double precision, that definitely surprises me.

• Chris Marshall - 2007-03-19

Logged In: YES
user_id=44920
Originator: NO

Well, since I thought perl floating point values were all
double precision, that definitely surprises me.