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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/C/atan.c
          +++ new/usr/src/lib/libm/common/C/atan.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 __atan = atan
  31   32  
  32      -/* INDENT OFF */
       33 +
  33   34  /*
  34   35   * atan(x)
  35   36   * Accurate Table look-up algorithm with polynomial approximation in
  36   37   * partially product form.
  37   38   *
  38   39   * -- K.C. Ng, October 17, 2004
  39   40   *
  40   41   * Algorithm
  41   42   *
  42   43   * (1). Purge off Inf and NaN and 0
↓ open down ↓ 2 lines elided ↑ open up ↑
  45   46   *      (2.1) if x < 2^(-prec/2), atan(x) = x  with inexact flag raised
  46   47   *      (2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x)
  47   48   *      (2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5)
  48   49   *      (2.4) Otherwise
  49   50   *              atan(x) = poly1(x) = x + A * B,
  50   51   *      where
  51   52   *              A = (p1*x*z) * (p2+z(p3+z))
  52   53   *              B = (p4+z)+z*z) * (p5+z(p6+z))
  53   54   *      Note: (i) domain of poly1 is [0, 1/8], (ii) remez relative
  54   55   *      approximation error of poly1 is bounded by
  55      - *              |(atan(x)-poly1(x))/x| <= 2^-57.61
       56 + *              |(atan(x)-poly1(x))/x| <= 2^-57.61
  56   57   * (4). For x >= 8 then
  57   58   *      (3.1) if x >= 2^prec,     atan(x) = atan(inf) - pio2lo
  58   59   *      (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x
  59   60   *      (3.3) if x <= 65,         atan(x) = atan(inf) - poly1(1/x)
  60   61   *      (3.4) otherwise           atan(x) = atan(inf) - poly2(1/x)
  61   62   *      where
  62   63   *              poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z),
  63   64   *      its domain is [0, 0.0154]; and its remez absolute
  64   65   *      approximation error is bounded by
  65   66   *              |atan(x)-poly2(x)|<= 2^-59.45
↓ open down ↓ 6 lines elided ↑ open up ↑
  72   73   *              atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j]))
  73   74   *      where y[j] are carefully chosen so that it matches x to around 4.5
  74   75   *      bits and at the same time atan(y[j]) is very close to an IEEE double
  75   76   *      floating point number. Calculation indicates that
  76   77   *              max|(x-y[j])/(1+x*y[j])| < 0.0154
  77   78   *              j,x
  78   79   *
  79   80   * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing
  80   81   * more than 10 million random arguments
  81   82   */
  82      -/* INDENT ON */
  83   83  
  84   84  #include "libm.h"
  85   85  #include "libm_protos.h"
  86   86  
  87   87  extern const double _TBL_atan[];
       88 +
  88   89  static const double g[] = {
  89      -/* one  = */  1.0,
  90      -/* p1   = */  8.02176624254765935351230154992663301527500152588e-0002,
  91      -/* p2   = */  1.27223421700559402580665846471674740314483642578e+0000,
  92      -/* p3   = */ -1.20606901800503640842521235754247754812240600586e+0000,
  93      -/* p4   = */ -2.36088967922325565496066701598465442657470703125e+0000,
  94      -/* p5   = */  1.38345799501389166152875986881554126739501953125e+0000,
  95      -/* p6   = */  1.06742368078953453469637224770849570631980895996e+0000,
       90 +/* one  = */
       91 +        1.0,
       92 +/* p1   = */8.02176624254765935351230154992663301527500152588e-0002,
       93 +/* p2   = */1.27223421700559402580665846471674740314483642578e+0000,
       94 +/* p3   = */-1.20606901800503640842521235754247754812240600586e+0000,
       95 +/* p4   = */-2.36088967922325565496066701598465442657470703125e+0000,
       96 +/* p5   = */1.38345799501389166152875986881554126739501953125e+0000,
       97 +/* p6   = */1.06742368078953453469637224770849570631980895996e+0000,
  96   98  /* q1   = */ -1.42796626333911796935538518482644576579332351685e-0001,
  97      -/* q2   = */  3.51427110447873227059810477159863497078605962912e+0000,
  98      -/* q3   = */  5.92129112708164262457444237952586263418197631836e-0001,
       99 +/* q2   = */ 3.51427110447873227059810477159863497078605962912e+0000,
      100 +/* q3   = */ 5.92129112708164262457444237952586263418197631836e-0001,
  99  101  /* q4   = */ -1.99272234785683144409063061175402253866195678711e+0000,
 100      -/* pio2hi */  1.570796326794896558e+00,
 101      -/* pio2lo */  6.123233995736765886e-17,
      102 +/* pio2hi */ 1.570796326794896558e+00,
      103 +/* pio2lo */ 6.123233995736765886e-17,
 102  104  /* t1   = */ -0.333333333333333333333333333333333,
 103      -/* t2   = */  0.2,
      105 +/* t2   = */ 0.2,
 104  106  /* t3   = */ -1.666666666666666666666666666666666,
 105  107  };
 106  108  
 107      -#define one g[0]
 108      -#define p1 g[1]
 109      -#define p2 g[2]
 110      -#define p3 g[3]
 111      -#define p4 g[4]
 112      -#define p5 g[5]
 113      -#define p6 g[6]
 114      -#define q1 g[7]
 115      -#define q2 g[8]
 116      -#define q3 g[9]
 117      -#define q4 g[10]
 118      -#define pio2hi g[11]
 119      -#define pio2lo g[12]
 120      -#define t1 g[13]
 121      -#define t2 g[14]
 122      -#define t3 g[15]
 123      -
      109 +#define one             g[0]
      110 +#define p1              g[1]
      111 +#define p2              g[2]
      112 +#define p3              g[3]
      113 +#define p4              g[4]
      114 +#define p5              g[5]
      115 +#define p6              g[6]
      116 +#define q1              g[7]
      117 +#define q2              g[8]
      118 +#define q3              g[9]
      119 +#define q4              g[10]
      120 +#define pio2hi          g[11]
      121 +#define pio2lo          g[12]
      122 +#define t1              g[13]
      123 +#define t2              g[14]
      124 +#define t3              g[15]
 124  125  
 125  126  double
 126      -atan(double x) {
      127 +atan(double x)
      128 +{
 127  129          double y, z, r, p, s;
 128  130          int ix, lx, hx, j;
 129  131  
 130      -        hx = ((int *) &x)[HIWORD];
 131      -        lx = ((int *) &x)[LOWORD];
      132 +        hx = ((int *)&x)[HIWORD];
      133 +        lx = ((int *)&x)[LOWORD];
 132  134          ix = hx & ~0x80000000;
 133  135          j = ix >> 20;
 134  136  
 135  137          /* for |x| < 1/8 */
 136  138          if (j < 0x3fc) {
 137      -                if (j < 0x3f5) {        /* when |x| < 2**(-prec/6-2) */
 138      -                        if (j < 0x3e3) {        /* if |x| < 2**(-prec/2-2) */
 139      -                                return ((int) x == 0 ? x : one);
 140      -                        }
      139 +                if (j < 0x3f5) {                /* when |x| < 2**(-prec/6-2) */
      140 +                        if (j < 0x3e3)          /* if |x| < 2**(-prec/2-2) */
      141 +                                return ((int)x == 0 ? x : one);
      142 +
 141  143                          if (j < 0x3f1) {        /* if |x| < 2**(-prec/4-1) */
 142  144                                  return (x + (x * t1) * (x * x));
 143      -                        } else {        /* if |x| < 2**(-prec/6-2) */
      145 +                        } else {                /* if |x| < 2**(-prec/6-2) */
 144  146                                  z = x * x;
 145  147                                  s = t2 * x;
 146  148                                  return (x + (t3 + z) * (s * z));
 147  149                          }
 148  150                  }
 149      -                z = x * x; s = p1 * x;
      151 +
      152 +                z = x * x;
      153 +                s = p1 * x;
 150  154                  return (x + ((s * z) * (p2 + z * (p3 + z))) *
 151      -                                (((p4 + z) + z * z) * (p5 + z * (p6 + z))));
      155 +                    (((p4 + z) + z * z) * (p5 + z * (p6 + z))));
 152  156          }
 153  157  
 154  158          /* for |x| >= 8.0 */
 155  159          if (j >= 0x402) {
 156  160                  if (j < 0x436) {
 157  161                          r = one / x;
      162 +
 158  163                          if (hx >= 0) {
 159      -                                y =  pio2hi; p =  pio2lo;
      164 +                                y = pio2hi;
      165 +                                p = pio2lo;
 160  166                          } else {
 161      -                                y = -pio2hi; p = -pio2lo;
      167 +                                y = -pio2hi;
      168 +                                p = -pio2lo;
 162  169                          }
      170 +
 163  171                          if (ix < 0x40504000) {  /* x <  65 */
 164  172                                  z = r * r;
 165  173                                  s = p1 * r;
 166      -                                return (y + ((p - r) - ((s * z) *
 167      -                                        (p2 + z * (p3 + z))) *
 168      -                                        (((p4 + z) + z * z) *
 169      -                                        (p5 + z * (p6 + z)))));
      174 +                                return (y + ((p - r) - ((s * z) * (p2 + z *
      175 +                                    (p3 + z))) * (((p4 + z) + z * z) *
      176 +                                    (p5 + z * (p6 + z)))));
 170  177                          } else if (j < 0x412) {
 171  178                                  z = r * r;
 172      -                                return (y + (p - ((q1 * r) * (q4 + z)) *
 173      -                                        (q2 + z * (q3 + z))));
 174      -                        } else
      179 +                                return (y + (p - ((q1 * r) * (q4 + z)) * (q2 +
      180 +                                    z * (q3 + z))));
      181 +                        } else {
 175  182                                  return (y + (p - r));
      183 +                        }
 176  184                  } else {
 177      -                        if (j >= 0x7ff) /* x is inf or NaN */
      185 +                        if (j >= 0x7ff) /* x is inf or NaN */
 178  186                                  if (((ix - 0x7ff00000) | lx) != 0)
 179  187  #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 180  188                                          return (ix >= 0x7ff80000 ? x : x - x);
 181      -                                        /* assumes sparc-like QNaN */
      189 +
      190 +                        /* assumes sparc-like QNaN */
 182  191  #else
 183  192                                          return (x - x);
 184  193  #endif
 185  194                          y = -pio2lo;
 186  195                          return (hx >= 0 ? pio2hi - y : y - pio2hi);
 187  196                  }
 188      -        } else {        /* now x is between 1/8 and 8 */
      197 +        } else {                        /* now x is between 1/8 and 8 */
 189  198                  double *w, w0, w1, s, z;
 190      -                w = (double *) _TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
 191      -                w0 = (hx >= 0)? w[0] : -w[0];
      199 +
      200 +                w = (double *)_TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
      201 +                w0 = (hx >= 0) ? w[0] : -w[0];
 192  202                  s = (x - w0) / (one + x * w0);
 193      -                w1 = (hx >= 0)? w[1] : -w[1];
      203 +                w1 = (hx >= 0) ? w[1] : -w[1];
 194  204                  z = s * s;
 195  205                  return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1);
 196  206          }
 197  207  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX