1 |
740 |
jeremybenn |
/* 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 */
|
65 |
|
|
static 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 |
|
|
|
157 |
|
|
static 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 |
|
|
|
170 |
|
|
|
171 |
|
|
__float128
|
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 (sign)
|
201 |
|
|
x = -x;
|
202 |
|
|
|
203 |
|
|
if (k >= 0x40024800) /* 10.25 */
|
204 |
|
|
{
|
205 |
|
|
k = 83;
|
206 |
|
|
t = -1.0/x;
|
207 |
|
|
}
|
208 |
|
|
else
|
209 |
|
|
{
|
210 |
|
|
/* Index of nearest table element.
|
211 |
|
|
Roundoff to integer is asymmetrical to avoid cancellation when t < 0
|
212 |
|
|
(cf. fdlibm). */
|
213 |
|
|
k = 8.0Q * x + 0.25Q;
|
214 |
|
|
u = 0.125Q * k;
|
215 |
|
|
/* Small arctan argument. */
|
216 |
|
|
t = (x - u) / (1.0 + x * u);
|
217 |
|
|
}
|
218 |
|
|
|
219 |
|
|
/* Arctan of small argument t. */
|
220 |
|
|
u = t * t;
|
221 |
|
|
p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
|
222 |
|
|
q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
|
223 |
|
|
u = t * u * p / q + t;
|
224 |
|
|
|
225 |
|
|
/* arctan x = arctan u + arctan t */
|
226 |
|
|
u = atantbl[k] + u;
|
227 |
|
|
if (sign)
|
228 |
|
|
return (-u);
|
229 |
|
|
else
|
230 |
|
|
return u;
|
231 |
|
|
}
|