1/* Software floating-point emulation. 2 Definitions for IEEE Extended Precision. 3 Copyright (C) 1999 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 Library General Public License as 9 published by the Free Software Foundation; either version 2 of the 10 License, or (at your option) any later version. 11 12 The GNU C Library is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 Library General Public License for more details. 16 17 You should have received a copy of the GNU Library General Public 18 License along with the GNU C Library; see the file COPYING.LIB. If 19 not, write to the Free Software Foundation, Inc., 20 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ 21 22 23#ifndef __MATH_EMU_EXTENDED_H__ 24#define __MATH_EMU_EXTENDED_H__ 25 26#if _FP_W_TYPE_SIZE < 32 27#error "Here's a nickel, kid. Go buy yourself a real computer." 28#endif 29 30#if _FP_W_TYPE_SIZE < 64 31#define _FP_FRACTBITS_E (4*_FP_W_TYPE_SIZE) 32#else 33#define _FP_FRACTBITS_E (2*_FP_W_TYPE_SIZE) 34#endif 35 36#define _FP_FRACBITS_E 64 37#define _FP_FRACXBITS_E (_FP_FRACTBITS_E - _FP_FRACBITS_E) 38#define _FP_WFRACBITS_E (_FP_WORKBITS + _FP_FRACBITS_E) 39#define _FP_WFRACXBITS_E (_FP_FRACTBITS_E - _FP_WFRACBITS_E) 40#define _FP_EXPBITS_E 15 41#define _FP_EXPBIAS_E 16383 42#define _FP_EXPMAX_E 32767 43 44#define _FP_QNANBIT_E \ 45 ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE) 46#define _FP_IMPLBIT_E \ 47 ((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE) 48#define _FP_OVERFLOW_E \ 49 ((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE)) 50 51#if _FP_W_TYPE_SIZE < 64 52 53union _FP_UNION_E 54{ 55 long double flt; 56 struct 57 { 58#if __BYTE_ORDER == __BIG_ENDIAN 59 unsigned long pad1 : _FP_W_TYPE_SIZE; 60 unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E); 61 unsigned long sign : 1; 62 unsigned long exp : _FP_EXPBITS_E; 63 unsigned long frac1 : _FP_W_TYPE_SIZE; 64 unsigned long frac0 : _FP_W_TYPE_SIZE; 65#else 66 unsigned long frac0 : _FP_W_TYPE_SIZE; 67 unsigned long frac1 : _FP_W_TYPE_SIZE; 68 unsigned exp : _FP_EXPBITS_E; 69 unsigned sign : 1; 70#endif /* not bigendian */ 71 } bits __attribute__((packed)); 72}; 73 74 75#define FP_DECL_E(X) _FP_DECL(4,X) 76 77#define FP_UNPACK_RAW_E(X, val) \ 78 do { \ 79 union _FP_UNION_E _flo; _flo.flt = (val); \ 80 \ 81 X##_f[2] = 0; X##_f[3] = 0; \ 82 X##_f[0] = _flo.bits.frac0; \ 83 X##_f[1] = _flo.bits.frac1; \ 84 X##_e = _flo.bits.exp; \ 85 X##_s = _flo.bits.sign; \ 86 if (!X##_e && (X##_f[1] || X##_f[0]) \ 87 && !(X##_f[1] & _FP_IMPLBIT_E)) \ 88 { \ 89 X##_e++; \ 90 FP_SET_EXCEPTION(FP_EX_DENORM); \ 91 } \ 92 } while (0) 93 94#define FP_UNPACK_RAW_EP(X, val) \ 95 do { \ 96 union _FP_UNION_E *_flo = \ 97 (union _FP_UNION_E *)(val); \ 98 \ 99 X##_f[2] = 0; X##_f[3] = 0; \ 100 X##_f[0] = _flo->bits.frac0; \ 101 X##_f[1] = _flo->bits.frac1; \ 102 X##_e = _flo->bits.exp; \ 103 X##_s = _flo->bits.sign; \ 104 if (!X##_e && (X##_f[1] || X##_f[0]) \ 105 && !(X##_f[1] & _FP_IMPLBIT_E)) \ 106 { \ 107 X##_e++; \ 108 FP_SET_EXCEPTION(FP_EX_DENORM); \ 109 } \ 110 } while (0) 111 112#define FP_PACK_RAW_E(val, X) \ 113 do { \ 114 union _FP_UNION_E _flo; \ 115 \ 116 if (X##_e) X##_f[1] |= _FP_IMPLBIT_E; \ 117 else X##_f[1] &= ~(_FP_IMPLBIT_E); \ 118 _flo.bits.frac0 = X##_f[0]; \ 119 _flo.bits.frac1 = X##_f[1]; \ 120 _flo.bits.exp = X##_e; \ 121 _flo.bits.sign = X##_s; \ 122 \ 123 (val) = _flo.flt; \ 124 } while (0) 125 126#define FP_PACK_RAW_EP(val, X) \ 127 do { \ 128 if (!FP_INHIBIT_RESULTS) \ 129 { \ 130 union _FP_UNION_E *_flo = \ 131 (union _FP_UNION_E *)(val); \ 132 \ 133 if (X##_e) X##_f[1] |= _FP_IMPLBIT_E; \ 134 else X##_f[1] &= ~(_FP_IMPLBIT_E); \ 135 _flo->bits.frac0 = X##_f[0]; \ 136 _flo->bits.frac1 = X##_f[1]; \ 137 _flo->bits.exp = X##_e; \ 138 _flo->bits.sign = X##_s; \ 139 } \ 140 } while (0) 141 142#define FP_UNPACK_E(X,val) \ 143 do { \ 144 FP_UNPACK_RAW_E(X,val); \ 145 _FP_UNPACK_CANONICAL(E,4,X); \ 146 } while (0) 147 148#define FP_UNPACK_EP(X,val) \ 149 do { \ 150 FP_UNPACK_RAW_2_P(X,val); \ 151 _FP_UNPACK_CANONICAL(E,4,X); \ 152 } while (0) 153 154#define FP_PACK_E(val,X) \ 155 do { \ 156 _FP_PACK_CANONICAL(E,4,X); \ 157 FP_PACK_RAW_E(val,X); \ 158 } while (0) 159 160#define FP_PACK_EP(val,X) \ 161 do { \ 162 _FP_PACK_CANONICAL(E,4,X); \ 163 FP_PACK_RAW_EP(val,X); \ 164 } while (0) 165 166#define FP_ISSIGNAN_E(X) _FP_ISSIGNAN(E,4,X) 167#define FP_NEG_E(R,X) _FP_NEG(E,4,R,X) 168#define FP_ADD_E(R,X,Y) _FP_ADD(E,4,R,X,Y) 169#define FP_SUB_E(R,X,Y) _FP_SUB(E,4,R,X,Y) 170#define FP_MUL_E(R,X,Y) _FP_MUL(E,4,R,X,Y) 171#define FP_DIV_E(R,X,Y) _FP_DIV(E,4,R,X,Y) 172#define FP_SQRT_E(R,X) _FP_SQRT(E,4,R,X) 173 174/* 175 * Square root algorithms: 176 * We have just one right now, maybe Newton approximation 177 * should be added for those machines where division is fast. 178 * This has special _E version because standard _4 square 179 * root would not work (it has to start normally with the 180 * second word and not the first), but as we have to do it 181 * anyway, we optimize it by doing most of the calculations 182 * in two UWtype registers instead of four. 183 */ 184 185#define _FP_SQRT_MEAT_E(R, S, T, X, q) \ 186 do { \ 187 q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \ 188 _FP_FRAC_SRL_4(X, (_FP_WORKBITS)); \ 189 while (q) \ 190 { \ 191 T##_f[1] = S##_f[1] + q; \ 192 if (T##_f[1] <= X##_f[1]) \ 193 { \ 194 S##_f[1] = T##_f[1] + q; \ 195 X##_f[1] -= T##_f[1]; \ 196 R##_f[1] += q; \ 197 } \ 198 _FP_FRAC_SLL_2(X, 1); \ 199 q >>= 1; \ 200 } \ 201 q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \ 202 while (q) \ 203 { \ 204 T##_f[0] = S##_f[0] + q; \ 205 T##_f[1] = S##_f[1]; \ 206 if (T##_f[1] < X##_f[1] || \ 207 (T##_f[1] == X##_f[1] && \ 208 T##_f[0] <= X##_f[0])) \ 209 { \ 210 S##_f[0] = T##_f[0] + q; \ 211 S##_f[1] += (T##_f[0] > S##_f[0]); \ 212 _FP_FRAC_DEC_2(X, T); \ 213 R##_f[0] += q; \ 214 } \ 215 _FP_FRAC_SLL_2(X, 1); \ 216 q >>= 1; \ 217 } \ 218 _FP_FRAC_SLL_4(R, (_FP_WORKBITS)); \ 219 if (X##_f[0] | X##_f[1]) \ 220 { \ 221 if (S##_f[1] < X##_f[1] || \ 222 (S##_f[1] == X##_f[1] && \ 223 S##_f[0] < X##_f[0])) \ 224 R##_f[0] |= _FP_WORK_ROUND; \ 225 R##_f[0] |= _FP_WORK_STICKY; \ 226 } \ 227 } while (0) 228 229#define FP_CMP_E(r,X,Y,un) _FP_CMP(E,4,r,X,Y,un) 230#define FP_CMP_EQ_E(r,X,Y) _FP_CMP_EQ(E,4,r,X,Y) 231 232#define FP_TO_INT_E(r,X,rsz,rsg) _FP_TO_INT(E,4,r,X,rsz,rsg) 233#define FP_TO_INT_ROUND_E(r,X,rsz,rsg) _FP_TO_INT_ROUND(E,4,r,X,rsz,rsg) 234#define FP_FROM_INT_E(X,r,rs,rt) _FP_FROM_INT(E,4,X,r,rs,rt) 235 236#define _FP_FRAC_HIGH_E(X) (X##_f[2]) 237#define _FP_FRAC_HIGH_RAW_E(X) (X##_f[1]) 238 239#else /* not _FP_W_TYPE_SIZE < 64 */ 240union _FP_UNION_E 241{ 242 long double flt /* __attribute__((mode(TF))) */ ; 243 struct { 244#if __BYTE_ORDER == __BIG_ENDIAN 245 unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E); 246 unsigned sign : 1; 247 unsigned exp : _FP_EXPBITS_E; 248 unsigned long frac : _FP_W_TYPE_SIZE; 249#else 250 unsigned long frac : _FP_W_TYPE_SIZE; 251 unsigned exp : _FP_EXPBITS_E; 252 unsigned sign : 1; 253#endif 254 } bits; 255}; 256 257#define FP_DECL_E(X) _FP_DECL(2,X) 258 259#define FP_UNPACK_RAW_E(X, val) \ 260 do { \ 261 union _FP_UNION_E _flo; _flo.flt = (val); \ 262 \ 263 X##_f0 = _flo.bits.frac; \ 264 X##_f1 = 0; \ 265 X##_e = _flo.bits.exp; \ 266 X##_s = _flo.bits.sign; \ 267 if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E)) \ 268 { \ 269 X##_e++; \ 270 FP_SET_EXCEPTION(FP_EX_DENORM); \ 271 } \ 272 } while (0) 273 274#define FP_UNPACK_RAW_EP(X, val) \ 275 do { \ 276 union _FP_UNION_E *_flo = \ 277 (union _FP_UNION_E *)(val); \ 278 \ 279 X##_f0 = _flo->bits.frac; \ 280 X##_f1 = 0; \ 281 X##_e = _flo->bits.exp; \ 282 X##_s = _flo->bits.sign; \ 283 if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E)) \ 284 { \ 285 X##_e++; \ 286 FP_SET_EXCEPTION(FP_EX_DENORM); \ 287 } \ 288 } while (0) 289 290#define FP_PACK_RAW_E(val, X) \ 291 do { \ 292 union _FP_UNION_E _flo; \ 293 \ 294 if (X##_e) X##_f0 |= _FP_IMPLBIT_E; \ 295 else X##_f0 &= ~(_FP_IMPLBIT_E); \ 296 _flo.bits.frac = X##_f0; \ 297 _flo.bits.exp = X##_e; \ 298 _flo.bits.sign = X##_s; \ 299 \ 300 (val) = _flo.flt; \ 301 } while (0) 302 303#define FP_PACK_RAW_EP(fs, val, X) \ 304 do { \ 305 if (!FP_INHIBIT_RESULTS) \ 306 { \ 307 union _FP_UNION_E *_flo = \ 308 (union _FP_UNION_E *)(val); \ 309 \ 310 if (X##_e) X##_f0 |= _FP_IMPLBIT_E; \ 311 else X##_f0 &= ~(_FP_IMPLBIT_E); \ 312 _flo->bits.frac = X##_f0; \ 313 _flo->bits.exp = X##_e; \ 314 _flo->bits.sign = X##_s; \ 315 } \ 316 } while (0) 317 318 319#define FP_UNPACK_E(X,val) \ 320 do { \ 321 FP_UNPACK_RAW_E(X,val); \ 322 _FP_UNPACK_CANONICAL(E,2,X); \ 323 } while (0) 324 325#define FP_UNPACK_EP(X,val) \ 326 do { \ 327 FP_UNPACK_RAW_EP(X,val); \ 328 _FP_UNPACK_CANONICAL(E,2,X); \ 329 } while (0) 330 331#define FP_PACK_E(val,X) \ 332 do { \ 333 _FP_PACK_CANONICAL(E,2,X); \ 334 FP_PACK_RAW_E(val,X); \ 335 } while (0) 336 337#define FP_PACK_EP(val,X) \ 338 do { \ 339 _FP_PACK_CANONICAL(E,2,X); \ 340 FP_PACK_RAW_EP(val,X); \ 341 } while (0) 342 343#define FP_ISSIGNAN_E(X) _FP_ISSIGNAN(E,2,X) 344#define FP_NEG_E(R,X) _FP_NEG(E,2,R,X) 345#define FP_ADD_E(R,X,Y) _FP_ADD(E,2,R,X,Y) 346#define FP_SUB_E(R,X,Y) _FP_SUB(E,2,R,X,Y) 347#define FP_MUL_E(R,X,Y) _FP_MUL(E,2,R,X,Y) 348#define FP_DIV_E(R,X,Y) _FP_DIV(E,2,R,X,Y) 349#define FP_SQRT_E(R,X) _FP_SQRT(E,2,R,X) 350 351/* 352 * Square root algorithms: 353 * We have just one right now, maybe Newton approximation 354 * should be added for those machines where division is fast. 355 * We optimize it by doing most of the calculations 356 * in one UWtype registers instead of two, although we don't 357 * have to. 358 */ 359#define _FP_SQRT_MEAT_E(R, S, T, X, q) \ 360 do { \ 361 q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \ 362 _FP_FRAC_SRL_2(X, (_FP_WORKBITS)); \ 363 while (q) \ 364 { \ 365 T##_f0 = S##_f0 + q; \ 366 if (T##_f0 <= X##_f0) \ 367 { \ 368 S##_f0 = T##_f0 + q; \ 369 X##_f0 -= T##_f0; \ 370 R##_f0 += q; \ 371 } \ 372 _FP_FRAC_SLL_1(X, 1); \ 373 q >>= 1; \ 374 } \ 375 _FP_FRAC_SLL_2(R, (_FP_WORKBITS)); \ 376 if (X##_f0) \ 377 { \ 378 if (S##_f0 < X##_f0) \ 379 R##_f0 |= _FP_WORK_ROUND; \ 380 R##_f0 |= _FP_WORK_STICKY; \ 381 } \ 382 } while (0) 383 384#define FP_CMP_E(r,X,Y,un) _FP_CMP(E,2,r,X,Y,un) 385#define FP_CMP_EQ_E(r,X,Y) _FP_CMP_EQ(E,2,r,X,Y) 386 387#define FP_TO_INT_E(r,X,rsz,rsg) _FP_TO_INT(E,2,r,X,rsz,rsg) 388#define FP_TO_INT_ROUND_E(r,X,rsz,rsg) _FP_TO_INT_ROUND(E,2,r,X,rsz,rsg) 389#define FP_FROM_INT_E(X,r,rs,rt) _FP_FROM_INT(E,2,X,r,rs,rt) 390 391#define _FP_FRAC_HIGH_E(X) (X##_f1) 392#define _FP_FRAC_HIGH_RAW_E(X) (X##_f0) 393 394#endif /* not _FP_W_TYPE_SIZE < 64 */ 395 396#endif /* __MATH_EMU_EXTENDED_H__ */ 397