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 atanl = __atanl
  31 
  32 /*
  33  * atanl(x)
  34  * Table look-up algorithm
  35  * By K.C. Ng, March 9, 1989
  36  *
  37  * Algorithm.
  38  *
  39  * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)).
  40  * We use poly1(x) to approximate atan(x) for x in [0,1/8] with
  41  * error (relative)
  42  *      |(atan(x)-poly1(x))/x|<= 2^-115.94   long double
  43  *      |(atan(x)-poly1(x))/x|<= 2^-58.85    double
  44  *      |(atan(x)-poly1(x))/x|<= 2^-25.53    float
  45  * and use poly2(x) to approximate atan(x) for x in [0,1/65] with
  46  * error (absolute)
  47  *      |atan(x)-poly2(x)|<= 2^-122.15               long double
  48  *      |atan(x)-poly2(x)|<= 2^-64.79                double
  49  *      |atan(x)-poly2(x)|<= 2^-35.36                float
  50  * Here poly1 and poly2 are odd polynomial with the following form:
  51  *              x + x^3*(a1+x^2*(a2+...))
  52  *
  53  * (0). Purge off Inf and NaN and 0
  54  * (1). Reduce x to positive by atan(x) = -atan(-x).
  55  * (2). For x <= 1/8, use
  56  *      (2.1) if x < 2^(-prec/2-2), atan(x) =  x  with inexact
  57  *      (2.2) Otherwise
  58  *              atan(x) = poly1(x)
  59  * (3). For x >= 8 then
  60  *      (3.1) if x >= 2^(prec+2),   atan(x) = atan(inf) - pio2lo
  61  *      (3.2) if x >= 2^(prec/3+2), atan(x) = atan(inf) - 1/x
  62  *      (3.3) if x >  65,           atan(x) = atan(inf) - poly2(1/x)
  63  *      (3.4) Otherwise,            atan(x) = atan(inf) - poly1(1/x)
  64  *
  65  * (4). Now x is in (0.125, 8)
  66  *      Find y that match x to 4.5 bit after binary (easy).
  67  *      If iy is the high word of y, then
  68  *              single : j = (iy - 0x3e000000) >> 19
  69  *              double : j = (iy - 0x3fc00000) >> 16
  70  *              quad   : j = (iy - 0x3ffc0000) >> 12
  71  *
  72  *      Let s = (x-y)/(1+x*y). Then
  73  *      atan(x) = atan(y) + poly1(s)
  74  *              = _TBL_atanl_hi[j] + (_TBL_atanl_lo[j] + poly2(s) )
  75  *
  76  *      Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125
  77  *
  78  */
  79 
  80 #include "libm.h"
  81 
  82 extern const long double _TBL_atanl_hi[], _TBL_atanl_lo[];
  83 static const long double
  84         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         long double y, z, r, p, s;
 114         int *px = (int *) &x, *py = (int *) &y;
 115         int ix, iy, sign, j;
 116 
 117         ix = px[i0];
 118         sign = ix & 0x80000000;
 119         ix ^= sign;
 120 
 121         /* for |x| < 1/8 */
 122         if (ix < 0x3ffc0000) {
 123                 if (ix < 0x3feb0000) {       /* when |x| < 2**(-prec/6-2) */
 124                         if (ix < 0x3fc50000) {       /* if |x| < 2**(-prec/2-2) */
 125                                 s = one;
 126                                 *(3 - i0 + (int *) &s) = -1;        /* s = 1-ulp */
 127                                 *(1 + (int *) &s) = -1;
 128                                 *(2 + (int *) &s) = -1;
 129                                 *(i0 + (int *) &s) -= 1;
 130                                 if ((int) (s * x) < 1)
 131                                         return (x);     /* raise inexact */
 132                         }
 133                         z = x * x;
 134                         if (ix < 0x3fe20000) {       /* if |x| < 2**(-prec/4-1) */
 135                                 return (x + (x * z) * p1);
 136                         } else {        /* if |x| < 2**(-prec/6-2) */
 137                                 return (x + (x * z) * (p1 + z * p2));
 138                         }
 139                 }
 140                 z = x * x;
 141                 return (x + (x * z) * (p1 + z * (p2 + z * (p3 + z * (p4 +
 142                         z * (p5 + z * (p6 + z * (p7 + z * (p8 + z * (p9 +
 143                         z * (p10 + z * (p11 + z * (p12 + z * p13)))))))))))));
 144         }
 145 
 146         /* for |x| >= 8.0 */
 147         if (ix >= 0x40020000) {
 148                 px[i0] = ix;
 149                 if (ix < 0x40050400) {       /* x <  65 */
 150                         r = one / x;
 151                         z = r * r;
 152                         /*
 153                          * poly1
 154                          */
 155                         y = r * (one + z * (p1 + z * (p2 + z * (p3 +
 156                                 z * (p4 + z * (p5 + z * (p6 + z * (p7 +
 157                                 z * (p8 + z * (p9 + z * (p10 + z * (p11 +
 158                                 z * (p12 + z * p13)))))))))))));
 159                         y -= pio2lo;
 160                 } else if (ix < 0x40260000) {        /* x <  2**(prec/3+2) */
 161                         r = one / x;
 162                         z = r * r;
 163                         /*
 164                          * poly2
 165                          */
 166                         y = r * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
 167                                 z * (q5 + z * (q6 + z * q7)))))));
 168                         y -= pio2lo;
 169                 } else if (ix < 0x40720000) {        /* x <  2**(prec+2) */
 170                         y = one / x - pio2lo;
 171                 } else if (ix < 0x7fff0000) {        /* x <  inf */
 172                         y = -pio2lo;
 173                 } else {                /* x is inf or NaN */
 174                         if (((ix - 0x7fff0000) | px[1] | px[2] | px[i1]) != 0)
 175                                 return (x - x);
 176                         y = -pio2lo;
 177                 }
 178 
 179                 if (sign == 0)
 180                         return (pio2hi - y);
 181                 else
 182                         return (y - pio2hi);
 183         }
 184 
 185         /* now x is between 1/8 and 8 */
 186         px[i0] = ix;
 187         iy = (ix + 0x00000800) & 0x7ffff000;
 188         py[i0] = iy;
 189         py[1] = py[2] = py[i1] = 0;
 190         j = (iy - 0x3ffc0000) >> 12;
 191 
 192         if (sign == 0)
 193                 s = (x - y) / (one + x * y);
 194         else
 195                 s = (y - x) / (one + x * y);
 196         z = s * s;
 197         if (ix == iy)
 198                 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * q4))));
 199         else
 200                 p = s * (one + z * (q1 + z * (q2 + z * (q3 + z * (q4 +
 201                         z * (q5 + z * (q6 + z * q7)))))));
 202         if (sign == 0) {
 203                 r = p + _TBL_atanl_lo[j];
 204                 return (r + _TBL_atanl_hi[j]);
 205         } else {
 206                 r = p - _TBL_atanl_lo[j];
 207                 return (r - _TBL_atanl_hi[j]);
 208         }
 209 }