Skip to content

Instantly share code, notes, and snippets.

@Const-me
Created July 25, 2016 16:17
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Const-me/c582805223e79cc11f2ecb4935ee148e to your computer and use it in GitHub Desktop.
Save Const-me/c582805223e79cc11f2ecb4935ee148e to your computer and use it in GitHub Desktop.
A benchmark comparing three methods of sine + cosine computation: stdlib, DirectX which is 10-11 degree polynomial approximation, and linearly-interpolated lookup table.
#include "stdafx.h"
using std::vector;
typedef std::chrono::high_resolution_clock stopwatch;
using namespace DirectX;
using std::array;
struct StdLib
{
inline void sinCos( float val, float& s, float& c ) const
{
s = sinf( val );
c = cosf( val );
}
};
struct DX
{
inline void sinCos( float val, float& s, float& c ) const
{
XMScalarSinCos( &s, &c, val );
}
};
class LinInt
{
static const size_t size = 256;
// Two tables for sin and cos, from -PI/2 to +PI/2, interleaved for cache friendliness
array<float, ( size + 1 ) * 2> lookupTable;
const float indexMul;
// value should be from -PI/2 to +PI/2
inline void lookup( float value, float& sin, float& cos ) const
{
// Calculate index + coefficients for linear interpolation
value *= indexMul;
int i1 = int( floor( value ) );
value -= i1;
i1 += ( size / 2 );
const float b = 1.0f - value;
// Interpolate both sin + cos using same coefficients.
const float* entries = lookupTable.data() + ( i1 << 1 );
sin = entries[ 0 ] * b + entries[ 2 ] * value;
cos = entries[ 1 ] * b + entries[ 3 ] * value;
}
public:
LinInt() : indexMul( float( size / M_PI ) )
{
for( int i = 0; i <= size; i++ )
{
double val = M_PI * ( ( i - ( size / 2 ) ) / double( size ) );
lookupTable[ i * 2 ] = float( sin( val ) );
lookupTable[ i * 2 + 1 ] = float( cos( val ) );
}
}
inline void sinCos( float Value, float& sin, float& cos ) const
{
// Map Value to y in [-pi,pi], x = 2*pi*quotient + remainder.
float quotient = XM_1DIV2PI * Value;
if( Value >= 0.0f )
{
quotient = (float)( (int)( quotient + 0.5f ) );
}
else
{
quotient = (float)( (int)( quotient - 0.5f ) );
}
float y = Value - XM_2PI * quotient;
// Map y to [-pi/2,pi/2] with sin(y) = sin(Value).
float sign;
if( y > XM_PIDIV2 )
{
y = XM_PI - y;
sign = -1.0f;
}
else if( y < -XM_PIDIV2 )
{
y = -XM_PI - y;
sign = -1.0f;
}
else
{
sign = +1.0f;
}
// Interpolate both
float cosNoSign;
lookup( y, sin, cosNoSign );
// Restore cos sign
cos = cosNoSign * sign;
}
};
template<class Algo>
void testSinCos( const vector<float>& src, const Algo& algo )
{
float rs = 0, rc = 0;
float s, c;
auto start = stopwatch::now();
for( float i : src )
{
algo.sinCos( i, s, c );
rs += s;
rc += c;
}
auto stop = stopwatch::now();
std::chrono::duration<double> secs = stop - start;
const double avgSin = double( rs ) / double( src.size() );
const double avgCos = double( rc ) / double( src.size() );
printf( "%f seconds, the averages are %f / %f\n",
secs.count(), avgSin, avgCos );
}
int main()
{
srand( GetTickCount() );
// static const size_t testSize = 10 * 1000;
static const size_t testSize = 10 * 1000 * 1000;
vector<float> src;
src.resize( testSize );
for( float& f : src )
f = float( 2.0 * rand() * M_PI / RAND_MAX );
// const StdLib algo;
// const DX algo;
const LinInt algo;
testSinCos( src, algo );
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment