1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.c */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 /* 13 * See comments in atan.c. 14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. 15 */ 16 17 #include "libm.h" 18 19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 20 long double atanl(long double x) 21 { 22 return atan(x); 23 } 24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 25 26 #if LDBL_MANT_DIG == 64 27 #define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | (u.i.m>>55 & 0xff)) 28 29 static const long double atanhi[] = { 30 4.63647609000806116202e-01L, 31 7.85398163397448309628e-01L, 32 9.82793723247329067960e-01L, 33 1.57079632679489661926e+00L, 34 }; 35 36 static const long double atanlo[] = { 37 1.18469937025062860669e-20L, 38 -1.25413940316708300586e-20L, 39 2.55232234165405176172e-20L, 40 -2.50827880633416601173e-20L, 41 }; 42 43 static const long double aT[] = { 44 3.33333333333333333017e-01L, 45 -1.99999999999999632011e-01L, 46 1.42857142857046531280e-01L, 47 -1.11111111100562372733e-01L, 48 9.09090902935647302252e-02L, 49 -7.69230552476207730353e-02L, 50 6.66661718042406260546e-02L, 51 -5.88158892835030888692e-02L, 52 5.25499891539726639379e-02L, 53 -4.70119845393155721494e-02L, 54 4.03539201366454414072e-02L, 55 -2.91303858419364158725e-02L, 56 1.24822046299269234080e-02L, 57 }; 58 59 static long double T_even(long double x) 60 { 61 return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + 62 x * (aT[8] + x * (aT[10] + x * aT[12]))))); 63 } 64 65 static long double T_odd(long double x) 66 { 67 return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + 68 x * (aT[9] + x * aT[11])))); 69 } 70 #elif LDBL_MANT_DIG == 113 71 #define EXPMAN(u) ((u.i.se & 0x7fff)<<8 | u.i.top>>8) 72 73 static const long double atanhi[] = { 74 4.63647609000806116214256231461214397e-01L, 75 7.85398163397448309615660845819875699e-01L, 76 9.82793723247329067985710611014666038e-01L, 77 1.57079632679489661923132169163975140e+00L, 78 }; 79 80 static const long double atanlo[] = { 81 4.89509642257333492668618435220297706e-36L, 82 2.16795253253094525619926100651083806e-35L, 83 -2.31288434538183565909319952098066272e-35L, 84 4.33590506506189051239852201302167613e-35L, 85 }; 86 87 static const long double aT[] = { 88 3.33333333333333333333333333333333125e-01L, 89 -1.99999999999999999999999999999180430e-01L, 90 1.42857142857142857142857142125269827e-01L, 91 -1.11111111111111111111110834490810169e-01L, 92 9.09090909090909090908522355708623681e-02L, 93 -7.69230769230769230696553844935357021e-02L, 94 6.66666666666666660390096773046256096e-02L, 95 -5.88235294117646671706582985209643694e-02L, 96 5.26315789473666478515847092020327506e-02L, 97 -4.76190476189855517021024424991436144e-02L, 98 4.34782608678695085948531993458097026e-02L, 99 -3.99999999632663469330634215991142368e-02L, 100 3.70370363987423702891250829918659723e-02L, 101 -3.44827496515048090726669907612335954e-02L, 102 3.22579620681420149871973710852268528e-02L, 103 -3.03020767654269261041647570626778067e-02L, 104 2.85641979882534783223403715930946138e-02L, 105 -2.69824879726738568189929461383741323e-02L, 106 2.54194698498808542954187110873675769e-02L, 107 -2.35083879708189059926183138130183215e-02L, 108 2.04832358998165364349957325067131428e-02L, 109 -1.54489555488544397858507248612362957e-02L, 110 8.64492360989278761493037861575248038e-03L, 111 -2.58521121597609872727919154569765469e-03L, 112 }; 113 114 static long double T_even(long double x) 115 { 116 return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * (aT[8] + 117 x * (aT[10] + x * (aT[12] + x * (aT[14] + x * (aT[16] + 118 x * (aT[18] + x * (aT[20] + x * aT[22]))))))))))); 119 } 120 121 static long double T_odd(long double x) 122 { 123 return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * (aT[9] + 124 x * (aT[11] + x * (aT[13] + x * (aT[15] + x * (aT[17] + 125 x * (aT[19] + x * (aT[21] + x * aT[23]))))))))))); 126 } 127 #endif 128 129 long double atanl(long double x) 130 { 131 union ldshape u = {x}; 132 long double w, s1, s2, z; 133 int id; 134 unsigned e = u.i.se & 0x7fff; 135 unsigned sign = u.i.se >> 15; 136 unsigned expman; 137 138 if (e >= 0x3fff + LDBL_MANT_DIG + 1) { /* if |x| is large, atan(x)~=pi/2 */ 139 if (isnan(x)) 140 return x; 141 return sign ? -atanhi[3] : atanhi[3]; 142 } 143 /* Extract the exponent and the first few bits of the mantissa. */ 144 expman = EXPMAN(u); 145 if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */ 146 if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { /* if |x| is small, atanl(x)~=x */ 147 /* raise underflow if subnormal */ 148 if (e == 0) 149 FORCE_EVAL((float)x); 150 return x; 151 } 152 id = -1; 153 } else { 154 x = fabsl(x); 155 if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */ 156 if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */ 157 id = 0; 158 x = (2.0*x-1.0)/(2.0+x); 159 } else { /* 11/16 <= |x| < 19/16 */ 160 id = 1; 161 x = (x-1.0)/(x+1.0); 162 } 163 } else { 164 if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */ 165 id = 2; 166 x = (x-1.5)/(1.0+1.5*x); 167 } else { /* 2.4375 <= |x| */ 168 id = 3; 169 x = -1.0/x; 170 } 171 } 172 } 173 /* end of argument reduction */ 174 z = x*x; 175 w = z*z; 176 /* break sum aT[i]z**(i+1) into odd and even poly */ 177 s1 = z*T_even(w); 178 s2 = w*T_odd(w); 179 if (id < 0) 180 return x - x*(s1+s2); 181 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x); 182 return sign ? -z : z; 183 } 184 #endif 185