usr/src/lib/libm/common/R/atanf.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 __atanf = atanf
  31
  32 /* INDENT OFF */
  33 /*
  34  * float atanf(float x);
  35  * Table look-up algorithm
  36  * By K.C. Ng, March 9, 1989
  37  *
  38  * Algorithm.
  39  *
  40  * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
  41  * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
  42  * error (relative)
  43  *      |(atan(x)-poly1(x))/x|<= 2^-115.94      long double
  44  *      |(atan(x)-poly1(x))/x|<= 2^-58.85       double
  45  *      |(atan(x)-poly1(x))/x|<= 2^-25.53       float
  46  * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
  47  * error (absolute)
  48  *      |atan(x)-poly2(x)|<= 2^-122.15          long double
  49  *      |atan(x)-poly2(x)|<= 2^-64.79           double
  50  *      |atan(x)-poly2(x)|<= 2^-35.36           float
  51  * and use poly3(x) to approximate atan(x) for x in [1/8,7/16] with
  52  * error (relative, on for single precision)
  53  *      |(atan(x)-poly1(x))/x|<= 2^-25.53       float
  54  *
  55  * Here poly1-3 are odd polynomial with the following form:
  56  *              x + x^3*(a1+x^2*(a2+...))
  57  *
  58  * (0). Purge off Inf and NaN and 0
  59  * (1). Reduce x to positive by atan(x) = -atan(-x).
  60  * (2). For x <= 1/8, use
  61  *      (2.1) if x < 2^(-prec/2-2), atan(x) =  x  with inexact
  62  *      (2.2) Otherwise
  63  *              atan(x) = poly1(x)
  64  * (3). For x >= 8 then
  65  *      (3.1) if x >= 2^(prec+2),   atan(x) = atan(inf) - pio2lo
  66  *      (3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
  67  *      (3.3) if x >  65,           atan(x) = atan(inf) - poly2(1/x)
  68  *      (3.4) Otherwise,            atan(x) = atan(inf) - poly1(1/x)
  69  *
  70  * (4). Now x is in (0.125, 8)
  71  *      Find y that match x to 4.5 bit after binary (easy).
  72  *      If iy is the high word of y, then
  73  *              single : j = (iy - 0x3e000000) >> 19
  74  *              (single is modified to (iy-0x3f000000)>>19)
  75  *              double : j = (iy - 0x3fc00000) >> 16
  76  *              quad   : j = (iy - 0x3ffc0000) >> 12
  77  *
  78  *      Let s = (x-y)/(1+x*y). Then
  79  *      atan(x) = atan(y) + poly1(s)
  80  *              = _TBL_r_atan_hi[j] + (_TBL_r_atan_lo[j] + poly2(s) )
  81  *
  82  *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
  83  *
  84  */
  85
  86 #include "libm.h"
  87
  88 extern const float _TBL_r_atan_hi[], _TBL_r_atan_lo[];
  89 static const float
  90         big     =   1.0e37F,
  91         one     =   1.0F,
  92         p1      =  -3.333185951111688247225368498733544672172e-0001F,
  93         p2      =   1.969352894213455405211341983203180636021e-0001F,
  94         q1      =  -3.332921964095646819563419704110132937456e-0001F,
  95         a1      =  -3.333323465223893614063523351509338934592e-0001F,
  96         a2      =   1.999425625935277805494082274808174062403e-0001F,
  97         a3      =  -1.417547090509737780085769846290301788559e-0001F,
  98         a4      =   1.016250813871991983097273733227432685084e-0001F,
  99         a5      =  -5.137023693688358515753093811791755221805e-0002F,
 100         pio2hi  =   1.570796371e+0000F,
 101         pio2lo  =  -4.371139000e-0008F;
 102 /* INDENT ON */
 103
 104 float
 105 atanf(float xx) {
 106         float x, y, z, r, p, s;
 107         volatile double dummy;
 108         int ix, iy, sign, j;
 109
 110         x = xx;
 111         ix = *(int *) &x;
 112         sign = ix & 0x80000000;
 113         ix ^= sign;
 114
 115         /* for |x| < 1/8 */
 116         if (ix < 0x3e000000) {
 117                 if (ix < 0x38800000) {  /* if |x| < 2**(-prec/2-2) */
 118                         dummy = big + x;        /* get inexact flag if x != 0 */
 119 #ifdef lint
 120                         dummy = dummy;
 121 #endif
 122                         return (x);
 123                 }
 124                 z = x * x;
 125                 if (ix < 0x3c000000) {  /* if |x| < 2**(-prec/4-1) */
 126                         x = x + (x * z) * p1;
 127                         return (x);
 128                 } else {
 129                         x = x + (x * z) * (p1 + z * p2);
 130                         return (x);
 131                 }
 132         }
 133
 134         /* for |x| >= 8.0 */
 135         if (ix >= 0x41000000) {
 136                 *(int *) &x = ix;
 137                 if (ix < 0x42820000) {  /* x <  65 */
 138                         r = one / x;
 139                         z = r * r;
 140                         y = r * (one + z * (p1 + z * p2));      /* poly1 */
 141                         y -= pio2lo;
 142                 } else if (ix < 0x44800000) {   /* x <  2**(prec/3+2) */
 143                         r = one / x;
 144                         z = r * r;
 145                         y = r * (one + z * q1); /* poly2 */
 146                         y -= pio2lo;
 147                 } else if (ix < 0x4c800000) {   /* x <  2**(prec+2) */
 148                         y = one / x - pio2lo;
 149                 } else if (ix < 0x7f800000) {   /* x <  inf */
 150                         y = -pio2lo;
 151                 } else {                /* x is inf or NaN */
 152                         if (ix > 0x7f800000) {
 153                                 return (x * x); /* - -> * for Cheetah */
 154                         }
 155                         y = -pio2lo;
 156                 }
 157
 158                 if (sign == 0)
 159                         x = pio2hi - y;
 160                 else
 161                         x = y - pio2hi;
 162                 return (x);
 163         }
 164
 165
 166         /* now x is between 1/8 and 8 */
 167         if (ix < 0x3f000000) {  /* between 1/8 and 1/2 */
 168                 z = x * x;
 169                 x = x + (x * z) * (a1 + z * (a2 + z * (a3 + z * (a4 +
 170                         z * a5))));
 171                 return (x);
 172         }
 173         *(int *) &x = ix;
 174         iy = (ix + 0x00040000) & 0x7ff80000;
 175         *(int *) &y = iy;
 176         j = (iy - 0x3f000000) >> 19;
 177
 178         if (ix == iy)
 179                 p = x - y;      /* p=0.0 */
 180         else {
 181                 if (sign == 0)
 182                         s = (x - y) / (one + x * y);
 183                 else
 184                         s = (y - x) / (one + x * y);
 185                 z = s * s;
 186                 p = s * (one + z * q1);
 187         }
 188         if (sign == 0) {
 189                 r = p + _TBL_r_atan_lo[j];
 190                 x = r + _TBL_r_atan_hi[j];
 191         } else {
 192                 r = p - _TBL_r_atan_lo[j];
 193                 x = r - _TBL_r_atan_hi[j];
 194         }
 195         return (x);
 196 }