function [F,Fc] = floatRep(M_bit,e_bit)
%floatRep(M_bit,e_bit)
% take a look at binary floating point representations for a floating point
% number F = s*M*B^e for M_bit mantissa and e_bit exponent

% (c) 2007 Curt A. L. Szuberla
% University of Alaska Fairbanks, all rights reserved

s = [-1 1]'; % not really useful here
Mv = (0:2^M_bit-1)/2^M_bit;
B = 2;
ev = (-2^(e_bit-1)+1:2^(e_bit-1));

Bev = B.^ev;

bits = 1+M_bit+e_bit; % total # of bits in F

disp(['M = ' rats(Mv)])
disp(['e = ' rats(ev)])
disp(['B^e = ' rats(Bev)])

M = cols(Mv',length(Bev));
Be = rows(Bev,length(Mv));

Fc = M.*Be;
Fc = [Fc;-Fc]; % this step does s*M*B^e
Fc = sort(Fc(:));
F = unique(Fc);

disp(['F (complete{' num2str(length(Fc)) '}) = ' rats(Fc')])
disp(['F (unique{' num2str(length(F)) '}) = ' rats(F')])

% figure differences from 0-point
dF = diff(F((length(F)+1)/2:end));
dF = [flipud(dF);0;dF];

clf
plot(F,log2(dF),'bo',F,0,'rx')
set(gca,'ylim',[-bits 1])
xlabel('F = s \cdot M \cdot B^e')
ylabel('log_2(\delta F)')
title(['F_' num2str(bits) ' (s_1, M_' num2str(M_bit) ...
    ', B=2, e_' num2str(e_bit) ')' ])
powerpoint

return

function lg2 = log2(x);
%log(2) = 0.693147180559945;
lg2 = log(x)/0.693147180559945;
return