1 /* Software floating-point emulation. 2 Definitions for IEEE Extended Precision. 3 Copyright (C) 1999-2014 Free Software Foundation, Inc. 4 This file is part of the GNU C Library. 5 Contributed by Jakub Jelinek (jj@ultra.linux.cz). 6 7 The GNU C Library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public 9 License as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 12 In addition to the permissions in the GNU Lesser General Public 13 License, the Free Software Foundation gives you unlimited 14 permission to link the compiled version of this file into 15 combinations with other programs, and to distribute those 16 combinations without any restriction coming from the use of this 17 file. (The Lesser General Public License restrictions do apply in 18 other respects; for example, they cover modification of the file, 19 and distribution when not linked into a combine executable.) 20 21 The GNU C Library is distributed in the hope that it will be useful, 22 but WITHOUT ANY WARRANTY; without even the implied warranty of 23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 24 Lesser General Public License for more details. 25 26 You should have received a copy of the GNU Lesser General Public 27 License along with the GNU C Library; if not, see 28 <http://www.gnu.org/licenses/>. */ 29 30 #if _FP_W_TYPE_SIZE < 32 31 # error "Here's a nickel, kid. Go buy yourself a real computer." 32 #endif 33 34 #if _FP_W_TYPE_SIZE < 64 35 # define _FP_FRACTBITS_E (4*_FP_W_TYPE_SIZE) 36 # define _FP_FRACTBITS_DW_E (8*_FP_W_TYPE_SIZE) 37 #else 38 # define _FP_FRACTBITS_E (2*_FP_W_TYPE_SIZE) 39 # define _FP_FRACTBITS_DW_E (4*_FP_W_TYPE_SIZE) 40 #endif 41 42 #define _FP_FRACBITS_E 64 43 #define _FP_FRACXBITS_E (_FP_FRACTBITS_E - _FP_FRACBITS_E) 44 #define _FP_WFRACBITS_E (_FP_WORKBITS + _FP_FRACBITS_E) 45 #define _FP_WFRACXBITS_E (_FP_FRACTBITS_E - _FP_WFRACBITS_E) 46 #define _FP_EXPBITS_E 15 47 #define _FP_EXPBIAS_E 16383 48 #define _FP_EXPMAX_E 32767 49 50 #define _FP_QNANBIT_E \ 51 ((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE) 52 #define _FP_QNANBIT_SH_E \ 53 ((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-2+_FP_WORKBITS) % _FP_W_TYPE_SIZE) 54 #define _FP_IMPLBIT_E \ 55 ((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE) 56 #define _FP_IMPLBIT_SH_E \ 57 ((_FP_W_TYPE) 1 << (_FP_FRACBITS_E-1+_FP_WORKBITS) % _FP_W_TYPE_SIZE) 58 #define _FP_OVERFLOW_E \ 59 ((_FP_W_TYPE) 1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE)) 60 61 #define _FP_WFRACBITS_DW_E (2 * _FP_WFRACBITS_E) 62 #define _FP_WFRACXBITS_DW_E (_FP_FRACTBITS_DW_E - _FP_WFRACBITS_DW_E) 63 #define _FP_HIGHBIT_DW_E \ 64 ((_FP_W_TYPE) 1 << (_FP_WFRACBITS_DW_E - 1) % _FP_W_TYPE_SIZE) 65 66 typedef float XFtype __attribute__ ((mode (XF))); 67 68 #if _FP_W_TYPE_SIZE < 64 69 70 union _FP_UNION_E 71 { 72 XFtype flt; 73 struct _FP_STRUCT_LAYOUT 74 { 75 # if __BYTE_ORDER == __BIG_ENDIAN 76 unsigned long pad1 : _FP_W_TYPE_SIZE; 77 unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E); 78 unsigned long sign : 1; 79 unsigned long exp : _FP_EXPBITS_E; 80 unsigned long frac1 : _FP_W_TYPE_SIZE; 81 unsigned long frac0 : _FP_W_TYPE_SIZE; 82 # else 83 unsigned long frac0 : _FP_W_TYPE_SIZE; 84 unsigned long frac1 : _FP_W_TYPE_SIZE; 85 unsigned exp : _FP_EXPBITS_E; 86 unsigned sign : 1; 87 # endif /* not bigendian */ 88 } bits __attribute__ ((packed)); 89 }; 90 91 92 # define FP_DECL_E(X) _FP_DECL (4, X) 93 94 # define FP_UNPACK_RAW_E(X, val) \ 95 do \ 96 { \ 97 union _FP_UNION_E FP_UNPACK_RAW_E_flo; \ 98 FP_UNPACK_RAW_E_flo.flt = (val); \ 99 \ 100 X##_f[2] = 0; \ 101 X##_f[3] = 0; \ 102 X##_f[0] = FP_UNPACK_RAW_E_flo.bits.frac0; \ 103 X##_f[1] = FP_UNPACK_RAW_E_flo.bits.frac1; \ 104 X##_e = FP_UNPACK_RAW_E_flo.bits.exp; \ 105 X##_s = FP_UNPACK_RAW_E_flo.bits.sign; \ 106 } \ 107 while (0) 108 109 # define FP_UNPACK_RAW_EP(X, val) \ 110 do \ 111 { \ 112 union _FP_UNION_E *FP_UNPACK_RAW_EP_flo \ 113 = (union _FP_UNION_E *) (val); \ 114 \ 115 X##_f[2] = 0; \ 116 X##_f[3] = 0; \ 117 X##_f[0] = FP_UNPACK_RAW_EP_flo->bits.frac0; \ 118 X##_f[1] = FP_UNPACK_RAW_EP_flo->bits.frac1; \ 119 X##_e = FP_UNPACK_RAW_EP_flo->bits.exp; \ 120 X##_s = FP_UNPACK_RAW_EP_flo->bits.sign; \ 121 } \ 122 while (0) 123 124 # define FP_PACK_RAW_E(val, X) \ 125 do \ 126 { \ 127 union _FP_UNION_E FP_PACK_RAW_E_flo; \ 128 \ 129 if (X##_e) \ 130 X##_f[1] |= _FP_IMPLBIT_E; \ 131 else \ 132 X##_f[1] &= ~(_FP_IMPLBIT_E); \ 133 FP_PACK_RAW_E_flo.bits.frac0 = X##_f[0]; \ 134 FP_PACK_RAW_E_flo.bits.frac1 = X##_f[1]; \ 135 FP_PACK_RAW_E_flo.bits.exp = X##_e; \ 136 FP_PACK_RAW_E_flo.bits.sign = X##_s; \ 137 \ 138 (val) = FP_PACK_RAW_E_flo.flt; \ 139 } \ 140 while (0) 141 142 # define FP_PACK_RAW_EP(val, X) \ 143 do \ 144 { \ 145 if (!FP_INHIBIT_RESULTS) \ 146 { \ 147 union _FP_UNION_E *FP_PACK_RAW_EP_flo \ 148 = (union _FP_UNION_E *) (val); \ 149 \ 150 if (X##_e) \ 151 X##_f[1] |= _FP_IMPLBIT_E; \ 152 else \ 153 X##_f[1] &= ~(_FP_IMPLBIT_E); \ 154 FP_PACK_RAW_EP_flo->bits.frac0 = X##_f[0]; \ 155 FP_PACK_RAW_EP_flo->bits.frac1 = X##_f[1]; \ 156 FP_PACK_RAW_EP_flo->bits.exp = X##_e; \ 157 FP_PACK_RAW_EP_flo->bits.sign = X##_s; \ 158 } \ 159 } \ 160 while (0) 161 162 # define FP_UNPACK_E(X, val) \ 163 do \ 164 { \ 165 FP_UNPACK_RAW_E (X, (val)); \ 166 _FP_UNPACK_CANONICAL (E, 4, X); \ 167 } \ 168 while (0) 169 170 # define FP_UNPACK_EP(X, val) \ 171 do \ 172 { \ 173 FP_UNPACK_RAW_EP (X, (val)); \ 174 _FP_UNPACK_CANONICAL (E, 4, X); \ 175 } \ 176 while (0) 177 178 # define FP_UNPACK_SEMIRAW_E(X, val) \ 179 do \ 180 { \ 181 FP_UNPACK_RAW_E (X, (val)); \ 182 _FP_UNPACK_SEMIRAW (E, 4, X); \ 183 } \ 184 while (0) 185 186 # define FP_UNPACK_SEMIRAW_EP(X, val) \ 187 do \ 188 { \ 189 FP_UNPACK_RAW_EP (X, (val)); \ 190 _FP_UNPACK_SEMIRAW (E, 4, X); \ 191 } \ 192 while (0) 193 194 # define FP_PACK_E(val, X) \ 195 do \ 196 { \ 197 _FP_PACK_CANONICAL (E, 4, X); \ 198 FP_PACK_RAW_E ((val), X); \ 199 } \ 200 while (0) 201 202 # define FP_PACK_EP(val, X) \ 203 do \ 204 { \ 205 _FP_PACK_CANONICAL (E, 4, X); \ 206 FP_PACK_RAW_EP ((val), X); \ 207 } \ 208 while (0) 209 210 # define FP_PACK_SEMIRAW_E(val, X) \ 211 do \ 212 { \ 213 _FP_PACK_SEMIRAW (E, 4, X); \ 214 FP_PACK_RAW_E ((val), X); \ 215 } \ 216 while (0) 217 218 # define FP_PACK_SEMIRAW_EP(val, X) \ 219 do \ 220 { \ 221 _FP_PACK_SEMIRAW (E, 4, X); \ 222 FP_PACK_RAW_EP ((val), X); \ 223 } \ 224 while (0) 225 226 # define FP_ISSIGNAN_E(X) _FP_ISSIGNAN (E, 4, X) 227 # define FP_NEG_E(R, X) _FP_NEG (E, 4, R, X) 228 # define FP_ADD_E(R, X, Y) _FP_ADD (E, 4, R, X, Y) 229 # define FP_SUB_E(R, X, Y) _FP_SUB (E, 4, R, X, Y) 230 # define FP_MUL_E(R, X, Y) _FP_MUL (E, 4, R, X, Y) 231 # define FP_DIV_E(R, X, Y) _FP_DIV (E, 4, R, X, Y) 232 # define FP_SQRT_E(R, X) _FP_SQRT (E, 4, R, X) 233 # define FP_FMA_E(R, X, Y, Z) _FP_FMA (E, 4, 8, R, X, Y, Z) 234 235 /* Square root algorithms: 236 We have just one right now, maybe Newton approximation 237 should be added for those machines where division is fast. 238 This has special _E version because standard _4 square 239 root would not work (it has to start normally with the 240 second word and not the first), but as we have to do it 241 anyway, we optimize it by doing most of the calculations 242 in two UWtype registers instead of four. */ 243 244 # define _FP_SQRT_MEAT_E(R, S, T, X, q) \ 245 do \ 246 { \ 247 (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1); \ 248 _FP_FRAC_SRL_4 (X, (_FP_WORKBITS)); \ 249 while (q) \ 250 { \ 251 T##_f[1] = S##_f[1] + (q); \ 252 if (T##_f[1] <= X##_f[1]) \ 253 { \ 254 S##_f[1] = T##_f[1] + (q); \ 255 X##_f[1] -= T##_f[1]; \ 256 R##_f[1] += (q); \ 257 } \ 258 _FP_FRAC_SLL_2 (X, 1); \ 259 (q) >>= 1; \ 260 } \ 261 (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1); \ 262 while (q) \ 263 { \ 264 T##_f[0] = S##_f[0] + (q); \ 265 T##_f[1] = S##_f[1]; \ 266 if (T##_f[1] < X##_f[1] \ 267 || (T##_f[1] == X##_f[1] \ 268 && T##_f[0] <= X##_f[0])) \ 269 { \ 270 S##_f[0] = T##_f[0] + (q); \ 271 S##_f[1] += (T##_f[0] > S##_f[0]); \ 272 _FP_FRAC_DEC_2 (X, T); \ 273 R##_f[0] += (q); \ 274 } \ 275 _FP_FRAC_SLL_2 (X, 1); \ 276 (q) >>= 1; \ 277 } \ 278 _FP_FRAC_SLL_4 (R, (_FP_WORKBITS)); \ 279 if (X##_f[0] | X##_f[1]) \ 280 { \ 281 if (S##_f[1] < X##_f[1] \ 282 || (S##_f[1] == X##_f[1] \ 283 && S##_f[0] < X##_f[0])) \ 284 R##_f[0] |= _FP_WORK_ROUND; \ 285 R##_f[0] |= _FP_WORK_STICKY; \ 286 } \ 287 } \ 288 while (0) 289 290 # define FP_CMP_E(r, X, Y, un, ex) _FP_CMP (E, 4, (r), X, Y, (un), (ex)) 291 # define FP_CMP_EQ_E(r, X, Y, ex) _FP_CMP_EQ (E, 4, (r), X, Y, (ex)) 292 # define FP_CMP_UNORD_E(r, X, Y, ex) _FP_CMP_UNORD (E, 4, (r), X, Y, (ex)) 293 294 # define FP_TO_INT_E(r, X, rsz, rsg) _FP_TO_INT (E, 4, (r), X, (rsz), (rsg)) 295 # define FP_FROM_INT_E(X, r, rs, rt) _FP_FROM_INT (E, 4, X, (r), (rs), rt) 296 297 # define _FP_FRAC_HIGH_E(X) (X##_f[2]) 298 # define _FP_FRAC_HIGH_RAW_E(X) (X##_f[1]) 299 300 # define _FP_FRAC_HIGH_DW_E(X) (X##_f[4]) 301 302 #else /* not _FP_W_TYPE_SIZE < 64 */ 303 union _FP_UNION_E 304 { 305 XFtype flt; 306 struct _FP_STRUCT_LAYOUT 307 { 308 # if __BYTE_ORDER == __BIG_ENDIAN 309 _FP_W_TYPE pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E); 310 unsigned sign : 1; 311 unsigned exp : _FP_EXPBITS_E; 312 _FP_W_TYPE frac : _FP_W_TYPE_SIZE; 313 # else 314 _FP_W_TYPE frac : _FP_W_TYPE_SIZE; 315 unsigned exp : _FP_EXPBITS_E; 316 unsigned sign : 1; 317 # endif 318 } bits; 319 }; 320 321 # define FP_DECL_E(X) _FP_DECL (2, X) 322 323 # define FP_UNPACK_RAW_E(X, val) \ 324 do \ 325 { \ 326 union _FP_UNION_E FP_UNPACK_RAW_E_flo; \ 327 FP_UNPACK_RAW_E_flo.flt = (val); \ 328 \ 329 X##_f0 = FP_UNPACK_RAW_E_flo.bits.frac; \ 330 X##_f1 = 0; \ 331 X##_e = FP_UNPACK_RAW_E_flo.bits.exp; \ 332 X##_s = FP_UNPACK_RAW_E_flo.bits.sign; \ 333 } \ 334 while (0) 335 336 # define FP_UNPACK_RAW_EP(X, val) \ 337 do \ 338 { \ 339 union _FP_UNION_E *FP_UNPACK_RAW_EP_flo \ 340 = (union _FP_UNION_E *) (val); \ 341 \ 342 X##_f0 = FP_UNPACK_RAW_EP_flo->bits.frac; \ 343 X##_f1 = 0; \ 344 X##_e = FP_UNPACK_RAW_EP_flo->bits.exp; \ 345 X##_s = FP_UNPACK_RAW_EP_flo->bits.sign; \ 346 } \ 347 while (0) 348 349 # define FP_PACK_RAW_E(val, X) \ 350 do \ 351 { \ 352 union _FP_UNION_E FP_PACK_RAW_E_flo; \ 353 \ 354 if (X##_e) \ 355 X##_f0 |= _FP_IMPLBIT_E; \ 356 else \ 357 X##_f0 &= ~(_FP_IMPLBIT_E); \ 358 FP_PACK_RAW_E_flo.bits.frac = X##_f0; \ 359 FP_PACK_RAW_E_flo.bits.exp = X##_e; \ 360 FP_PACK_RAW_E_flo.bits.sign = X##_s; \ 361 \ 362 (val) = FP_PACK_RAW_E_flo.flt; \ 363 } \ 364 while (0) 365 366 # define FP_PACK_RAW_EP(fs, val, X) \ 367 do \ 368 { \ 369 if (!FP_INHIBIT_RESULTS) \ 370 { \ 371 union _FP_UNION_E *FP_PACK_RAW_EP_flo \ 372 = (union _FP_UNION_E *) (val); \ 373 \ 374 if (X##_e) \ 375 X##_f0 |= _FP_IMPLBIT_E; \ 376 else \ 377 X##_f0 &= ~(_FP_IMPLBIT_E); \ 378 FP_PACK_RAW_EP_flo->bits.frac = X##_f0; \ 379 FP_PACK_RAW_EP_flo->bits.exp = X##_e; \ 380 FP_PACK_RAW_EP_flo->bits.sign = X##_s; \ 381 } \ 382 } \ 383 while (0) 384 385 386 # define FP_UNPACK_E(X, val) \ 387 do \ 388 { \ 389 FP_UNPACK_RAW_E (X, (val)); \ 390 _FP_UNPACK_CANONICAL (E, 2, X); \ 391 } \ 392 while (0) 393 394 # define FP_UNPACK_EP(X, val) \ 395 do \ 396 { \ 397 FP_UNPACK_RAW_EP (X, (val)); \ 398 _FP_UNPACK_CANONICAL (E, 2, X); \ 399 } \ 400 while (0) 401 402 # define FP_UNPACK_SEMIRAW_E(X, val) \ 403 do \ 404 { \ 405 FP_UNPACK_RAW_E (X, (val)); \ 406 _FP_UNPACK_SEMIRAW (E, 2, X); \ 407 } \ 408 while (0) 409 410 # define FP_UNPACK_SEMIRAW_EP(X, val) \ 411 do \ 412 { \ 413 FP_UNPACK_RAW_EP (X, (val)); \ 414 _FP_UNPACK_SEMIRAW (E, 2, X); \ 415 } \ 416 while (0) 417 418 # define FP_PACK_E(val, X) \ 419 do \ 420 { \ 421 _FP_PACK_CANONICAL (E, 2, X); \ 422 FP_PACK_RAW_E ((val), X); \ 423 } \ 424 while (0) 425 426 # define FP_PACK_EP(val, X) \ 427 do \ 428 { \ 429 _FP_PACK_CANONICAL (E, 2, X); \ 430 FP_PACK_RAW_EP ((val), X); \ 431 } \ 432 while (0) 433 434 # define FP_PACK_SEMIRAW_E(val, X) \ 435 do \ 436 { \ 437 _FP_PACK_SEMIRAW (E, 2, X); \ 438 FP_PACK_RAW_E ((val), X); \ 439 } \ 440 while (0) 441 442 # define FP_PACK_SEMIRAW_EP(val, X) \ 443 do \ 444 { \ 445 _FP_PACK_SEMIRAW (E, 2, X); \ 446 FP_PACK_RAW_EP ((val), X); \ 447 } \ 448 while (0) 449 450 # define FP_ISSIGNAN_E(X) _FP_ISSIGNAN (E, 2, X) 451 # define FP_NEG_E(R, X) _FP_NEG (E, 2, R, X) 452 # define FP_ADD_E(R, X, Y) _FP_ADD (E, 2, R, X, Y) 453 # define FP_SUB_E(R, X, Y) _FP_SUB (E, 2, R, X, Y) 454 # define FP_MUL_E(R, X, Y) _FP_MUL (E, 2, R, X, Y) 455 # define FP_DIV_E(R, X, Y) _FP_DIV (E, 2, R, X, Y) 456 # define FP_SQRT_E(R, X) _FP_SQRT (E, 2, R, X) 457 # define FP_FMA_E(R, X, Y, Z) _FP_FMA (E, 2, 4, R, X, Y, Z) 458 459 /* Square root algorithms: 460 We have just one right now, maybe Newton approximation 461 should be added for those machines where division is fast. 462 We optimize it by doing most of the calculations 463 in one UWtype registers instead of two, although we don't 464 have to. */ 465 # define _FP_SQRT_MEAT_E(R, S, T, X, q) \ 466 do \ 467 { \ 468 (q) = (_FP_W_TYPE) 1 << (_FP_W_TYPE_SIZE - 1); \ 469 _FP_FRAC_SRL_2 (X, (_FP_WORKBITS)); \ 470 while (q) \ 471 { \ 472 T##_f0 = S##_f0 + (q); \ 473 if (T##_f0 <= X##_f0) \ 474 { \ 475 S##_f0 = T##_f0 + (q); \ 476 X##_f0 -= T##_f0; \ 477 R##_f0 += (q); \ 478 } \ 479 _FP_FRAC_SLL_1 (X, 1); \ 480 (q) >>= 1; \ 481 } \ 482 _FP_FRAC_SLL_2 (R, (_FP_WORKBITS)); \ 483 if (X##_f0) \ 484 { \ 485 if (S##_f0 < X##_f0) \ 486 R##_f0 |= _FP_WORK_ROUND; \ 487 R##_f0 |= _FP_WORK_STICKY; \ 488 } \ 489 } \ 490 while (0) 491 492 # define FP_CMP_E(r, X, Y, un, ex) _FP_CMP (E, 2, (r), X, Y, (un), (ex)) 493 # define FP_CMP_EQ_E(r, X, Y, ex) _FP_CMP_EQ (E, 2, (r), X, Y, (ex)) 494 # define FP_CMP_UNORD_E(r, X, Y, ex) _FP_CMP_UNORD (E, 2, (r), X, Y, (ex)) 495 496 # define FP_TO_INT_E(r, X, rsz, rsg) _FP_TO_INT (E, 2, (r), X, (rsz), (rsg)) 497 # define FP_FROM_INT_E(X, r, rs, rt) _FP_FROM_INT (E, 2, X, (r), (rs), rt) 498 499 # define _FP_FRAC_HIGH_E(X) (X##_f1) 500 # define _FP_FRAC_HIGH_RAW_E(X) (X##_f0) 501 502 # define _FP_FRAC_HIGH_DW_E(X) (X##_f[2]) 503 504 #endif /* not _FP_W_TYPE_SIZE < 64 */ 505