% function dumpfp(d)
% function [sgnb_bin, n, man_bin, expn_bin, expn] = dumpfp(d)
%
% Dump the binary representation of a IEEE floating point double on i386
% architectures.
%
% Without output arguments the number is dumped to the screen.
%
% Input:
%   d           number to dump, scalar
% Output:
%   sgnb_bin    sign bit, '0' = positive, '1' = negative, char array size 1
%   n           MSB, char array which is '1.' for normal numbers, '0.' for
%               denormalized numbers and 'x.' for +/- infinity and not-a-number
%   man_bin     mantissa or significand bits without the MSB, char array size 52
%   expn_bin    exponent bits, char array size 11
%   expn        exponent as decimal number, corrected with offset 1023
%
% Examples:
%   >> dumpfp(-3E-324)
%   -4.9406564584124654e-324 -0.0000000000000000000000000000000000000000000000000001 2^00000000000 (denorm)
%   >> dumpfp(1)
%    1                       +1.0000000000000000000000000000000000000000000000000000 2^01111111111 (2^0)
%   >> dumpfp(eps)
%    2.2204460492503131e-16  +1.0000000000000000000000000000000000000000000000000000 2^01111001011 (2^-52)
%   >> dumpfp(1+eps)
%    1.0000000000000002      +1.0000000000000000000000000000000000000000000000000001 2^01111111111 (2^0)
%
% Thu Oct  2 15:03:31 CEST 2003, <dirk.nuyens@cs.kuleuven.ac.be>
function [varargout] = dumpfp(d)

old_matlab = all(sprintf('%bx', 1+eps) == '010000000000f03f');
if old_matlab
    s = sprintf('%bx', d);
else
    % new Matlab does give this hex string in reverse order,
    % also changes in char() function for bytereversion
    s = fliplr(sprintf('%bx', d));
end;
ss = zeros(size(s));
ss([1:2:end-1]+1) = s([1:2:end-1]); % swap bytes
ss([1:2:end-1]) = s([1:2:end-1]+1); %
if old_matlab, s = char(ss); end;

man1_hex = s([1:8]);
man0_hex = s([1:5]+8);
expn_sgnb_hex = s([1:3]+13);

P = fliplr(eye(4));
man1_bin = fliplr(char(double(dec2bin(hex2dec(man1_hex), 32)) * kron(eye(32/4), P)));
man0_bin = fliplr(char(double(dec2bin(hex2dec(man0_hex), 20)) * kron(eye(20/4), P)));
expn_sgnb_bin = fliplr(char(double(dec2bin(hex2dec(expn_sgnb_hex), 12)) * kron(eye(12/4), P)));
expn_bin = expn_sgnb_bin(2:end);
sgnb_bin = expn_sgnb_bin(1);

expn = bin2dec(expn_bin)-1023;
n = '1.';
if expn == -1023, expn_str = 'denorm'; n = '0.';
elseif expn == 1024, expn_str = 'Inf/NaN'; n = 'x.';
else expn_str = sprintf('2^%d', expn);
end;
if sgnb_bin == '1', sgnb = '-';
else sgnb = '+';
end;
man_bin = [man0_bin man1_bin];

if nargout == 0
    fprintf('% -24.17g %s%s%s 2^%s (%s)\n', d, sgnb, n, man_bin, ...
            expn_bin, expn_str);
else
    varargout(1) = {sgnb_bin};
    varargout(2) = {n};
    varargout(3) = {man_bin};
    varargout(4) = {expn_bin};
    varargout(5) = {expn};
end;