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 static const long double huge = 1.0e300;
  27
  28 long double atanl(long double x)
  29 {
  30         union IEEEl2bits u;
  31         long double w,s1,s2,z;
  32         int id;
  33         int16_t expsign, expt;
  34         int32_t expman;
  35
  36         u.e = x;
  37         expsign = u.xbits.expsign;
  38         expt = expsign & 0x7fff;
  39         if (expt >= ATAN_CONST) { /* if |x| is large, atan(x)~=pi/2 */
  40                 if (expt == BIAS + LDBL_MAX_EXP &&
  41                     ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0)  /* NaN */
  42                         return x+x;
  43                 if (expsign > 0)
  44                         return  atanhi[3]+atanlo[3];
  45                 else
  46                         return -atanhi[3]-atanlo[3];
  47         }
  48         /* Extract the exponent and the first few bits of the mantissa. */
  49         /* XXX There should be a more convenient way to do this. */
  50         expman = (expt << 8) | ((u.bits.manh >> (MANH_SIZE - 9)) & 0xff);
  51         if (expman < ((BIAS - 2) << 8) + 0xc0) {  /* |x| < 0.4375 */
  52                 if (expt < ATAN_LINEAR) {   /* if |x| is small, atanl(x)~=x */
  53                         /* raise inexact */
  54                         if (huge+x > 1.0)
  55                                 return x;
  56                 }
  57                 id = -1;
  58         } else {
  59                 x = fabsl(x);
  60                 if (expman < (BIAS << 8) + 0x30) {  /* |x| < 1.1875 */
  61                         if (expman < ((BIAS - 1) << 8) + 0x60) { /*  7/16 <= |x| < 11/16 */
  62                                 id = 0;
  63                                 x = (2.0*x-1.0)/(2.0+x);
  64                         } else {                                 /* 11/16 <= |x| < 19/16 */
  65                                 id = 1;
  66                                 x = (x-1.0)/(x+1.0);
  67                         }
  68                 } else {
  69                         if (expman < ((BIAS + 1) << 8) + 0x38) { /* |x| < 2.4375 */
  70                                 id = 2;
  71                                 x = (x-1.5)/(1.0+1.5*x);
  72                         } else {                                 /* 2.4375 <= |x| < 2^ATAN_CONST */
  73                                 id = 3;
  74                                 x = -1.0/x;
  75                         }
  76                 }
  77         }
  78         /* end of argument reduction */
  79         z = x*x;
  80         w = z*z;
  81         /* break sum aT[i]z**(i+1) into odd and even poly */
  82         s1 = z*T_even(w);
  83         s2 = w*T_odd(w);
  84         if (id < 0)
  85                 return x - x*(s1+s2);
  86         z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
  87         return expsign < 0 ? -z : z;
  88 }
  89 #endif