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 } 210