in reply to Re^3: Determining the minimum representable increment/decrement possible?
in thread Determining the minimum representable increment/decrement possible?

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

[download]

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

[download]

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 );
}

[download]
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.
"Science is about questioning the status quo. Questioning authority". I knew I was on the right track :)
In the absence of evidence, opinion is indistinguishable from prejudice. Not understood.
.

Replies are listed 'Best First'.
Re^5: Determining the minimum representable increment/decrement possible?
by BrowserUk (Patriarch) on Jul 04, 2016 at 15:01 UTC

    I was asked off-line to explain what I meant by:

    he uses convoluted and complicated numerical methods rather than simple and pragmatic bitwise methods

    Which would you rather maintain or debug? My routine for formating a double into a hexadecimal floating point constant:

    use constant {
    SIGNBIT => 0x8000_0000_0000_0000,
    EXPBITS => 0x7FF0_0000_0000_0000,
    SIGBITS =>    0xF_FFFF_FFFF_FFFF,
    INDBITS =>    0x8_0000_0000_0000,
};

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;
}

    [download]

    Or same thing from Data::Float:

    my %float_hex_defaults = (
    infinite_string => "inf",
    nan_string => "nan",
    exp_neg_sign => "-",
    exp_pos_sign => "+",
    pos_sign => "+",
    neg_sign => "-",
    hex_prefix_string => "0x",
    subnormal_strategy => "SUBNORMAL",
    zero_strategy => "STRING=0.0",
    frac_digits => 0,
    frac_digits_bits_mod => "ATLEAST",
    frac_digits_value_mod => "ATLEAST",
    exp_digits => 0,
    exp_digits_range_mod => "IGNORE",
);

sub _float_hex_option($$) {
    my($options, $name) = @_;
    my $val = defined($options) ? $options->{$name} : undef;
    return defined($val) ? $val : $float_hex_defaults{$name};
}

use constant exp_digits_range => do {
    my $minexp = min_normal_exp - significand_bits;
    my $maxexp = max_finite_exp + 1;
    my $len_minexp = length(-$minexp);
    my $len_maxexp = length($maxexp);
    $len_minexp > $len_maxexp ? $len_minexp : $len_maxexp;
};
use constant frac_digits_bits => (significand_bits + 3) >> 2;
use constant frac_sections => do { use integer; (frac_digits_bits + 6)
+ / 7; };

sub float_hex($;$) {
    my($val, $options) = @_;
    return _float_hex_option($options, "nan_string") if $val != $val;
    if(have_infinite) {
        my $inf_sign;
        if($val == $pos_infinity) {
            $inf_sign = _float_hex_option($options, "pos_sign");
            EMIT_INFINITY:
            return $inf_sign.
                _float_hex_option($options, "infinite_string");
        } elsif($val == $neg_infinity) {
            $inf_sign = _float_hex_option($options, "neg_sign");
            goto EMIT_INFINITY;
        }
    }
    my($sign, $exp, $sgnf);
    if($val == 0.0) {
        $sign = float_sign($val);
        my $strat = _float_hex_option($options, "zero_strategy");
        if($strat =~ /\ASTRING=(.*)\z/s) {
            my $string = $1;
            return _float_hex_option($options,
                    $sign eq "-" ? "neg_sign" : "pos_sign").
                $string;
        } elsif($strat eq "SUBNORMAL") {
            $exp = min_normal_exp;
        } elsif($strat =~ /\AEXPONENT=([-+]?[0-9]+)\z/) {
            $exp = $1;
        } else {
            croak "unrecognised zero strategy `$strat'";
        }
        $sgnf = 0.0;
    } else {
        ($sign, $exp, $sgnf) = float_parts($val);
    }
    my $digits = int($sgnf);
    if($digits eq "0" && $sgnf != 0.0) {
        my $strat = _float_hex_option($options, "subnormal_strategy");
        if($strat eq "NORMAL") {
            my $add_exp;
            (undef, $add_exp, $sgnf) = float_parts($sgnf);
            $exp += $add_exp;
            $digits = "1";
        } elsif($strat eq "SUBNORMAL") {
            # do nothing extra
        } else {
            croak "unrecognised subnormal strategy `$strat'";
        }
    }
    $sgnf -= $digits;
    for(my $i = frac_sections; $i--; ) {
        $sgnf *= 268435456.0;
        my $section = int($sgnf);
        $digits .= sprintf("%07x", $section);
        $sgnf -= $section;
    }
    $digits =~ s/(.)0+\z/$1/;
    my $ndigits = 1 + _float_hex_option($options, "frac_digits");
    croak "negative number of digits requested" if $ndigits <= 0;
    my $mindigits = 1;
    my $maxdigits = $ndigits + frac_digits_bits;
    foreach my $constraint (["frac_digits_bits_mod", 1+frac_digits_bit
+s],
                ["frac_digits_value_mod", length($digits)]) {
        my($optname, $number) = @$constraint;
        my $mod = _float_hex_option($options, $optname);
        if($mod =~ /\A(?:ATLEAST|EXACTLY)\z/) {
            $mindigits = $number if $mindigits < $number;
        }
        if($mod =~ /\A(?:ATMOST|EXACTLY)\z/) {
            $maxdigits = $number if $maxdigits > $number;
        }
        croak "unrecognised length modification setting `$mod'"
            unless $mod =~ /\A(?:AT(?:MO|LEA)ST|EXACTLY|IGNORE)\z/;
    }
    croak "incompatible length constraints" if $maxdigits < $mindigits
+;
    $ndigits = $ndigits < $mindigits ? $mindigits :
            $ndigits > $maxdigits ? $maxdigits : $ndigits;
    if($ndigits > length($digits)) {
        $digits .= "0" x ($ndigits - length($digits));
    } elsif($ndigits < length($digits)) {
        my $chopped = substr($digits, $ndigits, length($digits), "");
        if($chopped =~ /\A[89abcdef]/ &&
                !($chopped =~ /\A80*\z/ &&
                  $digits =~ /[02468ace]\z/)) {
            for(my $i = length($digits) - 1; ; ) {
                my $d = substr($digits, $i, 1);
                $d =~ tr/0-9a-f/1-9a-f0/;
                substr($digits, $i, 1, $d);
                last unless $d eq "0";
            }
            if($digits =~ /\A2/) {
                $exp++;
                substr($digits, 0, 1, "1");
            }
        }
    }
    my $nexpdigits = _float_hex_option($options, "exp_digits");
    my $mod = _float_hex_option($options, "exp_digits_range_mod");
    if($mod eq "ATLEAST") {
        $nexpdigits = exp_digits_range
            if $nexpdigits < exp_digits_range;
    } elsif($mod ne "IGNORE") {
        croak "unrecognised exponent length ".
            "modification setting `$mod'";
    }
    $digits =~ s/\A(.)(.)/$1.$2/;
    return sprintf("%s%s%sp%s%0*d",
        _float_hex_option($options,
            $sign eq "-" ? "neg_sign" : "pos_sign"),
        _float_hex_option($options, "hex_prefix_string"),
        $digits,
        _float_hex_option($options,
            $exp < 0 ? "exp_neg_sign" : "exp_pos_sign"),
        $nexpdigits, abs($exp));
}

    [download]
    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.
    "Science is about questioning the status quo. Questioning authority". I knew I was on the right track :)
    In the absence of evidence, opinion is indistinguishable from prejudice. Not understood.
[reply]
[d/l]
[select]

In Section Seekers of Perl Wisdom