Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

#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 = {
'BADVAL_PER_PDL' => '0',
'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_POSIX_THREADS' => '1',
'HIDE_TRYLINK' => '1',
'WITH_HDF' => undef,
'HDF_INC' => undef,
'OPENGL_LIBS' => undef,
'WITH_BADVAL' => '0',
'WITH_GD' => undef,
'FITS_LEGACY' => '1',
'WITH_SLATEC' => '1',
'BADVAL_USENAN' => '0',
'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:
osname=linux, osvers=2.6.15.7, archname=i486-linux-gnu-thread-multi
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
usethreads=define use5005threads=undef useithreads=define usemultiplicity=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
Linker and Libraries:
ld='cc', ldflags =' -L/usr/local/lib'
libpth=/usr/local/lib /lib /usr/lib
libs=-lgdbm -lgdbm_compat -ldb -ldl -lm -lpthread -lc -lcrypt
perllibs=-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'
Dynamic Linking:
dlsrc=dl_dlopen.xs, dlext=so, d_dlsymun=undef, ccdlflags='-Wl,-E'
cccdlflags='-fPIC', lddlflags='-shared -L/usr/local/lib'
===========================================

Discussion

  • John Harker
    John Harker
    2007-03-19

    • summary: rint, ceil, floor should return their expected result --> rint, ceil, floor sometimes return the wrong value?
     
  • John Harker
    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
    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
    Craig DeForest
    2007-03-19

    • assigned_to: nobody --> zowie
    • status: open --> closed-works-for-me
     
  • Chris Marshall
    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
    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
    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
    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
    Craig DeForest
    2007-03-19

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

    • priority: 2 --> 5
    • labels: 101696 -->
    • status: open --> closed-works-for-me
     
  • Chris Marshall
    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
    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.