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
66/* arctan(k/8), k = 0, ..., 82 */
67static const _Float128 atantbl[84] = {
68 L(0.0000000000000000000000000000000000000000E0),
69 L(1.2435499454676143503135484916387102557317E-1), /* arctan(0.125) */
70 L(2.4497866312686415417208248121127581091414E-1),
71 L(3.5877067027057222039592006392646049977698E-1),
72 L(4.6364760900080611621425623146121440202854E-1),
73 L(5.5859931534356243597150821640166127034645E-1),
74 L(6.4350110879328438680280922871732263804151E-1),
75 L(7.1882999962162450541701415152590465395142E-1),
76 L(7.8539816339744830961566084581987572104929E-1),
77 L(8.4415398611317100251784414827164750652594E-1),
78 L(8.9605538457134395617480071802993782702458E-1),
79 L(9.4200004037946366473793717053459358607166E-1),
80 L(9.8279372324732906798571061101466601449688E-1),
81 L(1.0191413442663497346383429170230636487744E0),
82 L(1.0516502125483736674598673120862998296302E0),
83 L(1.0808390005411683108871567292171998202703E0),
84 L(1.1071487177940905030170654601785370400700E0),
85 L(1.1309537439791604464709335155363278047493E0),
86 L(1.1525719972156675180401498626127513797495E0),
87 L(1.1722738811284763866005949441337046149712E0),
88 L(1.1902899496825317329277337748293183376012E0),
89 L(1.2068173702852525303955115800565576303133E0),
90 L(1.2220253232109896370417417439225704908830E0),
91 L(1.2360594894780819419094519711090786987027E0),
92 L(1.2490457723982544258299170772810901230778E0),
93 L(1.2610933822524404193139408812473357720101E0),
94 L(1.2722973952087173412961937498224804940684E0),
95 L(1.2827408797442707473628852511364955306249E0),
96 L(1.2924966677897852679030914214070816845853E0),
97 L(1.3016288340091961438047858503666855921414E0),
98 L(1.3101939350475556342564376891719053122733E0),
99 L(1.3182420510168370498593302023271362531155E0),
100 L(1.3258176636680324650592392104284756311844E0),
101 L(1.3329603993374458675538498697331558093700E0),
102 L(1.3397056595989995393283037525895557411039E0),
103 L(1.3460851583802539310489409282517796256512E0),
104 L(1.3521273809209546571891479413898128509842E0),
105 L(1.3578579772154994751124898859640585287459E0),
106 L(1.3633001003596939542892985278250991189943E0),
107 L(1.3684746984165928776366381936948529556191E0),
108 L(1.3734007669450158608612719264449611486510E0),
109 L(1.3780955681325110444536609641291551522494E0),
110 L(1.3825748214901258580599674177685685125566E0),
111 L(1.3868528702577214543289381097042486034883E0),
112 L(1.3909428270024183486427686943836432060856E0),
113 L(1.3948567013423687823948122092044222644895E0),
114 L(1.3986055122719575950126700816114282335732E0),
115 L(1.4021993871854670105330304794336492676944E0),
116 L(1.4056476493802697809521934019958079881002E0),
117 L(1.4089588955564736949699075250792569287156E0),
118 L(1.4121410646084952153676136718584891599630E0),
119 L(1.4152014988178669079462550975833894394929E0),
120 L(1.4181469983996314594038603039700989523716E0),
121 L(1.4209838702219992566633046424614466661176E0),
122 L(1.4237179714064941189018190466107297503086E0),
123 L(1.4263547484202526397918060597281265695725E0),
124 L(1.4288992721907326964184700745371983590908E0),
125 L(1.4313562697035588982240194668401779312122E0),
126 L(1.4337301524847089866404719096698873648610E0),
127 L(1.4360250423171655234964275337155008780675E0),
128 L(1.4382447944982225979614042479354815855386E0),
129 L(1.4403930189057632173997301031392126865694E0),
130 L(1.4424730991091018200252920599377292525125E0),
131 L(1.4444882097316563655148453598508037025938E0),
132 L(1.4464413322481351841999668424758804165254E0),
133 L(1.4483352693775551917970437843145232637695E0),
134 L(1.4501726582147939000905940595923466567576E0),
135 L(1.4519559822271314199339700039142990228105E0),
136 L(1.4536875822280323362423034480994649820285E0),
137 L(1.4553696664279718992423082296859928222270E0),
138 L(1.4570043196511885530074841089245667532358E0),
139 L(1.4585935117976422128825857356750737658039E0),
140 L(1.4601391056210009726721818194296893361233E0),
141 L(1.4616428638860188872060496086383008594310E0),
142 L(1.4631064559620759326975975316301202111560E0),
143 L(1.4645314639038178118428450961503371619177E0),
144 L(1.4659193880646627234129855241049975398470E0),
145 L(1.4672716522843522691530527207287398276197E0),
146 L(1.4685896086876430842559640450619880951144E0),
147 L(1.4698745421276027686510391411132998919794E0),
148 L(1.4711276743037345918528755717617308518553E0),
149 L(1.4723501675822635384916444186631899205983E0),
150 L(1.4735431285433308455179928682541563973416E0), /* arctan(10.25) */
151 L(1.5707963267948966192313216916397514420986E0) /* pi/2 */
152};
153
154
155/* arctan t = t + t^3 p(t^2) / q(t^2)
156 |t| <= 0.09375
157 peak relative error 5.3e-37 */
158
159static const _Float128
160 p0 = L(-4.283708356338736809269381409828726405572E1),
161 p1 = L(-8.636132499244548540964557273544599863825E1),
162 p2 = L(-5.713554848244551350855604111031839613216E1),
163 p3 = L(-1.371405711877433266573835355036413750118E1),
164 p4 = L(-8.638214309119210906997318946650189640184E-1),
165 q0 = L(1.285112506901621042780814422948906537959E2),
166 q1 = L(3.361907253914337187957855834229672347089E2),
167 q2 = L(3.180448303864130128268191635189365331680E2),
168 q3 = L(1.307244136980865800160844625025280344686E2),
169 q4 = L(2.173623741810414221251136181221172551416E1);
170 /* q5 = 1.000000000000000000000000000000000000000E0 */
171
172static const _Float128 huge = L(1.0e4930);
173
174_Float128
175__atanl (_Float128 x)
176{
177 int k, sign;
178 _Float128 t, u, p, q;
179 ieee854_long_double_shape_type s;
180
181 s.value = x;
182 k = s.parts32.w0;
183 if (k & 0x80000000)
184 sign = 1;
185 else
186 sign = 0;
187
188 /* Check for IEEE special cases. */
189 k &= 0x7fffffff;
190 if (k >= 0x7fff0000)
191 {
192 /* NaN. */
193 if ((k & 0xffff) | s.parts32.w1 | s.parts32.w2 | s.parts32.w3)
194 return (x + x);
195
196 /* Infinity. */
197 if (sign)
198 return -atantbl[83];
199 else
200 return atantbl[83];
201 }
202
203 if (k <= 0x3fc50000) /* |x| < 2**-58 */
204 {
205 math_check_force_underflow (x);
206 /* Raise inexact. */
207 if (huge + x > 0.0)
208 return x;
209 }
210
211 if (k >= 0x40720000) /* |x| > 2**115 */
212 {
213 /* Saturate result to {-,+}pi/2 */
214 if (sign)
215 return -atantbl[83];
216 else
217 return atantbl[83];
218 }
219
220 if (sign)
221 x = -x;
222
223 if (k >= 0x40024800) /* 10.25 */
224 {
225 k = 83;
226 t = -1.0/x;
227 }
228 else
229 {
230 /* Index of nearest table element.
231 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
232 (cf. fdlibm). */
233 k = 8.0 * x + 0.25;
234 u = L(0.125) * k;
235 /* Small arctan argument. */
236 t = (x - u) / (1.0 + x * u);
237 }
238
239 /* Arctan of small argument t. */
240 u = t * t;
241 p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
242 q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
243 u = t * u * p / q + t;
244
245 /* arctan x = arctan u + arctan t */
246 u = atantbl[k] + u;
247 if (sign)
248 return (-u);
249 else
250 return u;
251}
252
253weak_alias (__atanl, atanl)
254