[337] | 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 | const long double atanhi[] = {
|
---|
| 74 | 4.63647609000806116214256231461214397e-01L,
|
---|
| 75 | 7.85398163397448309615660845819875699e-01L,
|
---|
| 76 | 9.82793723247329067985710611014666038e-01L,
|
---|
| 77 | 1.57079632679489661923132169163975140e+00L,
|
---|
| 78 | };
|
---|
| 79 |
|
---|
| 80 | const long double atanlo[] = {
|
---|
| 81 | 4.89509642257333492668618435220297706e-36L,
|
---|
| 82 | 2.16795253253094525619926100651083806e-35L,
|
---|
| 83 | -2.31288434538183565909319952098066272e-35L,
|
---|
| 84 | 4.33590506506189051239852201302167613e-35L,
|
---|
| 85 | };
|
---|
| 86 |
|
---|
| 87 | 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
|
---|