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 __unused; 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 } 197