use warnings;
use strict;

# http://www.perlmonks.org/bare/?node_id=703222

sub double_from_hex { unpack 'd', scalar reverse pack 'H*', $_[0] }

use constant POS_INF => double_from_hex '7FF0000000000000';
use constant NEG_INF => double_from_hex 'FFF0000000000000';
use constant NaN     => double_from_hex '7FF8000000000000';

my @nums = (
'b6 a6 ac 29 6b ca ea 3f',
'd0 06 d8 4a 48 37 eb 3f',
'd8 a5 51 20 82 63 e7 bf',
'10 44 70 ec d2 28 c0 3f',
'00 00 00 00 00 00 f0 7f',
);

for (@nums) {
    my @bytes = map {chr (hex $_)} split / /;
    my $bytes = join '', @bytes;
    print "", parse_double ($bytes), "\n";
}
exit;

=head2 parse_double

    my $double = parse_double ($bytes);

This converts the MAT-File double (eight bytes, the IEEE 754 "double
64" format) into a Perl floating point number.

=cut

sub parse_double
{
    my ($bytes) = @_;
    my ($bottom, $top) = unpack ("LL", $bytes);
    # Reference:
    # http://en.wikipedia.org/wiki/Double_precision_floating-point_format

    # Eight zero bytes represents 0.0.
    if ($bottom == 0) {
        if ($top == 0) {
            return 0;
        }
        elsif ($top == 0x80000000) {
            return -0;
        }
        elsif ($top == 0x7ff00000) {
            return POS_INF;
        }
        elsif ($top == 0xfff00000) {
            return NEG_INF;
        }
    }
    elsif ($top == 0x7ff00000) {
        return NaN;
    }
    my $sign = $top >> 31;
#    print "$sign\n";
    my $exponent = (($top >> 20) & 0x7FF) - 1023;
#    print "$exponent\n";
    my $e = ($top >> 20) & 0x7FF;
    my $t = $top & 0xFFFFF;
#    printf ("--> !%011b%020b \n--> %032b\n", $e, $t, $top);
    my $mantissa = ($bottom + ($t*(2**32))) / 2**52 + 1;
#    print "$mantissa\n";
    my $double = (-1)**$sign * 2**$exponent * $mantissa;
    return $double;
}