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