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
20long 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
29static const long double atanhi[] = {
30	 4.63647609000806116202e-01L,
31	 7.85398163397448309628e-01L,
32	 9.82793723247329067960e-01L,
33	 1.57079632679489661926e+00L,
34};
35
36static const long double atanlo[] = {
37	 1.18469937025062860669e-20L,
38	-1.25413940316708300586e-20L,
39	 2.55232234165405176172e-20L,
40	-2.50827880633416601173e-20L,
41};
42
43static 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
59static 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
65static 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
73static const long double atanhi[] = {
74	 4.63647609000806116214256231461214397e-01L,
75	 7.85398163397448309615660845819875699e-01L,
76	 9.82793723247329067985710611014666038e-01L,
77	 1.57079632679489661923132169163975140e+00L,
78};
79
80static const long double atanlo[] = {
81	 4.89509642257333492668618435220297706e-36L,
82	 2.16795253253094525619926100651083806e-35L,
83	-2.31288434538183565909319952098066272e-35L,
84	 4.33590506506189051239852201302167613e-35L,
85};
86
87static 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
114static 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
121static 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
129long 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
185