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
__atanl(_Float128 x)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)
256