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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __atanl = atanl
  32 
  33 /*
  34  * atanl(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  * Here poly1 and poly2 are odd polynomial with the following form:
  52  *              x + x^3*(a1+x^2*(a2+...))
  53  *
  54  * (0). Purge off Inf and NaN and 0
  55  * (1). Reduce x to positive by atan(x) = -atan(-x).
  56  * (2). For x <= 1/8, use
  57  *      (2.1) if x < 2^(-prec/2-2), atan(x) =  x  with inexact
  58  *      (2.2) Otherwise
  59  *              atan(x) = poly1(x)
  60  * (3). For x >= 8 then
  61  *      (3.1) if x >= 2^(prec+2),   atan(x) = atan(inf) - pio2lo
  62  *      (3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
  63  *      (3.3) if x >  65,           atan(x) = atan(inf) - poly2(1/x)
  64  *      (3.4) Otherwise,            atan(x) = atan(inf) - poly1(1/x)
  65  *
  66  * (4). Now x is in (0.125, 8)
  67  *      Find y that match x to 4.5 bit after binary (easy).
  68  *      If iy is the high word of y, then
  69  *              single : j = (iy - 0x3e000000) >> 19
  70  *              double : j = (iy - 0x3fc00000) >> 16
  71  *              quad   : j = (iy - 0x3ffc0000) >> 12
  72  *
  73  *      Let s = (x-y)/(1+x*y). Then
  74  *      atan(x) = atan(y) + poly1(s)
  75  *              = _TBL_atanl_hi[j] + (_TBL_atanl_lo[j] + poly2(s) )
  76  *
  77  *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
  78  *
  79  */
  80 
  81 #include "libm.h"
  82 
  83 extern const long double _TBL_atanl_hi[], _TBL_atanl_lo[];
  84 static const long double one = 1.0L,
  85         p1 = -3.333333333333333333333333333331344526118e-0001L,
  86         p2 = 1.999999999999999999999999989931277668570e-0001L,
  87         p3 = -1.428571428571428571428553606221309530901e-0001L,
  88         p4 = 1.111111111111111111095219842737139747418e-0001L,
  89         p5 = -9.090909090909090825503603835248061123323e-0002L,
  90         p6 = 7.692307692307664052130743214708925258904e-0002L,
  91         p7 = -6.666666666660213835187713228363717388266e-0002L,
  92         p8 = 5.882352940152439399097283359608661949504e-0002L,
  93         p9 = -5.263157780447533993046614040509529668487e-0002L,
  94         p10 = 4.761895816878184933175855990886788439447e-0002L,
  95         p11 = -4.347345005832274022681019724553538135922e-0002L,
  96         p12 = 3.983031914579635037502589204647752042736e-0002L,
  97         p13 = -3.348206704469830575196657749413894897554e-0002L,
  98         q1 = -3.333333333333333333333333333195273650186e-0001L,
  99         q2 = 1.999999999999999999999988146114392615808e-0001L,
 100         q3 = -1.428571428571428571057630319435467111253e-0001L,
 101         q4 = 1.111111111111105373263048208994541544098e-0001L,
 102         q5 = -9.090909090421834209167373258681021816441e-0002L,
 103         q6 = 7.692305377813692706850171767150701644539e-0002L,
 104         q7 = -6.660896644393861499914731734305717901330e-0002L,
 105         pio2hi = 1.570796326794896619231321691639751398740e+0000L,
 106         pio2lo = 4.335905065061890512398522013021675984381e-0035L;
 107 
 108 #define i0      0
 109 #define i1      3
 110 
 111 long double
 112 atanl(long double x)
 113 {
 114         long double y, z, r, p, s;
 115         int *px = (int *)&x, *py = (int *)&y;
 116         int ix, iy, sign, j;
 117 
 118         ix = px[i0];
 119         sign = ix & 0x80000000;
 120         ix ^= sign;
 121 
 122         /* for |x| < 1/8 */
 123         if (ix < 0x3ffc0000) {
 124                 if (ix < 0x3feb0000) {               /* when |x| < 2**(-prec/6-2) */
 125                         if (ix < 0x3fc50000) {       /* if |x| < 2**(-prec/2-2) */
 126                                 s = one;
 127                                 *(3 - i0 + (int *)&s) = -1; /* s = 1-ulp */
 128                                 *(1 + (int *)&s) = -1;
 129                                 *(2 + (int *)&s) = -1;
 130                                 *(i0 + (int *)&s) -= 1;
 131 
 132                                 if ((int)(s * x) < 1)
 133                                         return (x);     /* raise inexact */
 134                         }
 135 
 136                         z = x * x;
 137 
 138                         if (ix < 0x3fe20000) /* if |x| < 2**(-prec/4-1) */
 139                                 return (x + (x * z) * p1);
 140                         else                    /* if |x| < 2**(-prec/6-2) */
 141                                 return (x + (x * z) * (p1 + z * p2));
 142                 }
 143 
 144                 z = x * x;
 145                 return (x + (x * z) * (p1 + z * (p2 + z * (p3 + z * (p4 + z *
 146                     (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 + z *
 147                     (p10 + z * (p11 + z * (p12 + z * p13)))))))))))));
 148         }
 149 
 150         /* for |x| >= 8.0 */
 151         if (ix >= 0x40020000) {
 152                 px[i0] = ix;
 153 
 154                 if (ix < 0x40050400) {       /* x <  65 */
 155                         r = one / x;
 156                         z = r * r;
 157 
 158                         /*
 159                          * poly1
 160                          */
 161                         y = r * (one + z * (p1 + z * (p2 + z * (p3 + z * (p4 +
 162                             z * (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 +
 163                             z * (p10 + z * (p11 + z * (p12 + z *
 164                             p13)))))))))))));
 165                         y -= pio2lo;
 166                 } else if (ix < 0x40260000) {        /* x <  2**(prec/3+2) */
 167                         r = one / x;
 168                         z = r * r;
 169 
 170                         /*
 171                          * poly2
 172                          */
 173                         y = r * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
 174                             z * (q5 + z * (q6 + z * q7)))))));
 175                         y -= pio2lo;
 176                 } else if (ix < 0x40720000) {        /* x <  2**(prec+2) */
 177                         y = one / x - pio2lo;
 178                 } else if (ix < 0x7fff0000) {        /* x <  inf */
 179                         y = -pio2lo;
 180                 } else {                        /* x is inf or NaN */
 181                         if (((ix - 0x7fff0000) | px[1] | px[2] | px[i1]) != 0)
 182                                 return (x - x);
 183 
 184                         y = -pio2lo;
 185                 }
 186 
 187                 if (sign == 0)
 188                         return (pio2hi - y);
 189                 else
 190                         return (y - pio2hi);
 191         }
 192 
 193         /* now x is between 1/8 and 8 */
 194         px[i0] = ix;
 195         iy = (ix + 0x00000800) & 0x7ffff000;
 196         py[i0] = iy;
 197         py[1] = py[2] = py[i1] = 0;
 198         j = (iy - 0x3ffc0000) >> 12;
 199 
 200         if (sign == 0)
 201                 s = (x - y) / (one + x * y);
 202         else
 203                 s = (y - x) / (one + x * y);
 204 
 205         z = s * s;
 206 
 207         if (ix == iy)
 208                 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * q4))));
 209         else
 210                 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 + z * (q5 +
 211                     z * (q6 + z * q7)))))));
 212 
 213         if (sign == 0) {
 214                 r = p + _TBL_atanl_lo[j];
 215                 return (r + _TBL_atanl_hi[j]);
 216         } else {
 217                 r = p - _TBL_atanl_lo[j];
 218                 return (r - _TBL_atanl_hi[j]);
 219         }
 220 }