I was googling to fix a bug, and was led to Data::Float, which actually has most of the tools I was trying to implement in my module.
You should not let that stop you. That module gets (mostly) the right answers, but boy does it ever go about it in a convoluted way; which makes it very slow even for pure perl code.
For example, this compares his method of nextup() and nextdown() with a couple of trivial routines I threw together when exploring this earlier:
C:\test>FPstuff.pl
Took:9.723010063 seconds ## his
Took:0.949742079 seconds ## mine
And this compares his float_hex() routine with my asHex(): C:\test>FPstuff.pl
"my" variable $start masks earlier declaration in same scope at C:\tes
+t\FPstuff.pl line 108.
Took:9.740586996 seconds ## his
Took:0.699841022 seconds ## mine
His are more than a magnitude slower in both cases -- mostly because he uses convoluted and complicated numerical methods rather than simple and pragmatic bitwise methods -- and certainly for my OP problem of trying to tune some floating point optimisations, that would be a show stopper.
Here is the code that produce the above timings -- is just a mish-mash of bits and pieces from other scripts thrown together to test and validate that module. If any of it is useful to you; feel free to use it: #! perl -slw
use Config;
use Inline C => Config => BUILD_NOISY => 1;
use Inline C => <<'END_C', NAME => 'nextafter', CLEAN_AFTER_BUILD =>0
+;
#include <math.h>
double cNextAfter( double x, double y ) {
return _nextafter( x, y );
}
SV *cAsHex( double d ) {
SV * ret = newSVpv( "123456789012345678901234", 24 );
SvCUR( ret ) = sprintf( SvPVX(ret), "% -24.13a", d );
return ret;
}
END_C
use constant DEBUG => $DEBUG;
use strict;
no warnings "portable";
use Data::Float qw[ nextafter nextup nextdown float_hex ];
use Time::HiRes qw[ time ];
use constant {
SIGNBIT => 0x8000_0000_0000_0000,
EXPBITS => 0x7FF0_0000_0000_0000,
SIGBITS => 0xF_FFFF_FFFF_FFFF,
INDBITS => 0x8_0000_0000_0000,
};
use constant {
POS_QNAN_LST => 0b0_11111111111_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
POS_QNAN_1ST => 0b0_11111111111_1000_00000000_00000000_00000000
+_00000000_00000000_00000000,
POS_SNAN_LST => 0b0_11111111111_0111_11111111_11111111_11111111
+_11111111_11111111_11111111,
POS_SNAN_1ST => 0b0_11111111111_0000_00000000_00000000_00000000
+_00000000_00000000_00000001,
POS_INF => 0b0_11111111111_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
POS_NORM_LST => 0b0_11111111110_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
POS_NORM_1ST => 0b0_00000000001_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
POS_DENORM_LST => 0b0_00000000000_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
POS_DENORM_1ST => 0b0_00000000000_0000_00000000_00000000_00000000
+_00000000_00000000_00000001,
POS_ZERO => 0b0_00000000000_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
NEG_ZERO => 0b1_00000000000_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
NEG_DENORM_1ST => 0b1_00000000000_0000_00000000_00000000_00000000
+_00000000_00000000_00000001,
NEG_DENORM_LST => 0b1_00000000000_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
NEG_NORM_1ST => 0b1_00000000001_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
NEG_NORM_LST => 0b1_11111111110_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
NEG_INF => 0b1_11111111111_0000_00000000_00000000_00000000
+_00000000_00000000_00000000,
NEG_SNAN_1ST => 0b1_11111111111_0000_00000000_00000000_00000000
+_00000000_00000000_00000001,
NEG_SNAN_LST => 0b1_11111111111_0111_11111111_11111111_11111111
+_11111111_11111111_11111111,
NEG_IND => 0b1_11111111111_1000_00000000_00000000_00000000
+_00000000_00000000_00000000,
NEG_QNAN_1ST => 0b1_11111111111_1000_00000000_00000000_00000000
+_00000000_00000000_00000001,
NEG_QNAN_LST => 0b1_11111111111_1111_11111111_11111111_11111111
+_11111111_11111111_11111111,
};
use enum qw[ X Y ];
sub assert{
return unless DEBUG;
die sprintf "%s - %s(%d): Assertion failed.\n", caller() unless $_
+[0];
}
sub asBits{ join '_', unpack 'A1 A11 A4 (A8)6', unpack 'B64', scalar r
+everse pack 'd', $_[X] }
sub asHex{
my $bin = unpack 'Q', pack 'd', $_[X];
return "-0x1.#IND000000000p+0" if $bin == NEG_IND;
my $sign = ( $bin & SIGNBIT ) ? '-' : ' ';
my $exp = ( ( $bin & EXPBITS ) >> 52 );
my $mant = $bin & SIGBITS;
if( $exp == 2047 ) {
return $sign . "0x1.#INF000000000p+0" unless $mant;
return $sign . ( $mant < INDBITS ? "0x1.#SNAN00000000p+0" : "0
+x1.#QNAN00000000p+0" );
}
$exp -= 1023;
my $hid = $exp == -1023 ? ( ++$exp, 0 ) : 1;
sprintf "%s0x%1u.%013xp%+d", $sign, $hid, $mant, $exp;
}
sub asDouble{ unpack 'd', pack 'Q', $_[X] }
sub nextUp{
return $_[X] if $_[X] != $_[X];
my $i = unpack 'Q', pack 'd', $_[X];
return asDouble( NEG_NORM_LST ) if $i == NEG_INF;
return $_[X] if ( $i & POS_INF ) == POS_INF;
return asDouble( POS_DENORM_1ST ) if $i == NEG_ZERO;
return asDouble( $_[X] < 0 ? --$i : ++$i );
}
sub nextDown { - nextUp( - $_[X] ) }
sub nextAfter{
return $_[Y] if $_[X] == $_[Y];
return $_[X] < $_[Y] ? nextUp( $_[X] ) : nextDown( $_[X] );
}
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $hex = float_hex( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $hex = asHex( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
__END__
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $up = nextup( $rv );
my $dn = nextdown( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $up = nextUp( $rv );
my $dn = nextDown( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
__END__
for my $first (
2e-300, 2e-200, 2e-100, 2e-10, 2e-1, 2e0, 2e1, 2e10, 2e100, 2e
+200, 2e300,
map asDouble( $_ ),
NEG_SNAN_1ST, NEG_SNAN_LST, NEG_IND, NEG_QNAN_1ST, NEG_QNAN_LST,
NEG_INF, NEG_NORM_LST, NEG_NORM_1ST, NEG_DENORM_LST, NEG_DENORM_1S
+T,
NEG_ZERO,
POS_ZERO,
POS_DENORM_1ST, POS_DENORM_LST, POS_NORM_1ST, POS_NORM_LST, POS_IN
+F,
POS_SNAN_1ST, POS_SNAN_LST, POS_QNAN_1ST, POS_QNAN_LST,
) {
printf "%32.18g\t%-24s\t%-24s\n", $first,
asHex( nextAfter( $first, asDouble( POS_INF ) ) ),
asHex( nextafter( $first, asDouble( POS_INF ) ) );
printf "%32.18g\t%-24s\t%-24s\n", $first,
asHex( nextAfter( $first, asDouble( NEG_INF ) ) ),
asHex( nextafter( $first, asDouble( NEG_INF ) ) );
}
__END__
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $hex = float_hex( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
srand( 1 );
my $start = time;
for( 1 .. 100000 ) {
my $rv = unpack 'd', pack 'VV', int( rand 65536 ), int( rand 65536
+ );
my $hex = asHex( $rv );
}
printf "Took:%.9f seconds\n", time() - $start;
__END__
for my $first ( map asDouble( $_ ),
NEG_INF, NEG_NORM_LST, NEG_NORM_1ST, NEG_DENORM_LST, NEG_DENORM_1S
+T, NEG_ZERO,
POS_ZERO, POS_DENORM_1ST, POS_DENORM_LST, POS_NORM_1ST, POS_NORM_L
+ST, POS_INF,
) {
printf "nextDown(%s) % 32.17e\n (%s) % 32.17e\n nextUp(%s)
+ % 32.17e\n\n",
asBits( nextDown( $first ) ), nextDown( $first ),
asBits( $first ), $first,
asBits( nextUp( $first ) ), nextUp( $first );
}
for my $first (
2e-300, 2e-200, 2e-100, 2e-10, 2e-1, 2e0, 2e1, 2e10, 2e100, 2e
+200, 2e300,
map asDouble( $_ ),
NEG_SNAN_1ST, NEG_SNAN_LST, NEG_IND, NEG_QNAN_1ST, NEG_QNAN_LST,
NEG_INF, NEG_NORM_LST, NEG_NORM_1ST, NEG_DENORM_LST, NEG_DENORM_1S
+T,
NEG_ZERO,
POS_ZERO,
POS_DENORM_1ST, POS_DENORM_LST, POS_NORM_1ST, POS_NORM_LST, POS_IN
+F,
POS_SNAN_1ST, POS_SNAN_LST, POS_QNAN_1ST, POS_QNAN_LST,
) {
printf "%s\t%32.18g\t%s\t%s\n", asBits( $first ), $first, cAsHex(
+$first ), asHex( $first );
}
With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
In the absence of evidence, opinion is indistinguishable from prejudice. Not understood. |