Arc tangent function in SSE version 2 Algorithm Explanation: The function atan() is evaluated by reducing the argument domain to the range [0,1] using the identities: atan(x) = PI/2 - atan(1/x) atan(x) = -atan(-x) and then using a rational polynomial approximation to evaluate the function over the reduced domain: a1*x + a2*x^3 + a3*x^5 atan(x) ~ -------------------------- b1 + b2*x^2 + b3*x^
| 366 | // It is possible to further reduce the argument from [0,1] to [0, pi/12] |
| 367 | // which would provide accurate results up to 20 bits of mantissa. |
| 368 | inline __m128 sseAtan(const __m128 x) |
| 369 | { |
| 370 | // Rational polynomial coefficients for the arc tangent approximation. |
| 371 | // Original source: http://www.ganssle.com/approx/approx.pdf |
| 372 | static const __m128 PN_ATAN_A1 = _mm_set1_ps((float) 48.70107004404898384); |
| 373 | static const __m128 PN_ATAN_A2 = _mm_set1_ps((float) 49.5326263772254345); |
| 374 | static const __m128 PN_ATAN_A3 = _mm_set1_ps((float) 9.40604244231624); |
| 375 | static const __m128 PN_ATAN_B1 = _mm_set1_ps((float) 48.70107004404996166); |
| 376 | static const __m128 PN_ATAN_B2 = _mm_set1_ps((float) 65.7663163908956299); |
| 377 | static const __m128 PN_ATAN_B3 = _mm_set1_ps((float) 21.587934067020262); |
| 378 | |
| 379 | // Apply identity atan(x) = -atan(-x) to reduce domain to [0, Inf) |
| 380 | __m128 sign_x = _mm_and_ps(x, ESIGN_MASK); |
| 381 | __m128 abs_x = _mm_and_ps(x, EABS_MASK); |
| 382 | |
| 383 | // Apply identity atan(x) = PI/2 - atan(1/x) to reduce domain to [0,1] |
| 384 | __m128 inv_mask = _mm_cmpgt_ps(abs_x, EONE); |
| 385 | __m128 inv_abs_x = _mm_div_ps(EONE, abs_x); |
| 386 | __m128 norm_x = sseSelect(inv_mask, inv_abs_x, abs_x); |
| 387 | |
| 388 | // compute atan using a normalized input |
| 389 | __m128 norm_x2 = _mm_mul_ps(norm_x, norm_x); |
| 390 | |
| 391 | __m128 num = |
| 392 | _mm_mul_ps( |
| 393 | _mm_add_ps( |
| 394 | _mm_mul_ps( |
| 395 | _mm_add_ps( |
| 396 | _mm_mul_ps(norm_x2, PN_ATAN_A3), |
| 397 | PN_ATAN_A2), |
| 398 | norm_x2), |
| 399 | PN_ATAN_A1), |
| 400 | norm_x); |
| 401 | |
| 402 | __m128 denom = |
| 403 | _mm_add_ps( |
| 404 | _mm_mul_ps( |
| 405 | _mm_add_ps( |
| 406 | _mm_mul_ps( |
| 407 | _mm_add_ps(norm_x2, PN_ATAN_B3), |
| 408 | norm_x2), |
| 409 | PN_ATAN_B2), |
| 410 | norm_x2), |
| 411 | PN_ATAN_B1); |
| 412 | |
| 413 | __m128 res = _mm_div_ps(num, denom); |
| 414 | |
| 415 | // If the input was inverted during domain reduction, |
| 416 | // correct the result by subtracting it from PI/2. |
| 417 | res = sseSelect(inv_mask, _mm_sub_ps(E_PI_2, res), res); |
| 418 | |
| 419 | // If the input was negated during domain reduction, |
| 420 | // correct the result by negating it again. |
| 421 | return _mm_or_ps(sign_x, res); |
| 422 | } |
| 423 | |
| 424 | // Scalar version of sseAtan |
| 425 | inline float sseAtan(const float v) |