12b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This file is part of Eigen, a lightweight C++ template library 22b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for linear algebra. 32b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com) 52b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 62b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This Source Code Form is subject to the terms of the Mozilla 72b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Public License v. 2.0. If a copy of the MPL was not distributed 82b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 92b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifndef EIGEN_MATH_FUNCTIONS_AVX_H 112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#define EIGEN_MATH_FUNCTIONS_AVX_H 122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang/* The sin, cos, exp, and log functions of this file are loosely derived from 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/ 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace Eigen { 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wanginline Packet8i pshiftleft(Packet8i v, int n) 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef EIGEN_VECTORIZE_AVX2 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_slli_epi32(v, n); 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#else 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i lo = _mm_slli_epi32(_mm256_extractf128_si256(v, 0), n); 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i hi = _mm_slli_epi32(_mm256_extractf128_si256(v, 1), n); 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1); 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wanginline Packet8f pshiftright(Packet8f v, int n) 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef EIGEN_VECTORIZE_AVX2 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_cvtepi32_ps(_mm256_srli_epi32(_mm256_castps_si256(v), n)); 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#else 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i lo = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 0), n); 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i hi = _mm_srli_epi32(_mm256_extractf128_si256(_mm256_castps_si256(v), 1), n); 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_cvtepi32_ps(_mm256_insertf128_si256(_mm256_castsi128_si256(lo), (hi), 1)); 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Sine function 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Computes sin(x) by wrapping x to the interval [-Pi/4,3*Pi/4] and 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// evaluating interpolants in [-Pi/4,Pi/4] or [Pi/4,3*Pi/4]. The interpolants 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// are (anti-)symmetric and thus have only odd/even coefficients 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpsin<Packet8f>(const Packet8f& _x) { 502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x = _x; 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Some useful values. 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8i(one, 1); 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f); 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(two, 2.0f); 562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(one_over_four, 0.25f); 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(one_over_pi, 3.183098861837907e-01f); 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(neg_pi_first, -3.140625000000000e+00f); 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(neg_pi_second, -9.670257568359375e-04f); 602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(neg_pi_third, -6.278329571784980e-07f); 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(four_over_pi, 1.273239544735163e+00f); 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Map x from [-Pi/4,3*Pi/4] to z in [-1,3] and subtract the shifted period. 642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f z = pmul(x, p8f_one_over_pi); 652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f shift = _mm256_floor_ps(padd(z, p8f_one_over_four)); 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmadd(shift, p8f_neg_pi_first, x); 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmadd(shift, p8f_neg_pi_second, x); 682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmadd(shift, p8f_neg_pi_third, x); 692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang z = pmul(x, p8f_four_over_pi); 702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Make a mask for the entries that need flipping, i.e. wherever the shift 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // is odd. 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8i shift_ints = _mm256_cvtps_epi32(shift); 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8i shift_isodd = _mm256_castps_si256(_mm256_and_ps(_mm256_castsi256_ps(shift_ints), _mm256_castsi256_ps(p8i_one))); 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8i sign_flip_mask = pshiftleft(shift_isodd, 31); 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Create a mask for which interpolant to use, i.e. if z > 1, then the mask 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // is set to ones for that entry. 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f ival_mask = _mm256_cmp_ps(z, p8f_one, _CMP_GT_OQ); 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Evaluate the polynomial for the interval [1,3] in z. 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_right_0, 9.999999724233232e-01f); 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_right_2, -3.084242535619928e-01f); 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_right_4, 1.584991525700324e-02f); 852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_right_6, -3.188805084631342e-04f); 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f z_minus_two = psub(z, p8f_two); 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f z_minus_two2 = pmul(z_minus_two, z_minus_two); 882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f right = pmadd(p8f_coeff_right_6, z_minus_two2, p8f_coeff_right_4); 892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang right = pmadd(right, z_minus_two2, p8f_coeff_right_2); 902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang right = pmadd(right, z_minus_two2, p8f_coeff_right_0); 912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Evaluate the polynomial for the interval [-1,1] in z. 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_left_1, 7.853981525427295e-01f); 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_left_3, -8.074536727092352e-02f); 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_left_5, 2.489871967827018e-03f); 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(coeff_left_7, -3.587725841214251e-05f); 972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f z2 = pmul(z, z); 982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f left = pmadd(p8f_coeff_left_7, z2, p8f_coeff_left_5); 992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang left = pmadd(left, z2, p8f_coeff_left_3); 1002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang left = pmadd(left, z2, p8f_coeff_left_1); 1012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang left = pmul(left, z); 1022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Assemble the results, i.e. select the left and right polynomials. 1042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang left = _mm256_andnot_ps(ival_mask, left); 1052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang right = _mm256_and_ps(ival_mask, right); 1062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f res = _mm256_or_ps(left, right); 1072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Flip the sign on the odd intervals and return the result. 1092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang res = _mm256_xor_ps(res, _mm256_castsi256_ps(sign_flip_mask)); 1102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return res; 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 1122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Natural logarithm 1142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Computes log(x) as log(2^e * m) = C*e + log(m), where the constant C =log(2) 1152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// and m is in the range [sqrt(1/2),sqrt(2)). In this range, the logarithm can 1162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// be easily approximated by a polynomial centered on m=1 for stability. 1172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// TODO(gonnet): Further reduce the interval allowing for lower-degree 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// polynomial interpolants -> ... -> profit! 1192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangplog<Packet8f>(const Packet8f& _x) { 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x = _x; 1232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f); 1242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f); 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(126f, 126.0f); 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inv_mant_mask, ~0x7f800000); 1282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // The smallest non denormalized float number. 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(min_norm_pos, 0x00800000); 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(minus_inf, 0xff800000); 1322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Polynomial coefficients. 1342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_SQRTHF, 0.707106781186547524f); 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p0, 7.0376836292E-2f); 1362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p1, -1.1514610310E-1f); 1372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p2, 1.1676998740E-1f); 1382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p3, -1.2420140846E-1f); 1392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p4, +1.4249322787E-1f); 1402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p5, -1.6668057665E-1f); 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p6, +2.0000714765E-1f); 1422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p7, -2.4999993993E-1f); 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_p8, +3.3333331174E-1f); 1442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q1, -2.12194440e-4f); 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_log_q2, 0.693359375f); 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f invalid_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_NGE_UQ); // not greater equal is true if x is NaN 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f iszero_mask = _mm256_cmp_ps(x, _mm256_setzero_ps(), _CMP_EQ_OQ); 1492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Truncate input values to the minimum positive normal. 1512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmax(x, p8f_min_norm_pos); 1522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f emm0 = pshiftright(x,23); 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f e = _mm256_sub_ps(emm0, p8f_126f); 1552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Set the exponents to -1, i.e. x are in the range [0.5,1). 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = _mm256_and_ps(x, p8f_inv_mant_mask); 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = _mm256_or_ps(x, p8f_half); 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // part2: Shift the inputs from the range [0.5,1) to [sqrt(1/2),sqrt(2)) 1612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // and shift by -1. The values are then centered around 0, which improves 1622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // the stability of the polynomial evaluation. 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // if( x < SQRTHF ) { 1642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // e -= 1; 1652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // x = x + x - 1.0; 1662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // } else { x = x - 1.0; } 1672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f mask = _mm256_cmp_ps(x, p8f_cephes_SQRTHF, _CMP_LT_OQ); 1682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f tmp = _mm256_and_ps(x, mask); 1692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = psub(x, p8f_1); 1702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang e = psub(e, _mm256_and_ps(p8f_1, mask)); 1712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = padd(x, tmp); 1722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x2 = pmul(x, x); 1742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x3 = pmul(x2, x); 1752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Evaluate the polynomial approximant of degree 8 in three parts, probably 1772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // to improve instruction-level parallelism. 1782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f y, y1, y2; 1792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(p8f_cephes_log_p0, x, p8f_cephes_log_p1); 1802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y1 = pmadd(p8f_cephes_log_p3, x, p8f_cephes_log_p4); 1812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y2 = pmadd(p8f_cephes_log_p6, x, p8f_cephes_log_p7); 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, x, p8f_cephes_log_p2); 1832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y1 = pmadd(y1, x, p8f_cephes_log_p5); 1842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y2 = pmadd(y2, x, p8f_cephes_log_p8); 1852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, x3, y1); 1862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, x3, y2); 1872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmul(y, x3); 1882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Add the logarithm of the exponent back to the result of the interpolation. 1902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y1 = pmul(e, p8f_cephes_log_q1); 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang tmp = pmul(x2, p8f_half); 1922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = padd(y, y1); 1932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = psub(x, tmp); 1942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y2 = pmul(e, p8f_cephes_log_q2); 1952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = padd(x, y); 1962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = padd(x, y2); 1972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Filter out invalid inputs, i.e. negative arg will be NAN, 0 will be -INF. 1992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_or_ps( 2002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _mm256_andnot_ps(iszero_mask, _mm256_or_ps(x, invalid_mask)), 2012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _mm256_and_ps(iszero_mask, p8f_minus_inf)); 2022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 2032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Exponential function. Works by writing "x = m*log(2) + r" where 2052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// "m = floor(x/log(2)+1/2)" and "r" is the remainder. The result is then 2062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1). 2072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 2082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 2092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpexp<Packet8f>(const Packet8f& _x) { 2102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(1, 1.0f); 2112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(half, 0.5f); 2122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(127, 127.0f); 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(exp_hi, 88.3762626647950f); 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(exp_lo, -88.3762626647949f); 2162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_LOG2EF, 1.44269504088896341f); 2182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p0, 1.9875691500E-4f); 2202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p1, 1.3981999507E-3f); 2212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p2, 8.3334519073E-3f); 2222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p3, 4.1665795894E-2f); 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p4, 1.6666665459E-1f); 2242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_p5, 5.0000001201E-1f); 2252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Clamp x. 2272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x = pmax(pmin(_x, p8f_exp_hi), p8f_exp_lo); 2282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Express exp(x) as exp(m*ln(2) + r), start by extracting 2302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // m = floor(x/ln(2) + 0.5). 2312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f m = _mm256_floor_ps(pmadd(x, p8f_cephes_LOG2EF, p8f_half)); 2322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Get r = x - m*ln(2). If no FMA instructions are available, m*ln(2) is 2342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// subtracted out in two parts, m*C1+m*C2 = m*ln(2), to avoid accumulating 2352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// truncation errors. Note that we don't use the "pmadd" function here to 2362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// ensure that a precision-preserving FMA instruction is used. 2372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifdef EIGEN_VECTORIZE_FMA 2382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(nln2, -0.6931471805599453f); 2392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f r = _mm256_fmadd_ps(m, p8f_nln2, x); 2402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#else 2412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C1, 0.693359375f); 2422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(cephes_exp_C2, -2.12194440e-4f); 2432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f r = psub(x, pmul(m, p8f_cephes_exp_C1)); 2442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang r = psub(r, pmul(m, p8f_cephes_exp_C2)); 2452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 2462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f r2 = pmul(r, r); 2482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // TODO(gonnet): Split into odd/even polynomials and try to exploit 2502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // instruction-level parallelism. 2512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f y = p8f_cephes_exp_p0; 2522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r, p8f_cephes_exp_p1); 2532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r, p8f_cephes_exp_p2); 2542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r, p8f_cephes_exp_p3); 2552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r, p8f_cephes_exp_p4); 2562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r, p8f_cephes_exp_p5); 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = pmadd(y, r2, r); 2582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang y = padd(y, p8f_1); 2592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Build emm0 = 2^m. 2612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8i emm0 = _mm256_cvttps_epi32(padd(m, p8f_127)); 2622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang emm0 = pshiftleft(emm0, 23); 2632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Return 2^m * exp(r). 2652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return pmax(pmul(y, _mm256_castsi256_ps(emm0)), _x); 2662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 2672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Hyperbolic Tangent function. 2692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 2702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 2712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangptanh<Packet8f>(const Packet8f& x) { 2722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return internal::generic_fast_tanh_float(x); 2732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 2742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 2762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4d 2772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpexp<Packet4d>(const Packet4d& _x) { 2782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d x = _x; 2792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(1, 1.0); 2812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(2, 2.0); 2822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(half, 0.5); 2832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(exp_hi, 709.437); 2852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(exp_lo, -709.436139303); 2862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_LOG2EF, 1.4426950408889634073599); 2882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p0, 1.26177193074810590878e-4); 2902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p1, 3.02994407707441961300e-2); 2912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_p2, 9.99999999999999999910e-1); 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q0, 3.00198505138664455042e-6); 2942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q1, 2.52448340349684104192e-3); 2952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q2, 2.27265548208155028766e-1); 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_q3, 2.00000000000000000009e0); 2972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C1, 0.693145751953125); 2992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(cephes_exp_C2, 1.42860682030941723212e-6); 3002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4i(1023, 1023); 3012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d tmp, fx; 3032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // clamp x 3052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmax(pmin(x, p4d_exp_hi), p4d_exp_lo); 3062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Express exp(x) as exp(g + n*log(2)). 3072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang fx = pmadd(p4d_cephes_LOG2EF, x, p4d_half); 3082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Get the integer modulus of log(2), i.e. the "n" described above. 3102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang fx = _mm256_floor_pd(fx); 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Get the remainder modulo log(2), i.e. the "g" described above. Subtract 3132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // n*log(2) out in two steps, i.e. n*C1 + n*C2, C1+C2=log2 to get the last 3142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // digits right. 3152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang tmp = pmul(fx, p4d_cephes_exp_C1); 3162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d z = pmul(fx, p4d_cephes_exp_C2); 3172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = psub(x, tmp); 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = psub(x, z); 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d x2 = pmul(x, x); 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Evaluate the numerator polynomial of the rational interpolant. 3232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d px = p4d_cephes_exp_p0; 3242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang px = pmadd(px, x2, p4d_cephes_exp_p1); 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang px = pmadd(px, x2, p4d_cephes_exp_p2); 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang px = pmul(px, x); 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Evaluate the denominator polynomial of the rational interpolant. 3292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet4d qx = p4d_cephes_exp_q0; 3302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang qx = pmadd(qx, x2, p4d_cephes_exp_q1); 3312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang qx = pmadd(qx, x2, p4d_cephes_exp_q2); 3322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang qx = pmadd(qx, x2, p4d_cephes_exp_q3); 3332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // I don't really get this bit, copied from the SSE2 routines, so... 3352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // TODO(gonnet): Figure out what is going on here, perhaps find a better 3362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // rational interpolant? 3372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = _mm256_div_pd(px, psub(qx, px)); 3382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmadd(p4d_2, x, p4d_1); 3392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Build e=2^n by constructing the exponents in a 128-bit vector and 3412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // shifting them to where they belong in double-precision values. 3422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i emm0 = _mm256_cvtpd_epi32(fx); 3432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang emm0 = _mm_add_epi32(emm0, p4i_1023); 3442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(3, 1, 2, 0)); 3452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i lo = _mm_slli_epi64(emm0, 52); 3462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m128i hi = _mm_slli_epi64(_mm_srli_epi64(emm0, 32), 52); 3472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang __m256i e = _mm256_insertf128_si256(_mm256_setzero_si256(), lo, 0); 3482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang e = _mm256_insertf128_si256(e, hi, 1); 3492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Construct the result 2^n * exp(g) = e * x. The max is used to catch 3512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // non-finite values in the input. 3522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return pmax(pmul(x, _mm256_castsi256_pd(e)), _x); 3532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 3542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Functions for sqrt. 3562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step 3572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// of Newton's method, at a cost of 1-2 bits of precision as opposed to the 3582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// exact solution. It does not handle +inf, or denormalized numbers correctly. 3592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// The main advantage of this approach is not just speed, but also the fact that 3602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// it can be inlined and pipelined with other computations, further reducing its 3612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// effective latency. This is similar to Quake3's fast inverse square root. 3622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// For detail see here: http://www.beyond3d.com/content/articles/8/ 3632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#if EIGEN_FAST_MATH 3642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> 3652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f 3662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangpsqrt<Packet8f>(const Packet8f& _x) { 3672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f half = pmul(_x, pset1<Packet8f>(.5f)); 3682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f denormal_mask = _mm256_and_ps( 3692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _mm256_cmp_ps(_x, pset1<Packet8f>((std::numeric_limits<float>::min)()), 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _CMP_LT_OQ), 3712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_GE_OQ)); 3722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Compute approximate reciprocal sqrt. 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x = _mm256_rsqrt_ps(_x); 3752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Do a single step of Newton's iteration. 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmul(x, psub(pset1<Packet8f>(1.5f), pmul(half, pmul(x,x)))); 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Flush results for denormals to zero. 3782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_andnot_ps(denormal_mask, pmul(_x,x)); 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#else 3812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 3822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangPacket8f psqrt<Packet8f>(const Packet8f& x) { 3832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_sqrt_ps(x); 3842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 3852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 3862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangPacket4d psqrt<Packet4d>(const Packet4d& x) { 3882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_sqrt_pd(x); 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 3902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#if EIGEN_FAST_MATH 3912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 3932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangPacket8f prsqrt<Packet8f>(const Packet8f& _x) { 3942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(inf, 0x7f800000); 3952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(nan, 0x7fc00000); 3962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f); 3972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f); 3982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000); 3992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f neg_half = pmul(_x, p8f_minus_half); 4012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // select only the inverse sqrt of positive normal inputs (denormals are 4032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // flushed to zero and cause infs as well). 4042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f le_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_LT_OQ); 4052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f x = _mm256_andnot_ps(le_zero_mask, _mm256_rsqrt_ps(_x)); 4062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Fill in NaNs and Infs for the negative/zero entries. 4082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f neg_mask = _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_LT_OQ); 4092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f zero_mask = _mm256_andnot_ps(neg_mask, le_zero_mask); 4102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Packet8f infs_and_nans = _mm256_or_ps(_mm256_and_ps(neg_mask, p8f_nan), 4112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _mm256_and_ps(zero_mask, p8f_inf)); 4122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Do a single step of Newton's iteration. 4142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five)); 4152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Insert NaNs and Infs in all the right places. 4172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_or_ps(x, infs_and_nans); 4182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 4192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#else 4212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 4222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangPacket8f prsqrt<Packet8f>(const Packet8f& x) { 4232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet8f(one, 1.0f); 4242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_div_ps(p8f_one, _mm256_sqrt_ps(x)); 4252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 4262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif 4272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangPacket4d prsqrt<Packet4d>(const Packet4d& x) { 4302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang _EIGEN_DECLARE_CONST_Packet4d(one, 1.0); 4312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return _mm256_div_pd(p4d_one, _mm256_sqrt_pd(x)); 4322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 4332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 4362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace Eigen 4382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif // EIGEN_MATH_FUNCTIONS_AVX_H 440