1## Copyright 2014-2016 Oliver Heimlich
2##
3## This program is free software; you can redistribute it and/or modify
4## it under the terms of the GNU General Public License as published by
5## the Free Software Foundation; either version 3 of the License, or
6## (at your option) any later version.
7##
8## This program is distributed in the hope that it will be useful,
9## but WITHOUT ANY WARRANTY; without even the implied warranty of
10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11## GNU General Public License for more details.
12##
13## You should have received a copy of the GNU General Public License
14## along with this program; if not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @documentencoding UTF-8
18## @deftypefun {[@var{X}, @var{ISEXACT}] =} hex2double (@var{S}, @var{DIRECTION})
19##
20## Convert a hexadecimal floating point number @var{S} to double precision with
21## directed rounding.
22##
23## The input number format is [+-]0xh[,.]h[[pP][+-]d].
24##
25## @seealso{str2decimal}
26## @end deftypefun
27
28## Author: Oliver Heimlich
29## Created: 2014-10-20
30
31function [binary, isexact] = hex2double (string, direction)
32
33  ## Strip Sign
34  if (isempty (string))
35    hex.s = false;
36  else
37    hex.s = (string (1) == "-");
38    if (strfind("+-", string(1)))
39      string = string (2:end);
40    endif
41  endif
42
43  ## Strip hex indicator
44  if (strncmpi ("0x", string, 2))
45    string = string (3:end);
46  else
47    error ("interval:InvalidOperand", ...
48           ["invalid hex number does not start with 0x: " string]);
49  endif
50
51  ## Split mantissa & exponent
52  [hex.m, hex.e] = strtok (string, "pP");
53
54  ## Convert exponent from string to number
55  if (isempty (hex.e))
56    hex.e = int64(0); # 2 ^ 0 = 1
57  else
58    if (strfind (hex.e, ".") || strfind (hex.e, ","))
59      error ("interval:InvalidOperand", ...
60             ["invalid hex number with rational exponent: " string]);
61    endif
62    hex.e = str2double (hex.e(2:end)); # remove “p” and convert
63    if (isnan (hex.e) || ... # greater than realmax or illegal format
64        abs(hex.e) >= pow2 (53)) # possibly lost precision
65      error ("interval:InvalidOperand", ...
66             ["invalid hex number with big/invalid exponent: " string]);
67    endif
68    hex.e = int64 (hex.e);
69  endif
70
71  ## Split Mantissa at point
72  hex.m = strsplit (hex.m, {".",","});
73  switch length(hex.m)
74    case 0
75      hex.m = {"", ""};
76    case 1
77      hex.m{2} = "";
78    case 2
79                                # nothing to do
80    otherwise
81      error ("interval:InvalidOperand", ...
82             ["invalid hex number with multiple points: " string]);
83  endswitch
84
85  ## Normalize mantissa string: move point to the right
86  if (hex.e - length (hex.m{2}) * 4 <= intmin (class (hex.e)))
87    error ("interval:InvalidOperand", ...
88           ["exponent overflow during normalization: " string ]);
89  endif
90  hex.e -= length (hex.m{2}) * 4;
91  hex.m = strcat (hex.m{1}, hex.m{2});
92
93  hexvalues = rem (uint8 (hex.m) - 47, 39); # 1 .. 16
94  lookup = [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1;...
95            0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1;...
96            0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1;...
97            0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
98  hex.m = (logical (lookup (:, hexvalues))) (:)';
99
100  ## Normalize mantissa: remove leading zeroes
101  firstnonzerodigit = find (hex.m, 1, "first");
102  if (firstnonzerodigit > 1)
103    hex.m = hex.m (firstnonzerodigit:end);
104  elseif (isempty (firstnonzerodigit))
105    ## All digits are zero
106    isexact = true ();
107    if (hex.s)
108      binary = -0;
109    else
110      binary = 0;
111    endif
112    return
113  endif
114
115  ## Move point to the left, right after the first significand binary digit
116  if (length (hex.m) > 1)
117    if (hex.e + (length (hex.m) - 1) >= intmax (class (hex.e)))
118      error ("interval:InvalidOperand", ...
119             ["exponent overflow during normalization: " string ]);
120    endif
121    hex.e += length (hex.m) - 1;
122  endif
123
124  ## Overflow
125  if (hex.e > 1023)
126    isexact = false ();
127    if (hex.s && direction < 0)
128      binary = -inf;
129    elseif (hex.s)
130      binary = -realmax ();
131    elseif (direction > 0)
132      binary = inf;
133    else
134      binary = realmax ();
135    endif
136    return
137  endif
138
139  ## Underflow
140  if (hex.e < -1074)
141    isexact = false ();
142    if (hex.s && direction < 0)
143      binary = -pow2 (-1074);
144    elseif (hex.s)
145      binary = -0;
146    elseif (direction > 0)
147      binary = pow2 (-1074);
148    else
149      binary = 0;
150    endif
151    return
152  endif
153
154  if (hex.e < -1022)
155    ## Subnormal numbers
156    hex.m = [zeros(1, -1022 - hex.e), hex.m];
157    ieee754exponent = -1023;
158  else
159    ## Normal numbers
160    ieee754exponent = hex.e;
161  endif
162
163  ## Only the most significand 53 bits are relevant.
164  significand = postpad (hex.m, 53, 0, 2);
165  isexact = length (hex.m) <= length (significand) || ...
166            isempty (find (hex.m (54:end), 1));
167
168  ## The first bit can be omitted (normalization).
169  significand (1) = [];
170
171  ## The remaining 52 bits can be converted to 13 hex digits
172  ieee754mantissa = dec2hex (pow2 (51 : -1 : 0) * significand', 13);
173
174  if (not (hex.s))
175    ieee754signandexponent = ieee754exponent + 1023;
176  else
177    ieee754signandexponent = ieee754exponent + 1023 + pow2 (11);
178  endif
179
180  ieee754double = strcat (dec2hex (ieee754signandexponent, 3), ieee754mantissa);
181  binary = hex2num (ieee754double);
182
183  ## Last, consider the rounding direction
184  if (isexact || ...
185      (direction < 0 && not (hex.s)) ||
186      (direction > 0 && hex.s))
187    ## The number is exact or the truncation of digits above resulted in
188    ## correct rounding direction
189  else
190    delta = pow2 (-1074);
191    binary = mpfr_function_d ('plus', direction, binary, delta);
192  endif
193
194endfunction
195