1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 /* 22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #include "libm.h" /* __k_atan2 */ 30 #include "complex_wrapper.h" 31 32 /* 33 * double __k_atan2(double y, double x, double *e) 34 * 35 * Compute atan2 with error terms. 36 * 37 * Important formula: 38 * 3 5 39 * x x 40 * atan(x) = x - ----- + ----- - ... (for x <= 1) 41 * 3 5 42 * 43 * pi 1 1 44 * = --- - --- + --- - ... (for x > 1) 45 * 3 46 * 2 x 3x 47 * 48 * Arg(x + y i) = sign(y) * atan2(|y|, x) 49 * = sign(y) * atan(|y|/x) (for x > 0) 50 * sign(y) * (PI - atan(|y|/|x|)) (for x < 0) 51 * Thus if x >> y (IEEE double: EXP(x) - EXP(y) >= 60): 52 * 1. (x > 0): atan2(y,x) ~ y/x 53 * 2. (x < 0): atan2(y,x) ~ sign(y) (PI - |y/x|)) 54 * Otherwise if x << y: 55 * atan2(y,x) ~ sign(y)*PI/2 - x/y 56 * 57 * __k_atan2 call static functions mx_poly, mx_atan 58 */ 59 60 /* 61 * (void) mx_poly (double *z, double *a, double *e, int n) 62 * return 63 * e = a + z*(a + z*(a + ... z*(a + e)...)) 64 * 0 2 4 2n 65 * Note: 66 * 1. e and coefficient ai are represented by two double numbers. 67 * For e, the first one contain the leading 24 bits rounded, and the 68 * second one contain the remaining 53 bits (total 77 bits accuracy). 69 * For ai, the first one contian the leading 53 bits rounded, and the 70 * second is the remaining 53 bits (total 106 bits accuracy). 71 * 2. z is an array of three doubles. 72 * z[0] : the rounded value of Z (the intended value of z) 73 * z[1] : the leading 24 bits of Z rounded 74 * z[2] : the remaining 53 bits of Z 75 * Note that z[0] = z[1]+z[2] rounded. 76 * 77 */ 78 79 static void 80 mx_poly(const double *z, const double *a, double *e, int n) { 81 double r, s, t, p_h, p_l, z_h, z_l, p; 82 int i; 83 84 n = n + n; 85 p = e[0] + a[n]; 86 p_l = a[n + 1]; 87 p_h = (double) ((float) p); 88 p = a[n - 2] + z[0] * p; 89 z_h = z[1]; z_l = z[2]; 90 p_l += e[0] - (p_h - a[n]); 91 92 for (i = n - 2; i >= 2; i -= 2) { 93 /* compute p = ai + z * p */ 94 t = z_h * p_h; 95 s = z[0] * p_l + p_h * z_l; 96 p_h = (double) ((float) p); 97 s += a[i + 1]; 98 r = t - (p_h - a[i]); 99 p = a[i - 2] + z[0] * p; 100 p_l = r + s; 101 } 102 e[0] = (double)((float) p); 103 t = z_h * p_h; 104 s = z[0] * p_l + p_h * z_l; 105 r = t - (e[0] - a[0]); 106 e[1] = r + s; 107 } 108 109 /* 110 * Table of constants for atan from 0.125 to 8 111 * 0.125 -- 0x3fc00000 --- (increment at bit 16) 112 * 0x3fc10000 113 * 0x3fc20000 114 * ... ... 115 * 0x401f0000 116 * 8.000 -- 0x40200000 (total: 97) 117 * By K.C. Ng, March 9, 1989 118 */ 119 120 static const double TBL_atan_hi[] = { 121 1.243549945467614382e-01, 1.320397616146387620e-01, 1.397088742891636204e-01, 122 1.473614810886516302e-01, 1.549967419239409727e-01, 1.626138285979485676e-01, 123 1.702119252854744080e-01, 1.777902289926760471e-01, 1.853479499956947607e-01, 124 1.928843122579746439e-01, 2.003985538258785115e-01, 2.078899272022629863e-01, 125 2.153576996977380476e-01, 2.228011537593945213e-01, 2.302195872768437179e-01, 126 2.376123138654712419e-01, 2.449786631268641435e-01, 2.596296294082575118e-01, 127 2.741674511196587893e-01, 2.885873618940774099e-01, 3.028848683749714166e-01, 128 3.170557532091470287e-01, 3.310960767041321029e-01, 3.450021772071051318e-01, 129 3.587706702705721895e-01, 3.723984466767542023e-01, 3.858826693980737521e-01, 130 3.992207695752525431e-01, 4.124104415973872673e-01, 4.254496373700422662e-01, 131 4.383365598579578304e-01, 4.510696559885234436e-01, 4.636476090008060935e-01, 132 4.883339510564055352e-01, 5.123894603107377321e-01, 5.358112379604637043e-01, 133 5.585993153435624414e-01, 5.807563535676704136e-01, 6.022873461349641522e-01, 134 6.231993299340659043e-01, 6.435011087932843710e-01, 6.632029927060932861e-01, 135 6.823165548747480713e-01, 7.008544078844501923e-01, 7.188299996216245269e-01, 136 7.362574289814280970e-01, 7.531512809621944138e-01, 7.695264804056582975e-01, 137 7.853981633974482790e-01, 8.156919233162234217e-01, 8.441539861131710509e-01, 138 8.709034570756529758e-01, 8.960553845713439269e-01, 9.197196053504168578e-01, 139 9.420000403794636101e-01, 9.629943306809362058e-01, 9.827937232473290541e-01, 140 1.001483135694234639e+00, 1.019141344266349725e+00, 1.035841253008800145e+00, 141 1.051650212548373764e+00, 1.066630365315743623e+00, 1.080839000541168327e+00, 142 1.094328907321189925e+00, 1.107148717794090409e+00, 1.130953743979160375e+00, 143 1.152571997215667610e+00, 1.172273881128476303e+00, 1.190289949682531656e+00, 144 1.206817370285252489e+00, 1.222025323210989667e+00, 1.236059489478081863e+00, 145 1.249045772398254428e+00, 1.261093382252440387e+00, 1.272297395208717319e+00, 146 1.282740879744270757e+00, 1.292496667789785336e+00, 1.301628834009196156e+00, 147 1.310193935047555547e+00, 1.318242051016837113e+00, 1.325817663668032553e+00, 148 1.339705659598999565e+00, 1.352127380920954636e+00, 1.363300100359693845e+00, 149 1.373400766945015894e+00, 1.382574821490125894e+00, 1.390942827002418447e+00, 150 1.398605512271957618e+00, 1.405647649380269870e+00, 1.412141064608495311e+00, 151 1.418146998399631542e+00, 1.423717971406494032e+00, 1.428899272190732761e+00, 152 1.433730152484709031e+00, 1.438244794498222623e+00, 1.442473099109101931e+00, 153 1.446441332248135092e+00, 154 }; 155 156 static const double TBL_atan_lo[] = { 157 -3.125324142453938311e-18, -1.276925400709959526e-17, 2.479758919089733066e-17, 158 5.409599147666297957e-18, 9.585415594114323829e-18, 7.784470643106252464e-18, 159 -3.541164079802125137e-18, 2.372599351477449041e-17, 4.180692268843078977e-18, 160 2.034098543938166622e-17, 3.139954287184449286e-18, 7.333160666520898500e-18, 161 4.738160130078732886e-19, -5.498822172446843173e-18, 1.231340452914270316e-17, 162 1.058231431371112987e-17, 1.069875561873445139e-17, 1.923875492461530410e-17, 163 8.261353575163771936e-18, -1.428369957377257085e-17, -1.101082790300136900e-17, 164 -1.893928924292642146e-17, -7.952610375793798701e-18, -2.293880475557830393e-17, 165 3.088733564861919217e-17, 1.961231150484565340e-17, 2.378822732491940868e-17, 166 2.246598105617042065e-17, 3.963462895355093301e-17, 2.331553074189288466e-17, 167 -2.494277030626540909e-17, 3.280735600183735558e-17, 2.269877745296168709e-17, 168 -1.137323618932958456e-17, -2.546278147285580353e-17, -4.063795683482557497e-18, 169 -5.455630548591626394e-18, -1.441464378193066908e-17, 2.950430737228402307e-17, 170 2.672403885140095079e-17, 1.583478505144428617e-17, -3.076054864429649001e-17, 171 6.943223671560007740e-18, -1.987626234335816123e-17, -2.147838844445698302e-17, 172 3.473937648299456719e-17, -2.425693465918206812e-17, -3.704991905602721293e-17, 173 3.061616997868383018e-17, -1.071456562778743077e-17, -4.841337011934916763e-17, 174 -2.269823590747287052e-17, 2.923876285774304890e-17, -4.057439412852767923e-17, 175 5.460837485846687627e-17, -3.986660595210752445e-18, 1.390331103123099845e-17, 176 9.438308023545392000e-17, 1.000401886936679889e-17, 3.194313981784503706e-17, 177 -9.650564731467513515e-17, -5.956589637160374564e-17, -1.567632251135907253e-17, 178 -5.490676155022364226e-18, 9.404471373566379412e-17, 7.123833804538446299e-17, 179 -9.159738508900378819e-17, 8.385188614028674371e-17, 7.683333629842068806e-17, 180 4.172467638861439118e-17, -2.979162864892849274e-17, 7.879752739459421280e-17, 181 -2.196203799612310905e-18, 3.242139621534960503e-17, 2.245875015034507026e-17, 182 -9.283188754266129476e-18, -6.830804768926660334e-17, -1.236918499824626670e-17, 183 8.745413734780278834e-17, -6.319394031144676258e-17, -8.824429373951136321e-17, 184 -2.599011860304134377e-17, 2.147674250751150961e-17, 1.093246171526936217e-16, 185 -3.307710355769516504e-17, -3.561490438648230100e-17, -9.843712133488842595e-17, 186 -2.324061182591627982e-17, -8.922630138234492386e-17, -9.573807110557223276e-17, 187 -8.263883782511013632e-17, 8.721870922223967507e-17, -6.457134743238754385e-17, 188 -4.396204466767636187e-17, -2.493019910264565554e-17, -1.105119435430315713e-16, 189 9.211323971545051565e-17, 190 }; 191 192 /* 193 * mx_atan(x,err) 194 * Table look-up algorithm 195 * By K.C. Ng, March 9, 1989 196 * 197 * Algorithm. 198 * 199 * The algorithm is based on atan(x)=atan(y)+atan((x-y)/(1+x*y)). 200 * We use poly1(x) to approximate atan(x) for x in [0,1/8] with 201 * error (relative) 202 * |(atan(x)-poly1(x))/x|<= 2^-83.41 203 * 204 * and use poly2(x) to approximate atan(x) for x in [0,1/65] with 205 * error 206 * |atan(x)-poly2(x)|<= 2^-86.8 207 * 208 * Here poly1 and poly2 are odd polynomial with the following form: 209 * x + x^3*(a1+x^2*(a2+...)) 210 * 211 * (0). Purge off Inf and NaN and 0 212 * (1). Reduce x to positive by atan(x) = -atan(-x). 213 * (2). For x <= 1/8, use 214 * (2.1) if x < 2^(-prec/2), atan(x) = x with inexact flag raised 215 * (2.2) Otherwise 216 * atan(x) = poly1(x) 217 * (3). For x >= 8 then (prec = 78) 218 * (3.1) if x >= 2^prec, atan(x) = atan(inf) - pio2lo 219 * (3.2) if x >= 2^(prec/3), atan(x) = atan(inf) - 1/x 220 * (3.3) if x > 65, atan(x) = atan(inf) - poly2(1/x) 221 * (3.4) Otherwise, atan(x) = atan(inf) - poly1(1/x) 222 * 223 * (4). Now x is in (0.125, 8) 224 * Find y that match x to 4.5 bit after binary (easy). 225 * If iy is the high word of y, then 226 * single : j = (iy - 0x3e000000) >> 19 227 * double : j = (iy - 0x3fc00000) >> 16 228 * quad : j = (iy - 0x3ffc0000) >> 12 229 * 230 * Let s = (x-y)/(1+x*y). Then 231 * atan(x) = atan(y) + poly1(s) 232 * = _TBL_atan_hi[j] + (_TBL_atan_lo[j] + poly2(s) ) 233 * 234 * Note. |s| <= 1.5384615385e-02 = 1/65. Maxium occurs at x = 1.03125 235 * 236 */ 237 238 #define P1 p[2] 239 #define P4 p[8] 240 #define P5 p[9] 241 #define P6 p[10] 242 #define P7 p[11] 243 #define P8 p[12] 244 #define P9 p[13] 245 static const double p[] = { 246 1.0, 247 0.0, 248 -3.33333333333333314830e-01, /* p1 = BFD55555 55555555 */ 249 -1.85030852238476921863e-17, /* p1_l = BC755525 9783A49C */ 250 2.00000000000000011102e-01, /* p2 = 3FC99999 9999999A */ 251 -1.27263196576150347368e-17, /* p2_l = BC6D584B 0D874007 */ 252 -1.42857142857141405923e-01, /* p3 = BFC24924 9249245E */ 253 -1.34258204847170493327e-17, /* p3_l = BC6EF534 A112500D */ 254 1.11111111110486909803e-01, /* p4 = 3FBC71C7 1C71176A */ 255 -9.09090907557387889470e-02, /* p5 = BFB745D1 73B47A7D */ 256 7.69230541541713053189e-02, /* p6 = 3FB3B13A B1E68DE6 */ 257 -6.66645815401964159097e-02, /* p7 = BFB110EE 1584446A */ 258 5.87081768778560317279e-02, /* p8 = 3FAE0EFF 87657733 */ 259 -4.90818147456113240690e-02, /* p9 = BFA92140 6A524B5C */ 260 }; 261 #define Q1 q[2] 262 #define Q3 q[6] 263 #define Q4 q[7] 264 #define Q5 q[8] 265 static const double q[] = { 266 1.0, 267 0.0, 268 -3.33333333333333314830e-01, /* q1 = BFD55555 55555555 */ 269 -1.85022941571278638733e-17, /* q1_l = BC7554E9 D20EFA66 */ 270 1.99999999999999927836e-01, /* q2 = 3FC99999 99999997 */ 271 -1.28782564407438833398e-17, /* q2_l = BC6DB1FB 17217417 */ 272 -1.42857142855492280642e-01, /* q3 = BFC24924 92483C46 */ 273 1.11111097130183356096e-01, /* q4 = 3FBC71C6 E06595CC */ 274 -9.08553303569109294013e-02, /* q5 = BFB7424B 808CDA76 */ 275 }; 276 static const double 277 one = 1.0, 278 pio2hi = 1.570796326794896558e+00, 279 pio2lo = 6.123233995736765886e-17; 280 281 static double 282 mx_atan(double x, double *err) { 283 double y, z, r, s, t, w, s_h, s_l, x_h, x_l, zz[3], ee[2], z_h, 284 z_l, r_h, r_l, u, v; 285 int ix, iy, sign, j; 286 287 ix = ((int *) &x)[HIWORD]; 288 sign = ix & 0x80000000; 289 ix ^= sign; 290 291 /* for |x| < 1/8 */ 292 if (ix < 0x3fc00000) { 293 if (ix < 0x3f300000) { /* when |x| < 2**-12 */ 294 if (ix < 0x3d800000) { /* if |x| < 2**-39 */ 295 *err = (double) ((int) x); 296 return (x); 297 } 298 z = x * x; 299 t = x * z * (q[2] + z * (q[4] + z * q[6])); 300 r = x + t; 301 *err = t - (r - x); 302 return (r); 303 } 304 z = x * x; 305 306 /* use double precision at p4 and on */ 307 ee[0] = z * 308 (P4 + z * 309 (P5 + z * (P6 + z * (P7 + z * (P8 + z * P9))))); 310 311 x_h = (double) ((float) x); 312 z_h = (double) ((float) z); 313 x_l = x - x_h; 314 z_l = (x_h * x_h - z_h); 315 zz[0] = z; 316 zz[1] = z_h; 317 zz[2] = z_l + x_l * (x + x_h); 318 319 /* 320 * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call 321 * mx_poly 322 */ 323 324 mx_poly(zz, p, ee, 3); 325 326 /* finally x*(1+z*(p1+...)) */ 327 r = x_h * ee[0]; 328 t = x * ee[1] + x_l * ee[0]; 329 s = t + r; 330 *err = t - (s - r); 331 return (s); 332 } 333 /* for |x| >= 8.0 */ 334 if (ix >= 0x40200000) { /* x >= 8 */ 335 x = fabs(x); 336 if (ix >= 0x42600000) { /* x >= 2**39 */ 337 if (ix >= 0x44c00000) { /* x >= 2**77 */ 338 y = -pio2lo; 339 } else 340 y = one / x - pio2lo; 341 if (sign == 0) { 342 t = pio2hi - y; 343 *err = -(y - (pio2hi - t)); 344 } else { 345 t = y - pio2hi; 346 *err = y - (pio2hi + t); 347 } 348 return (t); 349 } else { 350 /* compute r = 1/x */ 351 r = one / x; 352 z = r * r; 353 if (ix < 0x40504000) { /* 8 < x < 65 */ 354 355 /* use double precision at p4 and on */ 356 ee[0] = z * 357 (P4 + z * 358 (P5 + z * 359 (P6 + z * (P7 + z * (P8 + z * P9))))); 360 x_h = (double) ((float) x); 361 r_h = (double) ((float) r); 362 z_h = (double) ((float) z); 363 r_l = r * ((x_h - x) * r_h - (x_h * r_h - one)); 364 z_l = (r_h * r_h - z_h); 365 zz[0] = z; 366 zz[1] = z_h; 367 zz[2] = z_l + r_l * (r + r_h); 368 /* 369 * compute (1+z*(p1+z*(p2+z*(p3+e)))) by call 370 * mx_poly 371 */ 372 mx_poly(zz, p, ee, 3); 373 } else { /* x < 65 < 2**39 */ 374 /* use double precision at q3 and on */ 375 ee[0] = z * (Q3 + z * (Q4 + z * Q5)); 376 x_h = (double) ((float) x); 377 r_h = (double) ((float) r); 378 z_h = (double) ((float) z); 379 r_l = r * ((x_h - x) * r_h - (x_h * r_h - one)); 380 z_l = (r_h * r_h - z_h); 381 zz[0] = z; 382 zz[1] = z_h; 383 zz[2] = z_l + r_l * (r + r_h); 384 /* 385 * compute (1+z*(q1+z*(q2+e))) by call 386 * mx_poly 387 */ 388 mx_poly(zz, q, ee, 2); 389 } 390 /* pio2 - r*(1+...) */ 391 v = r_h * ee[0]; 392 t = pio2lo - (r * ee[1] + r_l * ee[0]); 393 if (sign == 0) { 394 s = pio2hi - v; 395 t -= (v - (pio2hi - s)); 396 } else { 397 s = v - pio2hi; 398 t = -(t - (v - (s + pio2hi))); 399 } 400 w = s + t; 401 *err = t - (w - s); 402 return (w); 403 } 404 } 405 /* now x is between 1/8 and 8 */ 406 ((int *) &x)[HIWORD] = ix; 407 iy = (ix + 0x00008000) & 0x7fff0000; 408 ((int *) &y)[HIWORD] = iy; 409 ((int *) &y)[LOWORD] = 0; 410 j = (iy - 0x3fc00000) >> 16; 411 412 w = (x - y); 413 v = 1 / (one + x * y); 414 s = w * v; 415 z = s * s; 416 /* use double precision at q3 and on */ 417 ee[0] = z * (Q3 + z * (Q4 + z * Q5)); 418 s_h = (double) ((float) s); 419 z_h = (double) ((float) z); 420 x_h = (double) ((float) x); 421 t = (double) ((float) (one + x * y)); 422 r = -((x_h - x) * y - (x_h * y - (t - one))); 423 s_l = -v * (s_h * r - (w - s_h * t)); 424 z_l = (s_h * s_h - z_h); 425 zz[0] = z; 426 zz[1] = z_h; 427 zz[2] = z_l + s_l * (s + s_h); 428 /* compute (1+z*(q1+z*(q2+e))) by call mx_poly */ 429 mx_poly(zz, q, ee, 2); 430 v = s_h * ee[0]; 431 t = TBL_atan_lo[j] + (s * ee[1] + s_l * ee[0]); 432 u = TBL_atan_hi[j]; 433 s = u + v; 434 t += (v - (s - u)); 435 w = s + t; 436 *err = t - (w - s); 437 if (sign != 0) { 438 w = -w; 439 *err = -*err; 440 } 441 return (w); 442 } 443 444 static const double 445 twom768 = 6.441148769597133308e-232, /* 2^-768 */ 446 two768 = 1.552518092300708935e+231, /* 2^768 */ 447 pi = 3.1415926535897931159979634685, 448 pi_lo = 1.224646799147353177e-16, 449 pio2 = 1.570796326794896558e+00, 450 pio2_lo = 6.123233995736765886e-17, 451 pio4 = 0.78539816339744827899949, 452 pio4_lo = 3.061616997868382943e-17, 453 pi3o4 = 2.356194490192344836998, 454 pi3o4_lo = 9.184850993605148829195e-17; 455 456 double 457 __k_atan2(double y, double x, double *w) { 458 double t, xh, th, t1, t2, w1, w2; 459 int ix, iy, hx, hy, lx, ly; 460 461 hy = ((int *) &y)[HIWORD]; 462 ly = ((int *) &y)[LOWORD]; 463 iy = hy & ~0x80000000; 464 465 hx = ((int *) &x)[HIWORD]; 466 lx = ((int *) &x)[LOWORD]; 467 ix = hx & ~0x80000000; 468 469 *w = 0.0; 470 if (ix >= 0x7ff00000 || iy >= 0x7ff00000) { /* ignore inexact */ 471 if (isnan(x) || isnan(y)) 472 return (x * y); 473 else if (iy < 0x7ff00000) { 474 if (hx >= 0) { /* ATAN2(+-finite, +inf) is +-0 */ 475 *w *= y; 476 return (*w); 477 } else { /* ATAN2(+-finite, -inf) is +-pi */ 478 *w = copysign(pi_lo, y); 479 return (copysign(pi, y)); 480 } 481 } else if (ix < 0x7ff00000) { 482 /* ATAN2(+-inf, finite) is +-pi/2 */ 483 *w = (hy >= 0)? pio2_lo : -pio2_lo; 484 return ((hy >= 0)? pio2 : -pio2); 485 } else if (hx > 0) { /* ATAN2(+-INF,+INF) = +-pi/4 */ 486 *w = (hy >= 0)? pio4_lo : -pio4_lo; 487 return ((hy >= 0)? pio4 : -pio4); 488 } else { /* ATAN2(+-INF,-INF) = +-3pi/4 */ 489 *w = (hy >= 0)? pi3o4_lo : -pi3o4_lo; 490 return ((hy >= 0)? pi3o4 : -pi3o4); 491 } 492 } else if ((ix | lx) == 0 || (iy | ly) == 0) { 493 if ((iy | ly) == 0) { 494 if (hx >= 0) /* ATAN2(+-0, +(0 <= x <= inf)) is +-0 */ 495 return (y); 496 else { /* ATAN2(+-0, -(0 <= x <= inf)) is +-pi */ 497 *w = (hy >= 0)? pi_lo : -pi_lo; 498 return ((hy >= 0)? pi : -pi); 499 } 500 } else { /* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2 */ 501 *w = (hy >= 0)? pio2_lo : -pio2_lo; 502 return ((hy >= 0)? pio2 : -pio2); 503 } 504 } else if (iy - ix > 0x06400000) { /* |x/y| < 2 ** -100 */ 505 *w = (hy >= 0)? pio2_lo : -pio2_lo; 506 return ((hy >= 0)? pio2 : -pio2); 507 } else if (ix - iy > 0x06400000) { /* |y/x| < 2 ** -100 */ 508 if (hx < 0) { 509 *w = (hy >= 0)? pi_lo : -pi_lo; 510 return ((hy >= 0)? pi : -pi); 511 } else { 512 t = y / x; 513 th = t; 514 ((int *) &th)[LOWORD] &= 0xf8000000; 515 xh = x; 516 ((int *) &xh)[LOWORD] &= 0xf8000000; 517 t1 = (x - xh) * t + xh * (t - th); 518 t2 = y - xh * th; 519 *w = (t2 - t1) / x; 520 return (t); 521 } 522 } else { 523 if (ix >= 0x5f300000) { 524 x *= twom768; 525 y *= twom768; 526 } else if (ix < 0x23d00000) { 527 x *= two768; 528 y *= two768; 529 } 530 y = fabs(y); 531 x = fabs(x); 532 t = y / x; 533 th = t; 534 ((int *) &th)[LOWORD] &= 0xf8000000; 535 xh = x; 536 ((int *) &xh)[LOWORD] &= 0xf8000000; 537 t1 = (x - xh) * t + xh * (t - th); 538 t2 = y - xh * th; 539 w1 = mx_atan(t, &w2); 540 w2 += (t2 - t1) / (x + y * t); 541 if (hx < 0) { 542 t1 = pi - w1; 543 t2 = pi - t1; 544 w2 = (pi_lo - w2) - (w1 - t2); 545 w1 = t1; 546 } 547 *w = (hy >= 0)? w2 : -w2; 548 return ((hy >= 0)? w1 : -w1); 549 } 550 } 551