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