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


   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


  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 }


   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


  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 }