From 145c243c32a8f58f1bf9ebf15e214c7fe6bbb530 Mon Sep 17 00:00:00 2001 From: Paul Ferrand Date: Sun, 31 May 2020 12:59:52 +0200 Subject: [PATCH] Removed mathfuns (now unused) --- README.md | 2 - src/external/mathfuns/neon_mathfun.h | 301 ----------- src/external/mathfuns/sse_mathfun.h | 713 --------------------------- 3 files changed, 1016 deletions(-) delete mode 100644 src/external/mathfuns/neon_mathfun.h delete mode 100644 src/external/mathfuns/sse_mathfun.h diff --git a/README.md b/README.md index 60ab196c..4ee9292a 100644 --- a/README.md +++ b/README.md @@ -21,8 +21,6 @@ The sfizz library also uses in some subprojects: - [benchmark], licensed under the Apache License 2.0 - [LV2], licensed under the ISC license - [JACK], licensed under the GNU Lesser General Public License v2.1 -- `neon_mathfun.h` and `sse_mathfun.h` by Julien Pommier, - licensed under the zlib license [Abseil]: https://github.com/abseil/abseil-cpp [atomic_queue]: https://github.com/max0x7ba/atomic_queue diff --git a/src/external/mathfuns/neon_mathfun.h b/src/external/mathfuns/neon_mathfun.h deleted file mode 100644 index f51a61e1..00000000 --- a/src/external/mathfuns/neon_mathfun.h +++ /dev/null @@ -1,301 +0,0 @@ -/* NEON implementation of sin, cos, exp and log - - Inspired by Intel Approximate Math library, and based on the - corresponding algorithms of the cephes math library -*/ - -/* Copyright (C) 2011 Julien Pommier - - This software is provided 'as-is', without any express or implied - warranty. In no event will the authors be held liable for any damages - arising from the use of this software. - - Permission is granted to anyone to use this software for any purpose, - including commercial applications, and to alter it and redistribute it - freely, subject to the following restrictions: - - 1. The origin of this software must not be misrepresented; you must not - claim that you wrote the original software. If you use this software - in a product, an acknowledgment in the product documentation would be - appreciated but is not required. - 2. Altered source versions must be plainly marked as such, and must not be - misrepresented as being the original software. - 3. This notice may not be removed or altered from any source distribution. - - (this is the zlib license) -*/ - -#include - -typedef float32x4_t v4sf; // vector of 4 float -typedef uint32x4_t v4su; // vector of 4 uint32 -typedef int32x4_t v4si; // vector of 4 uint32 - -#define c_inv_mant_mask ~0x7f800000u -#define c_cephes_SQRTHF 0.707106781186547524 -#define c_cephes_log_p0 7.0376836292E-2 -#define c_cephes_log_p1 - 1.1514610310E-1 -#define c_cephes_log_p2 1.1676998740E-1 -#define c_cephes_log_p3 - 1.2420140846E-1 -#define c_cephes_log_p4 + 1.4249322787E-1 -#define c_cephes_log_p5 - 1.6668057665E-1 -#define c_cephes_log_p6 + 2.0000714765E-1 -#define c_cephes_log_p7 - 2.4999993993E-1 -#define c_cephes_log_p8 + 3.3333331174E-1 -#define c_cephes_log_q1 -2.12194440e-4 -#define c_cephes_log_q2 0.693359375 - -/* natural logarithm computed for 4 simultaneous float - return NaN for x <= 0 -*/ -v4sf log_ps(v4sf x) { - v4sf one = vdupq_n_f32(1); - - x = vmaxq_f32(x, vdupq_n_f32(0)); /* force flush to zero on denormal values */ - v4su invalid_mask = vcleq_f32(x, vdupq_n_f32(0)); - - v4si ux = vreinterpretq_s32_f32(x); - - v4si emm0 = vshrq_n_s32(ux, 23); - - /* keep only the fractional part */ - ux = vandq_s32(ux, vdupq_n_s32(c_inv_mant_mask)); - ux = vorrq_s32(ux, vreinterpretq_s32_f32(vdupq_n_f32(0.5f))); - x = vreinterpretq_f32_s32(ux); - - emm0 = vsubq_s32(emm0, vdupq_n_s32(0x7f)); - v4sf e = vcvtq_f32_s32(emm0); - - e = vaddq_f32(e, one); - - /* part2: - if( x < SQRTHF ) { - e -= 1; - x = x + x - 1.0; - } else { x = x - 1.0; } - */ - v4su mask = vcltq_f32(x, vdupq_n_f32(c_cephes_SQRTHF)); - v4sf tmp = vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(x), mask)); - x = vsubq_f32(x, one); - e = vsubq_f32(e, vreinterpretq_f32_u32(vandq_u32(vreinterpretq_u32_f32(one), mask))); - x = vaddq_f32(x, tmp); - - v4sf z = vmulq_f32(x,x); - - v4sf y = vdupq_n_f32(c_cephes_log_p0); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p1)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p2)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p3)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p4)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p5)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p6)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p7)); - y = vmulq_f32(y, x); - y = vaddq_f32(y, vdupq_n_f32(c_cephes_log_p8)); - y = vmulq_f32(y, x); - - y = vmulq_f32(y, z); - - - tmp = vmulq_f32(e, vdupq_n_f32(c_cephes_log_q1)); - y = vaddq_f32(y, tmp); - - - tmp = vmulq_f32(z, vdupq_n_f32(0.5f)); - y = vsubq_f32(y, tmp); - - tmp = vmulq_f32(e, vdupq_n_f32(c_cephes_log_q2)); - x = vaddq_f32(x, y); - x = vaddq_f32(x, tmp); - x = vreinterpretq_f32_u32(vorrq_u32(vreinterpretq_u32_f32(x), invalid_mask)); // negative arg will be NAN - return x; -} - -#define c_exp_hi 88.3762626647949f -#define c_exp_lo -88.3762626647949f - -#define c_cephes_LOG2EF 1.44269504088896341 -#define c_cephes_exp_C1 0.693359375 -#define c_cephes_exp_C2 -2.12194440e-4 - -#define c_cephes_exp_p0 1.9875691500E-4 -#define c_cephes_exp_p1 1.3981999507E-3 -#define c_cephes_exp_p2 8.3334519073E-3 -#define c_cephes_exp_p3 4.1665795894E-2 -#define c_cephes_exp_p4 1.6666665459E-1 -#define c_cephes_exp_p5 5.0000001201E-1 - -/* exp() computed for 4 float at once */ -v4sf exp_ps(v4sf x) { - v4sf tmp, fx; - - v4sf one = vdupq_n_f32(1); - x = vminq_f32(x, vdupq_n_f32(c_exp_hi)); - x = vmaxq_f32(x, vdupq_n_f32(c_exp_lo)); - - /* express exp(x) as exp(g + n*log(2)) */ - fx = vmlaq_f32(vdupq_n_f32(0.5f), x, vdupq_n_f32(c_cephes_LOG2EF)); - - /* perform a floorf */ - tmp = vcvtq_f32_s32(vcvtq_s32_f32(fx)); - - /* if greater, substract 1 */ - v4su mask = vcgtq_f32(tmp, fx); - mask = vandq_u32(mask, vreinterpretq_u32_f32(one)); - - - fx = vsubq_f32(tmp, vreinterpretq_f32_u32(mask)); - - tmp = vmulq_f32(fx, vdupq_n_f32(c_cephes_exp_C1)); - v4sf z = vmulq_f32(fx, vdupq_n_f32(c_cephes_exp_C2)); - x = vsubq_f32(x, tmp); - x = vsubq_f32(x, z); - - static const float cephes_exp_p[6] = { c_cephes_exp_p0, c_cephes_exp_p1, c_cephes_exp_p2, c_cephes_exp_p3, c_cephes_exp_p4, c_cephes_exp_p5 }; - v4sf y = vld1q_dup_f32(cephes_exp_p+0); - v4sf c1 = vld1q_dup_f32(cephes_exp_p+1); - v4sf c2 = vld1q_dup_f32(cephes_exp_p+2); - v4sf c3 = vld1q_dup_f32(cephes_exp_p+3); - v4sf c4 = vld1q_dup_f32(cephes_exp_p+4); - v4sf c5 = vld1q_dup_f32(cephes_exp_p+5); - - y = vmulq_f32(y, x); - z = vmulq_f32(x,x); - y = vaddq_f32(y, c1); - y = vmulq_f32(y, x); - y = vaddq_f32(y, c2); - y = vmulq_f32(y, x); - y = vaddq_f32(y, c3); - y = vmulq_f32(y, x); - y = vaddq_f32(y, c4); - y = vmulq_f32(y, x); - y = vaddq_f32(y, c5); - - y = vmulq_f32(y, z); - y = vaddq_f32(y, x); - y = vaddq_f32(y, one); - - /* build 2^n */ - int32x4_t mm; - mm = vcvtq_s32_f32(fx); - mm = vaddq_s32(mm, vdupq_n_s32(0x7f)); - mm = vshlq_n_s32(mm, 23); - v4sf pow2n = vreinterpretq_f32_s32(mm); - - y = vmulq_f32(y, pow2n); - return y; -} - -#define c_minus_cephes_DP1 -0.78515625 -#define c_minus_cephes_DP2 -2.4187564849853515625e-4 -#define c_minus_cephes_DP3 -3.77489497744594108e-8 -#define c_sincof_p0 -1.9515295891E-4 -#define c_sincof_p1 8.3321608736E-3 -#define c_sincof_p2 -1.6666654611E-1 -#define c_coscof_p0 2.443315711809948E-005 -#define c_coscof_p1 -1.388731625493765E-003 -#define c_coscof_p2 4.166664568298827E-002 -#define c_cephes_FOPI 1.27323954473516 // 4 / M_PI - -/* evaluation of 4 sines & cosines at once. - - The code is the exact rewriting of the cephes sinf function. - Precision is excellent as long as x < 8192 (I did not bother to - take into account the special handling they have for greater values - -- it does not return garbage for arguments over 8192, though, but - the extra precision is missing). - - Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the - surprising but correct result. - - Note also that when you compute sin(x), cos(x) is available at - almost no extra price so both sin_ps and cos_ps make use of - sincos_ps.. - */ -void sincos_ps(v4sf x, v4sf *ysin, v4sf *ycos) { // any x - v4sf xmm1, xmm2, xmm3, y; - - v4su emm2; - - v4su sign_mask_sin, sign_mask_cos; - sign_mask_sin = vcltq_f32(x, vdupq_n_f32(0)); - x = vabsq_f32(x); - - /* scale by 4/Pi */ - y = vmulq_f32(x, vdupq_n_f32(c_cephes_FOPI)); - - /* store the integer part of y in mm0 */ - emm2 = vcvtq_u32_f32(y); - /* j=(j+1) & (~1) (see the cephes sources) */ - emm2 = vaddq_u32(emm2, vdupq_n_u32(1)); - emm2 = vandq_u32(emm2, vdupq_n_u32(~1)); - y = vcvtq_f32_u32(emm2); - - /* get the polynom selection mask - there is one polynom for 0 <= x <= Pi/4 - and another one for Pi/4 - -/* yes I know, the top of this file is quite ugly */ - -#ifdef _MSC_VER /* visual c++ */ -# define ALIGN16_BEG __declspec(align(16)) -# define ALIGN16_END -#else /* gcc or icc */ -# define ALIGN16_BEG -# define ALIGN16_END __attribute__((aligned(16))) -#endif - -#define USE_SSE2 - -/* __m128 is ugly to write */ -typedef __m128 v4sf; // vector of 4 float (sse1) - -#ifdef USE_SSE2 -# include -typedef __m128i v4si; // vector of 4 int (sse2) -#else -typedef __m64 v2si; // vector of 2 int (mmx) -#endif - -/* declare some SSE constants -- why can't I figure a better way to do that? */ -#define _PS_CONST(Name, Val) \ - static const ALIGN16_BEG float _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val } -#define _PI32_CONST(Name, Val) \ - static const ALIGN16_BEG int _pi32_##Name[4] ALIGN16_END = { Val, Val, Val, Val } -#define _PS_CONST_TYPE(Name, Type, Val) \ - static const ALIGN16_BEG Type _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val } - -_PS_CONST(1 , 1.0f); -_PS_CONST(0p5, 0.5f); -/* the smallest non denormalized float number */ -_PS_CONST_TYPE(min_norm_pos, int, 0x00800000); -_PS_CONST_TYPE(mant_mask, int, 0x7f800000); -_PS_CONST_TYPE(inv_mant_mask, int, ~0x7f800000); - -_PS_CONST_TYPE(sign_mask, int, (int)0x80000000); -_PS_CONST_TYPE(inv_sign_mask, int, ~0x80000000); - -_PI32_CONST(1, 1); -_PI32_CONST(inv1, ~1); -_PI32_CONST(2, 2); -_PI32_CONST(4, 4); -_PI32_CONST(0x7f, 0x7f); - -_PS_CONST(cephes_SQRTHF, 0.707106781186547524f); -_PS_CONST(cephes_log_p0, 7.0376836292E-2f); -_PS_CONST(cephes_log_p1, - 1.1514610310E-1f); -_PS_CONST(cephes_log_p2, 1.1676998740E-1f); -_PS_CONST(cephes_log_p3, - 1.2420140846E-1f); -_PS_CONST(cephes_log_p4, + 1.4249322787E-1f); -_PS_CONST(cephes_log_p5, - 1.6668057665E-1f); -_PS_CONST(cephes_log_p6, + 2.0000714765E-1f); -_PS_CONST(cephes_log_p7, - 2.4999993993E-1f); -_PS_CONST(cephes_log_p8, + 3.3333331174E-1f); -_PS_CONST(cephes_log_q1, -2.12194440e-4f); -_PS_CONST(cephes_log_q2, 0.693359375f); - -#ifndef USE_SSE2 -typedef union xmm_mm_union { - __m128 xmm; - __m64 mm[2]; -} xmm_mm_union; - -#define COPY_XMM_TO_MM(xmm_, mm0_, mm1_) { \ - xmm_mm_union u; u.xmm = xmm_; \ - mm0_ = u.mm[0]; \ - mm1_ = u.mm[1]; \ -} - -#define COPY_MM_TO_XMM(mm0_, mm1_, xmm_) { \ - xmm_mm_union u; u.mm[0]=mm0_; u.mm[1]=mm1_; xmm_ = u.xmm; \ - } - -#endif // USE_SSE2 - -/* natural logarithm computed for 4 simultaneous float - return NaN for x <= 0 -*/ -v4sf log_ps(v4sf x) { -#ifdef USE_SSE2 - v4si emm0; -#else - v2si mm0, mm1; -#endif - v4sf one = *(v4sf*)_ps_1; - - v4sf invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps()); - - x = _mm_max_ps(x, *(v4sf*)_ps_min_norm_pos); /* cut off denormalized stuff */ - -#ifndef USE_SSE2 - /* part 1: x = frexpf(x, &e); */ - COPY_XMM_TO_MM(x, mm0, mm1); - mm0 = _mm_srli_pi32(mm0, 23); - mm1 = _mm_srli_pi32(mm1, 23); -#else - emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23); -#endif - /* keep only the fractional part */ - x = _mm_and_ps(x, *(v4sf*)_ps_inv_mant_mask); - x = _mm_or_ps(x, *(v4sf*)_ps_0p5); - -#ifndef USE_SSE2 - /* now e=mm0:mm1 contain the really base-2 exponent */ - mm0 = _mm_sub_pi32(mm0, *(v2si*)_pi32_0x7f); - mm1 = _mm_sub_pi32(mm1, *(v2si*)_pi32_0x7f); - v4sf e = _mm_cvtpi32x2_ps(mm0, mm1); - _mm_empty(); /* bye bye mmx */ -#else - emm0 = _mm_sub_epi32(emm0, *(v4si*)_pi32_0x7f); - v4sf e = _mm_cvtepi32_ps(emm0); -#endif - - e = _mm_add_ps(e, one); - - /* part2: - if( x < SQRTHF ) { - e -= 1; - x = x + x - 1.0; - } else { x = x - 1.0; } - */ - v4sf mask = _mm_cmplt_ps(x, *(v4sf*)_ps_cephes_SQRTHF); - v4sf tmp = _mm_and_ps(x, mask); - x = _mm_sub_ps(x, one); - e = _mm_sub_ps(e, _mm_and_ps(one, mask)); - x = _mm_add_ps(x, tmp); - - - v4sf z = _mm_mul_ps(x,x); - - v4sf y = *(v4sf*)_ps_cephes_log_p0; - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p1); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p2); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p3); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p4); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p5); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p6); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p7); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p8); - y = _mm_mul_ps(y, x); - - y = _mm_mul_ps(y, z); - - - tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q1); - y = _mm_add_ps(y, tmp); - - - tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5); - y = _mm_sub_ps(y, tmp); - - tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q2); - x = _mm_add_ps(x, y); - x = _mm_add_ps(x, tmp); - x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN - return x; -} - -_PS_CONST(exp_hi, 88.3762626647949f); -_PS_CONST(exp_lo, -88.3762626647949f); - -_PS_CONST(cephes_LOG2EF, 1.44269504088896341f); -_PS_CONST(cephes_exp_C1, 0.693359375f); -_PS_CONST(cephes_exp_C2, -2.12194440e-4f); - -_PS_CONST(cephes_exp_p0, 1.9875691500E-4f); -_PS_CONST(cephes_exp_p1, 1.3981999507E-3f); -_PS_CONST(cephes_exp_p2, 8.3334519073E-3f); -_PS_CONST(cephes_exp_p3, 4.1665795894E-2f); -_PS_CONST(cephes_exp_p4, 1.6666665459E-1f); -_PS_CONST(cephes_exp_p5, 5.0000001201E-1f); - -v4sf exp_ps(v4sf x) { - v4sf tmp = _mm_setzero_ps(), fx; -#ifdef USE_SSE2 - v4si emm0; -#else - v2si mm0, mm1; -#endif - v4sf one = *(v4sf*)_ps_1; - - x = _mm_min_ps(x, *(v4sf*)_ps_exp_hi); - x = _mm_max_ps(x, *(v4sf*)_ps_exp_lo); - - /* express exp(x) as exp(g + n*log(2)) */ - fx = _mm_mul_ps(x, *(v4sf*)_ps_cephes_LOG2EF); - fx = _mm_add_ps(fx, *(v4sf*)_ps_0p5); - - /* how to perform a floorf with SSE: just below */ -#ifndef USE_SSE2 - /* step 1 : cast to int */ - tmp = _mm_movehl_ps(tmp, fx); - mm0 = _mm_cvttps_pi32(fx); - mm1 = _mm_cvttps_pi32(tmp); - /* step 2 : cast back to float */ - tmp = _mm_cvtpi32x2_ps(mm0, mm1); -#else - emm0 = _mm_cvttps_epi32(fx); - tmp = _mm_cvtepi32_ps(emm0); -#endif - /* if greater, substract 1 */ - v4sf mask = _mm_cmpgt_ps(tmp, fx); - mask = _mm_and_ps(mask, one); - fx = _mm_sub_ps(tmp, mask); - - tmp = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C1); - v4sf z = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C2); - x = _mm_sub_ps(x, tmp); - x = _mm_sub_ps(x, z); - - z = _mm_mul_ps(x,x); - - v4sf y = *(v4sf*)_ps_cephes_exp_p0; - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p1); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p2); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p3); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p4); - y = _mm_mul_ps(y, x); - y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p5); - y = _mm_mul_ps(y, z); - y = _mm_add_ps(y, x); - y = _mm_add_ps(y, one); - - /* build 2^n */ -#ifndef USE_SSE2 - z = _mm_movehl_ps(z, fx); - mm0 = _mm_cvttps_pi32(fx); - mm1 = _mm_cvttps_pi32(z); - mm0 = _mm_add_pi32(mm0, *(v2si*)_pi32_0x7f); - mm1 = _mm_add_pi32(mm1, *(v2si*)_pi32_0x7f); - mm0 = _mm_slli_pi32(mm0, 23); - mm1 = _mm_slli_pi32(mm1, 23); - - v4sf pow2n; - COPY_MM_TO_XMM(mm0, mm1, pow2n); - _mm_empty(); -#else - emm0 = _mm_cvttps_epi32(fx); - emm0 = _mm_add_epi32(emm0, *(v4si*)_pi32_0x7f); - emm0 = _mm_slli_epi32(emm0, 23); - v4sf pow2n = _mm_castsi128_ps(emm0); -#endif - y = _mm_mul_ps(y, pow2n); - return y; -} - -_PS_CONST(minus_cephes_DP1, -0.78515625f); -_PS_CONST(minus_cephes_DP2, -2.4187564849853515625e-4f); -_PS_CONST(minus_cephes_DP3, -3.77489497744594108e-8f); -_PS_CONST(sincof_p0, -1.9515295891E-4f); -_PS_CONST(sincof_p1, 8.3321608736E-3f); -_PS_CONST(sincof_p2, -1.6666654611E-1f); -_PS_CONST(coscof_p0, 2.443315711809948E-005f); -_PS_CONST(coscof_p1, -1.388731625493765E-003f); -_PS_CONST(coscof_p2, 4.166664568298827E-002f); -_PS_CONST(cephes_FOPI, 1.27323954473516f); // 4 / M_PI - - -/* evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so - it runs also on old athlons XPs and the pentium III of your grand - mother. - - The code is the exact rewriting of the cephes sinf function. - Precision is excellent as long as x < 8192 (I did not bother to - take into account the special handling they have for greater values - -- it does not return garbage for arguments over 8192, though, but - the extra precision is missing). - - Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the - surprising but correct result. - - Performance is also surprisingly good, 1.33 times faster than the - macos vsinf SSE2 function, and 1.5 times faster than the - __vrs4_sinf of amd's ACML (which is only available in 64 bits). Not - too bad for an SSE1 function (with no special tuning) ! - However the latter libraries probably have a much better handling of NaN, - Inf, denormalized and other special arguments.. - - On my core 1 duo, the execution of this function takes approximately 95 cycles. - - From what I have observed on the experiments with Intel AMath lib, switching to an - SSE2 version would improve the perf by only 10%. - - Since it is based on SSE intrinsics, it has to be compiled at -O2 to - deliver full speed. -*/ -v4sf sin_ps(v4sf x) { // any x - v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y; - -#ifdef USE_SSE2 - v4si emm0, emm2; -#else - v2si mm0, mm1, mm2, mm3; -#endif - sign_bit = x; - /* take the absolute value */ - x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask); - /* extract the sign bit (upper one) */ - sign_bit = _mm_and_ps(sign_bit, *(v4sf*)_ps_sign_mask); - - /* scale by 4/Pi */ - y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI); - -#ifdef USE_SSE2 - /* store the integer part of y in mm0 */ - emm2 = _mm_cvttps_epi32(y); - /* j=(j+1) & (~1) (see the cephes sources) */ - emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1); - emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1); - y = _mm_cvtepi32_ps(emm2); - - /* get the swap sign flag */ - emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4); - emm0 = _mm_slli_epi32(emm0, 29); - /* get the polynom selection mask - there is one polynom for 0 <= x <= Pi/4 - and another one for Pi/4