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 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #ifdef __RESTRICT 31 #define restrict _Restrict 32 #else 33 #define restrict 34 #endif 35 36 extern const double __vlibm_TBL_atan1[]; 37 38 static const double 39 pio4 = 7.8539816339744827900e-01, 40 pio2 = 1.5707963267948965580e+00, 41 pi = 3.1415926535897931160e+00; 42 43 static const float 44 zero = 0.0f, 45 one = 1.0f, 46 q1 = -3.3333333333296428046e-01f, 47 q2 = 1.9999999186853752618e-01f, 48 twop24 = 16777216.0f; 49 50 void 51 __vatan2f(int n, float * restrict y, int stridey, float * restrict x, 52 int stridex, float * restrict z, int stridez) 53 { 54 float x0, x1, x2, y0, y1, y2, *pz0 = 0, *pz1, *pz2; 55 double ah0, ah1, ah2; 56 double t0, t1, t2; 57 double sx0, sx1, sx2; 58 double sign0, sign1, sign2; 59 int i, k0 = 0, k1, k2, hx, sx, sy; 60 int hy0, hy1, hy2; 61 float base0 = 0.0, base1, base2; 62 double num0, num1, num2; 63 double den0, den1, den2; 64 double dx0, dx1, dx2; 65 double dy0, dy1, dy2; 66 double db0, db1, db2; 67 68 do 69 { 70 loop0: 71 hy0 = *(int*)y; 72 hx = *(int*)x; 73 sign0 = one; 74 sy = hy0 & 0x80000000; 75 hy0 &= ~0x80000000; 76 77 sx = hx & 0x80000000; 78 hx &= ~0x80000000; 79 80 if (hy0 > hx) 81 { 82 x0 = *y; 83 y0 = *x; 84 i = hx; 85 hx = hy0; 86 hy0 = i; 87 if (sy) 88 { 89 x0 = -x0; 90 sign0 = -sign0; 91 } 92 if (sx) 93 { 94 y0 = -y0; 95 ah0 = pio2; 96 } 97 else 98 { 99 ah0 = -pio2; 100 sign0 = -sign0; 101 } 102 } 103 else 104 { 105 y0 = *y; 106 x0 = *x; 107 if (sy) 108 { 109 y0 = -y0; 110 sign0 = -sign0; 111 } 112 if (sx) 113 { 114 x0 = -x0; 115 ah0 = -pi; 116 sign0 = -sign0; 117 } 118 else 119 ah0 = zero; 120 } 121 122 if (hx >= 0x7f800000 || hx - hy0 >= 0x0c800000) 123 { 124 if (hx >= 0x7f800000) 125 { 126 if (hx ^ 0x7f800000) /* nan */ 127 ah0 = x0 + y0; 128 else if (hy0 >= 0x7f800000) 129 ah0 += pio4; 130 } 131 else if ((int) ah0 == 0) 132 ah0 = y0 / x0; 133 *z = (sign0 == one) ? ah0 : -ah0; 134 /* sign0*ah0 would change nan behavior relative to previous release */ 135 x += stridex; 136 y += stridey; 137 z += stridez; 138 i = 0; 139 if (--n <= 0) 140 break; 141 goto loop0; 142 } 143 if (hy0 < 0x00800000) { 144 if (hy0 == 0) 145 { 146 *z = sign0 * (float) ah0; 147 x += stridex; 148 y += stridey; 149 z += stridez; 150 i = 0; 151 if (--n <= 0) 152 break; 153 goto loop0; 154 } 155 y0 *= twop24; /* scale subnormal y */ 156 x0 *= twop24; /* scale possibly subnormal x */ 157 hy0 = *(int*)&y0; 158 hx = *(int*)&x0; 159 } 160 pz0 = z; 161 162 k0 = (hy0 - hx + 0x3f800000) & 0xfff80000; 163 if (k0 >= 0x3C800000) /* if |x| >= (1/64)... */ 164 { 165 *(int*)&base0 = k0; 166 k0 = (k0 - 0x3C800000) >> 18; /* (index >> 19) << 1) */ 167 k0 += 4; 168 /* skip over 0,0,pi/2,pi/2 */ 169 } 170 else /* |x| < 1/64 */ 171 { 172 k0 = 0; 173 base0 = zero; 174 } 175 176 x += stridex; 177 y += stridey; 178 z += stridez; 179 i = 1; 180 if (--n <= 0) 181 break; 182 183 184 loop1: 185 hy1 = *(int*)y; 186 hx = *(int*)x; 187 sign1 = one; 188 sy = hy1 & 0x80000000; 189 hy1 &= ~0x80000000; 190 191 sx = hx & 0x80000000; 192 hx &= ~0x80000000; 193 194 if (hy1 > hx) 195 { 196 x1 = *y; 197 y1 = *x; 198 i = hx; 199 hx = hy1; 200 hy1 = i; 201 if (sy) 202 { 203 x1 = -x1; 204 sign1 = -sign1; 205 } 206 if (sx) 207 { 208 y1 = -y1; 209 ah1 = pio2; 210 } 211 else 212 { 213 ah1 = -pio2; 214 sign1 = -sign1; 215 } 216 } 217 else 218 { 219 y1 = *y; 220 x1 = *x; 221 if (sy) 222 { 223 y1 = -y1; 224 sign1 = -sign1; 225 } 226 if (sx) 227 { 228 x1 = -x1; 229 ah1 = -pi; 230 sign1 = -sign1; 231 } 232 else 233 ah1 = zero; 234 } 235 236 if (hx >= 0x7f800000 || hx - hy1 >= 0x0c800000) 237 { 238 if (hx >= 0x7f800000) 239 { 240 if (hx ^ 0x7f800000) /* nan */ 241 ah1 = x1 + y1; 242 else if (hy1 >= 0x7f800000) 243 ah1 += pio4; 244 } 245 else if ((int) ah1 == 0) 246 ah1 = y1 / x1; 247 *z = (sign1 == one)? ah1 : -ah1; 248 x += stridex; 249 y += stridey; 250 z += stridez; 251 i = 1; 252 if (--n <= 0) 253 break; 254 goto loop1; 255 } 256 if (hy1 < 0x00800000) { 257 if (hy1 == 0) 258 { 259 *z = sign1 * (float) ah1; 260 x += stridex; 261 y += stridey; 262 z += stridez; 263 i = 1; 264 if (--n <= 0) 265 break; 266 goto loop1; 267 } 268 y1 *= twop24; /* scale subnormal y */ 269 x1 *= twop24; /* scale possibly subnormal x */ 270 hy1 = *(int*)&y1; 271 hx = *(int*)&x1; 272 } 273 pz1 = z; 274 275 k1 = (hy1 - hx + 0x3f800000) & 0xfff80000; 276 if (k1 >= 0x3C800000) /* if |x| >= (1/64)... */ 277 { 278 *(int*)&base1 = k1; 279 k1 = (k1 - 0x3C800000) >> 18; /* (index >> 19) << 1) */ 280 k1 += 4; 281 /* skip over 0,0,pi/2,pi/2 */ 282 } 283 else /* |x| < 1/64 */ 284 { 285 k1 = 0; 286 base1 = zero; 287 } 288 289 x += stridex; 290 y += stridey; 291 z += stridez; 292 i = 2; 293 if (--n <= 0) 294 break; 295 296 loop2: 297 hy2 = *(int*)y; 298 hx = *(int*)x; 299 sign2 = one; 300 sy = hy2 & 0x80000000; 301 hy2 &= ~0x80000000; 302 303 sx = hx & 0x80000000; 304 hx &= ~0x80000000; 305 306 if (hy2 > hx) 307 { 308 x2 = *y; 309 y2 = *x; 310 i = hx; 311 hx = hy2; 312 hy2 = i; 313 if (sy) 314 { 315 x2 = -x2; 316 sign2 = -sign2; 317 } 318 if (sx) 319 { 320 y2 = -y2; 321 ah2 = pio2; 322 } 323 else 324 { 325 ah2 = -pio2; 326 sign2 = -sign2; 327 } 328 } 329 else 330 { 331 y2 = *y; 332 x2 = *x; 333 if (sy) 334 { 335 y2 = -y2; 336 sign2 = -sign2; 337 } 338 if (sx) 339 { 340 x2 = -x2; 341 ah2 = -pi; 342 sign2 = -sign2; 343 } 344 else 345 ah2 = zero; 346 } 347 348 if (hx >= 0x7f800000 || hx - hy2 >= 0x0c800000) 349 { 350 if (hx >= 0x7f800000) 351 { 352 if (hx ^ 0x7f800000) /* nan */ 353 ah2 = x2 + y2; 354 else if (hy2 >= 0x7f800000) 355 ah2 += pio4; 356 } 357 else if ((int) ah2 == 0) 358 ah2 = y2 / x2; 359 *z = (sign2 == one)? ah2 : -ah2; 360 x += stridex; 361 y += stridey; 362 z += stridez; 363 i = 2; 364 if (--n <= 0) 365 break; 366 goto loop2; 367 } 368 if (hy2 < 0x00800000) { 369 if (hy2 == 0) 370 { 371 *z = sign2 * (float) ah2; 372 x += stridex; 373 y += stridey; 374 z += stridez; 375 i = 2; 376 if (--n <= 0) 377 break; 378 goto loop2; 379 } 380 y2 *= twop24; /* scale subnormal y */ 381 x2 *= twop24; /* scale possibly subnormal x */ 382 hy2 = *(int*)&y2; 383 hx = *(int*)&x2; 384 } 385 386 pz2 = z; 387 388 k2 = (hy2 - hx + 0x3f800000) & 0xfff80000; 389 if (k2 >= 0x3C800000) /* if |x| >= (1/64)... */ 390 { 391 *(int*)&base2 = k2; 392 k2 = (k2 - 0x3C800000) >> 18; /* (index >> 19) << 1) */ 393 k2 += 4; 394 /* skip over 0,0,pi/2,pi/2 */ 395 } 396 else /* |x| < 1/64 */ 397 { 398 k2 = 0; 399 base2 = zero; 400 } 401 402 goto endloop; 403 404 endloop: 405 406 ah2 += __vlibm_TBL_atan1[k2]; 407 ah1 += __vlibm_TBL_atan1[k1]; 408 ah0 += __vlibm_TBL_atan1[k0]; 409 410 db2 = base2; 411 db1 = base1; 412 db0 = base0; 413 dy2 = y2; 414 dy1 = y1; 415 dy0 = y0; 416 dx2 = x2; 417 dx1 = x1; 418 dx0 = x0; 419 420 num2 = dy2 - dx2 * db2; 421 den2 = dx2 + dy2 * db2; 422 423 num1 = dy1 - dx1 * db1; 424 den1 = dx1 + dy1 * db1; 425 426 num0 = dy0 - dx0 * db0; 427 den0 = dx0 + dy0 * db0; 428 429 t2 = num2 / den2; 430 t1 = num1 / den1; 431 t0 = num0 / den0; 432 433 sx2 = t2 * t2; 434 sx1 = t1 * t1; 435 sx0 = t0 * t0; 436 437 t2 += t2 * sx2 * (q1 + sx2 * q2); 438 t1 += t1 * sx1 * (q1 + sx1 * q2); 439 t0 += t0 * sx0 * (q1 + sx0 * q2); 440 441 t2 += ah2; 442 t1 += ah1; 443 t0 += ah0; 444 445 *pz2 = sign2 * t2; 446 *pz1 = sign1 * t1; 447 *pz0 = sign0 * t0; 448 449 x += stridex; 450 y += stridey; 451 z += stridez; 452 i = 0; 453 } while (--n > 0); 454 455 if (i > 1) 456 { 457 ah1 += __vlibm_TBL_atan1[k1]; 458 t1 = (y1 - x1 * (double)base1) / 459 (x1 + y1 * (double)base1); 460 sx1 = t1 * t1; 461 t1 += t1 * sx1 * (q1 + sx1 * q2); 462 t1 += ah1; 463 *pz1 = sign1 * t1; 464 } 465 466 if (i > 0) 467 { 468 ah0 += __vlibm_TBL_atan1[k0]; 469 t0 = (y0 - x0 * (double)base0) / 470 (x0 + y0 * (double)base0); 471 sx0 = t0 * t0; 472 t0 += t0 * sx0 * (q1 + sx0 * q2); 473 t0 += ah0; 474 *pz0 = sign0 * t0; 475 } 476 } 477