1/* s_atanl.c
2 *
3 * Inverse circular tangent for 128-bit long double precision
4 * (arctangent)
5 *
6 *
7 *
8 * SYNOPSIS:
9 *
10 * long double x, y, atanl();
11 *
12 * y = atanl( x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
19 *
20 * The function uses a rational approximation of the form
21 * t + t^3 P(t^2)/Q(t^2), optimized for |t| < 0.09375.
22 *
23 * The argument is reduced using the identity
24 * arctan x - arctan u = arctan ((x-u)/(1 + ux))
25 * and an 83-entry lookup table for arctan u, with u = 0, 1/8, ..., 10.25.
26 * Use of the table improves the execution speed of the routine.
27 *
28 *
29 *
30 * ACCURACY:
31 *
32 * Relative error:
33 * arithmetic domain # trials peak rms
34 * IEEE -19, 19 4e5 1.7e-34 5.4e-35
35 *
36 *
37 * WARNING:
38 *
39 * This program uses integer operations on bit fields of floating-point
40 * numbers. It does not work with data structures other than the
41 * structure assumed.
42 *
43 */
44
45/* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
46
47 This library is free software; you can redistribute it and/or
48 modify it under the terms of the GNU Lesser General Public
49 License as published by the Free Software Foundation; either
50 version 2.1 of the License, or (at your option) any later version.
51
52 This library is distributed in the hope that it will be useful,
53 but WITHOUT ANY WARRANTY; without even the implied warranty of
54 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
55 Lesser General Public License for more details.
56
57 You should have received a copy of the GNU Lesser General Public
58 License along with this library; if not, see
59 <http://www.gnu.org/licenses/>. */
60
61
62#include <float.h>
63#include <math.h>
64#include <math_private.h>
65#include <libm-alias-ldouble.h>
66
67/* arctan(k/8), k = 0, ..., 82 */
68static const _Float128 atantbl[84] = {
69 L(0.0000000000000000000000000000000000000000E0),
70 L(1.2435499454676143503135484916387102557317E-1), /* arctan(0.125) */
71 L(2.4497866312686415417208248121127581091414E-1),
72 L(3.5877067027057222039592006392646049977698E-1),
73 L(4.6364760900080611621425623146121440202854E-1),
74 L(5.5859931534356243597150821640166127034645E-1),
75 L(6.4350110879328438680280922871732263804151E-1),
76 L(7.1882999962162450541701415152590465395142E-1),
77 L(7.8539816339744830961566084581987572104929E-1),
78 L(8.4415398611317100251784414827164750652594E-1),
79 L(8.9605538457134395617480071802993782702458E-1),
80 L(9.4200004037946366473793717053459358607166E-1),
81 L(9.8279372324732906798571061101466601449688E-1),
82 L(1.0191413442663497346383429170230636487744E0),
83 L(1.0516502125483736674598673120862998296302E0),
84 L(1.0808390005411683108871567292171998202703E0),
85 L(1.1071487177940905030170654601785370400700E0),
86 L(1.1309537439791604464709335155363278047493E0),
87 L(1.1525719972156675180401498626127513797495E0),
88 L(1.1722738811284763866005949441337046149712E0),
89 L(1.1902899496825317329277337748293183376012E0),
90 L(1.2068173702852525303955115800565576303133E0),
91 L(1.2220253232109896370417417439225704908830E0),
92 L(1.2360594894780819419094519711090786987027E0),
93 L(1.2490457723982544258299170772810901230778E0),
94 L(1.2610933822524404193139408812473357720101E0),
95 L(1.2722973952087173412961937498224804940684E0),
96 L(1.2827408797442707473628852511364955306249E0),
97 L(1.2924966677897852679030914214070816845853E0),
98 L(1.3016288340091961438047858503666855921414E0),
99 L(1.3101939350475556342564376891719053122733E0),
100 L(1.3182420510168370498593302023271362531155E0),
101 L(1.3258176636680324650592392104284756311844E0),
102 L(1.3329603993374458675538498697331558093700E0),
103 L(1.3397056595989995393283037525895557411039E0),
104 L(1.3460851583802539310489409282517796256512E0),
105 L(1.3521273809209546571891479413898128509842E0),
106 L(1.3578579772154994751124898859640585287459E0),
107 L(1.3633001003596939542892985278250991189943E0),
108 L(1.3684746984165928776366381936948529556191E0),
109 L(1.3734007669450158608612719264449611486510E0),
110 L(1.3780955681325110444536609641291551522494E0),
111 L(1.3825748214901258580599674177685685125566E0),
112 L(1.3868528702577214543289381097042486034883E0),
113 L(1.3909428270024183486427686943836432060856E0),
114 L(1.3948567013423687823948122092044222644895E0),
115 L(1.3986055122719575950126700816114282335732E0),
116 L(1.4021993871854670105330304794336492676944E0),
117 L(1.4056476493802697809521934019958079881002E0),
118 L(1.4089588955564736949699075250792569287156E0),
119 L(1.4121410646084952153676136718584891599630E0),
120 L(1.4152014988178669079462550975833894394929E0),
121 L(1.4181469983996314594038603039700989523716E0),
122 L(1.4209838702219992566633046424614466661176E0),
123 L(1.4237179714064941189018190466107297503086E0),
124 L(1.4263547484202526397918060597281265695725E0),
125 L(1.4288992721907326964184700745371983590908E0),
126 L(1.4313562697035588982240194668401779312122E0),
127 L(1.4337301524847089866404719096698873648610E0),
128 L(1.4360250423171655234964275337155008780675E0),
129 L(1.4382447944982225979614042479354815855386E0),
130 L(1.4403930189057632173997301031392126865694E0),
131 L(1.4424730991091018200252920599377292525125E0),
132 L(1.4444882097316563655148453598508037025938E0),
133 L(1.4464413322481351841999668424758804165254E0),
134 L(1.4483352693775551917970437843145232637695E0),
135 L(1.4501726582147939000905940595923466567576E0),
136 L(1.4519559822271314199339700039142990228105E0),
137 L(1.4536875822280323362423034480994649820285E0),
138 L(1.4553696664279718992423082296859928222270E0),
139 L(1.4570043196511885530074841089245667532358E0),
140 L(1.4585935117976422128825857356750737658039E0),
141 L(1.4601391056210009726721818194296893361233E0),
142 L(1.4616428638860188872060496086383008594310E0),
143 L(1.4631064559620759326975975316301202111560E0),
144 L(1.4645314639038178118428450961503371619177E0),
145 L(1.4659193880646627234129855241049975398470E0),
146 L(1.4672716522843522691530527207287398276197E0),
147 L(1.4685896086876430842559640450619880951144E0),
148 L(1.4698745421276027686510391411132998919794E0),
149 L(1.4711276743037345918528755717617308518553E0),
150 L(1.4723501675822635384916444186631899205983E0),
151 L(1.4735431285433308455179928682541563973416E0), /* arctan(10.25) */
152 L(1.5707963267948966192313216916397514420986E0) /* pi/2 */
153};
154
155
156/* arctan t = t + t^3 p(t^2) / q(t^2)
157 |t| <= 0.09375
158 peak relative error 5.3e-37 */
159
160static const _Float128
161 p0 = L(-4.283708356338736809269381409828726405572E1),
162 p1 = L(-8.636132499244548540964557273544599863825E1),
163 p2 = L(-5.713554848244551350855604111031839613216E1),
164 p3 = L(-1.371405711877433266573835355036413750118E1),
165 p4 = L(-8.638214309119210906997318946650189640184E-1),
166 q0 = L(1.285112506901621042780814422948906537959E2),
167 q1 = L(3.361907253914337187957855834229672347089E2),
168 q2 = L(3.180448303864130128268191635189365331680E2),
169 q3 = L(1.307244136980865800160844625025280344686E2),
170 q4 = L(2.173623741810414221251136181221172551416E1);
171 /* q5 = 1.000000000000000000000000000000000000000E0 */
172
173static const _Float128 huge = L(1.0e4930);
174
175_Float128
176__atanl (_Float128 x)
177{
178 int k, sign;
179 _Float128 t, u, p, q;
180 ieee854_long_double_shape_type s;
181
182 s.value = x;
183 k = s.parts32.w0;
184 if (k & 0x80000000)
185 sign = 1;
186 else
187 sign = 0;
188
189 /* Check for IEEE special cases. */
190 k &= 0x7fffffff;
191 if (k >= 0x7fff0000)
192 {
193 /* NaN. */
194 if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3)
195 return (x + x);
196
197 /* Infinity. */
198 if (sign)
199 return -atantbl[83];
200 else
201 return atantbl[83];
202 }
203
204 if (k <= 0x3fc50000) /* |x| < 2**-58 */
205 {
206 math_check_force_underflow (x);
207 /* Raise inexact. */
208 if (huge + x > 0.0)
209 return x;
210 }
211
212 if (k >= 0x40720000) /* |x| > 2**115 */
213 {
214 /* Saturate result to {-,+}pi/2 */
215 if (sign)
216 return -atantbl[83];
217 else
218 return atantbl[83];
219 }
220
221 if (sign)
222 x = -x;
223
224 if (k >= 0x40024800) /* 10.25 */
225 {
226 k = 83;
227 t = -1.0/x;
228 }
229 else
230 {
231 /* Index of nearest table element.
232 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
233 (cf. fdlibm). */
234 k = 8.0 * x + 0.25;
235 u = L(0.125) * k;
236 /* Small arctan argument. */
237 t = (x - u) / (1.0 + x * u);
238 }
239
240 /* Arctan of small argument t. */
241 u = t * t;
242 p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
243 q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
244 u = t * u * p / q + t;
245
246 /* arctan x = arctan u + arctan t */
247 u = atantbl[k] + u;
248 if (sign)
249 return (-u);
250 else
251 return u;
252}
253
254libm_alias_ldouble (__atan, atan)
255