1 /* origin: FreeBSD /usr/src/lib/msun/src/s_atanl.c */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
10 * ====================================================
13 * See comments in atan.c.
14 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
20 long double atanl(long double x)
24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
26 static const long double atanhi[] = {
27 4.63647609000806116202e-01L,
28 7.85398163397448309628e-01L,
29 9.82793723247329067960e-01L,
30 1.57079632679489661926e+00L,
33 static const long double atanlo[] = {
34 1.18469937025062860669e-20L,
35 -1.25413940316708300586e-20L,
36 2.55232234165405176172e-20L,
37 -2.50827880633416601173e-20L,
40 static const long double aT[] = {
41 3.33333333333333333017e-01L,
42 -1.99999999999999632011e-01L,
43 1.42857142857046531280e-01L,
44 -1.11111111100562372733e-01L,
45 9.09090902935647302252e-02L,
46 -7.69230552476207730353e-02L,
47 6.66661718042406260546e-02L,
48 -5.88158892835030888692e-02L,
49 5.25499891539726639379e-02L,
50 -4.70119845393155721494e-02L,
51 4.03539201366454414072e-02L,
52 -2.91303858419364158725e-02L,
53 1.24822046299269234080e-02L,
56 static long double T_even(long double x)
58 return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] +
59 x * (aT[8] + x * (aT[10] + x * aT[12])))));
62 static long double T_odd(long double x)
64 return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] +
65 x * (aT[9] + x * aT[11]))));
68 long double atanl(long double x)
71 long double w,s1,s2,z;
73 uint16_t expsign, expt;
77 expsign = u.xbits.expsign;
78 expt = expsign & 0x7fff;
79 if (expt >= 0x3fff + 65) { /* if |x| is large, atan(x)~=pi/2 */
81 ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)!=0) /* NaN */
83 z = atanhi[3] + 0x1p-120f;
84 return expsign>>15 ? -z : z;
86 /* Extract the exponent and the first few bits of the mantissa. */
87 /* XXX There should be a more convenient way to do this. */
88 expman = (expt << 8) | ((u.bits.manh >> (LDBL_MANH_SIZE - 9)) & 0xff);
89 if (expman < ((0x3fff - 2) << 8) + 0xc0) { /* |x| < 0.4375 */
90 if (expt < 0x3fff - 32) { /* if |x| is small, atanl(x)~=x */
91 /* raise underflow if subnormal */
99 if (expman < (0x3fff << 8) + 0x30) { /* |x| < 1.1875 */
100 if (expman < ((0x3fff - 1) << 8) + 0x60) { /* 7/16 <= |x| < 11/16 */
102 x = (2.0*x-1.0)/(2.0+x);
103 } else { /* 11/16 <= |x| < 19/16 */
108 if (expman < ((0x3fff + 1) << 8) + 0x38) { /* |x| < 2.4375 */
110 x = (x-1.5)/(1.0+1.5*x);
111 } else { /* 2.4375 <= |x| < 2^ATAN_CONST */
117 /* end of argument reduction */
120 /* break sum aT[i]z**(i+1) into odd and even poly */
124 return x - x*(s1+s2);
125 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
126 return expsign>>15 ? -z : z;