On Nov 18, 2015, at 3:47 PM, Chris Marshall <devel.chm...@gmail.com
<mailto:devel.chm...@gmail.com>> wrote:
I agree that the this is a surprising result---even
if correct by our current implementation. Maybe
the default for NV perl scalars should be type
PDL_Double. Maybe a perl scalar promotion
setting that could set a default value, double
or other specific type, smallest size, or ...
--Chris
On Wed, Nov 18, 2015 at 4:58 PM, Craig DeForest
<defor...@boulder.swri.edu> wrote:
Yes, I believe this is a Bug in the conversion code logic.
Currently we use the smallest possible type that can contain a
value. Values that fit in a float get promoted to float and those
that require a double get promoted to double. Here’s a four-atan
trial, showing the autopromotion using a Perl scalar constructed
to be similar to 2479/1024, that cannot fit in a float.
sub try {
my $val = shift;
printf(“\nval: %.15f\nfloat: %.15f\ndefault:
%.15f\ndouble: %.15f\n”,
$val, atan(float($val)), atan($val), atan(double($val))
);
}
foreach( float(2379.0)/1024.0, double(2379.0)/1024.0,
2379.0/1024.0, 2479.00000001/1024.0 ) { try($_) }
yields
val: 2.323242187500000
float: 1.164332985877990
default: 1.164332985877990
double: 1.164332932171300
val: 2.323242187500000
float: 1.164332985877990
default: 1.164332932171300
double: 1.164332932171300
val: 2.323242187500000
float: 1.164332985877990
default: 1.164332985877990
double: 1.164332932171300
val: 2.420898437509766
float: 1.179073929786680
default: 1.179073913765370
double: 1.179073913765370
where the top is the floating-point value, second is the double
value (and both work as desired), third is the “bug” value (perl
scalar is promoted to floating point), and fourth is the “antibug”
value (perl scalar is promoted to double).
Because a float can represent 2379/1024 exactly (after all, it’s
just the integer 2379, shifted 10 bits to the right), you get a
float. But 1024 is a pretty special number, and highly unusual.
Changing the ratio even slightly yields promotion to double by
default.
pdl> try(2479/1025)
val: 2.418536585365854
float: 1.178729414939880
default: 1.178729370918130
double: 1.178729370918130
Since the vast majority of numeric Perl scalars cannot be
represented exactly in a float, it’s easy to get into the mindset
that Perl scalars are autopromoted to double PDLs. By the
Principle of Least Surprise, we should therefore always
autopromote to the most common case (double), and require explicit
cast to get to floating-point values.
Sorry for the wall of text here, this is a pretty cool bug.
On Nov 18, 2015, at 1:42 PM, Chris Marshall
<devel.chm...@gmail.com> wrote:
I think we may need something like
an option to enable/disable forced
conversion of perl scalars to double
piddles.
--Chris
On Wed, Nov 18, 2015 at 3:03 PM, Chris Marshall
<devel.chm...@gmail.com> wrote:
Hi Doug-
I think the problem is only when
a perl scalar is the input to atan().
The top two cases are the expected
result. Maybe the logic in the type
conversion is off. You don't say
what version of PDL but I see the
same output from PDL-2.014_02.
--Chris
On Wed, Nov 18, 2015 at 1:29 PM, Douglas Hunt
<dh...@ucar.edu> wrote:
Hi all: I just noticed a thing which caused me some trouble.
The PDL 'atan' function defaults to 'float' precision
instead of the
more natural 'double':
use PDL;
my $a1 = atan(double(2379.0)/double(1024.0));
my $a2 = atan(float(2379.0)/float(1024.0));
my $a3 = atan(2379.0/1024.0);
print "$a1, $a2, $a3\n";
This prints:
1.1643329321713, 1.1643328666687, 1.1643328666687
This feature causes a problem when porting to/from C as
the C atan
function operates on doubles.
I think this traces down to Basic::Math in the math.pd file:
my (@ufuncs1) = qw(acos asin atan cosh sinh tan tanh); #
F,D only
...
my $func;
foreach $func (@ufuncs1) {
pp_def($func,
HandleBad => 1,
NoBadifNaN => 1,
GenericTypes => ['F','D'],
Pars => 'a(); [o]b();',
Inplace => 1,
Doc => inplace_doc( $func ),
Code => code_ufunc($func),
BadCode => badcode_ufunc($func),
);
}
Note that GenericTypes isF, D instead of D, F.
I would go so far to say that this is a bug, as the
default pdl type is
double as in:
my $pdl = pdl(4,5,6);
p $pdl->info
PDL: Double D [3]
Regards,
Doug Hunt
------------------------------------------------------------------------------
_______________________________________________
pdl-devel mailing list
pdl-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-devel
------------------------------------------------------------------------------
_______________________________________________
pdl-devel mailing list
pdl-devel@lists.sourceforge.net
<mailto:pdl-devel@lists.sourceforge.net>
https://lists.sourceforge.net/lists/listinfo/pdl-devel
