Bugs item #1683903, was opened at 2007-03-19 15:27
Message generated for change (Comment added) made by marshallch
You can respond by visiting:
https://sourceforge.net/tracker/?func=detail&atid=100612&aid=1683903&group_id=612
Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: None
Group: None
Status: Closed
Resolution: Works For Me
Priority: 5
Private: No
Submitted By: John Harker (jharker)
Assigned to: Craig DeForest (zowie)
Summary: rint, ceil, floor sometimes return the wrong value?
Initial Comment:
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'
===========================================
----------------------------------------------------------------------
>Comment By: Chris Marshall (marshallch)
Date: 2007-03-19 17:35
Message:
Logged In: YES
user_id=44920
Originator: NO
Well, since I thought perl floating point values were all
double precision, that definitely surprises me.
----------------------------------------------------------------------
Comment By: Chris Marshall (marshallch)
Date: 2007-03-19 17:35
Message:
Logged In: YES
user_id=44920
Originator: NO
Well, since I thought perl floating point values were all
double precision, that definitely surprises me.
----------------------------------------------------------------------
Comment By: Craig DeForest (zowie)
Date: 2007-03-19 17:18
Message:
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...
:-)
----------------------------------------------------------------------
Comment By: John Harker (jharker)
Date: 2007-03-19 17:09
Message:
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
----------------------------------------------------------------------
Comment By: John Harker (jharker)
Date: 2007-03-19 17:03
Message:
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?
----------------------------------------------------------------------
Comment By: Chris Marshall (marshallch)
Date: 2007-03-19 17:00
Message:
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.
----------------------------------------------------------------------
Comment By: Craig DeForest (zowie)
Date: 2007-03-19 16:27
Message:
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
----------------------------------------------------------------------
Comment By: John Harker (jharker)
Date: 2007-03-19 15:52
Message:
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.
----------------------------------------------------------------------
You can respond by visiting:
https://sourceforge.net/tracker/?func=detail&atid=100612&aid=1683903&group_id=612
_______________________________________________
PDL-porters mailing list
PDL-porters@jach.hawaii.edu
http://mailman.jach.hawaii.edu/mailman/listinfo/pdl-porters
