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