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> */
51 /* arctan(k/8), k = 0, ..., 82 */
52 static const long double atantbl[84] = {
53 0.0000000000000000000000000000000000000000E0L,
54 1.2435499454676143503135484916387102557317E-1L, /* arctan(0.125) */
55 2.4497866312686415417208248121127581091414E-1L,
56 3.5877067027057222039592006392646049977698E-1L,
57 4.6364760900080611621425623146121440202854E-1L,
58 5.5859931534356243597150821640166127034645E-1L,
59 6.4350110879328438680280922871732263804151E-1L,
60 7.1882999962162450541701415152590465395142E-1L,
61 7.8539816339744830961566084581987572104929E-1L,
62 8.4415398611317100251784414827164750652594E-1L,
63 8.9605538457134395617480071802993782702458E-1L,
64 9.4200004037946366473793717053459358607166E-1L,
65 9.8279372324732906798571061101466601449688E-1L,
66 1.0191413442663497346383429170230636487744E0L,
67 1.0516502125483736674598673120862998296302E0L,
68 1.0808390005411683108871567292171998202703E0L,
69 1.1071487177940905030170654601785370400700E0L,
70 1.1309537439791604464709335155363278047493E0L,
71 1.1525719972156675180401498626127513797495E0L,
72 1.1722738811284763866005949441337046149712E0L,
73 1.1902899496825317329277337748293183376012E0L,
74 1.2068173702852525303955115800565576303133E0L,
75 1.2220253232109896370417417439225704908830E0L,
76 1.2360594894780819419094519711090786987027E0L,
77 1.2490457723982544258299170772810901230778E0L,
78 1.2610933822524404193139408812473357720101E0L,
79 1.2722973952087173412961937498224804940684E0L,
80 1.2827408797442707473628852511364955306249E0L,
81 1.2924966677897852679030914214070816845853E0L,
82 1.3016288340091961438047858503666855921414E0L,
83 1.3101939350475556342564376891719053122733E0L,
84 1.3182420510168370498593302023271362531155E0L,
85 1.3258176636680324650592392104284756311844E0L,
86 1.3329603993374458675538498697331558093700E0L,
87 1.3397056595989995393283037525895557411039E0L,
88 1.3460851583802539310489409282517796256512E0L,
89 1.3521273809209546571891479413898128509842E0L,
90 1.3578579772154994751124898859640585287459E0L,
91 1.3633001003596939542892985278250991189943E0L,
92 1.3684746984165928776366381936948529556191E0L,
93 1.3734007669450158608612719264449611486510E0L,
94 1.3780955681325110444536609641291551522494E0L,
95 1.3825748214901258580599674177685685125566E0L,
96 1.3868528702577214543289381097042486034883E0L,
97 1.3909428270024183486427686943836432060856E0L,
98 1.3948567013423687823948122092044222644895E0L,
99 1.3986055122719575950126700816114282335732E0L,
100 1.4021993871854670105330304794336492676944E0L,
101 1.4056476493802697809521934019958079881002E0L,
102 1.4089588955564736949699075250792569287156E0L,
103 1.4121410646084952153676136718584891599630E0L,
104 1.4152014988178669079462550975833894394929E0L,
105 1.4181469983996314594038603039700989523716E0L,
106 1.4209838702219992566633046424614466661176E0L,
107 1.4237179714064941189018190466107297503086E0L,
108 1.4263547484202526397918060597281265695725E0L,
109 1.4288992721907326964184700745371983590908E0L,
110 1.4313562697035588982240194668401779312122E0L,
111 1.4337301524847089866404719096698873648610E0L,
112 1.4360250423171655234964275337155008780675E0L,
113 1.4382447944982225979614042479354815855386E0L,
114 1.4403930189057632173997301031392126865694E0L,
115 1.4424730991091018200252920599377292525125E0L,
116 1.4444882097316563655148453598508037025938E0L,
117 1.4464413322481351841999668424758804165254E0L,
118 1.4483352693775551917970437843145232637695E0L,
119 1.4501726582147939000905940595923466567576E0L,
120 1.4519559822271314199339700039142990228105E0L,
121 1.4536875822280323362423034480994649820285E0L,
122 1.4553696664279718992423082296859928222270E0L,
123 1.4570043196511885530074841089245667532358E0L,
124 1.4585935117976422128825857356750737658039E0L,
125 1.4601391056210009726721818194296893361233E0L,
126 1.4616428638860188872060496086383008594310E0L,
127 1.4631064559620759326975975316301202111560E0L,
128 1.4645314639038178118428450961503371619177E0L,
129 1.4659193880646627234129855241049975398470E0L,
130 1.4672716522843522691530527207287398276197E0L,
131 1.4685896086876430842559640450619880951144E0L,
132 1.4698745421276027686510391411132998919794E0L,
133 1.4711276743037345918528755717617308518553E0L,
134 1.4723501675822635384916444186631899205983E0L,
135 1.4735431285433308455179928682541563973416E0L, /* arctan(10.25) */
136 1.5707963267948966192313216916397514420986E0L /* pi/2 */
140 /* arctan t = t + t^3 p(t^2) / q(t^2)
142 peak relative error 5.3e-37 */
144 static const long double
145 p0 = -4.283708356338736809269381409828726405572E1L,
146 p1 = -8.636132499244548540964557273544599863825E1L,
147 p2 = -5.713554848244551350855604111031839613216E1L,
148 p3 = -1.371405711877433266573835355036413750118E1L,
149 p4 = -8.638214309119210906997318946650189640184E-1L,
150 q0 = 1.285112506901621042780814422948906537959E2L,
151 q1 = 3.361907253914337187957855834229672347089E2L,
152 q2 = 3.180448303864130128268191635189365331680E2L,
153 q3 = 1.307244136980865800160844625025280344686E2L,
154 q4 = 2.173623741810414221251136181221172551416E1L;
155 /* q5 = 1.000000000000000000000000000000000000000E0 */
159 atanl (long double x)
162 long double t, u, p, q;
164 sign = (x < 0) ? -1 : 1;
166 /* Check for zero or NaN. */
167 if (x != x || x == 0.0)
170 /* Check for infinity. */
172 return sign * atantbl[83];
176 if (x >= 10.25) /* 10.25 */
183 /* Index of nearest table element.
184 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
188 /* Small arctan argument. */
189 t = (x - u) / (1.0 + x * u);
192 /* Arctan of small argument t. */
194 p = ((((p4 * u) + p3) * u + p2) * u + p1) * u + p0;
195 q = ((((u + q4) * u + q3) * u + q2) * u + q1) * u + q0;
196 u = t * u * p / q + t;
198 /* arctan x = arctan u + arctan t */
199 return sign * (atantbl[k] + u);