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