Print this page
5262 libm needs to be carefully unifdef'd
5268 libm doesn't need to hide symbols which are already local

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/C/sqrt.c
          +++ new/usr/src/lib/libm/common/C/sqrt.c
↓ open down ↓ 22 lines elided ↑ open up ↑
  23   23   */
  24   24  /*
  25   25   * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26   26   * Use is subject to license terms.
  27   27   */
  28   28  
  29   29  #pragma weak sqrt = __sqrt
  30   30  
  31   31  #include "libm.h"
  32   32  
  33      -#ifdef __INLINE
  34   33  
  35   34  extern double __inline_sqrt(double);
  36   35  
  37   36  double
  38   37  sqrt(double x) {
  39   38          double  z = __inline_sqrt(x);
  40   39  
  41   40          if (isnan(x))
  42   41                  return (z);
  43   42          return ((x < 0.0)? _SVID_libm_err(x, x, 26) : z);
  44   43  }
  45   44  
  46      -#else   /* defined(__INLINE) */
  47      -
  48      -/*
  49      - * Warning: This correctly rounded sqrt is extremely slow because it computes
  50      - * the sqrt bit by bit using integer arithmetic.
  51      - */
  52      -
  53      -static const double big = 1.0e30, small = 1.0e-30;
  54      -
  55      -double
  56      -sqrt(double x)
  57      -{
  58      -        double          z;
  59      -        unsigned        r, t1, s1, ix1, q1;
  60      -        int             ix0, s0, j, q, m, n, t;
  61      -        int              *px = (int *)&x, *pz = (int *)&z;
  62      -
  63      -        ix0 = px[HIWORD];
  64      -        ix1 = px[LOWORD];
  65      -        if ((ix0 & 0x7ff00000) == 0x7ff00000) { /* x is inf or NaN */
  66      -                if (ix0 == 0xfff00000 && ix1 == 0)
  67      -                        return (_SVID_libm_err(x, x, 26));
  68      -                return (x + x);
  69      -        }
  70      -        if (((ix0 & 0x7fffffff) | ix1) == 0)    /* x is zero */
  71      -                return (x);
  72      -
  73      -        /* extract exponent and significand */
  74      -        m = ilogb(x);
  75      -        z = scalbn(x, -m);
  76      -        ix0 = (pz[HIWORD] & 0x000fffff) | 0x00100000;
  77      -        ix1 = pz[LOWORD];
  78      -        n = m >> 1;
  79      -        if (n + n != m) {
  80      -                ix0 = (ix0 << 1) | (ix1 >> 31);
  81      -                ix1 <<= 1;
  82      -                m -= 1;
  83      -        }
  84      -
  85      -        /* generate sqrt(x) bit by bit */
  86      -        ix0 = (ix0 << 1) | (ix1 >> 31);
  87      -        ix1 <<= 1;
  88      -        q = q1 = s0 = s1 = 0;
  89      -        r = 0x00200000;
  90      -
  91      -        for (j = 1; j <= 22; j++) {
  92      -                t = s0 + r;
  93      -                if (t <= ix0) {
  94      -                        s0 = t + r;
  95      -                        ix0 -= t;
  96      -                        q += r;
  97      -                }
  98      -                ix0 = (ix0 << 1) | (ix1 >> 31);
  99      -                ix1 <<= 1;
 100      -                r >>= 1;
 101      -        }
 102      -
 103      -        r = 0x80000000;
 104      -        for (j = 1; j <= 32; j++) {
 105      -                t1 = s1 + r;
 106      -                t = s0;
 107      -                if (t < ix0 || (t == ix0 && t1 <= ix1)) {
 108      -                        s1 = t1 + r;
 109      -                        if ((t1 & 0x80000000) == 0x80000000 &&
 110      -                            (s1 & 0x80000000) == 0)
 111      -                                s0 += 1;
 112      -                        ix0 -= t;
 113      -                        if (ix1 < t1)
 114      -                                ix0 -= 1;
 115      -                        ix1 -= t1;
 116      -                        q1 += r;
 117      -                }
 118      -                ix0 = (ix0 << 1) | (ix1 >> 31);
 119      -                ix1 <<= 1;
 120      -                r >>= 1;
 121      -        }
 122      -
 123      -        /* round */
 124      -        if ((ix0 | ix1) == 0)
 125      -                goto done;
 126      -        z = big - small;        /* trigger inexact flag */
 127      -        if (z < big)
 128      -                goto done;
 129      -        if (q1 == 0xffffffff) {
 130      -                q1 = 0;
 131      -                q += 1;
 132      -                goto done;
 133      -        }
 134      -        z = big + small;
 135      -        if (z > big) {
 136      -                if (q1 == 0xfffffffe)
 137      -                        q += 1;
 138      -                q1 += 2;
 139      -                goto done;
 140      -        }
 141      -        q1 += (q1 & 1);
 142      -done:
 143      -        pz[HIWORD] = (q >> 1) + 0x3fe00000;
 144      -        pz[LOWORD] = q1 >> 1;
 145      -        if ((q & 1) == 1)
 146      -                pz[LOWORD] |= 0x80000000;
 147      -        return (scalbn(z, n));
 148      -}
 149      -
 150      -#endif  /* defined(__INLINE) */
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX