3 * Inverse circular tangent for 128-bit long double precision
10 * long double x, y, atanl();
18 * Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
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.
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.
33 * arithmetic domain # trials peak rms
34 * IEEE -19, 19 4e5 1.7e-34 5.4e-35
39 * This program uses integer operations on bit fields of floating-point
40 * numbers. It does not work with data structures other than the
45 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
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.
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.
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/>. */
62 #include <math_private.h>
64 /* arctan(k/8), k = 0, ..., 82 */
65 static const long double atantbl
[84] = {
66 0.0000000000000000000000000000000000000000E0L
,
67 1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125) */
68 2.4497866312686415417208248121127581091414E-1L,
69 3.5877067027057222039592006392646049977698E-1L,
70 4.6364760900080611621425623146121440202854E-1L,
71 5.5859931534356243597150821640166127034645E-1L,
72 6.4350110879328438680280922871732263804151E-1L,
73 7.1882999962162450541701415152590465395142E-1L,
74 7.8539816339744830961566084581987572104929E-1L,
75 8.4415398611317100251784414827164750652594E-1L,
76 8.9605538457134395617480071802993782702458E-1L,
77 9.4200004037946366473793717053459358607166E-1L,
78 9.8279372324732906798571061101466601449688E-1L,
79 1.0191413442663497346383429170230636487744E0L
,
80 1.0516502125483736674598673120862998296302E0L
,
81 1.0808390005411683108871567292171998202703E0L
,
82 1.1071487177940905030170654601785370400700E0L
,
83 1.1309537439791604464709335155363278047493E0L
,
84 1.1525719972156675180401498626127513797495E0L
,
85 1.1722738811284763866005949441337046149712E0L
,
86 1.1902899496825317329277337748293183376012E0L
,
87 1.2068173702852525303955115800565576303133E0L
,
88 1.2220253232109896370417417439225704908830E0L
,
89 1.2360594894780819419094519711090786987027E0L
,
90 1.2490457723982544258299170772810901230778E0L
,
91 1.2610933822524404193139408812473357720101E0L
,
92 1.2722973952087173412961937498224804940684E0L
,
93 1.2827408797442707473628852511364955306249E0L
,
94 1.2924966677897852679030914214070816845853E0L
,
95 1.3016288340091961438047858503666855921414E0L
,
96 1.3101939350475556342564376891719053122733E0L
,
97 1.3182420510168370498593302023271362531155E0L
,
98 1.3258176636680324650592392104284756311844E0L
,
99 1.3329603993374458675538498697331558093700E0L
,
100 1.3397056595989995393283037525895557411039E0L
,
101 1.3460851583802539310489409282517796256512E0L
,
102 1.3521273809209546571891479413898128509842E0L
,
103 1.3578579772154994751124898859640585287459E0L
,
104 1.3633001003596939542892985278250991189943E0L
,
105 1.3684746984165928776366381936948529556191E0L
,
106 1.3734007669450158608612719264449611486510E0L
,
107 1.3780955681325110444536609641291551522494E0L
,
108 1.3825748214901258580599674177685685125566E0L
,
109 1.3868528702577214543289381097042486034883E0L
,
110 1.3909428270024183486427686943836432060856E0L
,
111 1.3948567013423687823948122092044222644895E0L
,
112 1.3986055122719575950126700816114282335732E0L
,
113 1.4021993871854670105330304794336492676944E0L
,
114 1.4056476493802697809521934019958079881002E0L
,
115 1.4089588955564736949699075250792569287156E0L
,
116 1.4121410646084952153676136718584891599630E0L
,
117 1.4152014988178669079462550975833894394929E0L
,
118 1.4181469983996314594038603039700989523716E0L
,
119 1.4209838702219992566633046424614466661176E0L
,
120 1.4237179714064941189018190466107297503086E0L
,
121 1.4263547484202526397918060597281265695725E0L
,
122 1.4288992721907326964184700745371983590908E0L
,
123 1.4313562697035588982240194668401779312122E0L
,
124 1.4337301524847089866404719096698873648610E0L
,
125 1.4360250423171655234964275337155008780675E0L
,
126 1.4382447944982225979614042479354815855386E0L
,
127 1.4403930189057632173997301031392126865694E0L
,
128 1.4424730991091018200252920599377292525125E0L
,
129 1.4444882097316563655148453598508037025938E0L
,
130 1.4464413322481351841999668424758804165254E0L
,
131 1.4483352693775551917970437843145232637695E0L
,
132 1.4501726582147939000905940595923466567576E0L
,
133 1.4519559822271314199339700039142990228105E0L
,
134 1.4536875822280323362423034480994649820285E0L
,
135 1.4553696664279718992423082296859928222270E0L
,
136 1.4570043196511885530074841089245667532358E0L
,
137 1.4585935117976422128825857356750737658039E0L
,
138 1.4601391056210009726721818194296893361233E0L
,
139 1.4616428638860188872060496086383008594310E0L
,
140 1.4631064559620759326975975316301202111560E0L
,
141 1.4645314639038178118428450961503371619177E0L
,
142 1.4659193880646627234129855241049975398470E0L
,
143 1.4672716522843522691530527207287398276197E0L
,
144 1.4685896086876430842559640450619880951144E0L
,
145 1.4698745421276027686510391411132998919794E0L
,
146 1.4711276743037345918528755717617308518553E0L
,
147 1.4723501675822635384916444186631899205983E0L
,
148 1.4735431285433308455179928682541563973416E0L
, /* arctan(10.25) */
149 1.5707963267948966192313216916397514420986E0L
/* pi/2 */
153 /* arctan t = t + t^3 p(t^2) / q(t^2)
155 peak relative error 5.3e-37 */
157 static const long double
158 p0
= -4.283708356338736809269381409828726405572E1L
,
159 p1
= -8.636132499244548540964557273544599863825E1L
,
160 p2
= -5.713554848244551350855604111031839613216E1L
,
161 p3
= -1.371405711877433266573835355036413750118E1L
,
162 p4
= -8.638214309119210906997318946650189640184E-1L,
163 q0
= 1.285112506901621042780814422948906537959E2L
,
164 q1
= 3.361907253914337187957855834229672347089E2L
,
165 q2
= 3.180448303864130128268191635189365331680E2L
,
166 q3
= 1.307244136980865800160844625025280344686E2L
,
167 q4
= 2.173623741810414221251136181221172551416E1L
;
168 /* q5 = 1.000000000000000000000000000000000000000E0 */
170 static const long double huge
= 1.0e4930L
;
173 __atanl (long double x
)
176 long double t
, u
, p
, q
;
177 ieee854_long_double_shape_type s
;
186 /* Check for IEEE special cases. */
191 if ((k
& 0xffff) | s
.parts32
.w1
| s
.parts32
.w2
| s
.parts32
.w3
)
201 if (k
<= 0x3fc50000) /* |x| < 2**-58 */
208 if (k
>= 0x40720000) /* |x| > 2**115 */
210 /* Saturate result to {-,+}pi/2 */
220 if (k
>= 0x40024800) /* 10.25 */
227 /* Index of nearest table element.
228 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
232 /* Small arctan argument. */
233 t
= (x
- u
) / (1.0 + x
* u
);
236 /* Arctan of small argument 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
;
242 /* arctan x = arctan u + arctan t */
250 weak_alias (__atanl
, atanl
)