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
> <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
> <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
------------------------------------------------------------------------------
_______________________________________________
pdl-devel mailing list
pdl-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/pdl-devel
