Hi Craig: Thanks for exploring this further. I'm glad it's a cool bug
and not a boring bug or a nerdy bug!
--Doug
On 11/18/15 2:58 PM, Craig DeForest 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
<mailto: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 <mailto: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
<mailto: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
<mailto: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
------------------------------------------------------------------------------
_______________________________________________
pdl-devel mailing list
pdl-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-devel
