soft-fp/extended.h

   1 /* Software floating-point emulation.
   2    Definitions for IEEE Extended Precision.
   3    Copyright (C) 1999-2013 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     union _FP_UNION_E _flo; _flo.flt = (val);           \
  97                                                         \
  98     X##_f[2] = 0; X##_f[3] = 0;                         \
  99     X##_f[0] = _flo.bits.frac0;                         \
 100     X##_f[1] = _flo.bits.frac1;                         \
 101     X##_e  = _flo.bits.exp;                             \
 102     X##_s  = _flo.bits.sign;                            \
 103   } while (0)
 104
 105 #define FP_UNPACK_RAW_EP(X, val)                        \
 106   do {                                                  \
 107     union _FP_UNION_E *_flo =                           \
 108     (union _FP_UNION_E *)(val);                         \
 109                                                         \
 110     X##_f[2] = 0; X##_f[3] = 0;                         \
 111     X##_f[0] = _flo->bits.frac0;                        \
 112     X##_f[1] = _flo->bits.frac1;                        \
 113     X##_e  = _flo->bits.exp;                            \
 114     X##_s  = _flo->bits.sign;                           \
 115   } while (0)
 116
 117 #define FP_PACK_RAW_E(val, X)                           \
 118   do {                                                  \
 119     union _FP_UNION_E _flo;                             \
 120                                                         \
 121     if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;               \
 122     else X##_f[1] &= ~(_FP_IMPLBIT_E);                  \
 123     _flo.bits.frac0 = X##_f[0];                         \
 124     _flo.bits.frac1 = X##_f[1];                         \
 125     _flo.bits.exp   = X##_e;                            \
 126     _flo.bits.sign  = X##_s;                            \
 127                                                         \
 128     (val) = _flo.flt;                                   \
 129   } while (0)
 130
 131 #define FP_PACK_RAW_EP(val, X)                          \
 132   do {                                                  \
 133     if (!FP_INHIBIT_RESULTS)                            \
 134       {                                                 \
 135         union _FP_UNION_E *_flo =                       \
 136           (union _FP_UNION_E *)(val);                   \
 137                                                         \
 138         if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;           \
 139         else X##_f[1] &= ~(_FP_IMPLBIT_E);              \
 140         _flo->bits.frac0 = X##_f[0];                    \
 141         _flo->bits.frac1 = X##_f[1];                    \
 142         _flo->bits.exp   = X##_e;                       \
 143         _flo->bits.sign  = X##_s;                       \
 144       }                                                 \
 145   } while (0)
 146
 147 #define FP_UNPACK_E(X,val)              \
 148   do {                                  \
 149     FP_UNPACK_RAW_E(X,val);             \
 150     _FP_UNPACK_CANONICAL(E,4,X);        \
 151   } while (0)
 152
 153 #define FP_UNPACK_EP(X,val)             \
 154   do {                                  \
 155     FP_UNPACK_RAW_EP(X,val);            \
 156     _FP_UNPACK_CANONICAL(E,4,X);        \
 157   } while (0)
 158
 159 #define FP_UNPACK_SEMIRAW_E(X,val)      \
 160   do {                                  \
 161     FP_UNPACK_RAW_E(X,val);             \
 162     _FP_UNPACK_SEMIRAW(E,4,X);          \
 163   } while (0)
 164
 165 #define FP_UNPACK_SEMIRAW_EP(X,val)     \
 166   do {                                  \
 167     FP_UNPACK_RAW_EP(X,val);            \
 168     _FP_UNPACK_SEMIRAW(E,4,X);          \
 169   } while (0)
 170
 171 #define FP_PACK_E(val,X)                \
 172   do {                                  \
 173     _FP_PACK_CANONICAL(E,4,X);          \
 174     FP_PACK_RAW_E(val,X);               \
 175   } while (0)
 176
 177 #define FP_PACK_EP(val,X)               \
 178   do {                                  \
 179     _FP_PACK_CANONICAL(E,4,X);          \
 180     FP_PACK_RAW_EP(val,X);              \
 181   } while (0)
 182
 183 #define FP_PACK_SEMIRAW_E(val,X)        \
 184   do {                                  \
 185     _FP_PACK_SEMIRAW(E,4,X);            \
 186     FP_PACK_RAW_E(val,X);               \
 187   } while (0)
 188
 189 #define FP_PACK_SEMIRAW_EP(val,X)       \
 190   do {                                  \
 191     _FP_PACK_SEMIRAW(E,4,X);            \
 192     FP_PACK_RAW_EP(val,X);              \
 193   } while (0)
 194
 195 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,4,X)
 196 #define FP_NEG_E(R,X)           _FP_NEG(E,4,R,X)
 197 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,4,R,X,Y)
 198 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,4,R,X,Y)
 199 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,4,R,X,Y)
 200 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,4,R,X,Y)
 201 #define FP_SQRT_E(R,X)          _FP_SQRT(E,4,R,X)
 202 #define FP_FMA_E(R,X,Y,Z)       _FP_FMA(E,4,8,R,X,Y,Z)
 203
 204 /*
 205  * Square root algorithms:
 206  * We have just one right now, maybe Newton approximation
 207  * should be added for those machines where division is fast.
 208  * This has special _E version because standard _4 square
 209  * root would not work (it has to start normally with the
 210  * second word and not the first), but as we have to do it
 211  * anyway, we optimize it by doing most of the calculations
 212  * in two UWtype registers instead of four.
 213  */
 214
 215 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 216   do {                                                  \
 217     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 218     _FP_FRAC_SRL_4(X, (_FP_WORKBITS));                  \
 219     while (q)                                           \
 220       {                                                 \
 221         T##_f[1] = S##_f[1] + q;                        \
 222         if (T##_f[1] <= X##_f[1])                       \
 223           {                                             \
 224             S##_f[1] = T##_f[1] + q;                    \
 225             X##_f[1] -= T##_f[1];                       \
 226             R##_f[1] += q;                              \
 227           }                                             \
 228         _FP_FRAC_SLL_2(X, 1);                           \
 229         q >>= 1;                                        \
 230       }                                                 \
 231     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 232     while (q)                                           \
 233       {                                                 \
 234         T##_f[0] = S##_f[0] + q;                        \
 235         T##_f[1] = S##_f[1];                            \
 236         if (T##_f[1] < X##_f[1] ||                      \
 237             (T##_f[1] == X##_f[1] &&                    \
 238              T##_f[0] <= X##_f[0]))                     \
 239           {                                             \
 240             S##_f[0] = T##_f[0] + q;                    \
 241             S##_f[1] += (T##_f[0] > S##_f[0]);          \
 242             _FP_FRAC_DEC_2(X, T);                       \
 243             R##_f[0] += q;                              \
 244           }                                             \
 245         _FP_FRAC_SLL_2(X, 1);                           \
 246         q >>= 1;                                        \
 247       }                                                 \
 248     _FP_FRAC_SLL_4(R, (_FP_WORKBITS));                  \
 249     if (X##_f[0] | X##_f[1])                            \
 250       {                                                 \
 251         if (S##_f[1] < X##_f[1] ||                      \
 252             (S##_f[1] == X##_f[1] &&                    \
 253              S##_f[0] < X##_f[0]))                      \
 254           R##_f[0] |= _FP_WORK_ROUND;                   \
 255         R##_f[0] |= _FP_WORK_STICKY;                    \
 256       }                                                 \
 257   } while (0)
 258
 259 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,4,r,X,Y,un)
 260 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,4,r,X,Y)
 261 #define FP_CMP_UNORD_E(r,X,Y)   _FP_CMP_UNORD(E,4,r,X,Y)
 262
 263 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,4,r,X,rsz,rsg)
 264 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,4,X,r,rs,rt)
 265
 266 #define _FP_FRAC_HIGH_E(X)      (X##_f[2])
 267 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f[1])
 268
 269 #define _FP_FRAC_HIGH_DW_E(X)   (X##_f[4])
 270
 271 #else   /* not _FP_W_TYPE_SIZE < 64 */
 272 union _FP_UNION_E
 273 {
 274   XFtype flt;
 275   struct _FP_STRUCT_LAYOUT {
 276 #if __BYTE_ORDER == __BIG_ENDIAN
 277     _FP_W_TYPE pad  : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
 278     unsigned sign   : 1;
 279     unsigned exp    : _FP_EXPBITS_E;
 280     _FP_W_TYPE frac : _FP_W_TYPE_SIZE;
 281 #else
 282     _FP_W_TYPE frac : _FP_W_TYPE_SIZE;
 283     unsigned exp    : _FP_EXPBITS_E;
 284     unsigned sign   : 1;
 285 #endif
 286   } bits;
 287 };
 288
 289 #define FP_DECL_E(X)            _FP_DECL(2,X)
 290
 291 #define FP_UNPACK_RAW_E(X, val)                                 \
 292   do {                                                          \
 293     union _FP_UNION_E _flo; _flo.flt = (val);                   \
 294                                                                 \
 295     X##_f0 = _flo.bits.frac;                                    \
 296     X##_f1 = 0;                                                 \
 297     X##_e = _flo.bits.exp;                                      \
 298     X##_s = _flo.bits.sign;                                     \
 299   } while (0)
 300
 301 #define FP_UNPACK_RAW_EP(X, val)                                \
 302   do {                                                          \
 303     union _FP_UNION_E *_flo =                                   \
 304       (union _FP_UNION_E *)(val);                               \
 305                                                                 \
 306     X##_f0 = _flo->bits.frac;                                   \
 307     X##_f1 = 0;                                                 \
 308     X##_e = _flo->bits.exp;                                     \
 309     X##_s = _flo->bits.sign;                                    \
 310   } while (0)
 311
 312 #define FP_PACK_RAW_E(val, X)                                   \
 313   do {                                                          \
 314     union _FP_UNION_E _flo;                                     \
 315                                                                 \
 316     if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                         \
 317     else X##_f0 &= ~(_FP_IMPLBIT_E);                            \
 318     _flo.bits.frac = X##_f0;                                    \
 319     _flo.bits.exp  = X##_e;                                     \
 320     _flo.bits.sign = X##_s;                                     \
 321                                                                 \
 322     (val) = _flo.flt;                                           \
 323   } while (0)
 324
 325 #define FP_PACK_RAW_EP(fs, val, X)                              \
 326   do {                                                          \
 327     if (!FP_INHIBIT_RESULTS)                                    \
 328       {                                                         \
 329         union _FP_UNION_E *_flo =                               \
 330           (union _FP_UNION_E *)(val);                           \
 331                                                                 \
 332         if (X##_e) X##_f0 |= _FP_IMPLBIT_E;                     \
 333         else X##_f0 &= ~(_FP_IMPLBIT_E);                        \
 334         _flo->bits.frac = X##_f0;                               \
 335         _flo->bits.exp  = X##_e;                                \
 336         _flo->bits.sign = X##_s;                                \
 337       }                                                         \
 338   } while (0)
 339
 340
 341 #define FP_UNPACK_E(X,val)              \
 342   do {                                  \
 343     FP_UNPACK_RAW_E(X,val);             \
 344     _FP_UNPACK_CANONICAL(E,2,X);        \
 345   } while (0)
 346
 347 #define FP_UNPACK_EP(X,val)             \
 348   do {                                  \
 349     FP_UNPACK_RAW_EP(X,val);            \
 350     _FP_UNPACK_CANONICAL(E,2,X);        \
 351   } while (0)
 352
 353 #define FP_UNPACK_SEMIRAW_E(X,val)      \
 354   do {                                  \
 355     FP_UNPACK_RAW_E(X,val);             \
 356     _FP_UNPACK_SEMIRAW(E,2,X);          \
 357   } while (0)
 358
 359 #define FP_UNPACK_SEMIRAW_EP(X,val)     \
 360   do {                                  \
 361     FP_UNPACK_RAW_EP(X,val);            \
 362     _FP_UNPACK_SEMIRAW(E,2,X);          \
 363   } while (0)
 364
 365 #define FP_PACK_E(val,X)                \
 366   do {                                  \
 367     _FP_PACK_CANONICAL(E,2,X);          \
 368     FP_PACK_RAW_E(val,X);               \
 369   } while (0)
 370
 371 #define FP_PACK_EP(val,X)               \
 372   do {                                  \
 373     _FP_PACK_CANONICAL(E,2,X);          \
 374     FP_PACK_RAW_EP(val,X);              \
 375   } while (0)
 376
 377 #define FP_PACK_SEMIRAW_E(val,X)        \
 378   do {                                  \
 379     _FP_PACK_SEMIRAW(E,2,X);            \
 380     FP_PACK_RAW_E(val,X);               \
 381   } while (0)
 382
 383 #define FP_PACK_SEMIRAW_EP(val,X)       \
 384   do {                                  \
 385     _FP_PACK_SEMIRAW(E,2,X);            \
 386     FP_PACK_RAW_EP(val,X);              \
 387   } while (0)
 388
 389 #define FP_ISSIGNAN_E(X)        _FP_ISSIGNAN(E,2,X)
 390 #define FP_NEG_E(R,X)           _FP_NEG(E,2,R,X)
 391 #define FP_ADD_E(R,X,Y)         _FP_ADD(E,2,R,X,Y)
 392 #define FP_SUB_E(R,X,Y)         _FP_SUB(E,2,R,X,Y)
 393 #define FP_MUL_E(R,X,Y)         _FP_MUL(E,2,R,X,Y)
 394 #define FP_DIV_E(R,X,Y)         _FP_DIV(E,2,R,X,Y)
 395 #define FP_SQRT_E(R,X)          _FP_SQRT(E,2,R,X)
 396 #define FP_FMA_E(R,X,Y,Z)       _FP_FMA(E,2,4,R,X,Y,Z)
 397
 398 /*
 399  * Square root algorithms:
 400  * We have just one right now, maybe Newton approximation
 401  * should be added for those machines where division is fast.
 402  * We optimize it by doing most of the calculations
 403  * in one UWtype registers instead of two, although we don't
 404  * have to.
 405  */
 406 #define _FP_SQRT_MEAT_E(R, S, T, X, q)                  \
 407   do {                                                  \
 408     q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);         \
 409     _FP_FRAC_SRL_2(X, (_FP_WORKBITS));                  \
 410     while (q)                                           \
 411       {                                                 \
 412         T##_f0 = S##_f0 + q;                            \
 413         if (T##_f0 <= X##_f0)                           \
 414           {                                             \
 415             S##_f0 = T##_f0 + q;                        \
 416             X##_f0 -= T##_f0;                           \
 417             R##_f0 += q;                                \
 418           }                                             \
 419         _FP_FRAC_SLL_1(X, 1);                           \
 420         q >>= 1;                                        \
 421       }                                                 \
 422     _FP_FRAC_SLL_2(R, (_FP_WORKBITS));                  \
 423     if (X##_f0)                                         \
 424       {                                                 \
 425         if (S##_f0 < X##_f0)                            \
 426           R##_f0 |= _FP_WORK_ROUND;                     \
 427         R##_f0 |= _FP_WORK_STICKY;                      \
 428       }                                                 \
 429   } while (0)
 430
 431 #define FP_CMP_E(r,X,Y,un)      _FP_CMP(E,2,r,X,Y,un)
 432 #define FP_CMP_EQ_E(r,X,Y)      _FP_CMP_EQ(E,2,r,X,Y)
 433 #define FP_CMP_UNORD_E(r,X,Y)   _FP_CMP_UNORD(E,2,r,X,Y)
 434
 435 #define FP_TO_INT_E(r,X,rsz,rsg)        _FP_TO_INT(E,2,r,X,rsz,rsg)
 436 #define FP_FROM_INT_E(X,r,rs,rt)        _FP_FROM_INT(E,2,X,r,rs,rt)
 437
 438 #define _FP_FRAC_HIGH_E(X)      (X##_f1)
 439 #define _FP_FRAC_HIGH_RAW_E(X)  (X##_f0)
 440
 441 #define _FP_FRAC_HIGH_DW_E(X)   (X##_f[2])
 442
 443 #endif /* not _FP_W_TYPE_SIZE < 64 */