usr/src/lib/libm/common/C/atan.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak __atan = atan
  31
  32 /* INDENT OFF */
  33 /*
  34  * atan(x)
  35  * Accurate Table look-up algorithm with polynomial approximation in
  36  * partially product form.
  37  *
  38  * -- K.C. Ng, October 17, 2004
  39  *
  40  * Algorithm
  41  *
  42  * (1). Purge off Inf and NaN and 0
  43  * (2). Reduce x to positive by atan(x) = -atan(-x).
  44  * (3). For x <= 1/8 and let z = x*x, return
  45  *      (2.1) if x < 2^(-prec/2), atan(x) = x  with inexact flag raised
  46  *      (2.2) if x < 2^(-prec/4-1), atan(x) = x+(x/3)(x*x)
  47  *      (2.3) if x < 2^(-prec/6-2), atan(x) = x+(z-5/3)(z*x/5)
  48  *      (2.4) Otherwise
  49  *              atan(x) = poly1(x) = x + A * B,
  50  *      where
  51  *              A = (p1*x*z) * (p2+z(p3+z))
  52  *              B = (p4+z)+z*z) * (p5+z(p6+z))
  53  *      Note: (i) domain of poly1 is [0, 1/8], (ii) remez relative
  54  *      approximation error of poly1 is bounded by
  55  *              |(atan(x)-poly1(x))/x| <= 2^-57.61
  56  * (4). For x >= 8 then
  57  *      (3.1) if x >= 2^prec,     atan(x) = atan(inf) - pio2lo
  58  *      (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x
  59  *      (3.3) if x <= 65,         atan(x) = atan(inf) - poly1(1/x)
  60  *      (3.4) otherwise           atan(x) = atan(inf) - poly2(1/x)
  61  *      where
  62  *              poly2(r) = (q1*r) * (q2+z(q3+z)) * (q4+z),
  63  *      its domain is [0, 0.0154]; and its remez absolute
  64  *      approximation error is bounded by
  65  *              |atan(x)-poly2(x)|<= 2^-59.45
  66  *
  67  * (5). Now x is in (0.125, 8).
  68  *      Recall identity
  69  *              atan(x) = atan(y) + atan((x-y)/(1+x*y)).
  70  *      Let j = (ix - 0x3fc00000) >> 16, 0 <= j < 96, where ix is the high
  71  *      part of x in IEEE double format. Then
  72  *              atan(x) = atan(y[j]) + poly2((x-y[j])/(1+x*y[j]))
  73  *      where y[j] are carefully chosen so that it matches x to around 4.5
  74  *      bits and at the same time atan(y[j]) is very close to an IEEE double
  75  *      floating point number. Calculation indicates that
  76  *              max|(x-y[j])/(1+x*y[j])| < 0.0154
  77  *              j,x
  78  *
  79  * Accuracy: Maximum error observed is bounded by 0.6 ulp after testing
  80  * more than 10 million random arguments
  81  */
  82 /* INDENT ON */
  83
  84 #include "libm.h"
  85 #include "libm_protos.h"
  86
  87 extern const double _TBL_atan[];
  88 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,
  96 /* q1   = */ -1.42796626333911796935538518482644576579332351685e-0001,
  97 /* q2   = */  3.51427110447873227059810477159863497078605962912e+0000,
  98 /* q3   = */  5.92129112708164262457444237952586263418197631836e-0001,
  99 /* q4   = */ -1.99272234785683144409063061175402253866195678711e+0000,
 100 /* pio2hi */  1.570796326794896558e+00,
 101 /* pio2lo */  6.123233995736765886e-17,
 102 /* t1   = */ -0.333333333333333333333333333333333,
 103 /* t2   = */  0.2,
 104 /* t3   = */ -1.666666666666666666666666666666666,
 105 };
 106
 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
 124
 125 double
 126 atan(double x) {
 127         double y, z, r, p, s;
 128         int ix, lx, hx, j;
 129
 130         hx = ((int *) &x)[HIWORD];
 131         lx = ((int *) &x)[LOWORD];
 132         ix = hx & ~0x80000000;
 133         j = ix >> 20;
 134
 135         /* for |x| < 1/8 */
 136         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                         }
 141                         if (j < 0x3f1) {        /* if |x| < 2**(-prec/4-1) */
 142                                 return (x + (x * t1) * (x * x));
 143                         } else {        /* if |x| < 2**(-prec/6-2) */
 144                                 z = x * x;
 145                                 s = t2 * x;
 146                                 return (x + (t3 + z) * (s * z));
 147                         }
 148                 }
 149                 z = x * x; s = p1 * x;
 150                 return (x + ((s * z) * (p2 + z * (p3 + z))) *
 151                                 (((p4 + z) + z * z) * (p5 + z * (p6 + z))));
 152         }
 153
 154         /* for |x| >= 8.0 */
 155         if (j >= 0x402) {
 156                 if (j < 0x436) {
 157                         r = one / x;
 158                         if (hx >= 0) {
 159                                 y =  pio2hi; p =  pio2lo;
 160                         } else {
 161                                 y = -pio2hi; p = -pio2lo;
 162                         }
 163                         if (ix < 0x40504000) {  /* x <  65 */
 164                                 z = r * r;
 165                                 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)))));
 170                         } else if (j < 0x412) {
 171                                 z = r * r;
 172                                 return (y + (p - ((q1 * r) * (q4 + z)) *
 173                                         (q2 + z * (q3 + z))));
 174                         } else
 175                                 return (y + (p - r));
 176                 } else {
 177                         if (j >= 0x7ff) /* x is inf or NaN */
 178                                 if (((ix - 0x7ff00000) | lx) != 0)
 179 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 180                                         return (ix >= 0x7ff80000 ? x : x - x);
 181                                         /* assumes sparc-like QNaN */
 182 #else
 183                                         return (x - x);
 184 #endif
 185                         y = -pio2lo;
 186                         return (hx >= 0 ? pio2hi - y : y - pio2hi);
 187                 }
 188         } else {        /* now x is between 1/8 and 8 */
 189                 double *w, w0, w1, s, z;
 190                 w = (double *) _TBL_atan + (((ix - 0x3fc00000) >> 16) << 1);
 191                 w0 = (hx >= 0)? w[0] : -w[0];
 192                 s = (x - w0) / (one + x * w0);
 193                 w1 = (hx >= 0)? w[1] : -w[1];
 194                 z = s * s;
 195                 return (((q1 * s) * (q4 + z)) * (q2 + z * (q3 + z)) + w1);
 196         }
 197 }