569 lines
20 KiB
C
Executable File
569 lines
20 KiB
C
Executable File
/* -----------------------------------------------------------------------------
|
|
*
|
|
* Copyright (c) 2008-2016 Alexis Naveros.
|
|
*
|
|
*
|
|
* The SIMD trigonometry functions are Copyright (C) 2007 Julien Pommier
|
|
* See copyright notice for simd4f_sin_ps(), simd4f_cos_ps(), simd4f_sincos_ps()
|
|
*
|
|
*
|
|
* Some functions are Copyright (C) 2008 José Fonseca
|
|
* See copyright notice for simd4f_exp2_ps(), simd4f_log2_ps(), simd4f_pow_ps()
|
|
*
|
|
*
|
|
* Portions developed under contract to the SURVICE Engineering Company.
|
|
*
|
|
* 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.
|
|
*
|
|
* -----------------------------------------------------------------------------
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <stddef.h>
|
|
#include <stdint.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <float.h>
|
|
#include <limits.h>
|
|
|
|
|
|
#include "cpusimd.h"
|
|
|
|
|
|
////
|
|
|
|
|
|
#if CPU_SSE_SUPPORT
|
|
|
|
const uint32_t CPU_ALIGN16 simd4fSignMask[4] = { 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
|
|
const uint32_t CPU_ALIGN16 simd4fSignMaskInv[4] = { 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff };
|
|
const float CPU_ALIGN16 simd4fHalf[4] = { 0.5, 0.5, 0.5, 0.5 };
|
|
const float CPU_ALIGN16 simd4fOne[4] = { 1.0, 1.0, 1.0, 1.0 };
|
|
const float CPU_ALIGN16 simd4fTwo[4] = { 2.0, 2.0, 2.0, 2.0 };
|
|
const float CPU_ALIGN16 simd4fThree[4] = { 3.0, 3.0, 3.0, 3.0 };
|
|
const uint32_t CPU_ALIGN16 simd4uOne[4] = { 1, 1, 1, 1 };
|
|
const uint32_t CPU_ALIGN16 simd4uOneInv[4] = { ~1, ~1, ~1, ~1 };
|
|
const uint32_t CPU_ALIGN16 simd4uTwo[4] = { 2, 2, 2, 2 };
|
|
const uint32_t CPU_ALIGN16 simd4uFour[4] = { 4, 4, 4, 4 };
|
|
const float CPU_ALIGN16 simd4fQuarter[4] = { 0.25, 0.25, 0.25, 0.25 };
|
|
// Avoid double to float conversion warnings when using M_PI --> use our own PI
|
|
#define PI 3.14159265358979323846f
|
|
const float CPU_ALIGN16 simd4fPi[4] = { PI, PI, PI, PI };
|
|
const float CPU_ALIGN16 simd4fZeroOneTwoThree[4] = { 0.0, 1.0, 2.0, 3.0 };
|
|
const uint32_t CPU_ALIGN16 simd4fAlphaMask[4] = { 0x00000000, 0x00000000, 0x00000000, 0xffffffff };
|
|
const float CPU_ALIGN16 simd4f255[4] = { 255.0f, 255.0f, 255.0f, 255.0f };
|
|
const float CPU_ALIGN16 simd4f255Inv[4] = { 1.0f/255.0f, 1.0f/255.0f, 1.0f/255.0f, 1.0f/255.0f };
|
|
|
|
#endif
|
|
|
|
|
|
////
|
|
|
|
|
|
#if CPU_SSE2_SUPPORT
|
|
|
|
|
|
/* Copyright (C) 2007 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)
|
|
*/
|
|
|
|
static const float CPU_ALIGN16 simd4f_cephes_FOPI[4] = { 1.27323954473516f, 1.27323954473516f, 1.27323954473516f, 1.27323954473516f };
|
|
static const float CPU_ALIGN16 simd4f_minus_cephes_DP1[4] = { -0.78515625, -0.78515625, -0.78515625, -0.78515625 };
|
|
static const float CPU_ALIGN16 simd4f_minus_cephes_DP2[4] = { -2.4187564849853515625e-4, -2.4187564849853515625e-4, -2.4187564849853515625e-4, -2.4187564849853515625e-4 };
|
|
static const float CPU_ALIGN16 simd4f_minus_cephes_DP3[4] = { -3.77489497744594108e-8f, -3.77489497744594108e-8f, -3.77489497744594108e-8f, -3.77489497744594108e-8f };
|
|
static const float CPU_ALIGN16 simd4f_sincof_p0[4] = { -1.9515295891E-4f, -1.9515295891E-4f, -1.9515295891E-4f, -1.9515295891E-4f };
|
|
static const float CPU_ALIGN16 simd4f_sincof_p1[4] = { 8.3321608736E-3f, 8.3321608736E-3f, 8.3321608736E-3f, 8.3321608736E-3f };
|
|
static const float CPU_ALIGN16 simd4f_sincof_p2[4] = { -1.6666654611E-1f, -1.6666654611E-1f, -1.6666654611E-1f, -1.6666654611E-1f };
|
|
static const float CPU_ALIGN16 simd4f_coscof_p0[4] = { 2.443315711809948E-005f, 2.443315711809948E-005f, 2.443315711809948E-005f, 2.443315711809948E-005f };
|
|
static const float CPU_ALIGN16 simd4f_coscof_p1[4] = { -1.388731625493765E-003f, -1.388731625493765E-003f, -1.388731625493765E-003f, -1.388731625493765E-003f };
|
|
static const float CPU_ALIGN16 simd4f_coscof_p2[4] = { 4.166664568298827E-002f, 4.166664568298827E-002f, 4.166664568298827E-002f, 4.166664568298827E-002f };
|
|
|
|
__m128 simd4f_sin_ps( __m128 x )
|
|
{
|
|
__m128 xmm1, xmm2, xmm3, sign_bit, y;
|
|
__m128i emm0, emm2;
|
|
|
|
xmm2 = _mm_setzero_ps();
|
|
|
|
sign_bit = x;
|
|
/* take the absolute value */
|
|
x = _mm_and_ps( x, *(__m128 *)simd4fSignMaskInv );
|
|
/* extract the sign bit (upper one) */
|
|
sign_bit = _mm_and_ps(sign_bit, *(__m128 *)simd4fSignMask);
|
|
|
|
/* scale by 4/Pi */
|
|
y = _mm_mul_ps(x, *(__m128 *)simd4f_cephes_FOPI);
|
|
|
|
/* 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, *(__m128i*)simd4uOne);
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uOneInv);
|
|
y = _mm_cvtepi32_ps(emm2);
|
|
|
|
/* get the swap sign flag */
|
|
emm0 = _mm_and_si128(emm2, *(__m128i*)simd4uFour);
|
|
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<x<=Pi/2
|
|
Both branches will be computed.
|
|
*/
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uTwo);
|
|
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
|
|
|
|
__m128 swap_sign_bit = _mm_castsi128_ps(emm0);
|
|
__m128 poly_mask = _mm_castsi128_ps(emm2);
|
|
sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
|
|
|
|
/* The magic pass: "Extended precision modular arithmetic"
|
|
x = ((x - y * DP1) - y * DP2) - y * DP3; */
|
|
xmm1 = *(__m128 *)simd4f_minus_cephes_DP1;
|
|
xmm2 = *(__m128 *)simd4f_minus_cephes_DP2;
|
|
xmm3 = *(__m128 *)simd4f_minus_cephes_DP3;
|
|
xmm1 = _mm_mul_ps(y, xmm1);
|
|
xmm2 = _mm_mul_ps(y, xmm2);
|
|
xmm3 = _mm_mul_ps(y, xmm3);
|
|
x = _mm_add_ps(x, xmm1);
|
|
x = _mm_add_ps(x, xmm2);
|
|
x = _mm_add_ps(x, xmm3);
|
|
|
|
/* Evaluate the first polynom (0 <= x <= Pi/4) */
|
|
y = *(__m128 *)simd4f_coscof_p0;
|
|
__m128 z = _mm_mul_ps(x,x);
|
|
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128 *)simd4f_coscof_p1);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128 *)simd4f_coscof_p2);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_mul_ps(y, z);
|
|
__m128 tmp = _mm_mul_ps(z, *(__m128 *)simd4fHalf);
|
|
y = _mm_sub_ps(y, tmp);
|
|
y = _mm_add_ps(y, *(__m128 *)simd4fOne);
|
|
|
|
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
|
|
|
|
__m128 y2 = *(__m128 *)simd4f_sincof_p0;
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128 *)simd4f_sincof_p1);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128 *)simd4f_sincof_p2);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_mul_ps(y2, x);
|
|
y2 = _mm_add_ps(y2, x);
|
|
|
|
/* select the correct result from the two polynoms */
|
|
xmm3 = poly_mask;
|
|
y2 = _mm_and_ps(xmm3, y2);
|
|
y = _mm_andnot_ps(xmm3, y);
|
|
y = _mm_add_ps(y,y2);
|
|
/* update the sign */
|
|
y = _mm_xor_ps(y, sign_bit);
|
|
|
|
return y;
|
|
}
|
|
|
|
|
|
|
|
/* almost the same as sin_ps */
|
|
__m128 simd4f_cos_ps( __m128 x )
|
|
{
|
|
__m128 xmm1, xmm2, xmm3, y;
|
|
__m128i emm0, emm2;
|
|
|
|
xmm2 = _mm_setzero_ps();
|
|
|
|
/* take the absolute value */
|
|
x = _mm_and_ps(x, *(__m128*)simd4fSignMaskInv);
|
|
|
|
/* scale by 4/Pi */
|
|
y = _mm_mul_ps(x, *(__m128*)simd4f_cephes_FOPI);
|
|
|
|
/* 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, *(__m128i*)simd4uOne);
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uOneInv);
|
|
y = _mm_cvtepi32_ps(emm2);
|
|
|
|
emm2 = _mm_sub_epi32(emm2, *(__m128i*)simd4uTwo);
|
|
|
|
/* get the swap sign flag */
|
|
emm0 = _mm_andnot_si128(emm2, *(__m128i*)simd4uFour);
|
|
emm0 = _mm_slli_epi32(emm0, 29);
|
|
/* get the polynom selection mask */
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uTwo);
|
|
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
|
|
|
|
__m128 sign_bit = _mm_castsi128_ps(emm0);
|
|
__m128 poly_mask = _mm_castsi128_ps(emm2);
|
|
/* The magic pass: "Extended precision modular arithmetic"
|
|
x = ((x - y * DP1) - y * DP2) - y * DP3; */
|
|
xmm1 = *(__m128*)simd4f_minus_cephes_DP1;
|
|
xmm2 = *(__m128*)simd4f_minus_cephes_DP2;
|
|
xmm3 = *(__m128*)simd4f_minus_cephes_DP3;
|
|
xmm1 = _mm_mul_ps(y, xmm1);
|
|
xmm2 = _mm_mul_ps(y, xmm2);
|
|
xmm3 = _mm_mul_ps(y, xmm3);
|
|
x = _mm_add_ps(x, xmm1);
|
|
x = _mm_add_ps(x, xmm2);
|
|
x = _mm_add_ps(x, xmm3);
|
|
|
|
/* Evaluate the first polynom (0 <= x <= Pi/4) */
|
|
y = *(__m128*)simd4f_coscof_p0;
|
|
__m128 z = _mm_mul_ps(x,x);
|
|
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128*)simd4f_coscof_p1);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128*)simd4f_coscof_p2);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_mul_ps(y, z);
|
|
__m128 tmp = _mm_mul_ps(z, *(__m128*)simd4fHalf);
|
|
y = _mm_sub_ps(y, tmp);
|
|
y = _mm_add_ps(y, *(__m128*)simd4fOne);
|
|
|
|
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
|
|
|
|
__m128 y2 = *(__m128*)simd4f_sincof_p0;
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128*)simd4f_sincof_p1);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128*)simd4f_sincof_p2);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_mul_ps(y2, x);
|
|
y2 = _mm_add_ps(y2, x);
|
|
|
|
/* select the correct result from the two polynoms */
|
|
xmm3 = poly_mask;
|
|
y2 = _mm_and_ps(xmm3, y2); //, xmm3);
|
|
y = _mm_andnot_ps(xmm3, y);
|
|
y = _mm_add_ps(y,y2);
|
|
/* update the sign */
|
|
y = _mm_xor_ps(y, sign_bit);
|
|
|
|
return y;
|
|
}
|
|
|
|
/* since sin_ps and cos_ps are almost identical, sincos_ps could replace both of them..
|
|
it is almost as fast, and gives you a free cosine with your sine */
|
|
void simd4f_sincos_ps( __m128 x, __m128 *s, __m128 *c )
|
|
{
|
|
__m128 xmm1, xmm2, xmm3, sign_bit_sin, y;
|
|
__m128i emm0, emm2, emm4;
|
|
|
|
xmm3 = _mm_setzero_ps();
|
|
|
|
sign_bit_sin = x;
|
|
/* take the absolute value */
|
|
x = _mm_and_ps(x, *(__m128*)simd4fSignMaskInv);
|
|
/* extract the sign bit (upper one) */
|
|
sign_bit_sin = _mm_and_ps(sign_bit_sin, *(__m128*)simd4fSignMask);
|
|
|
|
/* scale by 4/Pi */
|
|
y = _mm_mul_ps(x, *(__m128*)simd4f_cephes_FOPI);
|
|
|
|
/* store the integer part of y in emm2 */
|
|
emm2 = _mm_cvttps_epi32(y);
|
|
|
|
/* j=(j+1) & (~1) (see the cephes sources) */
|
|
emm2 = _mm_add_epi32(emm2, *(__m128i*)simd4uOne);
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uOneInv);
|
|
y = _mm_cvtepi32_ps(emm2);
|
|
|
|
emm4 = emm2;
|
|
|
|
/* get the swap sign flag for the sine */
|
|
emm0 = _mm_and_si128(emm2, *(__m128i*)simd4uFour);
|
|
emm0 = _mm_slli_epi32(emm0, 29);
|
|
__m128 swap_sign_bit_sin = _mm_castsi128_ps(emm0);
|
|
|
|
/* get the polynom selection mask for the sine*/
|
|
emm2 = _mm_and_si128(emm2, *(__m128i*)simd4uTwo);
|
|
emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
|
|
__m128 poly_mask = _mm_castsi128_ps(emm2);
|
|
|
|
/* The magic pass: "Extended precision modular arithmetic"
|
|
x = ((x - y * DP1) - y * DP2) - y * DP3; */
|
|
xmm1 = *(__m128*)simd4f_minus_cephes_DP1;
|
|
xmm2 = *(__m128*)simd4f_minus_cephes_DP2;
|
|
xmm3 = *(__m128*)simd4f_minus_cephes_DP3;
|
|
xmm1 = _mm_mul_ps(y, xmm1);
|
|
xmm2 = _mm_mul_ps(y, xmm2);
|
|
xmm3 = _mm_mul_ps(y, xmm3);
|
|
x = _mm_add_ps(x, xmm1);
|
|
x = _mm_add_ps(x, xmm2);
|
|
x = _mm_add_ps(x, xmm3);
|
|
|
|
emm4 = _mm_sub_epi32(emm4, *(__m128i*)simd4uTwo);
|
|
emm4 = _mm_andnot_si128(emm4, *(__m128i*)simd4uFour);
|
|
emm4 = _mm_slli_epi32(emm4, 29);
|
|
__m128 sign_bit_cos = _mm_castsi128_ps(emm4);
|
|
|
|
sign_bit_sin = _mm_xor_ps(sign_bit_sin, swap_sign_bit_sin);
|
|
|
|
/* Evaluate the first polynom (0 <= x <= Pi/4) */
|
|
__m128 z = _mm_mul_ps(x,x);
|
|
y = *(__m128*)simd4f_coscof_p0;
|
|
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128*)simd4f_coscof_p1);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_add_ps(y, *(__m128*)simd4f_coscof_p2);
|
|
y = _mm_mul_ps(y, z);
|
|
y = _mm_mul_ps(y, z);
|
|
__m128 tmp = _mm_mul_ps(z, *(__m128*)simd4fHalf);
|
|
y = _mm_sub_ps(y, tmp);
|
|
y = _mm_add_ps(y, *(__m128*)simd4fOne);
|
|
|
|
/* Evaluate the second polynom (Pi/4 <= x <= 0) */
|
|
|
|
__m128 y2 = *(__m128*)simd4f_sincof_p0;
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128*)simd4f_sincof_p1);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_add_ps(y2, *(__m128*)simd4f_sincof_p2);
|
|
y2 = _mm_mul_ps(y2, z);
|
|
y2 = _mm_mul_ps(y2, x);
|
|
y2 = _mm_add_ps(y2, x);
|
|
|
|
/* select the correct result from the two polynoms */
|
|
xmm3 = poly_mask;
|
|
__m128 ysin2 = _mm_and_ps(xmm3, y2);
|
|
__m128 ysin1 = _mm_andnot_ps(xmm3, y);
|
|
y2 = _mm_sub_ps(y2,ysin2);
|
|
y = _mm_sub_ps(y, ysin1);
|
|
|
|
xmm1 = _mm_add_ps(ysin1,ysin2);
|
|
xmm2 = _mm_add_ps(y,y2);
|
|
|
|
/* update the sign */
|
|
*s = _mm_xor_ps(xmm1, sign_bit_sin);
|
|
*c = _mm_xor_ps(xmm2, sign_bit_cos);
|
|
}
|
|
|
|
|
|
#endif
|
|
|
|
|
|
////
|
|
|
|
|
|
#if CPU_SSE2_SUPPORT
|
|
|
|
|
|
/* Copyright (C) 2008 José Fonseca
|
|
http://jrfonseca.blogspot.ca/2008/09/fast-sse2-pow-tables-or-polynomials.html
|
|
MIT license
|
|
|
|
Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
of this software and associated documentation files (the "Software"), to deal
|
|
in the Software without restriction, including without limitation the rights
|
|
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
copies of the Software, and to permit persons to whom the Software is
|
|
furnished to do so, subject to the following conditions:
|
|
|
|
The above copyright notice and this permission notice shall be included in all
|
|
copies or substantial portions of the Software.
|
|
|
|
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
SOFTWARE.
|
|
*/
|
|
|
|
#define POLY0(x,c0) _mm_set1_ps(c0)
|
|
#define POLY1(x,c0,c1) _mm_add_ps(_mm_mul_ps(POLY0(x, c1), x), _mm_set1_ps(c0))
|
|
#define POLY2(x,c0,c1,c2) _mm_add_ps(_mm_mul_ps(POLY1(x, c1, c2), x), _mm_set1_ps(c0))
|
|
#define POLY3(x,c0,c1,c2,c3) _mm_add_ps(_mm_mul_ps(POLY2(x, c1, c2, c3), x), _mm_set1_ps(c0))
|
|
#define POLY4(x,c0,c1,c2,c3,c4) _mm_add_ps(_mm_mul_ps(POLY3(x, c1, c2, c3, c4), x), _mm_set1_ps(c0))
|
|
#define POLY5(x,c0,c1,c2,c3,c4,c5) _mm_add_ps(_mm_mul_ps(POLY4(x, c1, c2, c3, c4, c5), x), _mm_set1_ps(c0))
|
|
|
|
#define EXP_POLY_DEGREE 3
|
|
#define LOG_POLY_DEGREE 5
|
|
|
|
__m128 simd4f_exp2_ps( __m128 x )
|
|
{
|
|
__m128i ipart;
|
|
__m128 fpart, expipart, expfpart;
|
|
|
|
x = _mm_min_ps( x, _mm_set1_ps( 129.00000f ) );
|
|
x = _mm_max_ps( x, _mm_set1_ps( -126.99999f ) );
|
|
/* ipart = int(x - 0.5) */
|
|
ipart = _mm_cvtps_epi32( _mm_sub_ps( x, _mm_set1_ps( 0.5f ) ) );
|
|
/* fpart = x - ipart */
|
|
fpart = _mm_sub_ps( x, _mm_cvtepi32_ps( ipart ) );
|
|
/* expipart = (float) (1 << ipart) */
|
|
expipart = _mm_castsi128_ps( _mm_slli_epi32( _mm_add_epi32( ipart, _mm_set1_epi32( 127 ) ), 23 ) );
|
|
/* minimax polynomial fit of 2**x, in range [-0.5, 0.5[ */
|
|
#if EXP_POLY_DEGREE == 5
|
|
expfpart = POLY5( fpart, 9.9999994e-1f, 6.9315308e-1f, 2.4015361e-1f, 5.5826318e-2f, 8.9893397e-3f, 1.8775767e-3f );
|
|
#elif EXP_POLY_DEGREE == 4
|
|
expfpart = POLY4( fpart, 1.0000026f, 6.9300383e-1f, 2.4144275e-1f, 5.2011464e-2f, 1.3534167e-2f );
|
|
#elif EXP_POLY_DEGREE == 3
|
|
expfpart = POLY3( fpart, 9.9992520e-1f, 6.9583356e-1f, 2.2606716e-1f, 7.8024521e-2f );
|
|
#elif EXP_POLY_DEGREE == 2
|
|
expfpart = POLY2( fpart, 1.0017247f, 6.5763628e-1f, 3.3718944e-1f );
|
|
#else
|
|
#error
|
|
#endif
|
|
return _mm_mul_ps(expipart, expfpart);
|
|
}
|
|
|
|
__m128 simd4f_log2_ps( __m128 x )
|
|
{
|
|
__m128i expmask, mantmask, i;
|
|
__m128 one, vexp, mant, logmant;
|
|
|
|
expmask = _mm_set1_epi32( 0x7f800000 );
|
|
mantmask = _mm_set1_epi32( 0x007fffff );
|
|
one = _mm_set1_ps( 1.0f );
|
|
i = _mm_castps_si128( x );
|
|
/* exp = (float) exponent(x) */
|
|
vexp = _mm_cvtepi32_ps( _mm_sub_epi32( _mm_srli_epi32( _mm_and_si128( i, expmask ), 23 ), _mm_set1_epi32( 127 ) ) );
|
|
/* mant = (float) mantissa(x) */
|
|
mant = _mm_or_ps( _mm_castsi128_ps( _mm_and_si128( i, mantmask ) ), one );
|
|
/* Minimax polynomial fit of log2(x)/(x - 1), for x in range [1, 2[
|
|
* These coefficients can be generate with
|
|
* http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals2/minimax.html
|
|
*/
|
|
#if LOG_POLY_DEGREE == 6
|
|
logmant = POLY5( mant, 3.11578814719469302614f, -3.32419399085241980044f, 2.59883907202499966007f, -1.23152682416275988241f, 0.318212422185251071475f, -0.0344359067839062357313f );
|
|
#elif LOG_POLY_DEGREE == 5
|
|
logmant = POLY4( mant, 2.8882704548164776201f, -2.52074962577807006663f, 1.48116647521213171641f, -0.465725644288844778798f, 0.0596515482674574969533f );
|
|
#elif LOG_POLY_DEGREE == 4
|
|
logmant = POLY3( mant, 2.61761038894603480148f, -1.75647175389045657003f, 0.688243882994381274313f, -0.107254423828329604454f );
|
|
#elif LOG_POLY_DEGREE == 3
|
|
logmant = POLY2( mant, 2.28330284476918490682f, -1.04913055217340124191f, 0.204446009836232697516f );
|
|
#else
|
|
#error
|
|
#endif
|
|
/* This effectively increases the polynomial degree by one, but ensures that log2(1) == 0*/
|
|
logmant = _mm_mul_ps( logmant, _mm_sub_ps(mant, one ) );
|
|
return _mm_add_ps( logmant, vexp );
|
|
}
|
|
|
|
|
|
__m128 simd4f_pow_ps( __m128 x, __m128 y )
|
|
{
|
|
return simd4f_exp2_ps( _mm_mul_ps( simd4f_log2_ps( x ), y ) );
|
|
}
|
|
|
|
|
|
#endif
|
|
|
|
|
|
////
|
|
|
|
|
|
#if CPU_SSE2_SUPPORT
|
|
|
|
|
|
/*
|
|
By Potatoswatter
|
|
http://stackoverflow.com/questions/6475373/optimizations-for-pow-with-const-non-integer-exponent
|
|
*/
|
|
|
|
#ifndef CC_ALWAYSINLINE
|
|
#if defined(__GNUC__) || defined(__INTEL_COMPILER)
|
|
#define CC_ALWAYSINLINE __attribute__((always_inline))
|
|
#else
|
|
#define CC_ALWAYSINLINE
|
|
#endif
|
|
#endif
|
|
|
|
static inline CC_ALWAYSINLINE __m128 simd4f_fastpow_ps( __m128 arg, uint32_t expnum, uint32_t expden, uint32_t coeffnum, uint32_t coeffden )
|
|
{
|
|
__m128 ret = arg;
|
|
float corrfactor, powfactor;
|
|
/* Apply a constant pre-correction factor. */
|
|
corrfactor = (float)(exp2( 127.0 * expden / expnum - 127.0 ) * pow( 1.0 * coeffnum / coeffden, 1.0 * expden / expnum ));
|
|
powfactor = 1.0f * expnum / expden;
|
|
ret = _mm_mul_ps( ret, _mm_set1_ps( corrfactor ) );
|
|
/* Reinterpret arg as integer to obtain logarithm. */
|
|
ret = _mm_cvtepi32_ps( _mm_castps_si128( ret ) );
|
|
/* Multiply logarithm by power. */
|
|
ret = _mm_mul_ps( ret, _mm_set1_ps( powfactor ) );
|
|
/* Convert back to "integer" to exponentiate. */
|
|
ret = _mm_castsi128_ps( _mm_cvtps_epi32( ret ) );
|
|
return ret;
|
|
}
|
|
|
|
__m128 simd4f_pow12d5_ps( __m128 arg )
|
|
{
|
|
/* Lower exponents provide lower initial error, but too low causes overflow. */
|
|
__m128 xf = simd4f_fastpow_ps( arg, 4, 5, (int)( 1.38316186f * (float)1e9 ), (int)1e9 );
|
|
/* Imprecise 4-cycle sqrt is still far better than fastpow, good enough. */
|
|
__m128 xfm4 = _mm_rsqrt_ps( xf );
|
|
__m128 xf4 = _mm_mul_ps( xf, xfm4 );
|
|
/* Precisely calculate x^2 and x^3 */
|
|
__m128 x2 = _mm_mul_ps( arg, arg );
|
|
__m128 x3 = _mm_mul_ps( x2, arg );
|
|
/* Overestimate of x^2 * x^0.4 */
|
|
x2 = _mm_mul_ps( x2, xf4 );
|
|
/* Get x^-0.2 from x^0.4, and square it for x^-0.4. Combine into x^-0.6. */
|
|
__m128 xfm2 = _mm_rsqrt_ps( xf4 );
|
|
x3 = _mm_mul_ps( x3, xfm4 );
|
|
x3 = _mm_mul_ps( x3, xfm2 );
|
|
return _mm_mul_ps( _mm_add_ps( x2, x3 ), _mm_set1_ps( 1.0f/1.960131704207789f * 0.9999f ) );
|
|
}
|
|
|
|
__m128 simd4f_pow5d12_ps( __m128 arg )
|
|
{
|
|
/* 5/12 is too small, so compute the 4th root of 20/12 instead. */
|
|
/* 20/12 = 5/3 = 1 + 2/3 = 2 - 1/3. 2/3 is a suitable argument for fastpow. */
|
|
/* weighting coefficient: a^-1/2 = 2 a; a = 2^-2/3 */
|
|
__m128 xf = simd4f_fastpow_ps( arg, 2, 3, (int)( 0.629960524947437f * (float)1e9 ), (int)1e9 );
|
|
__m128 xover = _mm_mul_ps( arg, xf );
|
|
__m128 xfm1 = _mm_rsqrt_ps( xf );
|
|
__m128 x2 = _mm_mul_ps( arg, arg );
|
|
__m128 xunder = _mm_mul_ps( x2, xfm1 );
|
|
/* sqrt2 * over + 2 * sqrt2 * under */
|
|
__m128 xavg = _mm_mul_ps( _mm_set1_ps( 1.0f/( 3.0f * 0.629960524947437f ) * 0.999852f ), _mm_add_ps( xover, xunder ) );
|
|
xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
|
|
xavg = _mm_mul_ps( xavg, _mm_rsqrt_ps( xavg ) );
|
|
return xavg;
|
|
}
|
|
|
|
#endif
|
|
|
|
|
|
////
|
|
|