Hi all: When testing the new PDL interface to plplot, we came across this
odd PDL behavior:
use PDL;
print "acos(-1.) = ", acos(-1.), "\n";
print "acos(-1 ) = ", acos(-1 ), "\n";
The acos(-1.) version incorrectly gives a single precision result, whereas
the acos(1) version gives the expected double precision result. This is
due to the fact that there is no integer version of acos and PP routines
default to 'double' instead of float in this case.
This stems from the fact that:
perldl> p PDL::get_datatype(-1.)
5 # $PDL_F
Which seems odd--most PDL operations default to double precision.
After a bit of digging with gdb, it turns out this is caused by:
From: pdlcore.c.PL
/* Check minimum, at least float, datatype required to represent number */
int pdl_whichdatatype_double (double nv) {
TESTTYPE(PDL_F,PDL_Float)
TESTTYPE(PDL_D,PDL_Double)
if( !finite(nv) ) { return PDL_D; }
croak("Something's gone wrong: %lf cannot be converted by
whichdatatype_double",
nv);
}
So, I would like to propose one of two remedies (attached).
The first (pdlcore.c.PL.diff) would be to get rid of
pdl_whichdatatype_double and just assume that scalars convert to doubles,
not floats. This would be a change to a long-standing rule, but I think
it would be sensible. It breaks one test in conv.t that was specifically
designed to verify this behavior--that test would have to be changed.
Another approach would be to leave the 'PDL::get_datatype(-1.) == $PDL_F'
problem alone and just fix 'acos' and the rest of the math functions that
can return either float or double values. The other attached patch
(math.pd.diff) just changes
GenericTypes => ['F', 'D']
to
GenericTypes => ['D', 'F']
in several pp_def calls. This has less impact, but solves the irksome
'acos(-1.) != acos(-1)' problem without changing the fact that '1.'
defaults to 'float' instead of 'double'.
I've tested both of these out and (with the exception of the conv.t test)
they both pass the test suite and solve the problem.
Any ideas on which approach to take?
Regards,
Doug
dh...@ucar.edu
Software Engineer IV
UCAR - COSMIC, Tel. (303) 497-2611
*** ./Basic/Core/pdlcore.c.PL.orig Wed Nov 12 21:13:53 2008
--- ./Basic/Core/pdlcore.c.PL Fri Dec 12 17:56:11 2008
***************
*** 175,192 ****
nv);
}
- /* Check minimum, at least float, datatype required to represent number */
-
- int pdl_whichdatatype_double (double nv) {
- TESTTYPE(PDL_F,PDL_Float)
- TESTTYPE(PDL_D,PDL_Double)
-
- if( !finite(nv) ) { return PDL_D; }
-
- croak("Something's gone wrong: %lf cannot be converted by
whichdatatype_double",
- nv);
- }
-
#ifdef _MSC_VER
#pragma optimize("", on)
#endif
--- 175,180 ----
***************
*** 258,293 ****
/* Figure datatype to use */
if ( !SvIOK(sv) && SvNOK(sv) && SvNIOK(sv) ) {/* Perl Double (e.g.
2.0) */
data = SvNV(sv);
!
! !NO!SUBS!
!
! # XXX HACK this may not be sensible (DJB 08/31/00)
! # - only relevant if BADVAL_USENAN is true in config file
! #
!
! if ( $bvalflag and $usenan ) {
! print OUT <<'!NO!SUBS!';
!
! /*
! * default NaN/Infs to double
! * XXX sensible ?
! */
! if ( finite(data) == 0 ) {
! datatype = PDL_D;
! } else {
! datatype = pdl_whichdatatype_double(data);
! }
!
! !NO!SUBS!
! } else {
!
! print OUT "\tdatatype = pdl_whichdatatype_double(data);\n";
!
! } # if: $bvalflag
!
! print OUT <<'!NO!SUBS!';
!
! }
else { /* Perl Int (e.g. 2) */
data = SvNV(sv);
datatype = pdl_whichdatatype(data);
--- 246,253 ----
/* Figure datatype to use */
if ( !SvIOK(sv) && SvNOK(sv) && SvNIOK(sv) ) {/* Perl Double (e.g.
2.0) */
data = SvNV(sv);
! datatype = PDL_D;
! }
else { /* Perl Int (e.g. 2) */
data = SvNV(sv);
datatype = pdl_whichdatatype(data);
*** ./Basic/Math/math.pd.orig Wed Nov 12 21:13:54 2008
--- ./Basic/Math/math.pd Fri Dec 12 17:49:11 2008
***************
*** 167,173 ****
pp_def($func,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['F','D'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $func ),
--- 167,173 ----
pp_def($func,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['D','F'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $func ),
***************
*** 205,211 ****
pp_def($func,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['F','D'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $func ),
--- 205,211 ----
pp_def($func,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['D', 'F'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $func ),
***************
*** 219,225 ****
pp_def($fname,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['F','D'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $fname ),
--- 219,225 ----
pp_def($fname,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['D', 'F'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $fname ),
***************
*** 233,239 ****
pp_def($fname,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['F','D'],
Pars => 'a(); int n(); [o]b();',
Inplace => [ 'a' ],
Doc => inplace_doc( $fname ),
--- 233,239 ----
pp_def($fname,
HandleBad => 1,
NoBadifNaN => 1,
! GenericTypes => ['D', 'F'],
Pars => 'a(); int n(); [o]b();',
Inplace => [ 'a' ],
Doc => inplace_doc( $fname ),
------------------------------------------------------------------------------
SF.Net email is Sponsored by MIX09, March 18-20, 2009 in Las Vegas, Nevada.
The future of the web can't happen without you. Join us at MIX09 to help
pave the way to the Next Web now. Learn more and register at
http://ad.doubleclick.net/clk;208669438;13503038;i?http://2009.visitmix.com/
_______________________________________________
Plplot-devel mailing list
Plplot-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/plplot-devel
