Print this page
11210 libm should be cstyle(1ONBLD) clean

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/m9x/fmal.c
          +++ new/usr/src/lib/libm/common/m9x/fmal.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
  21   21  
  22   22  /*
  23   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
       25 +
  25   26  /*
  26   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  #pragma weak fmal = __fmal
  31   32  
  32   33  #include "libm.h"
  33   34  #include "fma.h"
  34   35  #include "fenv_inlines.h"
  35   36  
  36   37  #if defined(__sparc)
  37      -
  38   38  static const union {
  39   39          unsigned i[2];
  40   40          double d;
  41   41  } C[] = {
  42   42          { 0x3fe00000u, 0 },
  43   43          { 0x40000000u, 0 },
  44   44          { 0x3ef00000u, 0 },
  45   45          { 0x3e700000u, 0 },
  46   46          { 0x41300000u, 0 },
  47   47          { 0x3e300000u, 0 },
↓ open down ↓ 2 lines elided ↑ open up ↑
  50   50          { 0x42300000u, 0 },
  51   51          { 0x3df00000u, 0 },
  52   52          { 0x7fe00000u, 0 },
  53   53          { 0x00100000u, 0 },
  54   54          { 0x00100001u, 0 },
  55   55          { 0, 0 },
  56   56          { 0x7ff00000u, 0 },
  57   57          { 0x7ff00001u, 0 }
  58   58  };
  59   59  
  60      -#define half    C[0].d
  61      -#define two     C[1].d
  62      -#define twom16  C[2].d
  63      -#define twom24  C[3].d
  64      -#define two20   C[4].d
  65      -#define twom28  C[5].d
  66      -#define twom76  C[6].d
  67      -#define twom124 C[7].d
  68      -#define two36   C[8].d
  69      -#define twom32  C[9].d
  70      -#define huge    C[10].d
  71      -#define tiny    C[11].d
  72      -#define tiny2   C[12].d
  73      -#define zero    C[13].d
  74      -#define inf     C[14].d
  75      -#define snan    C[15].d
       60 +#define half            C[0].d
       61 +#define two             C[1].d
       62 +#define twom16          C[2].d
       63 +#define twom24          C[3].d
       64 +#define two20           C[4].d
       65 +#define twom28          C[5].d
       66 +#define twom76          C[6].d
       67 +#define twom124         C[7].d
       68 +#define two36           C[8].d
       69 +#define twom32          C[9].d
       70 +#define huge            C[10].d
       71 +#define tiny            C[11].d
       72 +#define tiny2           C[12].d
       73 +#define zero            C[13].d
       74 +#define inf             C[14].d
       75 +#define snan            C[15].d
  76   76  
  77   77  static const unsigned int fsr_rm = 0xc0000000u;
  78   78  
  79   79  /*
  80   80   * fmal for SPARC: 128-bit quad precision, big-endian
  81   81   */
  82   82  long double
  83      -__fmal(long double x, long double y, long double z) {
       83 +__fmal(long double x, long double y, long double z)
       84 +{
  84   85          union {
  85   86                  unsigned int i[4];
  86   87                  long double q;
  87   88          } xx, yy, zz;
  88   89          union {
  89   90                  unsigned int i[2];
  90   91                  double d;
  91   92          } u;
       93 +
  92   94          double dx[5], dy[5], dxy[9], c, s;
  93   95          unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
  94   96          unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
  95   97          unsigned int rm, sticky;
  96   98          unsigned int fsr;
  97   99          int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit;
  98  100          int cx, cy, cz;
  99      -        volatile double dummy;
      101 +        volatile double dummy;
 100  102  
 101  103          /* extract the high order words of the arguments */
 102  104          xx.q = x;
 103  105          yy.q = y;
 104  106          zz.q = z;
 105  107          hx = xx.i[0] & ~0x80000000;
 106  108          hy = yy.i[0] & ~0x80000000;
 107  109          hz = zz.i[0] & ~0x80000000;
 108  110  
 109  111          /*
↓ open down ↓ 2 lines elided ↑ open up ↑
 112  114           */
 113  115          if (hx >= 0x7fff0000) {
 114  116                  if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) {
 115  117                          if (!(hx & 0x8000)) {
 116  118                                  /* signaling nan, raise invalid */
 117  119                                  dummy = snan;
 118  120                                  dummy += snan;
 119  121                                  xx.i[0] |= 0x8000;
 120  122                                  return (xx.q);
 121  123                          }
 122      -                        cx = 3; /* quiet nan */
 123      -                } else
 124      -                        cx = 2; /* inf */
      124 +
      125 +                        cx = 3;         /* quiet nan */
      126 +                } else {
      127 +                        cx = 2;         /* inf */
      128 +                }
 125  129          } else if (hx == 0) {
 126  130                  cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
 127      -                                /* subnormal or zero */
 128      -        } else
 129      -                cx = 1;         /* finite nonzero */
      131 +                /* subnormal or zero */
      132 +        } else {
      133 +                cx = 1;                 /* finite nonzero */
      134 +        }
 130  135  
 131  136          if (hy >= 0x7fff0000) {
 132  137                  if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
 133  138                          if (!(hy & 0x8000)) {
 134  139                                  dummy = snan;
 135  140                                  dummy += snan;
 136  141                                  yy.i[0] |= 0x8000;
 137  142                                  return (yy.q);
 138  143                          }
      144 +
 139  145                          cy = 3;
 140      -                } else
      146 +                } else {
 141  147                          cy = 2;
      148 +                }
 142  149          } else if (hy == 0) {
 143  150                  cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
 144      -        } else
      151 +        } else {
 145  152                  cy = 1;
      153 +        }
 146  154  
 147  155          if (hz >= 0x7fff0000) {
 148  156                  if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
 149  157                          if (!(hz & 0x8000)) {
 150  158                                  dummy = snan;
 151  159                                  dummy += snan;
 152  160                                  zz.i[0] |= 0x8000;
 153  161                                  return (zz.q);
 154  162                          }
      163 +
 155  164                          cz = 3;
 156      -                } else
      165 +                } else {
 157  166                          cz = 2;
      167 +                }
 158  168          } else if (hz == 0) {
 159  169                  cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
 160      -        } else
      170 +        } else {
 161  171                  cz = 1;
      172 +        }
 162  173  
 163  174          /* get the fsr and clear current exceptions */
 164  175          __fenv_getfsr32(&fsr);
 165  176          fsr &= ~FSR_CEXC;
 166  177  
 167  178          /* handle all other zero, inf, and nan cases */
 168  179          if (cx != 1 || cy != 1 || cz != 1) {
 169  180                  /* if x or y is a quiet nan, return it */
 170  181                  if (cx == 3) {
 171  182                          __fenv_setfsr32(&fsr);
 172  183                          return (x);
 173  184                  }
      185 +
 174  186                  if (cy == 3) {
 175  187                          __fenv_setfsr32(&fsr);
 176  188                          return (y);
 177  189                  }
 178  190  
 179  191                  /* if x*y is 0*inf, raise invalid and return the default nan */
 180  192                  if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) {
 181  193                          dummy = zero;
 182  194                          dummy *= inf;
 183  195                          zz.i[0] = 0x7fffffff;
↓ open down ↓ 10 lines elided ↑ open up ↑
 194  206                  /*
 195  207                   * now none of x, y, or z is nan; handle cases where x or y
 196  208                   * is inf
 197  209                   */
 198  210                  if (cx == 2 || cy == 2) {
 199  211                          /*
 200  212                           * if z is also inf, either we have inf-inf or
 201  213                           * the result is the same as z depending on signs
 202  214                           */
 203  215                          if (cz == 2) {
 204      -                                if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
      216 +                                if ((int)((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
 205  217                                          dummy = inf;
 206  218                                          dummy -= inf;
 207  219                                          zz.i[0] = 0x7fffffff;
 208  220                                          zz.i[1] = zz.i[2] = zz.i[3] =
 209      -                                                0xffffffff;
      221 +                                            0xffffffff;
 210  222                                          return (zz.q);
 211  223                                  }
      224 +
 212  225                                  __fenv_setfsr32(&fsr);
 213  226                                  return (z);
 214  227                          }
 215  228  
 216  229                          /* otherwise the result is inf with appropriate sign */
 217  230                          zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) |
 218      -                                0x7fff0000;
      231 +                            0x7fff0000;
 219  232                          zz.i[1] = zz.i[2] = zz.i[3] = 0;
 220  233                          __fenv_setfsr32(&fsr);
 221  234                          return (zz.q);
 222  235                  }
 223  236  
 224  237                  /* if z is inf, return it */
 225  238                  if (cz == 2) {
 226  239                          __fenv_setfsr32(&fsr);
 227  240                          return (z);
 228  241                  }
 229  242  
 230  243                  /*
 231  244                   * now x, y, and z are all finite; handle cases where x or y
 232  245                   * is zero
 233  246                   */
 234  247                  if (cx == 0 || cy == 0) {
 235  248                          /* either we have 0-0 or the result is the same as z */
 236      -                        if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
 237      -                                0) {
      249 +                        if (cz == 0 && (int)((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
      250 +                            0) {
 238  251                                  zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
 239      -                                        0;
      252 +                                    0;
 240  253                                  __fenv_setfsr32(&fsr);
 241  254                                  return (zz.q);
 242  255                          }
      256 +
 243  257                          __fenv_setfsr32(&fsr);
 244  258                          return (z);
 245  259                  }
 246  260  
 247  261                  /* if we get here, x and y are nonzero finite, z must be zero */
 248  262                  return (x * y);
 249  263          }
 250  264  
 251  265          /*
 252  266           * now x, y, and z are all finite and nonzero; set round-to-
↓ open down ↓ 1 lines elided ↑ open up ↑
 254  268           */
 255  269          __fenv_setfsr32(&fsr_rm);
 256  270  
 257  271          /*
 258  272           * get the signs and exponents and normalize the significands
 259  273           * of x and y
 260  274           */
 261  275          sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000;
 262  276          ex = hx >> 16;
 263  277          hx &= 0xffff;
      278 +
 264  279          if (!ex) {
 265  280                  if (hx | (xx.i[1] & 0xfffe0000)) {
 266  281                          ex = 1;
 267  282                  } else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) {
 268  283                          hx = xx.i[1];
 269  284                          xx.i[1] = xx.i[2];
 270  285                          xx.i[2] = xx.i[3];
 271  286                          xx.i[3] = 0;
 272  287                          ex = -31;
 273  288                  } else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) {
 274  289                          hx = xx.i[2];
 275  290                          xx.i[1] = xx.i[3];
 276  291                          xx.i[2] = xx.i[3] = 0;
 277  292                          ex = -63;
 278  293                  } else {
 279  294                          hx = xx.i[3];
 280  295                          xx.i[1] = xx.i[2] = xx.i[3] = 0;
 281  296                          ex = -95;
 282  297                  }
      298 +
 283  299                  while ((hx & 0x10000) == 0) {
 284  300                          hx = (hx << 1) | (xx.i[1] >> 31);
 285  301                          xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
 286  302                          xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
 287  303                          xx.i[3] <<= 1;
 288  304                          ex--;
 289  305                  }
 290      -        } else
      306 +        } else {
 291  307                  hx |= 0x10000;
      308 +        }
      309 +
 292  310          ey = hy >> 16;
 293  311          hy &= 0xffff;
      312 +
 294  313          if (!ey) {
 295  314                  if (hy | (yy.i[1] & 0xfffe0000)) {
 296  315                          ey = 1;
 297  316                  } else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
 298  317                          hy = yy.i[1];
 299  318                          yy.i[1] = yy.i[2];
 300  319                          yy.i[2] = yy.i[3];
 301  320                          yy.i[3] = 0;
 302  321                          ey = -31;
 303  322                  } else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) {
 304  323                          hy = yy.i[2];
 305  324                          yy.i[1] = yy.i[3];
 306  325                          yy.i[2] = yy.i[3] = 0;
 307  326                          ey = -63;
 308  327                  } else {
 309  328                          hy = yy.i[3];
 310  329                          yy.i[1] = yy.i[2] = yy.i[3] = 0;
 311  330                          ey = -95;
 312  331                  }
      332 +
 313  333                  while ((hy & 0x10000) == 0) {
 314  334                          hy = (hy << 1) | (yy.i[1] >> 31);
 315  335                          yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
 316  336                          yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
 317  337                          yy.i[3] <<= 1;
 318  338                          ey--;
 319  339                  }
 320      -        } else
      340 +        } else {
 321  341                  hy |= 0x10000;
      342 +        }
      343 +
 322  344          exy = ex + ey - 0x3fff;
 323  345  
 324  346          /* convert the significands of x and y to doubles */
 325  347          c = twom16;
 326      -        dx[0] = (double) ((int) hx) * c;
 327      -        dy[0] = (double) ((int) hy) * c;
      348 +        dx[0] = (double)((int)hx) * c;
      349 +        dy[0] = (double)((int)hy) * c;
 328  350  
 329  351          c *= twom24;
 330      -        dx[1] = (double) ((int) (xx.i[1] >> 8)) * c;
 331      -        dy[1] = (double) ((int) (yy.i[1] >> 8)) * c;
      352 +        dx[1] = (double)((int)(xx.i[1] >> 8)) * c;
      353 +        dy[1] = (double)((int)(yy.i[1] >> 8)) * c;
 332  354  
 333  355          c *= twom24;
 334      -        dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) &
      356 +        dx[2] = (double)((int)(((xx.i[1] << 16) | (xx.i[2] >> 16)) &
 335  357              0xffffff)) * c;
 336      -        dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) &
      358 +        dy[2] = (double)((int)(((yy.i[1] << 16) | (yy.i[2] >> 16)) &
 337  359              0xffffff)) * c;
 338  360  
 339  361          c *= twom24;
 340      -        dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) &
 341      -            0xffffff)) * c;
 342      -        dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) &
 343      -            0xffffff)) * c;
      362 +        dx[3] = (double)((int)(((xx.i[2] << 8) | (xx.i[3] >> 24)) & 0xffffff)) *
      363 +            c;
      364 +        dy[3] = (double)((int)(((yy.i[2] << 8) | (yy.i[3] >> 24)) & 0xffffff)) *
      365 +            c;
 344  366  
 345  367          c *= twom24;
 346      -        dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c;
 347      -        dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c;
      368 +        dx[4] = (double)((int)(xx.i[3] & 0xffffff)) * c;
      369 +        dy[4] = (double)((int)(yy.i[3] & 0xffffff)) * c;
 348  370  
 349  371          /* form the "digits" of the product */
 350  372          dxy[0] = dx[0] * dy[0];
 351  373          dxy[1] = dx[0] * dy[1] + dx[1] * dy[0];
 352  374          dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
 353      -        dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] +
 354      -            dx[3] * dy[0];
 355      -        dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] +
 356      -            dx[3] * dy[1] + dx[4] * dy[0];
 357      -        dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] +
 358      -            dx[4] * dy[1];
      375 +        dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] + dx[3] * dy[0];
      376 +        dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] + dx[3] * dy[1] +
      377 +            dx[4] * dy[0];
      378 +        dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] + dx[4] * dy[1];
 359  379          dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
 360  380          dxy[7] = dx[3] * dy[4] + dx[4] * dy[3];
 361  381          dxy[8] = dx[4] * dy[4];
 362  382  
 363  383          /* split odd-numbered terms and combine into even-numbered terms */
 364  384          c = (dxy[1] + two20) - two20;
 365  385          dxy[0] += c;
 366  386          dxy[1] -= c;
 367  387          c = (dxy[3] + twom28) - twom28;
 368  388          dxy[2] += c + dxy[1];
↓ open down ↓ 3 lines elided ↑ open up ↑
 372  392          dxy[5] -= c;
 373  393          c = (dxy[7] + twom124) - twom124;
 374  394          dxy[6] += c + dxy[5];
 375  395          dxy[8] += (dxy[7] - c);
 376  396  
 377  397          /* propagate carries, adjusting the exponent if need be */
 378  398          dxy[7] = dxy[6] + dxy[8];
 379  399          dxy[5] = dxy[4] + dxy[7];
 380  400          dxy[3] = dxy[2] + dxy[5];
 381  401          dxy[1] = dxy[0] + dxy[3];
      402 +
 382  403          if (dxy[1] >= two) {
 383  404                  dxy[0] *= half;
 384  405                  dxy[1] *= half;
 385  406                  dxy[2] *= half;
 386  407                  dxy[3] *= half;
 387  408                  dxy[4] *= half;
 388  409                  dxy[5] *= half;
 389  410                  dxy[6] *= half;
 390  411                  dxy[7] *= half;
 391  412                  dxy[8] *= half;
↓ open down ↓ 49 lines elided ↑ open up ↑
 441  462          dxy[8] -= c;
 442  463  
 443  464          s *= twom32;
 444  465          u.d = c = dxy[8] + s;
 445  466          xy7 = u.i[1];
 446  467  
 447  468          /* extract the sign, exponent, and significand of z */
 448  469          sz = zz.i[0] & 0x80000000;
 449  470          ez = hz >> 16;
 450  471          z0 = hz & 0xffff;
      472 +
 451  473          if (!ez) {
 452  474                  if (z0 | (zz.i[1] & 0xfffe0000)) {
 453  475                          z1 = zz.i[1];
 454  476                          z2 = zz.i[2];
 455  477                          z3 = zz.i[3];
 456  478                          ez = 1;
 457  479                  } else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) {
 458  480                          z0 = zz.i[1];
 459  481                          z1 = zz.i[2];
 460  482                          z2 = zz.i[3];
↓ open down ↓ 2 lines elided ↑ open up ↑
 463  485                  } else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) {
 464  486                          z0 = zz.i[2];
 465  487                          z1 = zz.i[3];
 466  488                          z2 = z3 = 0;
 467  489                          ez = -63;
 468  490                  } else {
 469  491                          z0 = zz.i[3];
 470  492                          z1 = z2 = z3 = 0;
 471  493                          ez = -95;
 472  494                  }
      495 +
 473  496                  while ((z0 & 0x10000) == 0) {
 474  497                          z0 = (z0 << 1) | (z1 >> 31);
 475  498                          z1 = (z1 << 1) | (z2 >> 31);
 476  499                          z2 = (z2 << 1) | (z3 >> 31);
 477  500                          z3 <<= 1;
 478  501                          ez--;
 479  502                  }
 480  503          } else {
 481  504                  z0 |= 0x10000;
 482  505                  z1 = zz.i[1];
 483  506                  z2 = zz.i[2];
 484  507                  z3 = zz.i[3];
 485  508          }
      509 +
 486  510          z4 = z5 = z6 = z7 = 0;
 487  511  
 488  512          /*
 489  513           * now x*y is represented by sxy, exy, and xy[0-7], and z is
 490  514           * represented likewise; swap if need be so |xy| <= |z|
 491  515           */
 492  516          if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
 493      -                (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 ||
 494      -                (xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
 495      -                e = sxy; sxy = sz; sz = e;
 496      -                e = exy; exy = ez; ez = e;
 497      -                e = xy0; xy0 = z0; z0 = e;
 498      -                e = xy1; xy1 = z1; z1 = e;
 499      -                e = xy2; xy2 = z2; z2 = e;
 500      -                e = xy3; xy3 = z3; z3 = e;
 501      -                z4 = xy4; xy4 = 0;
 502      -                z5 = xy5; xy5 = 0;
 503      -                z6 = xy6; xy6 = 0;
 504      -                z7 = xy7; xy7 = 0;
      517 +            (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 || (xy3 == z3 &&
      518 +            (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
      519 +                e = sxy;
      520 +                sxy = sz;
      521 +                sz = e;
      522 +                e = exy;
      523 +                exy = ez;
      524 +                ez = e;
      525 +                e = xy0;
      526 +                xy0 = z0;
      527 +                z0 = e;
      528 +                e = xy1;
      529 +                xy1 = z1;
      530 +                z1 = e;
      531 +                e = xy2;
      532 +                xy2 = z2;
      533 +                z2 = e;
      534 +                e = xy3;
      535 +                xy3 = z3;
      536 +                z3 = e;
      537 +                z4 = xy4;
      538 +                xy4 = 0;
      539 +                z5 = xy5;
      540 +                xy5 = 0;
      541 +                z6 = xy6;
      542 +                xy6 = 0;
      543 +                z7 = xy7;
      544 +                xy7 = 0;
 505  545          }
 506  546  
 507  547          /* shift the significand of xy keeping a sticky bit */
 508  548          e = ez - exy;
      549 +
 509  550          if (e > 236) {
 510  551                  xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
 511  552                  xy7 = 1;
 512  553          } else if (e >= 224) {
 513      -                sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 |
 514      -                        ((xy0 << 1) << (255 - e));
      554 +                sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 | ((xy0 <<
      555 +                    1) << (255 - e));
 515  556                  xy7 = xy0 >> (e - 224);
      557 +
 516  558                  if (sticky)
 517  559                          xy7 |= 1;
      560 +
 518  561                  xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
 519  562          } else if (e >= 192) {
 520      -                sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 |
 521      -                        ((xy1 << 1) << (223 - e));
      563 +                sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | ((xy1 << 1) <<
      564 +                    (223 - e));
 522  565                  xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
      566 +
 523  567                  if (sticky)
 524  568                          xy7 |= 1;
      569 +
 525  570                  xy6 = xy0 >> (e - 192);
 526  571                  xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0;
 527  572          } else if (e >= 160) {
 528      -                sticky = xy7 | xy6 | xy5 | xy4 | xy3 |
 529      -                        ((xy2 << 1) << (191 - e));
      573 +                sticky = xy7 | xy6 | xy5 | xy4 | xy3 | ((xy2 << 1) << (191 -
      574 +                    e));
 530  575                  xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
      576 +
 531  577                  if (sticky)
 532  578                          xy7 |= 1;
      579 +
 533  580                  xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
 534  581                  xy5 = xy0 >> (e - 160);
 535  582                  xy0 = xy1 = xy2 = xy3 = xy4 = 0;
 536  583          } else if (e >= 128) {
 537  584                  sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
 538  585                  xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
      586 +
 539  587                  if (sticky)
 540  588                          xy7 |= 1;
      589 +
 541  590                  xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
 542  591                  xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
 543  592                  xy4 = xy0 >> (e - 128);
 544  593                  xy0 = xy1 = xy2 = xy3 = 0;
 545  594          } else if (e >= 96) {
 546  595                  sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
 547  596                  xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
      597 +
 548  598                  if (sticky)
 549  599                          xy7 |= 1;
      600 +
 550  601                  xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
 551  602                  xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
 552  603                  xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
 553  604                  xy3 = xy0 >> (e - 96);
 554  605                  xy0 = xy1 = xy2 = 0;
 555  606          } else if (e >= 64) {
 556  607                  sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
 557  608                  xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
      609 +
 558  610                  if (sticky)
 559  611                          xy7 |= 1;
      612 +
 560  613                  xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
 561  614                  xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
 562  615                  xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
 563  616                  xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
 564  617                  xy2 = xy0 >> (e - 64);
 565  618                  xy0 = xy1 = 0;
 566  619          } else if (e >= 32) {
 567  620                  sticky = xy7 | ((xy6 << 1) << (63 - e));
 568  621                  xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
      622 +
 569  623                  if (sticky)
 570  624                          xy7 |= 1;
      625 +
 571  626                  xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
 572  627                  xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
 573  628                  xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
 574  629                  xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
 575  630                  xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
 576  631                  xy1 = xy0 >> (e - 32);
 577  632                  xy0 = 0;
 578  633          } else if (e) {
 579  634                  sticky = (xy7 << 1) << (31 - e);
 580  635                  xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
      636 +
 581  637                  if (sticky)
 582  638                          xy7 |= 1;
      639 +
 583  640                  xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
 584  641                  xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
 585  642                  xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
 586  643                  xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
 587  644                  xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
 588  645                  xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
 589  646                  xy0 >>= e;
 590  647          }
 591  648  
 592  649          /* if this is a magnitude subtract, negate the significand of xy */
 593  650          if (sxy ^ sz) {
 594  651                  xy0 = ~xy0;
 595  652                  xy1 = ~xy1;
 596  653                  xy2 = ~xy2;
 597  654                  xy3 = ~xy3;
 598  655                  xy4 = ~xy4;
 599  656                  xy5 = ~xy5;
 600  657                  xy6 = ~xy6;
 601  658                  xy7 = -xy7;
      659 +
 602  660                  if (xy7 == 0)
 603  661                          if (++xy6 == 0)
 604  662                                  if (++xy5 == 0)
 605  663                                          if (++xy4 == 0)
 606  664                                                  if (++xy3 == 0)
 607  665                                                          if (++xy2 == 0)
 608  666                                                                  if (++xy1 == 0)
 609  667                                                                          xy0++;
 610  668          }
 611  669  
 612  670          /* add, propagating carries */
 613  671          z7 += xy7;
 614  672          e = (z7 < xy7);
 615  673          z6 += xy6;
      674 +
 616  675          if (e) {
 617  676                  z6++;
 618  677                  e = (z6 <= xy6);
 619      -        } else
      678 +        } else {
 620  679                  e = (z6 < xy6);
      680 +        }
      681 +
 621  682          z5 += xy5;
      683 +
 622  684          if (e) {
 623  685                  z5++;
 624  686                  e = (z5 <= xy5);
 625      -        } else
      687 +        } else {
 626  688                  e = (z5 < xy5);
      689 +        }
      690 +
 627  691          z4 += xy4;
      692 +
 628  693          if (e) {
 629  694                  z4++;
 630  695                  e = (z4 <= xy4);
 631      -        } else
      696 +        } else {
 632  697                  e = (z4 < xy4);
      698 +        }
      699 +
 633  700          z3 += xy3;
      701 +
 634  702          if (e) {
 635  703                  z3++;
 636  704                  e = (z3 <= xy3);
 637      -        } else
      705 +        } else {
 638  706                  e = (z3 < xy3);
      707 +        }
      708 +
 639  709          z2 += xy2;
      710 +
 640  711          if (e) {
 641  712                  z2++;
 642  713                  e = (z2 <= xy2);
 643      -        } else
      714 +        } else {
 644  715                  e = (z2 < xy2);
      716 +        }
      717 +
 645  718          z1 += xy1;
      719 +
 646  720          if (e) {
 647  721                  z1++;
 648  722                  e = (z1 <= xy1);
 649      -        } else
      723 +        } else {
 650  724                  e = (z1 < xy1);
      725 +        }
      726 +
 651  727          z0 += xy0;
      728 +
 652  729          if (e)
 653  730                  z0++;
 654  731  
 655  732          /* postnormalize and collect rounding information into z4 */
 656  733          if (ez < 1) {
 657  734                  /* result is tiny; shift right until exponent is within range */
 658  735                  e = 1 - ez;
      736 +
 659  737                  if (e > 116) {
 660      -                        z4 = 1; /* result can't be exactly zero */
      738 +                        z4 = 1;         /* result can't be exactly zero */
 661  739                          z0 = z1 = z2 = z3 = 0;
 662  740                  } else if (e >= 96) {
 663      -                        sticky = z7 | z6 | z5 | z4 | z3 | z2 |
 664      -                                ((z1 << 1) << (127 - e));
      741 +                        sticky = z7 | z6 | z5 | z4 | z3 | z2 | ((z1 << 1) <<
      742 +                            (127 - e));
 665  743                          z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
      744 +
 666  745                          if (sticky)
 667  746                                  z4 |= 1;
      747 +
 668  748                          z3 = z0 >> (e - 96);
 669  749                          z0 = z1 = z2 = 0;
 670  750                  } else if (e >= 64) {
 671      -                        sticky = z7 | z6 | z5 | z4 | z3 |
 672      -                                ((z2 << 1) << (95 - e));
      751 +                        sticky = z7 | z6 | z5 | z4 | z3 | ((z2 << 1) << (95 -
      752 +                            e));
 673  753                          z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
      754 +
 674  755                          if (sticky)
 675  756                                  z4 |= 1;
      757 +
 676  758                          z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
 677  759                          z2 = z0 >> (e - 64);
 678  760                          z0 = z1 = 0;
 679  761                  } else if (e >= 32) {
 680  762                          sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
 681  763                          z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
      764 +
 682  765                          if (sticky)
 683  766                                  z4 |= 1;
      767 +
 684  768                          z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
 685  769                          z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
 686  770                          z1 = z0 >> (e - 32);
 687  771                          z0 = 0;
 688  772                  } else {
 689  773                          sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
 690  774                          z4 = (z4 >> e) | ((z3 << 1) << (31 - e));
      775 +
 691  776                          if (sticky)
 692  777                                  z4 |= 1;
      778 +
 693  779                          z3 = (z3 >> e) | ((z2 << 1) << (31 - e));
 694  780                          z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
 695  781                          z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
 696  782                          z0 >>= e;
 697  783                  }
      784 +
 698  785                  ez = 1;
 699  786          } else if (z0 >= 0x20000) {
 700  787                  /* carry out; shift right by one */
 701  788                  sticky = (z4 & 1) | z5 | z6 | z7;
 702  789                  z4 = (z4 >> 1) | (z3 << 31);
      790 +
 703  791                  if (sticky)
 704  792                          z4 |= 1;
      793 +
 705  794                  z3 = (z3 >> 1) | (z2 << 31);
 706  795                  z2 = (z2 >> 1) | (z1 << 31);
 707  796                  z1 = (z1 >> 1) | (z0 << 31);
 708  797                  z0 >>= 1;
 709  798                  ez++;
 710  799          } else {
 711      -                if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7)
 712      -                        != 0) {
      800 +                if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7) !=
      801 +                    0) {
 713  802                          /*
 714  803                           * borrow/cancellation; shift left as much as
 715  804                           * exponent allows
 716  805                           */
 717  806                          while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
 718  807                                  z0 = z1;
 719  808                                  z1 = z2;
 720  809                                  z2 = z3;
 721  810                                  z3 = z4;
 722  811                                  z4 = z5;
 723  812                                  z5 = z6;
 724  813                                  z6 = z7;
 725  814                                  z7 = 0;
 726  815                                  ez -= 32;
 727  816                          }
      817 +
 728  818                          while (z0 < 0x10000 && ez > 1) {
 729  819                                  z0 = (z0 << 1) | (z1 >> 31);
 730  820                                  z1 = (z1 << 1) | (z2 >> 31);
 731  821                                  z2 = (z2 << 1) | (z3 >> 31);
 732  822                                  z3 = (z3 << 1) | (z4 >> 31);
 733  823                                  z4 = (z4 << 1) | (z5 >> 31);
 734  824                                  z5 = (z5 << 1) | (z6 >> 31);
 735  825                                  z6 = (z6 << 1) | (z7 >> 31);
 736  826                                  z7 <<= 1;
 737  827                                  ez--;
 738  828                          }
 739  829                  }
      830 +
 740  831                  if (z5 | z6 | z7)
 741  832                          z4 |= 1;
 742  833          }
 743  834  
 744  835          /* get the rounding mode */
 745  836          rm = fsr >> 30;
 746  837  
 747  838          /* strip off the integer bit, if there is one */
 748  839          ibit = z0 & 0x10000;
 749      -        if (ibit)
      840 +
      841 +        if (ibit) {
 750  842                  z0 -= 0x10000;
 751      -        else {
      843 +        } else {
 752  844                  ez = 0;
 753      -                if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
      845 +
      846 +                if (!(z0 | z1 | z2 | z3 | z4)) {        /* exact zero */
 754  847                          zz.i[0] = rm == FSR_RM ? 0x80000000 : 0;
 755  848                          zz.i[1] = zz.i[2] = zz.i[3] = 0;
 756  849                          __fenv_setfsr32(&fsr);
 757  850                          return (zz.q);
 758  851                  }
 759  852          }
 760  853  
 761  854          /*
 762  855           * flip the sense of directed roundings if the result is negative;
 763  856           * the logic below applies to a positive result
 764  857           */
 765  858          if (sz)
 766  859                  rm ^= rm >> 1;
 767  860  
 768  861          /* round and raise exceptions */
 769  862          if (z4) {
 770  863                  fsr |= FSR_NXC;
 771  864  
 772  865                  /* decide whether to round the fraction up */
 773      -                if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u ||
 774      -                        (z4 == 0x80000000u && (z3 & 1))))) {
      866 +                if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u || (z4 ==
      867 +                    0x80000000u && (z3 & 1))))) {
 775  868                          /* round up and renormalize if necessary */
 776  869                          if (++z3 == 0)
 777  870                                  if (++z2 == 0)
 778  871                                          if (++z1 == 0)
 779  872                                                  if (++z0 == 0x10000) {
 780  873                                                          z0 = 0;
 781  874                                                          ez++;
 782  875                                                  }
 783  876                  }
 784  877          }
 785  878  
 786  879          /* check for under/overflow */
 787  880          if (ez >= 0x7fff) {
 788  881                  if (rm == FSR_RN || rm == FSR_RP) {
 789  882                          zz.i[0] = sz | 0x7fff0000;
 790  883                          zz.i[1] = zz.i[2] = zz.i[3] = 0;
 791  884                  } else {
 792  885                          zz.i[0] = sz | 0x7ffeffff;
 793  886                          zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
 794  887                  }
      888 +
 795  889                  fsr |= FSR_OFC | FSR_NXC;
 796  890          } else {
 797  891                  zz.i[0] = sz | (ez << 16) | z0;
 798  892                  zz.i[1] = z1;
 799  893                  zz.i[2] = z2;
 800  894                  zz.i[3] = z3;
 801  895  
 802  896                  /*
 803  897                   * !ibit => exact result was tiny before rounding,
 804  898                   * z4 nonzero => result delivered is inexact
↓ open down ↓ 2 lines elided ↑ open up ↑
 807  901                          if (z4)
 808  902                                  fsr |= FSR_UFC | FSR_NXC;
 809  903                          else if (fsr & FSR_UFM)
 810  904                                  fsr |= FSR_UFC;
 811  905                  }
 812  906          }
 813  907  
 814  908          /* restore the fsr and emulate exceptions as needed */
 815  909          if ((fsr & FSR_CEXC) & (fsr >> 23)) {
 816  910                  __fenv_setfsr32(&fsr);
      911 +
 817  912                  if (fsr & FSR_OFC) {
 818  913                          dummy = huge;
 819  914                          dummy *= huge;
 820  915                  } else if (fsr & FSR_UFC) {
 821  916                          dummy = tiny;
      917 +
 822  918                          if (fsr & FSR_NXC)
 823  919                                  dummy *= tiny;
 824  920                          else
 825  921                                  dummy -= tiny2;
 826  922                  } else {
 827  923                          dummy = huge;
 828  924                          dummy += tiny;
 829  925                  }
 830  926          } else {
 831  927                  fsr |= (fsr & 0x1f) << 5;
 832  928                  __fenv_setfsr32(&fsr);
 833  929          }
      930 +
 834  931          return (zz.q);
 835  932  }
 836      -
 837  933  #elif defined(__x86)
 838      -
 839  934  static const union {
 840  935          unsigned i[2];
 841  936          double d;
 842  937  } C[] = {
 843  938          { 0, 0x3fe00000u },
 844  939          { 0, 0x40000000u },
 845  940          { 0, 0x3df00000u },
 846  941          { 0, 0x3bf00000u },
 847  942          { 0, 0x41f00000u },
 848  943          { 0, 0x43e00000u },
 849  944          { 0, 0x7fe00000u },
 850  945          { 0, 0x00100000u },
 851  946          { 0, 0x00100001u }
 852  947  };
 853  948  
 854      -#define half    C[0].d
 855      -#define two     C[1].d
 856      -#define twom32  C[2].d
 857      -#define twom64  C[3].d
 858      -#define two32   C[4].d
 859      -#define two63   C[5].d
 860      -#define huge    C[6].d
 861      -#define tiny    C[7].d
 862      -#define tiny2   C[8].d
      949 +#define half            C[0].d
      950 +#define two             C[1].d
      951 +#define twom32          C[2].d
      952 +#define twom64          C[3].d
      953 +#define two32           C[4].d
      954 +#define two63           C[5].d
      955 +#define huge            C[6].d
      956 +#define tiny            C[7].d
      957 +#define tiny2           C[8].d
 863  958  
 864  959  #if defined(__amd64)
 865      -#define NI      4
      960 +#define NI              4
 866  961  #else
 867      -#define NI      3
      962 +#define NI              3
 868  963  #endif
 869  964  
 870  965  /*
 871  966   * fmal for x86: 80-bit extended double precision, little-endian
 872  967   */
 873  968  long double
 874      -__fmal(long double x, long double y, long double z) {
      969 +__fmal(long double x, long double y, long double z)
      970 +{
 875  971          union {
 876  972                  unsigned i[NI];
 877  973                  long double e;
 878  974          } xx, yy, zz;
      975 +
 879  976          long double xhi, yhi, xlo, ylo, t;
 880  977          unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
 881  978          unsigned oldcwsw, cwsw, rm, sticky, carry;
 882  979          int ex, ey, ez, exy, sxy, sz, e, tinyafter;
 883      -        volatile double dummy;
      980 +        volatile double dummy;
 884  981  
 885  982          /* extract the exponents of the arguments */
 886  983          xx.e = x;
 887  984          yy.e = y;
 888  985          zz.e = z;
 889  986          ex = xx.i[2] & 0x7fff;
 890  987          ey = yy.i[2] & 0x7fff;
 891  988          ez = zz.i[2] & 0x7fff;
 892  989  
 893  990          /* dispense with inf, nan, and zero cases */
 894  991          if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 ||
 895      -                (ey | yy.i[1] | yy.i[0]) == 0)  /* x or y is inf, nan, or 0 */
      992 +            (ey | yy.i[1] | yy.i[0]) == 0)      /* x or y is inf, nan, or 0 */
 896  993                  return (x * y + z);
 897  994  
 898      -        if (ez == 0x7fff)                       /* z is inf or nan */
      995 +        if (ez == 0x7fff)       /* z is inf or nan */
 899  996                  return (x + z); /* avoid spurious under/overflow in x * y */
 900  997  
 901  998          if ((ez | zz.i[1] | zz.i[0]) == 0)      /* z is zero */
 902  999                  /*
 903 1000                   * x * y isn't zero but could underflow to zero,
 904 1001                   * so don't add z, lest we perturb the sign
 905 1002                   */
 906 1003                  return (x * y);
 907 1004  
 908 1005          /*
 909 1006           * now x, y, and z are all finite and nonzero; extract signs and
 910 1007           * normalize the significands (this will raise the denormal operand
 911 1008           * exception if need be)
 912 1009           */
 913 1010          sxy = (xx.i[2] ^ yy.i[2]) & 0x8000;
 914 1011          sz = zz.i[2] & 0x8000;
     1012 +
 915 1013          if (!ex) {
 916 1014                  xx.e = x * two63;
 917 1015                  ex = (xx.i[2] & 0x7fff) - 63;
 918 1016          }
     1017 +
 919 1018          if (!ey) {
 920 1019                  yy.e = y * two63;
 921 1020                  ey = (yy.i[2] & 0x7fff) - 63;
 922 1021          }
     1022 +
 923 1023          if (!ez) {
 924 1024                  zz.e = z * two63;
 925 1025                  ez = (zz.i[2] & 0x7fff) - 63;
 926 1026          }
 927 1027  
 928 1028          /*
 929 1029           * save the control and status words, mask all exceptions, and
 930 1030           * set rounding to 64-bit precision and toward-zero
 931 1031           */
 932 1032          __fenv_getcwsw(&oldcwsw);
↓ open down ↓ 5 lines elided ↑ open up ↑
 938 1038          xx.i[2] = 0x3fff;
 939 1039          yy.i[2] = 0x3fff;
 940 1040          x = xx.e;
 941 1041          y = yy.e;
 942 1042          xhi = ((x + twom32) + two32) - two32;
 943 1043          yhi = ((y + twom32) + two32) - two32;
 944 1044          xlo = x - xhi;
 945 1045          ylo = y - yhi;
 946 1046          x *= y;
 947 1047          y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo;
     1048 +
 948 1049          if (x >= two) {
 949 1050                  x *= half;
 950 1051                  y *= half;
 951 1052                  exy++;
 952 1053          }
 953 1054  
 954 1055          /* extract the significands */
 955 1056          xx.e = x;
 956 1057          xy0 = xx.i[1];
 957 1058          xy1 = xx.i[0];
↓ open down ↓ 3 lines elided ↑ open up ↑
 961 1062          xy3 = yy.i[0];
 962 1063          xy4 = 0;
 963 1064          z0 = zz.i[1];
 964 1065          z1 = zz.i[0];
 965 1066          z2 = z3 = z4 = 0;
 966 1067  
 967 1068          /*
 968 1069           * now x*y is represented by sxy, exy, and xy[0-4], and z is
 969 1070           * represented likewise; swap if need be so |xy| <= |z|
 970 1071           */
 971      -        if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 &&
 972      -                (xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) {
 973      -                e = sxy; sxy = sz; sz = e;
 974      -                e = exy; exy = ez; ez = e;
 975      -                e = xy0; xy0 = z0; z0 = e;
 976      -                e = xy1; xy1 = z1; z1 = e;
 977      -                z2 = xy2; xy2 = 0;
 978      -                z3 = xy3; xy3 = 0;
     1072 +        if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
     1073 +            (xy1 == z1 && (xy2 | xy3) != 0)))))) {
     1074 +                e = sxy;
     1075 +                sxy = sz;
     1076 +                sz = e;
     1077 +                e = exy;
     1078 +                exy = ez;
     1079 +                ez = e;
     1080 +                e = xy0;
     1081 +                xy0 = z0;
     1082 +                z0 = e;
     1083 +                e = xy1;
     1084 +                xy1 = z1;
     1085 +                z1 = e;
     1086 +                z2 = xy2;
     1087 +                xy2 = 0;
     1088 +                z3 = xy3;
     1089 +                xy3 = 0;
 979 1090          }
 980 1091  
 981 1092          /* shift the significand of xy keeping a sticky bit */
 982 1093          e = ez - exy;
     1094 +
 983 1095          if (e > 130) {
 984 1096                  xy0 = xy1 = xy2 = xy3 = 0;
 985 1097                  xy4 = 1;
 986 1098          } else if (e >= 128) {
 987 1099                  sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e));
 988 1100                  xy4 = xy0 >> (e - 128);
     1101 +
 989 1102                  if (sticky)
 990 1103                          xy4 |= 1;
     1104 +
 991 1105                  xy0 = xy1 = xy2 = xy3 = 0;
 992 1106          } else if (e >= 96) {
 993 1107                  sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e));
 994 1108                  xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
     1109 +
 995 1110                  if (sticky)
 996 1111                          xy4 |= 1;
     1112 +
 997 1113                  xy3 = xy0 >> (e - 96);
 998 1114                  xy0 = xy1 = xy2 = 0;
 999 1115          } else if (e >= 64) {
1000 1116                  sticky = xy3 | ((xy2 << 1) << (95 - e));
1001 1117                  xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
     1118 +
1002 1119                  if (sticky)
1003 1120                          xy4 |= 1;
     1121 +
1004 1122                  xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
1005 1123                  xy2 = xy0 >> (e - 64);
1006 1124                  xy0 = xy1 = 0;
1007 1125          } else if (e >= 32) {
1008 1126                  sticky = (xy3 << 1) << (63 - e);
1009 1127                  xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
     1128 +
1010 1129                  if (sticky)
1011 1130                          xy4 |= 1;
     1131 +
1012 1132                  xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
1013 1133                  xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
1014 1134                  xy1 = xy0 >> (e - 32);
1015 1135                  xy0 = 0;
1016 1136          } else if (e) {
1017 1137                  xy4 = (xy3 << 1) << (31 - e);
1018 1138                  xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
1019 1139                  xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
1020 1140                  xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
1021 1141                  xy0 >>= e;
1022 1142          }
1023 1143  
1024 1144          /* if this is a magnitude subtract, negate the significand of xy */
1025 1145          if (sxy ^ sz) {
1026 1146                  xy0 = ~xy0;
1027 1147                  xy1 = ~xy1;
1028 1148                  xy2 = ~xy2;
1029 1149                  xy3 = ~xy3;
1030 1150                  xy4 = -xy4;
     1151 +
1031 1152                  if (xy4 == 0)
1032 1153                          if (++xy3 == 0)
1033 1154                                  if (++xy2 == 0)
1034 1155                                          if (++xy1 == 0)
1035 1156                                                  xy0++;
1036 1157          }
1037 1158  
1038 1159          /* add, propagating carries */
1039 1160          z4 += xy4;
1040 1161          carry = (z4 < xy4);
1041 1162          z3 += xy3;
     1163 +
1042 1164          if (carry) {
1043 1165                  z3++;
1044 1166                  carry = (z3 <= xy3);
1045      -        } else
     1167 +        } else {
1046 1168                  carry = (z3 < xy3);
     1169 +        }
     1170 +
1047 1171          z2 += xy2;
     1172 +
1048 1173          if (carry) {
1049 1174                  z2++;
1050 1175                  carry = (z2 <= xy2);
1051      -        } else
     1176 +        } else {
1052 1177                  carry = (z2 < xy2);
     1178 +        }
     1179 +
1053 1180          z1 += xy1;
     1181 +
1054 1182          if (carry) {
1055 1183                  z1++;
1056 1184                  carry = (z1 <= xy1);
1057      -        } else
     1185 +        } else {
1058 1186                  carry = (z1 < xy1);
     1187 +        }
     1188 +
1059 1189          z0 += xy0;
     1190 +
1060 1191          if (carry) {
1061 1192                  z0++;
1062 1193                  carry = (z0 <= xy0);
1063      -        } else
     1194 +        } else {
1064 1195                  carry = (z0 < xy0);
     1196 +        }
1065 1197  
1066 1198          /* for a magnitude subtract, ignore the last carry out */
1067 1199          if (sxy ^ sz)
1068 1200                  carry = 0;
1069 1201  
1070 1202          /* postnormalize and collect rounding information into z2 */
1071 1203          if (ez < 1) {
1072 1204                  /* result is tiny; shift right until exponent is within range */
1073 1205                  e = 1 - ez;
     1206 +
1074 1207                  if (e > 67) {
1075      -                        z2 = 1; /* result can't be exactly zero */
     1208 +                        z2 = 1;         /* result can't be exactly zero */
1076 1209                          z0 = z1 = 0;
1077 1210                  } else if (e >= 64) {
1078 1211                          sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e));
1079 1212                          z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e));
     1213 +
1080 1214                          if (sticky)
1081 1215                                  z2 |= 1;
     1216 +
1082 1217                          z1 = carry >> (e - 64);
1083 1218                          z0 = 0;
1084 1219                  } else if (e >= 32) {
1085 1220                          sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e));
1086 1221                          z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
     1222 +
1087 1223                          if (sticky)
1088 1224                                  z2 |= 1;
     1225 +
1089 1226                          z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e));
1090 1227                          z0 = carry >> (e - 32);
1091 1228                  } else {
1092 1229                          sticky = z4 | z3 | (z2 << 1) << (31 - e);
1093 1230                          z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
     1231 +
1094 1232                          if (sticky)
1095 1233                                  z2 |= 1;
     1234 +
1096 1235                          z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
1097 1236                          z0 = (z0 >> e) | ((carry << 1) << (31 - e));
1098 1237                  }
     1238 +
1099 1239                  ez = 1;
1100 1240          } else if (carry) {
1101 1241                  /* carry out; shift right by one */
1102 1242                  sticky = (z2 & 1) | z3 | z4;
1103 1243                  z2 = (z2 >> 1) | (z1 << 31);
     1244 +
1104 1245                  if (sticky)
1105 1246                          z2 |= 1;
     1247 +
1106 1248                  z1 = (z1 >> 1) | (z0 << 31);
1107 1249                  z0 = (z0 >> 1) | 0x80000000;
1108 1250                  ez++;
1109 1251          } else {
1110 1252                  if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) {
1111 1253                          /*
1112 1254                           * borrow/cancellation; shift left as much as
1113 1255                           * exponent allows
1114 1256                           */
1115 1257                          while (!z0 && ez >= 33) {
1116 1258                                  z0 = z1;
1117 1259                                  z1 = z2;
1118 1260                                  z2 = z3;
1119 1261                                  z3 = z4;
1120 1262                                  z4 = 0;
1121 1263                                  ez -= 32;
1122 1264                          }
     1265 +
1123 1266                          while (z0 < 0x80000000u && ez > 1) {
1124 1267                                  z0 = (z0 << 1) | (z1 >> 31);
1125 1268                                  z1 = (z1 << 1) | (z2 >> 31);
1126 1269                                  z2 = (z2 << 1) | (z3 >> 31);
1127 1270                                  z3 = (z3 << 1) | (z4 >> 31);
1128 1271                                  z4 <<= 1;
1129 1272                                  ez--;
1130 1273                          }
1131 1274                  }
     1275 +
1132 1276                  if (z3 | z4)
1133 1277                          z2 |= 1;
1134 1278          }
1135 1279  
1136 1280          /* get the rounding mode */
1137 1281          rm = oldcwsw & 0x0c000000;
1138 1282  
1139 1283          /* adjust exponent if result is subnormal */
1140 1284          tinyafter = 0;
     1285 +
1141 1286          if (!(z0 & 0x80000000)) {
1142 1287                  ez = 0;
1143 1288                  tinyafter = 1;
1144      -                if (!(z0 | z1 | z2)) { /* exact zero */
     1289 +
     1290 +                if (!(z0 | z1 | z2)) {  /* exact zero */
1145 1291                          zz.i[2] = rm == FCW_RM ? 0x8000 : 0;
1146 1292                          zz.i[1] = zz.i[0] = 0;
1147 1293                          __fenv_setcwsw(&oldcwsw);
1148 1294                          return (zz.e);
1149 1295                  }
1150 1296          }
1151 1297  
1152 1298          /*
1153 1299           * flip the sense of directed roundings if the result is negative;
1154 1300           * the logic below applies to a positive result
1155 1301           */
1156 1302          if (sz && (rm == FCW_RM || rm == FCW_RP))
1157 1303                  rm = (FCW_RM + FCW_RP) - rm;
1158 1304  
1159 1305          /* round */
1160 1306          if (z2) {
1161      -                if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u ||
1162      -                        (z2 == 0x80000000u && (z1 & 1))))) {
     1307 +                if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u || (z2 ==
     1308 +                    0x80000000u && (z1 & 1))))) {
1163 1309                          /* round up and renormalize if necessary */
1164 1310                          if (++z1 == 0) {
1165 1311                                  if (++z0 == 0) {
1166 1312                                          z0 = 0x80000000;
1167 1313                                          ez++;
1168 1314                                  } else if (z0 == 0x80000000) {
1169 1315                                          /* rounded up to smallest normal */
1170 1316                                          ez = 1;
     1317 +
1171 1318                                          if ((rm == FCW_RP && z2 >
1172      -                                                0x80000000u) || (rm == FCW_RN &&
1173      -                                                z2 >= 0xc0000000u))
     1319 +                                            0x80000000u) || (rm == FCW_RN &&
     1320 +                                            z2 >= 0xc0000000u))
1174 1321                                                  /*
1175 1322                                                   * would have rounded up to
1176 1323                                                   * smallest normal even with
1177 1324                                                   * unbounded range
1178 1325                                                   */
1179 1326                                                  tinyafter = 0;
1180 1327                                  }
1181 1328                          }
1182 1329                  }
1183 1330          }
1184 1331  
1185 1332          /* restore the control and status words, check for over/underflow */
1186 1333          __fenv_setcwsw(&oldcwsw);
     1334 +
1187 1335          if (ez >= 0x7fff) {
1188 1336                  if (rm == FCW_RN || rm == FCW_RP) {
1189 1337                          zz.i[2] = sz | 0x7fff;
1190 1338                          zz.i[1] = 0x80000000;
1191 1339                          zz.i[0] = 0;
1192 1340                  } else {
1193 1341                          zz.i[2] = sz | 0x7ffe;
1194 1342                          zz.i[1] = 0xffffffff;
1195 1343                          zz.i[0] = 0xffffffff;
1196 1344                  }
     1345 +
1197 1346                  dummy = huge;
1198 1347                  dummy *= huge;
1199 1348          } else {
1200 1349                  zz.i[2] = sz | ez;
1201 1350                  zz.i[1] = z0;
1202 1351                  zz.i[0] = z1;
1203 1352  
1204 1353                  /*
1205 1354                   * tinyafter => result rounded w/ unbounded range would be tiny,
1206 1355                   * z2 nonzero => result delivered is inexact
1207 1356                   */
1208 1357                  if (tinyafter) {
1209 1358                          dummy = tiny;
     1359 +
1210 1360                          if (z2)
1211 1361                                  dummy *= tiny;
1212 1362                          else
1213 1363                                  dummy -= tiny2;
1214 1364                  } else if (z2) {
1215 1365                          dummy = huge;
1216 1366                          dummy += tiny;
1217 1367                  }
1218 1368          }
1219 1369  
1220 1370          return (zz.e);
1221 1371  }
1222      -
1223 1372  #else
1224 1373  #error Unknown architecture
1225 1374  #endif
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX