function b = rnd24(a)
%   b = rnd24(a);
% This function rounds a to 24 bits for simulated IEEE single
% precision.  The input argument a is a scalar, vector or matrix.
% The return variable b is a scalar, matrix or vector of the
% same size, but rounded.
%
% For example, try:
%    >> format long
%    >> rnd24(-1:0.1:1)


% Save the signs in sg, and take absolute values.

sg = sign(a);
a1 = abs(a);

% This next statement is intended to prevent warnings of log of 0
% in the statement that follows.
a1 = a1 + (a1 == 0); 

% Take the log base 2 of a1.

r = log(a1) / log(2);

% Get the integer exponent in base 2.
e = floor(r);

% Recover the mantissa, between 1 and 2.
m = exp((r - e) * log(2));

% Round the mantissa to 24 bits by multiplying it by 2^23
% rounding, and dividing by 2^23.  Then multiply by the
% saved exponents, and by the signs.

b = round(m*2^23) / 2^23 .* exp(e * log(2)) .* sg;


