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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/Q/sqrtl.c
          +++ new/usr/src/lib/libm/common/Q/sqrtl.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 __sqrtl = sqrtl
  31   32  
  32   33  #include "libm.h"
  33   34  #include "longdouble.h"
  34   35  
↓ open down ↓ 2 lines elided ↑ open up ↑
  37   38  extern enum fp_direction_type __swapRD(enum fp_direction_type);
  38   39  
  39   40  /*
  40   41   * in struct longdouble, msw consists of
  41   42   *      unsigned short  sgn:1;
  42   43   *      unsigned short  exp:15;
  43   44   *      unsigned short  frac1:16;
  44   45   */
  45   46  
  46   47  #ifdef __LITTLE_ENDIAN
  47      -
  48   48  /* array indices used to access words within a double */
  49      -#define HIWORD  1
  50      -#define LOWORD  0
       49 +#define HIWORD          1
       50 +#define LOWORD          0
  51   51  
  52   52  /* structure used to access words within a quad */
  53   53  union longdouble {
  54   54          struct {
  55      -                unsigned int    frac4;
  56      -                unsigned int    frac3;
  57      -                unsigned int    frac2;
  58      -                unsigned int    msw;
       55 +                unsigned int frac4;
       56 +                unsigned int frac3;
       57 +                unsigned int frac2;
       58 +                unsigned int msw;
  59   59          } l;
  60      -        long double     d;
       60 +        long double d;
  61   61  };
  62   62  
  63   63  /* default NaN returned for sqrt(neg) */
  64      -static const union longdouble
  65      -        qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
       64 +static const union longdouble qnan = {
       65 +        0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff
       66 +};
  66   67  
  67   68  /* signalling NaN used to raise invalid */
  68   69  static const union {
  69   70          unsigned u[2];
  70   71          double d;
  71   72  } snan = { 0, 0x7ff00001 };
  72      -
  73   73  #else
  74      -
  75   74  /* array indices used to access words within a double */
  76      -#define HIWORD  0
  77      -#define LOWORD  1
       75 +#define HIWORD          0
       76 +#define LOWORD          1
  78   77  
  79   78  /* structure used to access words within a quad */
  80   79  union longdouble {
  81   80          struct {
  82      -                unsigned int    msw;
  83      -                unsigned int    frac2;
  84      -                unsigned int    frac3;
  85      -                unsigned int    frac4;
       81 +                unsigned int msw;
       82 +                unsigned int frac2;
       83 +                unsigned int frac3;
       84 +                unsigned int frac4;
  86   85          } l;
  87      -        long double     d;
       86 +        long double d;
  88   87  };
  89   88  
  90   89  /* default NaN returned for sqrt(neg) */
  91      -static const union longdouble
  92      -        qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
       90 +static const union longdouble qnan = {
       91 +        0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff
       92 +};
  93   93  
  94   94  /* signalling NaN used to raise invalid */
  95   95  static const union {
  96   96          unsigned u[2];
  97   97          double d;
  98   98  } snan = { 0x7ff00001, 0 };
  99      -
 100   99  #endif /* __LITTLE_ENDIAN */
 101  100  
 102      -
 103      -static const double
 104      -        zero = 0.0,
      101 +static const double zero = 0.0,
 105  102          half = 0.5,
 106  103          one = 1.0,
 107  104          huge = 1.0e300,
 108  105          tiny = 1.0e-300,
 109  106          two36 = 6.87194767360000000000e+10,
 110  107          two30 = 1.07374182400000000000e+09,
 111  108          two6 = 6.40000000000000000000e+01,
 112  109          two4 = 1.60000000000000000000e+01,
 113  110          twom18 = 3.81469726562500000000e-06,
 114  111          twom28 = 3.72529029846191406250e-09,
 115  112          twom42 = 2.27373675443232059479e-13,
 116  113          twom60 = 8.67361737988403547206e-19,
 117  114          twom62 = 2.16840434497100886801e-19,
 118  115          twom66 = 1.35525271560688054251e-20,
 119  116          twom90 = 8.07793566946316088742e-28,
 120  117          twom113 = 9.62964972193617926528e-35,
 121  118          twom124 = 4.70197740328915003187e-38;
 122  119  
 123      -
 124  120  /*
 125      -*       Extract the exponent and normalized significand (represented as
 126      -*       an array of five doubles) from a finite, nonzero quad.
 127      -*/
      121 + *      Extract the exponent and normalized significand (represented as
      122 + *      an array of five doubles) from a finite, nonzero quad.
      123 + */
 128  124  static int
 129  125  __q_unpack(const union longdouble *x, double *s)
 130  126  {
 131  127          union {
 132      -                double                  d;
 133      -                unsigned int    l[2];
      128 +                double d;
      129 +                unsigned int l[2];
 134  130          } u;
 135      -        double                  b;
 136      -        unsigned int    lx, w[3];
 137      -        int                             ex;
      131 +
      132 +        double b;
      133 +        unsigned int lx, w[3];
      134 +        int ex;
 138  135  
 139  136          /* get the normalized significand and exponent */
 140      -        ex = (int) ((x->l.msw & 0x7fffffff) >> 16);
      137 +        ex = (int)((x->l.msw & 0x7fffffff) >> 16);
 141  138          lx = x->l.msw & 0xffff;
 142      -        if (ex)
 143      -        {
      139 +
      140 +        if (ex) {
 144  141                  lx |= 0x10000;
 145  142                  w[0] = x->l.frac2;
 146  143                  w[1] = x->l.frac3;
 147  144                  w[2] = x->l.frac4;
 148      -        }
 149      -        else
 150      -        {
 151      -                if (lx | (x->l.frac2 & 0xfffe0000))
 152      -                {
      145 +        } else {
      146 +                if (lx | (x->l.frac2 & 0xfffe0000)) {
 153  147                          w[0] = x->l.frac2;
 154  148                          w[1] = x->l.frac3;
 155  149                          w[2] = x->l.frac4;
 156  150                          ex = 1;
 157      -                }
 158      -                else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000))
 159      -                {
      151 +                } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
 160  152                          lx = x->l.frac2;
 161  153                          w[0] = x->l.frac3;
 162  154                          w[1] = x->l.frac4;
 163  155                          w[2] = 0;
 164  156                          ex = -31;
 165      -                }
 166      -                else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000))
 167      -                {
      157 +                } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
 168  158                          lx = x->l.frac3;
 169  159                          w[0] = x->l.frac4;
 170  160                          w[1] = w[2] = 0;
 171  161                          ex = -63;
 172      -                }
 173      -                else
 174      -                {
      162 +                } else {
 175  163                          lx = x->l.frac4;
 176  164                          w[0] = w[1] = w[2] = 0;
 177  165                          ex = -95;
 178  166                  }
 179      -                while ((lx & 0x10000) == 0)
 180      -                {
      167 +
      168 +                while ((lx & 0x10000) == 0) {
 181  169                          lx = (lx << 1) | (w[0] >> 31);
 182  170                          w[0] = (w[0] << 1) | (w[1] >> 31);
 183  171                          w[1] = (w[1] << 1) | (w[2] >> 31);
 184  172                          w[2] <<= 1;
 185  173                          ex--;
 186  174                  }
 187  175          }
 188  176  
 189  177          /* extract the significand into five doubles */
 190  178          u.l[HIWORD] = 0x42300000;
↓ open down ↓ 21 lines elided ↑ open up ↑
 212  200          u.l[HIWORD] |= w[1] & 0xffff;
 213  201          u.l[LOWORD] = w[2] & 0xff000000;
 214  202          s[3] = u.d - b;
 215  203  
 216  204          u.l[HIWORD] = 0x3c300000;
 217  205          u.l[LOWORD] = 0;
 218  206          b = u.d;
 219  207          u.l[LOWORD] = w[2] & 0xffffff;
 220  208          s[4] = u.d - b;
 221  209  
 222      -        return ex - 0x3fff;
      210 +        return (ex - 0x3fff);
 223  211  }
 224  212  
 225      -
 226  213  /*
 227      -*       Pack an exponent and array of three doubles representing a finite,
 228      -*       nonzero number into a quad.  Assume the sign is already there and
 229      -*       the rounding mode has been fudged accordingly.
 230      -*/
      214 + *      Pack an exponent and array of three doubles representing a finite,
      215 + *      nonzero number into a quad.  Assume the sign is already there and
      216 + *      the rounding mode has been fudged accordingly.
      217 + */
 231  218  static void
 232  219  __q_pack(const double *z, int exp, enum fp_direction_type rm,
 233      -        union longdouble *x, int *inexact)
      220 +    union longdouble *x, int *inexact)
 234  221  {
 235  222          union {
 236      -                double                  d;
 237      -                unsigned int    l[2];
      223 +                double d;
      224 +                unsigned int l[2];
 238  225          } u;
 239      -        double                  s[3], t, t2;
 240      -        unsigned int    msw, frac2, frac3, frac4;
      226 +
      227 +        double s[3], t, t2;
      228 +        unsigned int msw, frac2, frac3, frac4;
 241  229  
 242  230          /* bias exponent and strip off integer bit */
 243  231          exp += 0x3fff;
 244  232          s[0] = z[0] - one;
 245  233          s[1] = z[1];
 246  234          s[2] = z[2];
 247  235  
 248  236          /*
 249  237           * chop the significand to obtain the fraction;
 250  238           * use round-to-minus-infinity to ensure chopping
↓ open down ↓ 24 lines elided ↑ open up ↑
 275  263          frac4 = u.l[LOWORD];
 276  264          s[0] -= (u.d - twom60);
 277  265  
 278  266          /*
 279  267           * keep track of what's left for rounding; note that
 280  268           * t2 will be non-negative due to rounding mode
 281  269           */
 282  270          t = s[0] + s[1];
 283  271          t2 = (s[0] - t) + s[1];
 284  272  
 285      -        if (t != zero)
 286      -        {
      273 +        if (t != zero) {
 287  274                  *inexact = 1;
 288  275  
 289  276                  /* decide whether to round the fraction up */
 290  277                  if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
 291      -                        (t == twom113 && (t2 != zero || frac4 & 1)))))
 292      -                {
      278 +                    (t == twom113 && (t2 != zero || frac4 & 1))))) {
 293  279                          /* round up and renormalize if necessary */
 294      -                        if (++frac4 == 0)
 295      -                                if (++frac3 == 0)
 296      -                                        if (++frac2 == 0)
 297      -                                                if (++msw == 0x10000)
 298      -                                                {
      280 +                        if (++frac4 == 0) {
      281 +                                if (++frac3 == 0) {
      282 +                                        if (++frac2 == 0) {
      283 +                                                if (++msw == 0x10000) {
 299  284                                                          msw = 0;
 300  285                                                          exp++;
 301  286                                                  }
      287 +                                        }
      288 +                                }
      289 +                        }
 302  290                  }
 303  291          }
 304  292  
 305  293          /* assemble the result */
 306  294          x->l.msw |= msw | (exp << 16);
 307  295          x->l.frac2 = frac2;
 308  296          x->l.frac3 = frac3;
 309  297          x->l.frac4 = frac4;
 310  298  }
 311  299  
 312      -
 313  300  /*
 314      -*       Compute the square root of x and place the TP result in s.
 315      -*/
      301 + *      Compute the square root of x and place the TP result in s.
      302 + */
 316  303  static void
 317  304  __q_tp_sqrt(const double *x, double *s)
 318  305  {
 319      -        double  c, rr, r[3], tt[3], t[5];
      306 +        double c, rr, r[3], tt[3], t[5];
 320  307  
 321  308          /* approximate the divisor for the Newton iteration */
 322  309          c = sqrt((x[0] + x[1]) + x[2]);
 323  310          rr = half / c;
 324  311  
 325  312          /* compute the first five "digits" of the square root */
 326  313          t[0] = (c + two30) - two30;
 327  314          tt[0] = t[0] + t[0];
 328  315          r[0] = ((x[0] - t[0] * t[0]) + x[1]) + x[2];
 329  316  
↓ open down ↓ 16 lines elided ↑ open up ↑
 346  333          t[3] = (rr * (r[0] + r[1]) + twom42) - twom42;
 347  334          r[0] = ((r[0] - tt[0] * t[3]) + r[1]) - tt[1] * t[3];
 348  335          r[1] = -tt[2] * t[3];
 349  336          c = (r[1] + twom90) - twom90;
 350  337          r[0] += c;
 351  338          r[1] = (r[1] - c) - t[3] * t[3];
 352  339  
 353  340          t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
 354  341  
 355  342          /* here we just need to get the sign of the remainder */
 356      -        c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1])
 357      -                - tt[2] * t[4]) - (t[3] + t[3]) * t[4]) - t[4] * t[4];
      343 +        c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1]) - tt[2] * t[4]) -
      344 +            (t[3] + t[3]) * t[4]) - t[4] * t[4];
 358  345  
 359  346          /* reduce to three doubles */
 360  347          t[0] += t[1];
 361  348          t[1] = t[2] + t[3];
 362  349          t[2] = t[4];
 363  350  
 364  351          /* if the third term might lie on a rounding boundary, perturb it */
 365      -        if (c != zero && t[2] == (twom62 + t[2]) - twom62)
 366      -        {
      352 +        if (c != zero && t[2] == (twom62 + t[2]) - twom62) {
 367  353                  if (c < zero)
 368  354                          t[2] -= twom124;
 369  355                  else
 370  356                          t[2] += twom124;
 371  357          }
 372  358  
 373  359          /* condense the square root */
 374  360          c = t[1] + t[2];
 375  361          t[2] += (t[1] - c);
 376  362          t[1] = c;
 377  363          c = t[0] + t[1];
 378  364          s[1] = t[1] + (t[0] - c);
 379  365          s[0] = c;
 380      -        if (s[1] == zero)
 381      -        {
      366 +
      367 +        if (s[1] == zero) {
 382  368                  c = s[0] + t[2];
 383  369                  s[1] = t[2] + (s[0] - c);
 384  370                  s[0] = c;
 385  371                  s[2] = zero;
 386      -        }
 387      -        else
 388      -        {
      372 +        } else {
 389  373                  c = s[1] + t[2];
 390  374                  s[2] = t[2] + (s[1] - c);
 391  375                  s[1] = c;
 392  376          }
 393  377  }
 394  378  
 395      -
 396  379  long double
 397  380  sqrtl(long double ldx)
 398  381  {
 399      -        union   longdouble              x;
 400      -        volatile double                 t;
 401      -        double                                  xx[5], zz[3];
 402      -        enum fp_direction_type  rm;
 403      -        int                             ex, inexact, exc, traps;
      382 +        union   longdouble x;
      383 +        volatile double t;
      384 +        double xx[5], zz[3];
      385 +        enum fp_direction_type rm;
      386 +        int ex, inexact, exc, traps;
 404  387  
 405  388          /* clear cexc */
 406  389          t = zero;
 407  390          t -= zero;
 408  391  
 409  392          /* check for zero operand */
 410  393          x.d = ldx;
      394 +
 411  395          if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
 412      -                return ldx;
      396 +                return (ldx);
 413  397  
 414  398          /* handle nan and inf cases */
 415      -        if ((x.l.msw & 0x7fffffff) >= 0x7fff0000)
 416      -        {
 417      -                if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4)
 418      -                {
 419      -                        if (!(x.l.msw & 0x8000))
 420      -                        {
      399 +        if ((x.l.msw & 0x7fffffff) >= 0x7fff0000) {
      400 +                if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4) {
      401 +                        if (!(x.l.msw & 0x8000)) {
 421  402                                  /* snan, signal invalid */
 422  403                                  t += snan.d;
 423  404                          }
      405 +
 424  406                          x.l.msw |= 0x8000;
 425      -                        return x.d;
      407 +                        return (x.d);
 426  408                  }
 427      -                if (x.l.msw & 0x80000000)
 428      -                {
      409 +
      410 +                if (x.l.msw & 0x80000000) {
 429  411                          /* sqrt(-inf), signal invalid */
 430  412                          t = -one;
 431  413                          t = sqrt(t);
 432      -                        return qnan.d;
      414 +                        return (qnan.d);
 433  415                  }
      416 +
 434  417                  /* sqrt(inf), return inf */
 435      -                return x.d;
      418 +                return (x.d);
 436  419          }
 437  420  
 438  421          /* handle negative numbers */
 439      -        if (x.l.msw & 0x80000000)
 440      -        {
      422 +        if (x.l.msw & 0x80000000) {
 441  423                  t = -one;
 442  424                  t = sqrt(t);
 443      -                return qnan.d;
      425 +                return (qnan.d);
 444  426          }
 445  427  
 446  428          /* now x is finite, positive */
 447  429  
 448  430          traps = __swapTE(0);
 449  431          exc = __swapEX(0);
 450  432          rm = __swapRD(fp_nearest);
 451  433  
 452  434          ex = __q_unpack(&x, xx);
 453      -        if (ex & 1)
 454      -        {
      435 +
      436 +        if (ex & 1) {
 455  437                  /* make exponent even */
 456  438                  xx[0] += xx[0];
 457  439                  xx[1] += xx[1];
 458  440                  xx[2] += xx[2];
 459  441                  xx[3] += xx[3];
 460  442                  xx[4] += xx[4];
 461  443                  ex--;
 462  444          }
      445 +
 463  446          __q_tp_sqrt(xx, zz);
 464  447  
 465  448          /* put everything together */
 466  449          x.l.msw = 0;
 467  450          inexact = 0;
 468  451          __q_pack(zz, ex >> 1, rm, &x, &inexact);
 469  452  
 470  453          (void) __swapRD(rm);
 471  454          (void) __swapEX(exc);
 472  455          (void) __swapTE(traps);
 473      -        if (inexact)
 474      -        {
      456 +
      457 +        if (inexact) {
 475  458                  t = huge;
 476  459                  t += tiny;
 477  460          }
 478      -        return x.d;
      461 +
      462 +        return (x.d);
 479  463  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX