% halfprecision converts IEEE 754 floating point to half precision IEEE 754r
%******************************************************************************
%
% MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
% Function: halfprecision
% Filename: halfprecision.c
% Programmer: James Tursa
% Version: 1.0
% Date: March 3, 2009
% Copyright: (c) 2009 by James Tursa, All Rights Reserved
%
% This code uses the BSD License:
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% halfprecision converts the input argument to/from a half precision floating
% point bit pattern corresponding to IEEE 754r. The bit pattern is stored in a
% uint16 class variable. Please note that halfprecision is *not* a class. That
% is, you cannot do any arithmetic with the half precision bit patterns.
% halfprecision is simply a function that converts the IEEE 754r half precision
% bit pattern to/from other numeric MATLAB variables. You can, however, take
% the half precision bit patterns, convert them to single or double, do the
% operation, and then convert the result back manually.
%
% 1 bit sign bit
% 5 bits exponent, biased by 15
% 10 bits mantissa, hidden leading bit, normalized to 1.0
%
% Special floating point bit patterns recognized and supported:
%
% All exponent bits zero:
% - If all mantissa bits are zero, then number is zero (possibly signed)
% - Otherwise, number is a denormalized bit pattern
%
% All exponent bits set to 1:
% - If all mantissa bits are zero, then number is +Infinity or -Infinity
% - Otherwise, number is NaN (Not a Number)
%
% Building:
%
% halfprecision requires that a mex routine be built (one time only). This
% process is typically self-building the first time you call the function
% as long as you have the files halfprecision.m and halfprecision.c in the
% same directory somewhere on the MATLAB path. If you need to manually build
% the mex function, here are the commands:
%
% >> mex -setup
% (then follow instructions to select a C / C++ compiler of your choice)
% >> mex halfprecision.c
%
% If you have an older version of MATLAB, you may need to use this command:
%
% >> mex -DDEFINEMWSIZE halfprecision.c
%
% Syntax
%
% B = halfprecision(A)
% C = halfprecision(B,S)
% halfprecision(B,'disp')
%
% Description
%
% A = a MATLAB numeric array, char array, or logical array.
%
% B = the variable A converted into half precision floating point bit pattern.
% The bit pattern will be returned as a uint16 class variable. The values
% displayed are simply the bit pattern interpreted as if it were an unsigned
% 16-bit integer. To see the halfprecision values, use the 'disp' option, which
% simply converts the bit patterns into a single class and then displays them.
%
% C = the half precision floating point bit pattern in B converted into class S.
% B must be a uint16 or int16 class variable.
%
% S = char string naming the desired class (e.g., 'single', 'int32', etc.)
% If S = 'disp', then the floating point bit values are simply displayed.
%
% Examples
%
% >> a = [-inf -1e30 -1.2 NaN 1.2 1e30 inf]
% a =
% 1.0e+030 *
% -Inf -1.0000 -0.0000 NaN 0.0000 1.0000 Inf
%
% >> b = halfprecision(a)
% b =
% 64512 64512 48333 65024 15565 31744 31744
%
% >> halfprecision(b,'disp')
% -Inf -Inf -1.2002 NaN 1.2002 Inf Inf
%
% >> halfprecision(b,'double')
% ans =
% -Inf -Inf -1.2002 NaN 1.2002 Inf Inf
%
% >> 2^(-24)
% ans =
% 5.9605e-008
%
% >> halfprecision(ans)
% ans =
% 1
%
% >> halfprecision(ans,'disp')
% 5.9605e-008
%
% >> 2^(-25)
% ans =
% 2.9802e-008
%
% >> halfprecision(ans)
% ans =
% 1
%
% >> halfprecision(ans,'disp')
% 5.9605e-008
%
% >> 2^(-26)
% ans =
% 1.4901e-008
%
% >> halfprecision(ans)
% ans =
% 0
%
% >> halfprecision(ans,'disp')
% 0
%
% Note that the special cases of -Inf, +Inf, and NaN are handled correctly.
% Also, note that the -1e30 and 1e30 values overflow the half precision format
% and are converted into half precision -Inf and +Inf values, and stay that
% way when they are converted back into doubles.
%
% For the denormalized cases, note that 2^(-24) is the smallest number that can
% be represented in half precision exactly. 2^(-25) will convert to 2^(-24)
% because of the rounding algorithm used, and 2^(-26) is too small and underflows
% to zero.
%
%**************************************************************************
function varargout = halfprecision(varargin)
disp(' ');
disp('You must build the mex routine before you can use halfprecision.');
disp('Attempting to do so now ...');
disp(' ');
mname = mfilename('fullpath');
cname = [mname '.c'];
if( isempty(dir(cname)) )
disp('Cannot find the file halfprecision.c in the same directory as the');
disp('file halfprecision.m. Please ensure that they are in the same');
disp('directory and try again. The following file was not found:');
disp(' ');
disp(cname);
disp(' ');
error('Unable to compile halprecision.c');
else
disp(['Found file halfprecision.c in ' cname]);
disp(' ');
disp('Now attempting to compile ...');
disp('(If prompted, please press the Enter key and then select any C/C++');
disp('compiler that is available, such as lcc.)');
disp(' ');
disp(['mex(''' cname ''')']);
disp(' ');
try
mex(cname);
disp('mex halfprecision.c build completed ... you may now use halfprecision.');
disp(' ');
catch
disp(' ');
disp('Well, *that* didn''t work ... now trying it with mwSize defined ...');
disp(' ');
try
disp(' ');
disp(['mex(''-DDEFINEMWSIZE'',''' cname ''')']);
disp(' ');
mex('-DDEFINEMWSIZE',cname);
disp('mex halfprecision.c build completed ... you may now use halfprecision.');
disp(' ');
catch
disp('Hmmm ... That didn''t work either.');
disp(' ');
disp('The mex command failed. This may be because you have already run');
disp('mex -setup and selected a non-C compiler, such as Fortran. If this');
disp('is the case, then rerun mex -setup and select a C/C++ compiler.');
disp(' ');
error('Unable to compile halprecision.c');
end
end
end
if false
varargout = varargin; % Get rid of the lint message
end
end