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