nsz Git - musl/blob - 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 #include "__invtrigl.h"
  26
  27 #define ATAN_CONST      (BIAS + 65)     /* 2**65 */
  28 #define ATAN_LINEAR     (BIAS - 32)     /* 2**-32 */
  29 static const long double huge = 1.0e300;
  30
  31 static const long double atanhi[] = {
  32          4.63647609000806116202e-01L,
  33          7.85398163397448309628e-01L,
  34          9.82793723247329067960e-01L,
  35          1.57079632679489661926e+00L,
  36 };
  37
  38 static const long double atanlo[] = {
  39          1.18469937025062860669e-20L,
  40         -1.25413940316708300586e-20L,
  41          2.55232234165405176172e-20L,
  42         -2.50827880633416601173e-20L,
  43 };
  44
  45 static const long double aT[] = {
  46          3.33333333333333333017e-01L,
  47         -1.99999999999999632011e-01L,
  48          1.42857142857046531280e-01L,
  49         -1.11111111100562372733e-01L,
  50          9.09090902935647302252e-02L,
  51         -7.69230552476207730353e-02L,
  52          6.66661718042406260546e-02L,
  53         -5.88158892835030888692e-02L,
  54          5.25499891539726639379e-02L,
  55         -4.70119845393155721494e-02L,
  56          4.03539201366454414072e-02L,
  57         -2.91303858419364158725e-02L,
  58          1.24822046299269234080e-02L,
  59 };
  60
  61 static long double T_even(long double x)
  62 {
  63         return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] +
  64                 x * (aT[8] + x * (aT[10] + x * aT[12])))));
  65 }
  66
  67 static long double T_odd(long double x)
  68 {
  69         return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] +
  70                 x * (aT[9] + x * aT[11]))));
  71 }
  72
  73 long double atanl(long double x)
  74 {
  75         union IEEEl2bits u;
  76         long double w,s1,s2,z;
  77         int id;
  78         int16_t expsign, expt;
  79         int32_t expman;
  80
  81         u.e = x;
  82         expsign = u.xbits.expsign;
  83         expt = expsign & 0x7fff;
  84         if (expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */
  85                 if (expt == BIAS + LDBL_MAX_EXP &&
  86                     ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0)  /* NaN */
  87                         return x+x;
  88                 if (expsign > 0)
  89                         return  atanhi[3]+atanlo[3];
  90                 else
  91                         return -atanhi[3]-atanlo[3];
  92         }
  93         /* Extract the exponent and the first few bits of the mantissa. */
  94         /* XXX There should be a more convenient way to do this. */
  95         expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff);
  96         if (expman < ((BIAS - 2) << 8) + 0xc0) {  /* |x| < 0.4375 */
  97                 if (expt < ATAN_LINEAR) {   /* if |x| is small, atanl(x)~=x */
  98                         /* raise inexact */
  99                         if (huge+x > 1.0)
 100                                 return x;
 101                 }
 102                 id = -1;
 103         } else {
 104                 x = fabsl(x);
 105                 if (expman < (BIAS << 8) + 0x30) {  /* |x| < 1.1875 */
 106                         if (expman < ((BIAS - 1) << 8) + 0x60) { /*  7/16 <= |x| < 11/16 */
 107                                 id = 0;
 108                                 x = (2.0*x-1.0)/(2.0+x);
 109                         } else {                                 /* 11/16 <= |x| < 19/16 */
 110                                 id = 1;
 111                                 x = (x-1.0)/(x+1.0);
 112                         }
 113                 } else {
 114                         if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
 115                                 id = 2;
 116                                 x = (x-1.5)/(1.0+1.5*x);
 117                         } else {                                 /* 2.4375 <= |x| < 2^ATAN_CONST */
 118                                 id = 3;
 119                                 x = -1.0/x;
 120                         }
 121                 }
 122         }
 123         /* end of argument reduction */
 124         z = x*x;
 125         w = z*z;
 126         /* break sum aT[i]z**(i+1) into odd and even poly */
 127         s1 = z*T_even(w);
 128         s2 = w*T_odd(w);
 129         if (id < 0)
 130                 return x - x*(s1+s2);
 131         z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
 132         return expsign < 0 ? -z : z;
 133 }
 134 #endif