1/* s_atanl.c
2 *
3 * Inverse circular tangent for 128-bit __float128 precision
4 * (arctangent)
5 *
6 *
7 *
8 * SYNOPSIS:
9 *
10 * __float128 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, write to the Free Software
59 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
60
61
62#include "quadmath-imp.h"
63
64/* arctan(k/8), k = 0, ..., 82 */
65static const __float128 atantbl[84] = {
66 0.0000000000000000000000000000000000000000E0Q,
67 1.2435499454676143503135484916387102557317E-1Q, /* arctan(0.125) */
68 2.4497866312686415417208248121127581091414E-1Q,
69 3.5877067027057222039592006392646049977698E-1Q,
70 4.6364760900080611621425623146121440202854E-1Q,
71 5.5859931534356243597150821640166127034645E-1Q,
72 6.4350110879328438680280922871732263804151E-1Q,
73 7.1882999962162450541701415152590465395142E-1Q,
74 7.8539816339744830961566084581987572104929E-1Q,
75 8.4415398611317100251784414827164750652594E-1Q,
76 8.9605538457134395617480071802993782702458E-1Q,
77 9.4200004037946366473793717053459358607166E-1Q,
78 9.8279372324732906798571061101466601449688E-1Q,
79 1.0191413442663497346383429170230636487744E0Q,
80 1.0516502125483736674598673120862998296302E0Q,
81 1.0808390005411683108871567292171998202703E0Q,
82 1.1071487177940905030170654601785370400700E0Q,
83 1.1309537439791604464709335155363278047493E0Q,
84 1.1525719972156675180401498626127513797495E0Q,
85 1.1722738811284763866005949441337046149712E0Q,
86 1.1902899496825317329277337748293183376012E0Q,
87 1.2068173702852525303955115800565576303133E0Q,
88 1.2220253232109896370417417439225704908830E0Q,
89 1.2360594894780819419094519711090786987027E0Q,
90 1.2490457723982544258299170772810901230778E0Q,
91 1.2610933822524404193139408812473357720101E0Q,
92 1.2722973952087173412961937498224804940684E0Q,
93 1.2827408797442707473628852511364955306249E0Q,
94 1.2924966677897852679030914214070816845853E0Q,
95 1.3016288340091961438047858503666855921414E0Q,
96 1.3101939350475556342564376891719053122733E0Q,
97 1.3182420510168370498593302023271362531155E0Q,
98 1.3258176636680324650592392104284756311844E0Q,
99 1.3329603993374458675538498697331558093700E0Q,
100 1.3397056595989995393283037525895557411039E0Q,
101 1.3460851583802539310489409282517796256512E0Q,
102 1.3521273809209546571891479413898128509842E0Q,
103 1.3578579772154994751124898859640585287459E0Q,
104 1.3633001003596939542892985278250991189943E0Q,
105 1.3684746984165928776366381936948529556191E0Q,
106 1.3734007669450158608612719264449611486510E0Q,
107 1.3780955681325110444536609641291551522494E0Q,
108 1.3825748214901258580599674177685685125566E0Q,
109 1.3868528702577214543289381097042486034883E0Q,
110 1.3909428270024183486427686943836432060856E0Q,
111 1.3948567013423687823948122092044222644895E0Q,
112 1.3986055122719575950126700816114282335732E0Q,
113 1.4021993871854670105330304794336492676944E0Q,
114 1.4056476493802697809521934019958079881002E0Q,
115 1.4089588955564736949699075250792569287156E0Q,
116 1.4121410646084952153676136718584891599630E0Q,
117 1.4152014988178669079462550975833894394929E0Q,
118 1.4181469983996314594038603039700989523716E0Q,
119 1.4209838702219992566633046424614466661176E0Q,
120 1.4237179714064941189018190466107297503086E0Q,
121 1.4263547484202526397918060597281265695725E0Q,
122 1.4288992721907326964184700745371983590908E0Q,
123 1.4313562697035588982240194668401779312122E0Q,
124 1.4337301524847089866404719096698873648610E0Q,
125 1.4360250423171655234964275337155008780675E0Q,
126 1.4382447944982225979614042479354815855386E0Q,
127 1.4403930189057632173997301031392126865694E0Q,
128 1.4424730991091018200252920599377292525125E0Q,
129 1.4444882097316563655148453598508037025938E0Q,
130 1.4464413322481351841999668424758804165254E0Q,
131 1.4483352693775551917970437843145232637695E0Q,
132 1.4501726582147939000905940595923466567576E0Q,
133 1.4519559822271314199339700039142990228105E0Q,
134 1.4536875822280323362423034480994649820285E0Q,
135 1.4553696664279718992423082296859928222270E0Q,
136 1.4570043196511885530074841089245667532358E0Q,
137 1.4585935117976422128825857356750737658039E0Q,
138 1.4601391056210009726721818194296893361233E0Q,
139 1.4616428638860188872060496086383008594310E0Q,
140 1.4631064559620759326975975316301202111560E0Q,
141 1.4645314639038178118428450961503371619177E0Q,
142 1.4659193880646627234129855241049975398470E0Q,
143 1.4672716522843522691530527207287398276197E0Q,
144 1.4685896086876430842559640450619880951144E0Q,
145 1.4698745421276027686510391411132998919794E0Q,
146 1.4711276743037345918528755717617308518553E0Q,
147 1.4723501675822635384916444186631899205983E0Q,
148 1.4735431285433308455179928682541563973416E0Q, /* arctan(10.25) */
149 1.5707963267948966192313216916397514420986E0Q /* pi/2 */
150};
151
152
153/* arctan t = t + t^3 p(t^2) / q(t^2)
154 |t| <= 0.09375
155 peak relative error 5.3e-37 */
156
157static const __float128
158 p0 = -4.283708356338736809269381409828726405572E1Q,
159 p1 = -8.636132499244548540964557273544599863825E1Q,
160 p2 = -5.713554848244551350855604111031839613216E1Q,
161 p3 = -1.371405711877433266573835355036413750118E1Q,
162 p4 = -8.638214309119210906997318946650189640184E-1Q,
163 q0 = 1.285112506901621042780814422948906537959E2Q,
164 q1 = 3.361907253914337187957855834229672347089E2Q,
165 q2 = 3.180448303864130128268191635189365331680E2Q,
166 q3 = 1.307244136980865800160844625025280344686E2Q,
167 q4 = 2.173623741810414221251136181221172551416E1Q;
168 /* q5 = 1.000000000000000000000000000000000000000E0 */
169
170static const __float128 huge = 1.0e4930Q;
171
172__float128
173atanq (__float128 x)
174{
175 int k, sign;
176 __float128 t, u, p, q;
177 ieee854_float128 s;
178
179 s.value = x;
180 k = s.words32.w0;
181 if (k & 0x80000000)
182 sign = 1;
183 else
184 sign = 0;
185
186 /* Check for IEEE special cases. */
187 k &= 0x7fffffff;
188 if (k >= 0x7fff0000)
189 {
190 /* NaN. */
191 if ((k & 0xffff) | s.words32.w1 | s.words32.w2 | s.words32.w3)
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 <= 0x3fc50000) /* |x| < 2**-58 */
202 {
203 math_check_force_underflow (x);
204 /* Raise inexact. */
205 if (huge + x > 0.0)
206 return x;
207 }
208
209 if (k >= 0x40720000) /* |x| > 2**115 */
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 >= 0x40024800) /* 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.0Q * x + 0.25Q;
232 u = 0.125Q * 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