1 /* s_atanl.c 2 * 3 * Inverse circular tangent for 128-bit long double precision 4 * (arctangent) 5 * 6 * 7 * 8 * SYNOPSIS: 9 * 10 * long double x, y, atanq(); 11 * 12 * y = atanq( 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, see 59 <http://www.gnu.org/licenses/>. */ 60 61 #include "quadmath-imp.h" 62 63 /* arctan(k/8), k = 0, ..., 82 */ 64 static const __float128 atantbl[84] = { 65 0.0000000000000000000000000000000000000000E0Q, 66 1.2435499454676143503135484916387102557317E-1Q, /* arctan(0.125) */ 67 2.4497866312686415417208248121127581091414E-1Q, 68 3.5877067027057222039592006392646049977698E-1Q, 69 4.6364760900080611621425623146121440202854E-1Q, 70 5.5859931534356243597150821640166127034645E-1Q, 71 6.4350110879328438680280922871732263804151E-1Q, 72 7.1882999962162450541701415152590465395142E-1Q, 73 7.8539816339744830961566084581987572104929E-1Q, 74 8.4415398611317100251784414827164750652594E-1Q, 75 8.9605538457134395617480071802993782702458E-1Q, 76 9.4200004037946366473793717053459358607166E-1Q, 77 9.8279372324732906798571061101466601449688E-1Q, 78 1.0191413442663497346383429170230636487744E0Q, 79 1.0516502125483736674598673120862998296302E0Q, 80 1.0808390005411683108871567292171998202703E0Q, 81 1.1071487177940905030170654601785370400700E0Q, 82 1.1309537439791604464709335155363278047493E0Q, 83 1.1525719972156675180401498626127513797495E0Q, 84 1.1722738811284763866005949441337046149712E0Q, 85 1.1902899496825317329277337748293183376012E0Q, 86 1.2068173702852525303955115800565576303133E0Q, 87 1.2220253232109896370417417439225704908830E0Q, 88 1.2360594894780819419094519711090786987027E0Q, 89 1.2490457723982544258299170772810901230778E0Q, 90 1.2610933822524404193139408812473357720101E0Q, 91 1.2722973952087173412961937498224804940684E0Q, 92 1.2827408797442707473628852511364955306249E0Q, 93 1.2924966677897852679030914214070816845853E0Q, 94 1.3016288340091961438047858503666855921414E0Q, 95 1.3101939350475556342564376891719053122733E0Q, 96 1.3182420510168370498593302023271362531155E0Q, 97 1.3258176636680324650592392104284756311844E0Q, 98 1.3329603993374458675538498697331558093700E0Q, 99 1.3397056595989995393283037525895557411039E0Q, 100 1.3460851583802539310489409282517796256512E0Q, 101 1.3521273809209546571891479413898128509842E0Q, 102 1.3578579772154994751124898859640585287459E0Q, 103 1.3633001003596939542892985278250991189943E0Q, 104 1.3684746984165928776366381936948529556191E0Q, 105 1.3734007669450158608612719264449611486510E0Q, 106 1.3780955681325110444536609641291551522494E0Q, 107 1.3825748214901258580599674177685685125566E0Q, 108 1.3868528702577214543289381097042486034883E0Q, 109 1.3909428270024183486427686943836432060856E0Q, 110 1.3948567013423687823948122092044222644895E0Q, 111 1.3986055122719575950126700816114282335732E0Q, 112 1.4021993871854670105330304794336492676944E0Q, 113 1.4056476493802697809521934019958079881002E0Q, 114 1.4089588955564736949699075250792569287156E0Q, 115 1.4121410646084952153676136718584891599630E0Q, 116 1.4152014988178669079462550975833894394929E0Q, 117 1.4181469983996314594038603039700989523716E0Q, 118 1.4209838702219992566633046424614466661176E0Q, 119 1.4237179714064941189018190466107297503086E0Q, 120 1.4263547484202526397918060597281265695725E0Q, 121 1.4288992721907326964184700745371983590908E0Q, 122 1.4313562697035588982240194668401779312122E0Q, 123 1.4337301524847089866404719096698873648610E0Q, 124 1.4360250423171655234964275337155008780675E0Q, 125 1.4382447944982225979614042479354815855386E0Q, 126 1.4403930189057632173997301031392126865694E0Q, 127 1.4424730991091018200252920599377292525125E0Q, 128 1.4444882097316563655148453598508037025938E0Q, 129 1.4464413322481351841999668424758804165254E0Q, 130 1.4483352693775551917970437843145232637695E0Q, 131 1.4501726582147939000905940595923466567576E0Q, 132 1.4519559822271314199339700039142990228105E0Q, 133 1.4536875822280323362423034480994649820285E0Q, 134 1.4553696664279718992423082296859928222270E0Q, 135 1.4570043196511885530074841089245667532358E0Q, 136 1.4585935117976422128825857356750737658039E0Q, 137 1.4601391056210009726721818194296893361233E0Q, 138 1.4616428638860188872060496086383008594310E0Q, 139 1.4631064559620759326975975316301202111560E0Q, 140 1.4645314639038178118428450961503371619177E0Q, 141 1.4659193880646627234129855241049975398470E0Q, 142 1.4672716522843522691530527207287398276197E0Q, 143 1.4685896086876430842559640450619880951144E0Q, 144 1.4698745421276027686510391411132998919794E0Q, 145 1.4711276743037345918528755717617308518553E0Q, 146 1.4723501675822635384916444186631899205983E0Q, 147 1.4735431285433308455179928682541563973416E0Q, /* arctan(10.25) */ 148 1.5707963267948966192313216916397514420986E0Q /* pi/2 */ 149 }; 150 151 152 /* arctan t = t + t^3 p(t^2) / q(t^2) 153 |t| <= 0.09375 154 peak relative error 5.3e-37 */ 155 156 static const __float128 157 p0 = -4.283708356338736809269381409828726405572E1Q, 158 p1 = -8.636132499244548540964557273544599863825E1Q, 159 p2 = -5.713554848244551350855604111031839613216E1Q, 160 p3 = -1.371405711877433266573835355036413750118E1Q, 161 p4 = -8.638214309119210906997318946650189640184E-1Q, 162 q0 = 1.285112506901621042780814422948906537959E2Q, 163 q1 = 3.361907253914337187957855834229672347089E2Q, 164 q2 = 3.180448303864130128268191635189365331680E2Q, 165 q3 = 1.307244136980865800160844625025280344686E2Q, 166 q4 = 2.173623741810414221251136181221172551416E1Q; 167 /* q5 = 1.000000000000000000000000000000000000000E0 */ 168 169 static const __float128 huge = 1.0e4930Q; 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 (k <= 0x3fc50000) /* |x| < 2**-58 */ 201 { 202 math_check_force_underflow (x); 203 /* Raise inexact. */ 204 if (huge + x > 0.0) 205 return x; 206 } 207 208 if (k >= 0x40720000) /* |x| > 2**115 */ 209 { 210 /* Saturate result to {-,+}pi/2 */ 211 if (sign) 212 return -atantbl[83]; 213 else 214 return atantbl[83]; 215 } 216 217 if (sign) 218 x = -x; 219 220 if (k >= 0x40024800) /* 10.25 */ 221 { 222 k = 83; 223 t = -1.0/x; 224 } 225 else 226 { 227 /* Index of nearest table element. 228 Roundoff to integer is asymmetrical to avoid cancellation when t < 0 229 (cf. fdlibm). */ 230 k = 8.0 * x + 0.25; 231 u = 0.125Q * k; 232 /* Small arctan argument. */ 233 t = (x - u) / (1.0 + x * u); 234 } 235 236 /* Arctan of small argument t. */ 237 u = t * 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; 241 242 /* arctan x = arctan u + arctan t */ 243 u = atantbl[k] + u; 244 if (sign) 245 return (-u); 246 else 247 return u; 248 } 249