src/math/atanl.c

   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