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 program is free software; you can redistribute it and/or modify
48 it under the terms of the GNU General Public License as published by
49 the Free Software Foundation; either version 2, or (at your option)
52 This program 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
55 GNU General Public License for more details.
57 You should have received a copy of the GNU General Public License
58 along with this program; see the file COPYING.
59 If not, write to the Free Software Foundation,
60 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
66 /* arctan(k/8), k = 0, ..., 82 */
67 static const long double atantbl[84] = {
68 0.0000000000000000000000000000000000000000E0L,
69 1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125) */
70 2.4497866312686415417208248121127581091414E-1L,
71 3.5877067027057222039592006392646049977698E-1L,
72 4.6364760900080611621425623146121440202854E-1L,
73 5.5859931534356243597150821640166127034645E-1L,
74 6.4350110879328438680280922871732263804151E-1L,
75 7.1882999962162450541701415152590465395142E-1L,
76 7.8539816339744830961566084581987572104929E-1L,
77 8.4415398611317100251784414827164750652594E-1L,
78 8.9605538457134395617480071802993782702458E-1L,
79 9.4200004037946366473793717053459358607166E-1L,
80 9.8279372324732906798571061101466601449688E-1L,
81 1.0191413442663497346383429170230636487744E0L,
82 1.0516502125483736674598673120862998296302E0L,
83 1.0808390005411683108871567292171998202703E0L,
84 1.1071487177940905030170654601785370400700E0L,
85 1.1309537439791604464709335155363278047493E0L,
86 1.1525719972156675180401498626127513797495E0L,
87 1.1722738811284763866005949441337046149712E0L,
88 1.1902899496825317329277337748293183376012E0L,
89 1.2068173702852525303955115800565576303133E0L,
90 1.2220253232109896370417417439225704908830E0L,
91 1.2360594894780819419094519711090786987027E0L,
92 1.2490457723982544258299170772810901230778E0L,
93 1.2610933822524404193139408812473357720101E0L,
94 1.2722973952087173412961937498224804940684E0L,
95 1.2827408797442707473628852511364955306249E0L,
96 1.2924966677897852679030914214070816845853E0L,
97 1.3016288340091961438047858503666855921414E0L,
98 1.3101939350475556342564376891719053122733E0L,
99 1.3182420510168370498593302023271362531155E0L,
100 1.3258176636680324650592392104284756311844E0L,
101 1.3329603993374458675538498697331558093700E0L,
102 1.3397056595989995393283037525895557411039E0L,
103 1.3460851583802539310489409282517796256512E0L,
104 1.3521273809209546571891479413898128509842E0L,
105 1.3578579772154994751124898859640585287459E0L,
106 1.3633001003596939542892985278250991189943E0L,
107 1.3684746984165928776366381936948529556191E0L,
108 1.3734007669450158608612719264449611486510E0L,
109 1.3780955681325110444536609641291551522494E0L,
110 1.3825748214901258580599674177685685125566E0L,
111 1.3868528702577214543289381097042486034883E0L,
112 1.3909428270024183486427686943836432060856E0L,
113 1.3948567013423687823948122092044222644895E0L,
114 1.3986055122719575950126700816114282335732E0L,
115 1.4021993871854670105330304794336492676944E0L,
116 1.4056476493802697809521934019958079881002E0L,
117 1.4089588955564736949699075250792569287156E0L,
118 1.4121410646084952153676136718584891599630E0L,
119 1.4152014988178669079462550975833894394929E0L,
120 1.4181469983996314594038603039700989523716E0L,
121 1.4209838702219992566633046424614466661176E0L,
122 1.4237179714064941189018190466107297503086E0L,
123 1.4263547484202526397918060597281265695725E0L,
124 1.4288992721907326964184700745371983590908E0L,
125 1.4313562697035588982240194668401779312122E0L,
126 1.4337301524847089866404719096698873648610E0L,
127 1.4360250423171655234964275337155008780675E0L,
128 1.4382447944982225979614042479354815855386E0L,
129 1.4403930189057632173997301031392126865694E0L,
130 1.4424730991091018200252920599377292525125E0L,
131 1.4444882097316563655148453598508037025938E0L,
132 1.4464413322481351841999668424758804165254E0L,
133 1.4483352693775551917970437843145232637695E0L,
134 1.4501726582147939000905940595923466567576E0L,
135 1.4519559822271314199339700039142990228105E0L,
136 1.4536875822280323362423034480994649820285E0L,
137 1.4553696664279718992423082296859928222270E0L,
138 1.4570043196511885530074841089245667532358E0L,
139 1.4585935117976422128825857356750737658039E0L,
140 1.4601391056210009726721818194296893361233E0L,
141 1.4616428638860188872060496086383008594310E0L,
142 1.4631064559620759326975975316301202111560E0L,
143 1.4645314639038178118428450961503371619177E0L,
144 1.4659193880646627234129855241049975398470E0L,
145 1.4672716522843522691530527207287398276197E0L,
146 1.4685896086876430842559640450619880951144E0L,
147 1.4698745421276027686510391411132998919794E0L,
148 1.4711276743037345918528755717617308518553E0L,
149 1.4723501675822635384916444186631899205983E0L,
150 1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
151 1.5707963267948966192313216916397514420986E0L /* pi/2 */
155 /* arctan t = t + t^3 p(t^2) / q(t^2)
157 peak relative error 5.3e-37 */
159 static const long double
160 p0 = -4.283708356338736809269381409828726405572E1L,
161 p1 = -8.636132499244548540964557273544599863825E1L,
162 p2 = -5.713554848244551350855604111031839613216E1L,
163 p3 = -1.371405711877433266573835355036413750118E1L,
164 p4 = -8.638214309119210906997318946650189640184E-1L,
165 q0 = 1.285112506901621042780814422948906537959E2L,
166 q1 = 3.361907253914337187957855834229672347089E2L,
167 q2 = 3.180448303864130128268191635189365331680E2L,
168 q3 = 1.307244136980865800160844625025280344686E2L,
169 q4 = 2.173623741810414221251136181221172551416E1L;
170 /* q5 = 1.000000000000000000000000000000000000000E0 */
174 atanl (long double x)
177 long double t, u, p, q;
181 /* Check for zero or NaN. */
182 if (x != x || x == 0.0)
204 /* Index of nearest table element.
205 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
209 /* Small arctan argument. */
210 t = (x - u) / (1.0 + x * u);
213 /* Arctan of small argument t. */
215 p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
216 q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
217 u = t * u * p / q + t;
219 /* arctan x = arctan u + arctan t */