#include "xmmLibm_prefix.h"
#include "math.h"
#include "fenv.h"
static inline xDouble _xlog( xDouble x ) ALWAYS_INLINE;
static const double A1[2] = {
0x1.5555555550286p-4, 0x1.999A0BC712416p-7 };
static const double A2[4] = {
0x1.55555555554E6p-4, 0x1.9999999BAC6D4p-7, 0x1.2492307F1519Fp-9, 0x1.C8034C85DFFF0p-12 };
static const double C_lead[129] = {
0x0, 0x1.FE02A6B100000p-8, 0x1.FC0A8B0FC0000p-7, 0x1.7B91B07D60000p-6, 0x1.F829B0E780000p-6, 0x1.39E87B9FE8000p-5, 0x1.77458F6330000p-5, 0x1.B42DD71198000p-5, 0x1.F0A30C0118000p-5, 0x1.16536EEA38000p-4, 0x1.341D7961BC000p-4, 0x1.51B073F060000p-4, 0x1.6F0D28AE58000p-4, 0x1.8C345D6318000p-4, 0x1.A926D3A4AC000p-4, 0x1.C5E548F5BC000p-4, 0x1.E27076E2B0000p-4, 0x1.FEC9131DC0000p-4, 0x1.0D77E7CD08000p-3, 0x1.1B72AD52F6000p-3, 0x1.29552F8200000p-3, 0x1.371FC201E8000p-3, 0x1.44D2B6CCB8000p-3, 0x1.526E5E3A1C000p-3, 0x1.5FF3070A7A000p-3, 0x1.6D60FE719E000p-3, 0x1.7AB890210E000p-3, 0x1.87FA06520C000p-3, 0x1.9525A9CF46000p-3, 0x1.A23BC1FE2C000p-3, 0x1.AF3C94E80C000p-3, 0x1.BC286742D8000p-3, 0x1.C8FF7C79AA000p-3, 0x1.D5C216B4FC000p-3, 0x1.E27076E2B0000p-3, 0x1.EF0ADCBDC6000p-3, 0x1.FB9186D5E4000p-3, 0x1.0402594B4D000p-2, 0x1.0A324E2739000p-2, 0x1.1058BF9AE5000p-2, 0x1.1675CABABA000p-2, 0x1.1C898C169A000p-2, 0x1.22941FBCF8000p-2, 0x1.2895A13DE8000p-2, 0x1.2E8E2BAE12000p-2, 0x1.347DD9A988000p-2, 0x1.3A64C55694000p-2, 0x1.404308686A000p-2, 0x1.4618BC21C6000p-2, 0x1.4BE5F95778000p-2, 0x1.51AAD872E0000p-2, 0x1.5767717456000p-2, 0x1.5D1BDBF581000p-2, 0x1.62C82F2B9C000p-2, 0x1.686C81E9B1000p-2, 0x1.6E08EAA2BA000p-2, 0x1.739D7F6BBD000p-2, 0x1.792A55FDD4000p-2, 0x1.7EAF83B82B000p-2, 0x1.842D1DA1E9000p-2, 0x1.89A3386C14000p-2, 0x1.8F11E87366000p-2, 0x1.947941C211000p-2, 0x1.99D958117E000p-2, 0x1.9F323ECBFA000p-2, 0x1.A484090E5C000p-2, 0x1.A9CEC9A9A1000p-2, 0x1.AF12932478000p-2, 0x1.B44F77BCC9000p-2, 0x1.B985896931000p-2, 0x1.BEB4D9DA72000p-2, 0x1.C3DD7A7CDB000p-2, 0x1.C8FF7C79AA000p-2, 0x1.CE1AF0B85F000p-2, 0x1.D32FE7E00F000p-2, 0x1.D83E7258A3000p-2, 0x1.DD46A04C1C000p-2, 0x1.E24881A7C7000p-2, 0x1.E744261D68000p-2, 0x1.EC399D2469000p-2, 0x1.F128F5FAF0000p-2, 0x1.F6123FA703000p-2, 0x1.FAF588F78F000p-2, 0x1.FFD2E0857F000p-2, 0x1.02552A5A5D000p-1, 0x1.04BDF9DA92800p-1, 0x1.0723E5C1CE000p-1, 0x1.0986F4F573800p-1, 0x1.0BE72E4252800p-1, 0x1.0E44985D1D000p-1, 0x1.109F39E2D5000p-1, 0x1.12F719593F000p-1, 0x1.154C3D2F4D800p-1, 0x1.179EABBD89800p-1, 0x1.19EE6B467C800p-1, 0x1.1C3B81F714000p-1, 0x1.1E85F5E704000p-1, 0x1.20CDCD192A800p-1, 0x1.23130D7BEC000p-1, 0x1.2555BCE98F800p-1, 0x1.2795E1289B000p-1, 0x1.29D37FEC2B000p-1, 0x1.2C0E9ED449000p-1, 0x1.2E47436E40000p-1, 0x1.307D7334F1000p-1, 0x1.32B1339122000p-1, 0x1.34E289D9CE000p-1, 0x1.37117B5474800p-1, 0x1.393E0D3562800p-1, 0x1.3B6844A000000p-1, 0x1.3D9026A715800p-1, 0x1.3FB5B84D17000p-1, 0x1.41D8FE8467000p-1, 0x1.43F9FE2F9D000p-1, 0x1.4618BC21C6000p-1, 0x1.48353D1EA8800p-1, 0x1.4A4F85DB04000p-1, 0x1.4C679AFCCF000p-1, 0x1.4E7D811B75800p-1, 0x1.50913CC016800p-1, 0x1.52A2D265BC800p-1, 0x1.54B2467999800p-1, 0x1.56BF9D5B3F000p-1, 0x1.58CADB5CD7800p-1, 0x1.5AD404C35A000p-1, 0x1.5CDB1DC6C1800p-1, 0x1.5EE02A9241800p-1, 0x1.60E32F4478800p-1, 0x1.62E42FEFA3800p-1, };
static const double C_trail[129] = {
0.0, 0x1.9E23F0DDA40E4p-46, 0x1.F1E7CF6D3A69Cp-50, -0x1.3B955B602ACE4p-44, 0x1.980267C7E09E4p-45, 0x1.EAFD480AD9015p-44, -0x1.181DCE586AF09p-44, -0x1.C827AE5D6704Cp-46, -0x1.D599E83368E91p-45, -0x1.47C5E768FA309p-46, 0x1.1D09299837610p-44, 0x1.83F69278E686Ap-44, -0x1.4B4641B664613p-44, 0x1.B20F5ACB42A66p-44, 0x1.563650BD22A9Cp-44, 0x1.D0C57585FBE06p-46, -0x1.A342C2AF0003Cp-45, -0x1.54555D1AE6607p-44, 0x1.CB2CD2EE2F482p-44, 0x1.E80A41811A396p-45, -0x1.5B967F4471DFCp-44, 0x1.EE8779B2D8ABCp-44, -0x1.70CC16135783Cp-46, -0x1.790BA37FC5238p-44, -0x1.8586F183BEBF2p-44, -0x1.BC6E557134767p-44, -0x1.BDB9072534A58p-45, 0x1.22120401202FCp-44, -0x1.297137D9F158Fp-44, -0x1.539CD91DC9F0Bp-44, -0x1.A4E633FCD9066p-52, 0x1.9AC53F39D121Cp-44, -0x1.7794F689F8434p-45, -0x1.1BA91BBCA681Bp-45, -0x1.A342C2AF0003Cp-44, -0x1.B26B79C86AF24p-45, -0x1.D572AAB993C87p-47, 0x1.036B89EF42D7Fp-48, 0x1.C6BEE7EF4030Ep-47, -0x1.4AB9D817D52CDp-44, 0x1.8380E731F55C4p-44, -0x1.81410E5C62AFFp-44, -0x1.A6976F5EB0963p-44, 0x1.A8D7AD24C13F0p-44, -0x1.67B1E99B72BD8p-45, -0x1.5594DD4C58092p-45, 0x1.7A71CBCD735D0p-44, 0x1.F8EF43049F7D3p-44, -0x1.3D82F484C84CCp-46, -0x1.D7C92CD9AD824p-44, -0x1.F4BD8DB0A7CC1p-44, -0x1.64EAD9524D7CAp-44, -0x1.8D6BDC9C7C238p-44, 0x1.E54BDBD7C8A98p-44, 0x1.2BB110AF84054p-44, 0x1.E38C139318D71p-46, 0x1.A7389314FEB50p-52, 0x1.E89F057691FEAp-44, -0x1.E4DA62D0C25ADp-49, -0x1.3A2DB13AE687Cp-44, 0x1.2D5AD38C40882p-45, 0x1.63BF0BB4EAB4Cp-45, 0x1.BEAE9337451F4p-44, 0x1.1597525DD88F0p-47, -0x1.ED03525CA2643p-44, -0x1.3D7500D6523C5p-44, -0x1.ED9CADEC02B43p-44, -0x1.E53BB31EED7A9p-44, -0x1.3AE68224AA2CEp-47, 0x1.F6B31F629F11Ep-47, -0x1.21021E78B2151p-44, -0x1.5946261F5A42Bp-45, -0x1.7794F689F8434p-44, 0x1.F5BDBE95E5568p-45, -0x1.0AA7884DCD050p-44, -0x1.835F5D48BA26Dp-47, 0x1.282FB989A9274p-44, -0x1.ECF1A1385D356p-45, 0x1.E1F8DF68DBCF3p-44, -0x1.9FF45188D6065p-45, 0x1.BB2CD720EC44Cp-44, -0x1.D4E7AEA4F0D25p-44, 0x1.8F6CD7D9F2754p-45, 0x1.261565F40D932p-44, 0x1.FD8D38D2BAFDDp-46, -0x1.2D9A033EFF74Ep-45, -0x1.7F6350D38EDDDp-46, -0x1.6FA37012B5806p-44, 0x1.415B4C4BDD99Fp-44, -0x1.BA048A8D10B4Bp-44, -0x1.B4810E09B27A4p-44, -0x1.0EB3FB7398E0Cp-47, -0x1.0B2B38662E34Dp-44, 0x1.A0BFC60E6FA08p-45, 0x1.6ECC5CBDD7782p-45, -0x1.EDA1B58389902p-44, 0x1.A07BD8B34BE7Cp-46, 0x1.B6C9A81E87BAEp-44, -0x1.7AFA4392F1BA7p-46, -0x1.A61FDE292977Ep-48, 0x1.1AEB783F3DB97p-45, 0x1.1590B9AD974BAp-46, -0x1.74468563CE45Dp-45, 0x1.34202A10C3491p-44, 0x1.7C3F6B2143EADp-46, -0x1.4766FD54A4C27p-44, 0x1.D316EB92D885Dp-45, -0x1.28E88BF6DEEC9p-47, 0x1.0CD4E221301B7p-44, -0x1.EEA838909F3D3p-44, -0x1.055BFBD9C2F53p-45, -0x1.7B4962C55F46Bp-46, 0x1.5732325E617A3p-44, -0x1.98858D84649F1p-45, -0x1.3D82F484C84CCp-45, 0x1.BEE7ABD176604p-46, -0x1.44FDD840B8591p-45, -0x1.C64E971322CE8p-45, 0x1.D84E584C2B22Cp-44, 0x1.AD2F2CE96C2D6p-47, -0x1.2A88C41BA8752p-44, -0x1.B42B755EBA5E1p-44, 0x1.CCA08E310B9B2p-44, 0x1.893092F25D931p-45, -0x1.A609ACAAB41FCp-46, -0x1.36E612387451Fp-46, -0x1.8A8F29F6A02DCp-45, 0x1.B194F912B416Ap-46, 0x1.EF35793C76730p-45, };
static const double T1 = 0x1.E0FABFBC702A3p-1; static const double T2 = 0x1.1082B577D34EEp+0; static const double tiny = 0x0.0000000000001p-1022; static const xDouble plusinf = { INFINITY, 0.0 };
static const xUInt64 smallestNormal = { 0x0010000000000000ULL, 0 };
static const double denormBias = 0.0;
static const xDouble one = { 1.0, 0.0 };
static const double d128 = 128.0;
static const double d0_0078125 = 0.0078125;
static const xSInt64 bias[2] = { {2045, 0}, {1023, 0} };
static const xSInt64 *biasp = bias + 1;
static const double mOne = -1.0;
static const double two = 2.0;
static const xSInt64 logNaN = { 0x7ff8000000000033LL, 0 };
static inline xDouble mantissa( xDouble x, xDouble *exponent ) ALWAYS_INLINE;
static inline xDouble mantissa( xDouble x, xDouble *exponent )
{
xDouble mantissa = _mm_andnot_pd( plusinf, x ); xDouble isDenormal = _mm_cmplt_sdm( x, (double*) &smallestNormal ); xDouble denormExp = _mm_and_pd( one, isDenormal ); int biasSelector = _mm_cvtsi128_si32( (__m128i) isDenormal ); x = _mm_or_pd( x, denormExp ); x = _mm_sub_sd( x, denormExp ); mantissa = _mm_andnot_pd( plusinf, x ); xUInt64 iExp = (__m128i) _mm_and_pd( x, plusinf ); iExp = _mm_srli_epi64( iExp, 52 ); iExp = _mm_sub_epi64( iExp, biasp[ biasSelector ] ); *exponent = _mm_cvtepi32_pd( iExp ); return _mm_or_pd( mantissa, one ); }
static inline xDouble _xlog( xDouble x )
{
static const double half = 0.5;
xDouble xIsNaN = _mm_cmpunord_sd( x, x );
xDouble safeX = _mm_andnot_pd( xIsNaN, x );
xDouble xGTzero = _mm_cmpgt_sdm( safeX, (double*) &minusZeroD );
xDouble xLTinf = _mm_cmplt_sdm( safeX, (double*) &plusinf );
xDouble Y, F, f, f1, L_lead, L_trail, R, u, u1, u2, v, q, g, m;
int j;
if( _mm_istrue_sd( _mm_and_pd( xGTzero, xLTinf ) ) ) {
xDouble xLET1 = _mm_cmple_sdm( x, &T1 );
xDouble xGET2 = _mm_cmpge_sdm( x, &T2 );
if( _mm_istrue_sd( _mm_or_pd( xLET1, xGET2 ) ) )
{
Y = mantissa( x, &m );
j = _mm_cvttsd_si32( _mm_add_sdm( _mm_mul_sdm( Y, &d128 ), &half) );
F = _mm_mul_sdm( _mm_cvtsi32_sd( Y, j) , &d0_0078125 );
j -= 128;
f = _mm_sub_sd( Y, F );
L_lead = _mm_add_sdm( _mm_mul_sdm( m, &C_lead[128] ), &C_lead[j] );
L_trail = _mm_add_sdm( _mm_mul_sdm( m, &C_trail[128] ), &C_trail[j] );
u = _mm_div_sd( _mm_add_sd( f, f), _mm_add_sd( Y, F ) );
v = _mm_mul_sd( u, u );
q = _mm_add_sdm( _mm_mul_sdm( v, &A1[1] ), &A1[0] );
q = _mm_mul_sd( q, v );
q = _mm_mul_sd( q, u );
R = _mm_add_sd( L_lead, _mm_add_sd( u , _mm_add_sd( q, L_trail ) ) );
}
else
{
xDouble xEQone = _mm_cmpeq_sdm( x, (double*) &one );
f = _mm_sub_sdm( x, (double*) &one );
g = _mm_div_sd( one, _mm_add_sdm( f, &two ) );
u = _mm_mul_sd( f, g ); u = _mm_add_sd( u, u ); v = _mm_mul_sd( u, u );
q = _mm_add_sdm( _mm_mul_sdm( v, &A2[3] ), &A2[2] );
q = _mm_add_sdm( _mm_mul_sd( q, v ), &A2[1] );
q = _mm_add_sdm( _mm_mul_sd( q, v ), &A2[0] );
q = _mm_mul_sd( _mm_mul_sd( q, v ), u );
u1 = _mm_cvtss_sd( u, _mm_cvtsd_ss( (xFloat) u, u ) );
f1 = _mm_cvtss_sd( f, _mm_cvtsd_ss( (xFloat) f, f ) );
xDouble t0 = _mm_sub_sd( f, u1 ); xDouble t1 = _mm_sub_sd( f, f1 ); xDouble t2 = _mm_mul_sd( u1, f1 ); t0 = _mm_add_sd( t0, t0 ); t1 = _mm_mul_sd( t1, u1 ); t0 = _mm_sub_sd( t0, t2 ); t0 = _mm_sub_sd( t0, t1 ); u2 = _mm_mul_sd( t0, g ); R = _mm_add_sd( u2, q );
R = _mm_add_sd( R, u1 );
R = _mm_andnot_pd( xEQone, R );
}
}
else
{ if( _mm_istrue_sd( _mm_cmpeq_sd( x, minusZeroD ) ) )
R = _mm_div_sd( _mm_load_sd( &two ), minusZeroD );
else
{
xDouble inf = plusinf;
xDouble isInf = _mm_cmpeq_sd( x, inf );
inf = _mm_andnot_pd( xIsNaN, inf );
inf = _mm_andnot_pd( isInf, inf );
xDouble xLTzero = _mm_cmplt_sdm( safeX, (double*) &minusZeroD );
R = _mm_sub_sd( x, inf ); R = _mm_add_sd( R, inf );
xLTzero = _mm_and_pd( xLTzero, (xDouble) logNaN );
R = _mm_add_sd( xLTzero, R ); }
}
return R;
}
static inline double _xlog2( double x ) ALWAYS_INLINE;
static inline double _xlog2( double x )
{
static const double conversion = 1.0 / 6.9314718055994530942E-1;
xDouble r = _xlog( DOUBLE_2_XDOUBLE( x ) );
r = _mm_mul_sdm( r, &conversion );
return XDOUBLE_2_DOUBLE( r );
}
#pragma mark -
#pragma mark log2
#if defined( BUILDING_FOR_CARBONCORE_LEGACY )
double log2( double x )
{
return _xlog2( x );
}
#else
float log2f( float x )
{
return (float) _xlog2( (double) x );
}
#pragma mark -
#pragma mark log
double log ( double x )
{
xDouble xd = DOUBLE_2_XDOUBLE( x );
xd = _xlog( xd );
return XDOUBLE_2_DOUBLE( xd );
}
float logf( float x )
{ xDouble xd = FLOAT_2_XDOUBLE( x );
xd = _xlog( xd );
return XDOUBLE_2_FLOAT( xd );
}
#pragma mark -
#pragma mark log1p
static inline xDouble _xlog1p( xDouble x ) ALWAYS_INLINE;
static const double T1p = -0x1.F0540438FD5C4p-5; static const double T2p = 0x1.082B577D34ED8p-4; static const double T4p = 0x1.0000000000000p-53;
static inline xDouble _xlog1p( xDouble x )
{
static const double mTwo = -2.0;
static const double x52 = 52.0;
static const double eight = 8.0;
static const double oneEighth = 0.125;
static const double twom1 = 0x1.0p-1;
static const double twom9 = 0x1.0p-9;
xDouble Y, F, f, f1, L_lead, L_trail, R, u, u1, u2, v, q, g, m;
xSInt64 im, imc;
int j;
xDouble xIsNaN = _mm_cmpunord_sd( x, x );
if( _mm_istrue_sd( xIsNaN ) )
return _mm_add_sd( x, x );
xDouble xGTminusOne = _mm_cmpgt_sdm( x, &mOne );
xDouble xLTinfinity = _mm_cmplt_sdm( x, (double*) &plusinf );
if( _mm_istrue_sd( _mm_and_pd( xGTminusOne, xLTinfinity ) ) ) {
xDouble xLEt1p = _mm_cmple_sdm( x, &T1p );
xDouble xGEt2p = _mm_cmpge_sdm( x, &T2p );
xDouble fabsX = _mm_andnot_pd( minusZeroD, x );
xDouble bt = _mm_or_pd( xLEt1p, xGEt2p );
if( _mm_istrue_sd( bt ) ) {
double fparts[4];
int index;
Y = _mm_sub_sdm( x, &mOne );
im = _mm_srli_epi64( (xSInt64) Y, 52 ); R = _mm_mul_sdm( Y, &twom9 ); imc = _mm_sub_epi64( im, bias[1] ); im = _mm_add_epi64( im, (xSInt64) bt ); Y = _mm_mul_sdm( Y, &twom1 ); R = _mm_and_pd( R, plusinf ); m = _mm_cvtepi32_pd( imc ); F = _mm_add_sd( Y, R ); x = _mm_mul_sdm( x, &twom1 ); F = (xDouble) _mm_srli_epi64( (xSInt64) F, 52-7 ); im = _mm_slli_epi64( im, 7 );
im = _mm_sub_epi64( (xSInt64) F, im );
F = (xDouble) _mm_slli_epi64( (xSInt64) F, 52-7 );
j = _mm_cvtsi128_si32( im );
fparts[0] = 0.0;
fparts[1] = 0.5;
_mm_store_sd( &fparts[2], x );
_mm_store_sd( &fparts[3], Y );
xDouble mLEmtwo = _mm_cmple_sdm( m, &mTwo );
xDouble mLE52 = _mm_cmple_sdm( m, &x52 );
index = _mm_cvtsi128_si32( _mm_add_epi64( (xSInt64) mLEmtwo, (xSInt64) mLE52 ) );
f = _mm_sub_sdm( F, fparts + 2 + index ); f = _mm_sub_sdm( f, fparts + 1 - index );
L_lead = _mm_add_sdm( _mm_mul_sdm( m, C_lead + 128 ), C_lead + j );
L_trail = _mm_add_sdm( _mm_mul_sdm( m, C_trail + 128 ), C_trail + j);
Y = _mm_add_sd( Y, F ); f = _mm_mul_sdm( f, &mTwo ); u = _mm_div_sd( f, Y ); v = _mm_mul_sd( u, u );
q = _mm_add_sdm( _mm_mul_sdm( v, A1 + 1 ), A1 );
v = _mm_mul_sd( v, u );
q = _mm_mul_sd( q, v );
R = _mm_add_sd( q, L_trail );
R = _mm_add_sd( R, u );
R = _mm_add_sd( R, L_lead );
}
else if( _mm_istrue_sd( _mm_cmple_sdm( fabsX, &T4p ) ) )
{
xDouble xEQzero = _mm_cmpeq_sdm( x, (double*) &minusZeroD );
xDouble ittyBit = _mm_load_sd( &tiny );
ittyBit = _mm_andnot_pd( xEQzero, ittyBit );
R = _mm_mul_sdm( x, &eight );
R = _mm_sub_sd( R, ittyBit );
R = _mm_mul_sdm( R, &oneEighth );
R = _mm_sel_pd( R, x, xEQzero ); }
else
{
g = _mm_div_sd( one, _mm_sub_sdm( x, &mTwo ) );
u = _mm_mul_sd( x, g );
u = _mm_add_sd( u, u );
v = _mm_mul_sd( u, u );
q = _mm_add_sdm( _mm_mul_sdm( v, A2 + 3 ), A2 + 2 );
q = _mm_add_sdm( _mm_mul_sd( q, v ), A2 + 1 );
q = _mm_add_sdm( _mm_mul_sd( q, v ), A2 + 0 );
v = _mm_mul_sd( v, u );
q = _mm_mul_sd( q, v );
u1 = _mm_cvtss_sd( u, _mm_cvtsd_ss( (xFloat) u, u ) );
f1 = _mm_cvtss_sd( x, _mm_cvtsd_ss( (xFloat) x, x ) );
xDouble t0 = _mm_sub_sd( x, u1 ); xDouble t1 = _mm_sub_sd( x, f1 ); xDouble t2 = _mm_mul_sd( u1, f1 ); t0 = _mm_add_sd( t0, t0 ); t1 = _mm_mul_sd( t1, u1 ); t0 = _mm_sub_sd( t0, t2 ); t0 = _mm_sub_sd( t0, t1 ); u2 = _mm_mul_sd( t0, g );
R = _mm_add_sd( u2, q );
R = _mm_add_sd( R, u1 );
}
}
else
{ if( _mm_istrue_sd( _mm_cmpeq_sdm( x, &mOne ) ) )
R = _mm_div_sd( _mm_load_sd( &two ), minusZeroD );
else
{
xDouble inf = plusinf;
xDouble isInf = _mm_cmpeq_sd( x, inf );
inf = _mm_andnot_pd( xIsNaN, inf );
inf = _mm_andnot_pd( isInf, inf );
xDouble xLTzero = _mm_cmplt_sdm( x, (double*) &minusZeroD );
R = _mm_sub_sd( x, inf ); R = _mm_add_sd( R, inf );
xLTzero = _mm_and_pd( xLTzero, (xDouble) logNaN );
R = _mm_add_sd( xLTzero, R ); }
}
return R;
}
double log1p ( double x )
{
xDouble xd = DOUBLE_2_XDOUBLE( x );
xd = _xlog1p( xd );
return XDOUBLE_2_DOUBLE( xd );
}
float log1pf( float x )
{ xDouble xd = FLOAT_2_XDOUBLE( x );
xd = _xlog1p( xd );
return XDOUBLE_2_FLOAT( xd );
}
#pragma mark -
#pragma mark log10
static inline double _xlog10( double ) ALWAYS_INLINE;
static xUInt64 LOGORITHMIC_NAN = { 0x7FF8000000000000ULL, 0 };
static inline double _xlog10( double x )
{
static const double half = 0.5;
static const double log10e = 0.434294481903251827651128918916605082; xDouble xx = DOUBLE_2_XDOUBLE( x );
xDouble xIsNaN = _mm_cmpunord_sd( xx, xx );
xDouble safeX = _mm_andnot_pd( xIsNaN, xx );
xDouble xGTzero = _mm_cmpgt_sdm( safeX, (double*) &minusZeroD );
xDouble xLTinf = _mm_cmplt_sdm( safeX, (double*) &plusinf );
xDouble R;
if( _mm_istrue_sd( _mm_and_pd( xGTzero, xLTinf ) ) ) {
if( (x <= T1) || x >= T2 )
{
xDouble m;
xDouble Y = mantissa( xx, &m );
int j = _mm_cvttsd_si32( _mm_add_sdm( _mm_mul_sdm( Y, &d128 ), &half ));
xDouble F = _mm_mul_sdm( _mm_cvtsi32_sd( Y, j) , &d0_0078125 );
j -= 128;
xDouble f = _mm_sub_sd( Y, F );
xDouble L_lead = _mm_add_sdm( _mm_mul_sdm( m, &C_lead[128] ), &C_lead[j] );
xDouble L_trail = _mm_add_sdm( _mm_mul_sdm( m, &C_trail[128] ), &C_trail[j] );
xDouble u = _mm_div_sd( _mm_add_sd( f, f), _mm_add_sd( Y, F ) );
xDouble v = _mm_mul_sd( u, u );
xDouble q = _mm_add_sdm( _mm_mul_sdm( v, &A1[1] ), &A1[0] );
q = _mm_mul_sd( q, v );
q = _mm_mul_sd( q, u );
R = _mm_add_sd( L_lead, _mm_add_sd( u , _mm_add_sd( q, L_trail ) ) );
R = _mm_mul_sdm( R, &log10e );
}
else
{
xDouble xEQone = _mm_cmpeq_sdm( xx, (double*) &one );
xDouble f = _mm_sub_sdm( xx, (double*) &one );
xDouble g = _mm_div_sd( one, _mm_add_sdm( f, &two ) );
xDouble u = _mm_mul_sd( f, g ); u = _mm_add_sd( u, u ); xDouble v = _mm_mul_sd( u, u );
xDouble q = _mm_add_sdm( _mm_mul_sdm( v, &A2[3] ), &A2[2] );
q = _mm_add_sdm( _mm_mul_sd( q, v ), &A2[1] );
q = _mm_add_sdm( _mm_mul_sd( q, v ), &A2[0] );
q = _mm_mul_sd( _mm_mul_sd( q, v ), u );
xDouble u1 = _mm_cvtss_sd( u, _mm_cvtsd_ss( (xFloat) u, u ) );
xDouble f1 = _mm_cvtss_sd( f, _mm_cvtsd_ss( (xFloat) f, f ) );
xDouble t0 = _mm_sub_sd( f, u1 ); xDouble t1 = _mm_sub_sd( f, f1 ); xDouble t2 = _mm_mul_sd( u1, f1 ); t0 = _mm_add_sd( t0, t0 ); t1 = _mm_mul_sd( t1, u1 ); t0 = _mm_sub_sd( t0, t2 ); t0 = _mm_sub_sd( t0, t1 ); xDouble u2 = _mm_mul_sd( t0, g ); R = _mm_add_sd( u2, q );
R = _mm_add_sd( R, u1 );
R = _mm_mul_sdm( R, &log10e );
R = _mm_andnot_pd( xEQone, R );
}
}
else
{ if( _mm_istrue_sd( _mm_cmpeq_sd( xx, minusZeroD ) ) )
R = _mm_div_sd( _mm_load_sd( &two ), minusZeroD );
else
{
xDouble inf = plusinf;
xDouble isInf = _mm_cmpeq_sd( xx, inf );
inf = _mm_andnot_pd( xIsNaN, inf );
inf = _mm_andnot_pd( isInf, inf );
xDouble xLTzero = _mm_cmplt_sdm( safeX, (double*) &minusZeroD );
R = _mm_sub_sd( xx, inf ); R = _mm_add_sd( R, inf );
xLTzero = _mm_and_pd( xLTzero, (xDouble) logNaN );
R = _mm_add_sd( xLTzero, R ); }
}
return XDOUBLE_2_DOUBLE( R );
}
double log10 ( double x )
{
return _xlog10( x );
}
float log10f( float x )
{ return (float) _xlog10( (double) x );
}
#endif