Whenever one needs a curve, a spline of some sort is usually involved. The most common forms are Catmull-Rom splines and Cubic Bezier Splines. Both of these splines have relatively simple polynomial forms, and are easily evaluated in the vertex shader.

A more complicated issue is the question of how to move along such a curve at a constant speed. Anyone looking for such a process is likely to stumble the root finding method, specified in the paper “Moving Along a Curve with Specified Speed”, by David Eberly. The process is straightforward – given a length, solve for the time that the curve reaches that length through a modified Newton’s method. This method involves repeated computations of length, which is a slow operation. Calculating more than a few dozen reparameterizations is not possible on high end hardware, so the problem becomes one of acceleration.

The first step is to simply generate better initial guesses for Newton’s method, so that the number of iterations is minimized. At a cost of 4 bytes per spline, storing the previous result and feeding it (plus a small constant) to the next attempt to find a root reduces the number of iterations by two or three. This saves several integrals, but is still too slow. A more aggressive caching policy involves a preprocessing step. Instead of solving for time, a single integral can generate lengths for every time step. This generates a set of unevenly distributed length, time pairs, which can be linearly interpolated to approximate a reparameterization.

The number of samples directly correlates to the error in the approximation, as shown in the image below. The numbers shown in the graph are derived from many random splines. A larger version can be found here

4 samples will keep the difference below 30%, and each additional doubling of the number of samples halves the error level. 16 samples has roughly 5% error, which is good enough for most graphical applications. 16 samples also happens to fit into 4 float4 registers in a shader, so the linear fit interpolation can be done in hardware without any additional CPU work beyond computing the initial fit.

There appear to be some major issues with the node application that runs this website. It crashes after some time, without a visible stack trace or anything. I'm not sure why this happens - the application is running under a daemon that should restart it if it dies, but this appears also not to work.

Either way, the golang application thats running on the same server is running just fine after several months. More investigation is needed, but I can't help but think that a better language wouldn't have these problems.

While reading Hacker News, I came upon an article that purported that Java (or other high level languages) should be used in preference to C++ due to the quality of modern interpreters and JIT compilers. I will not discuss the conclusion, which is shaky at best, but rather the microbenchmark that was used to display the quality of the Java virtual machine.

The benchmark used was a cumulative sum method:

void slowishSum(vector<int> & data) 
{
    if(data.size() == 0) 
        return;

    for (size_t i = 1; i < data.size(); ++i) 
    {
        data[i] += data[i - 1];
    }
}

On my machine, this does about 460~ million integers per second. This would be great, if it was all that the processor could do - but its not.

Many articles have been written on the benefits of cache coherency and data locality, with the idea that if you process your data in order, the CPU will load up the requested data (and a little more) into L1 cache, making further accesses to memory around there very fast. Cache awareness in programming can lead to significant performance improvements. Changing an algorithm from a AOS (Array of Structures) access pattern to a SOA (Structure of Arrays) access pattern can improve performance upwards of 200%.

This does not apply to the cumulative sum function. All memory accesses here are sequential, allowing the CPU to have optimal cache use. So how do you make this method faster? One can unroll the loop, through usage of -funroll-loops (the best gcc flag name ever), or manually:

void sum(vector<int> & data) 
{
    if(data.size() == 0) 
        return;

    size_t i = 1;
    for (; i < data.size() - 1; i+=2) 
    {
        data[i] += data[i - 1];
        data[i+1] += data[i ];
    }

    for (; i != data.size(); ++i) 
    {
        data[i] += data[i - 1];
    }
}

This does not help, however, computing a similar 470~ million integers per second.

To further optimize this method, a basic understanding of CPU architecture is required. A very basic view of CPU operation is that it goes through several distinct steps while processing an instruction - fetching the instruction from memory, decoding the instruction, execution of the instruction, and writing the results back into memory. Each of these steps are independent of each other, so a basic CPU can pipeline four instructions at once - one doing writeback, one doing the execution, one decoding and the last fetching from memory. Modern CPUs are more complicated, and their pipelines deeper, but the general concept remains the same. However, pipelining introduces data hazards. Consider a program that executes 'INC A; STORE B A'. Without hazard mitigation, the increment may execute, and then while it is in the write-back stage, the store instruction executes, fetching the unincremented value of A and storing it into B. Hazards can be mitigated by hand, inserting additional instructions into the pipeline so that this situation never happens, or by the compiler, doing a similar operation. At the CPU level, the CPU can detect data dependencies and intentionally add in bubbles to the pipeline, so that most hazards are avoided.

A naive cumulative sum implementation, as in the first code block above, is almost the worst case for the CPU pipeline. Each iteration of the loop needs to load up two memory locations, one of which is dependent upon the operation in the previous loop. Even considering the instructions that are used to do loop maintenance (incrementing the counter, comparing against the size, jump if zero), the CPU must still insert a bubble to delay the dependent load. With processors having rather long pipelines, this can cause a tremendous slowdown. Therefore, removing this bubble should show some performance improvement. Without thinking too much, it is completely possible to compute several elements of the cumulative sum at the same time, without having any dependent data loads within the body of the loop -

void avoidStalls(vector<int>& data)
{
	if (data.size() == 0) 
	{
		return;
	}

	int remain = data[0];
	size_t i = 1;
	size_t size = data.size() - 2;
	for (i = 1; i < size; i += 3) 
	{
		int a = data[i + 0];
		data[i + 0] = remain + a;

		int b = data[i + 1];
		data[i + 1] = remain + a + b;

		int c = data[i + 2];
		data[i + 2] = remain + a + b + c;

		remain = data[i + 2];
	}

	for (; i < data.size(); ++i) 
	{
		data[i] += data[i - 1];
	}
}

At first glance, this code looks less efficient, since there is a lot of repeated work. However, doing the computation this way removes most of the data dependencies from the loop. A quick performance comparision shows that this version can compute 1040 million integers per second or more, an improvement of over 100%. This benchmark was compiled on MSVC 2012, with default settings for the 'Release' configuration. Other people, using GCC, have reported numbers similar (and some much higher) than those I obtained using the stall avoiding algorithm. It may be that they have better computers, or that GCC is capable of optimizing this dependent load into a better-for-the-pipeline algorithm, but I cannot make any definitive remarks on that.

Practically, I would like to say that it is always important to consider the processor and its memory while writing performance oriented code, but this is honestly a micro-optimization, and should be looked over for more fundamental changes for performance. In general, it may be hard to find patterns that are so obviously stalling the pipeline, and harder to change the algorithm to avoid such stalls. In any case, knowing this property of processors is always helpful while writing performance sensitive code, so beware.

Just a quick aside for anyone that may be using luajit with luabind on x86 platforms - be careful about the way that luabind throws exceptions. Luajit on x86 platforms has no exception interopability, which means that throwing exceptions through lua frames is not supported, and may corrupt the stack. This usually is not a problem, since luabind generally only throws exceptions after an error, or some other confluence of events conspires so that the exception does not corrupt anything. However, the following code will create a buffer overflow in the lua internals -

void show() {}

// ... Setup omitted. Assume show is bound to lua.

lua_loadstring(L, "show(1);");
luabind::object obj(luabind::from_stack(L, -1));
obj();

A workaround, without having to remove exception support from luabind, is to use the raw lua API to invoke a chunk, instead of turning it into a luabind::object - that is, invoke lua_pcall by hand and handle errors manually.

According to the most recent Steam Hardware Survey, 99% of steam users have SSE2 and SSE3 instruction support. The most immediately obvious thing to do with SSE is to try to optimize 4-vector math. However, a naiive reimplementation of vector methods doesn't come close to the theoretical 4x improvement that one would expect from an SSE implementation

Lets start with the naiive implementation of a vector4, and the naiive implementation of operator +=:

struct Vector4
{
	//Assume that this has constructors
	//and other reasonable methods on it,
	//like dot, operator[], other operators, 
	//implemented in the most straightforward way.

	Vector4& operator+=(const Vector4& other)
	{
		for (size_t i = 0; i < 4; ++i)
		{
			elements[i] += other.elements[i];
		}

		return this;
	}

	float elements[4];	
};

Wonderful. Next, a test rig just to test this basic functionality:

Vector4 v = Vector4(r[0], r[1], r[2], r[3]);
Vector4 b = Vector4(r[4], r[5], r[6], r[7]);

LARGE_INTEGER begin;
QueryPerformanceCounter(&begin);

Vector4 c = Vector4::Zero;
for (size_t i = 0; i < 100000000; ++i)
{
	c += v + b;
	c -= b + v;
}

LARGE_INTEGER end;
QueryPerformanceCounter(&end);

LARGE_INTEGER freq;
QueryPerformanceFrequency(&freq);

double time = (double)(end.QuadPart - begin.QuadPart) / (freq.QuadPart);
std::cout << "Took: " << 1000 * time << "ms\n";

The expected output is [0, 0, 0, 0], and r is populated with random data (no denormals).

Using the naiive implementation, with /Ox (Full optimization) and no architecture flags, the generated code simply loads up every floating point value, converts to double, adds it, and then converts backwards. The core of the operation is shown below - a bunch of flds, fadds and fstp.

009714AB  fadd        dword ptr [esp+5Ch]  
009714AF  fstp        dword ptr [esp+98h]  
009714B6  fld         dword ptr [esp+8Ch]  
009714BD  fadd        dword ptr [esp+7Ch]  
009714C1  fstp        dword ptr [esp+7Ch]  
009714C5  fld         dword ptr [esp+90h]  
009714CC  fadd        dword ptr [esp+80h]  
009714D3  fstp        dword ptr [esp+80h]  
009714DA  fld         dword ptr [esp+94h]  
009714E1  fadd        dword ptr [esp+84h]  
009714E8  fstp        dword ptr [esp+84h]  
009714EF  fld         dword ptr [esp+98h]  
009714F6  fadd        dword ptr [esp+88h]  
009714FD  fstp        dword ptr [esp+88h]  

The problem is that the floating point model is set to precise, which forces 80-bit precision on the FPU. Changing the floating point model to fast does not net any gains in this department - the generated code is functionally identical between the versions.

The baseline performance on my i7 here is roughly 810 milliseconds with /fp:fast, and 1080 milliseconds with /fp:precise. This means that the target performance for this loop is 202.5 milliseconds, which will represent a 4x speedup.

The least labor intensive way of preforming SSE optimizations is to let the compiler do it for you. While Visual Studio 2010 doesn't allow targeting SSE3 architectures specifically, the SSE3 instructions do not really add very much - the only interesting instruction for short vector operations is hadd, the horizontal add. HADD can be used to implement SSE dot products (ignoring the SSE4 dot product intrinsic, which is actually pretty slow), but only shows significant gains if you do 4 dot products at once.

Either way, enabling /arch:SSE2 produces shockingly bad code. Basically, it loads every element into its own XMM register with movss, and then adds them one by one. Since there are not enough registers to load all 13 floating point values at one value per register, the compiler generates a load of movss into memory, which incurs even more slowdown.

;;; This loads up a bunch of elements into XMM registers 0-3.
0024139E  movss       xmm1,dword ptr [esp+80h]  
002413A7  movss       xmm2,dword ptr [esp+84h]  
002413B0  movss       xmm3,dword ptr [esp+88h]  
002413B9  addss       xmm1,dword ptr [esp+70h]  
002413BF  addss       xmm2,dword ptr [esp+74h]  
002413C5  addss       xmm3,dword ptr [esp+78h]  
002413CB  movss       dword ptr [esp+30h],xmm0  
002413D1  movss       xmm0,dword ptr [Math::Vector<float,4>::Zero+4 (24631Ch)]  
002413D9  movss       dword ptr [esp+34h],xmm0  
002413DF  movss       xmm0,dword ptr [Math::Vector<float,4>::Zero+8 (246320h)]  
002413E7  movss       dword ptr [esp+38h],xmm0  
002413ED  movss       xmm0,dword ptr [Math::Vector<float,4>::Zero+0Ch (246324h)]  
002413F5  movss       dword ptr [esp+3Ch],xmm0  
002413FB  movss       xmm0,dword ptr [esp+7Ch]  
00241401  addss       xmm0,dword ptr [esp+6Ch]  

;;; Integer conversion code omitted 

00241442  movss       xmm4,dword ptr [esp+10h]  
00241448  movaps      xmm5,xmm4  
0024144B  addss       xmm5,xmm6  
0024144F  addss       xmm5,dword ptr [esp+30h]  
00241455  movaps      xmm6,xmm4  
00241458  addss       xmm6,xmm7  
0024145C  movss       dword ptr [esp+54h],xmm6  
00241462  movaps      xmm6,xmm4  
00241465  movaps      xmm7,xmm2  
00241468  addss       xmm6,xmm7  
0024146C  movss       xmm7,dword ptr [esp+78h]  
00241472  addss       xmm6,dword ptr [esp+38h]  
00241478  addss       xmm7,xmm4  
0024147C  addss       xmm7,dword ptr [esp+3Ch]  
00241482  movss       dword ptr [esp+3Ch],xmm7  

After all of this, the actual performance is actually a bit better - 700 milliseconds for the same loop, a 1.15x improvement on the naiive implementation, but still far from the idealized 210 millisecond runtime that an SSE implementation should aim for.

Before actually implementing the SSE version, certain constraints need to be added here. SSE performs optimally when loading from 16 byte aligned addresses, which means that, by default, any heap allocation will likely be incompatible with the movps instructions. One can specify that any stack allocation should be aligned with a declspec, but this causes certain problems. In particular, the vector implementation that ships with Visual Studio contains at least this troublesome line:

void resize(size_type _Newsize, _Ty _Val)

When _Ty is an aligned type, the template will fail to compile, since it is not possible, in the general case, to push an aligned value onto the stack in preparation for a function call. Changing the container to use boost::container::vector, which uses a constant reference instead of a value, allows the container to compile.

However, actually using the container causes memory problems - movps expects an aligned memory address, but the container doesn't guarantee any such alignment. This will result in a runtime crash. A solution here is to store pointers to your data, but this design will invalidate your cache like mad. A proper solution is to use an aligned allocator, but this causes template type issues. Without support for template template aliases, the next best solution is to use some type computation to generate a container with the properly aligned allocator, but this has maintainability problems.

Another problem to consider here is that repeated calls to movps actually will do redundant loads if you're not careful. The easiest way to make sure that move instructions are generated only when they are required is to invoke a small bit of undefined behavior, and store both the __m128 register and the floating point representation together:

__declspec(align(16)) union storage
{
	__m128 sse
	float  fpu[4];
};

The declspec above is not strictly required, since that can be inferred from the __m128 value, but it doesn't harm the code generation at all, and it is always nice to be explicit. The actual implementation lives in the rmath repository, and not be posted here. The generated code here is not quite optimal, but it is reasonable:

00F71406  movaps      xmm2,xmmword ptr [esp+40h]  
00F7140B  movdqa      xmm1,xmmword ptr [Math::Vector<float,4>::Zero (0F766C0h)]  
00F71413  addps       xmm2,xmmword ptr [esp+10h]  

;;; Integer division and conversion ommitted ...

; Move (float)(i / 200) into xmm0, and shuffle it into all components, and 
; some more math operations 

00F71441  movss       xmm0,dword ptr [esp+0Ch]  
00F71447  fstp        dword ptr [esp+10h]  
00F7144B  shufps      xmm0,xmm0,0  
00F7144F  addps       xmm0,xmm2  
00F71452  addps       xmm0,xmm1  
00F71455  movss       xmm1,dword ptr [esp+10h]  
00F7145B  shufps      xmm1,xmm1,0  
00F7145F  addps       xmm1,xmm2  
00F71462  subps       xmm0,xmm1  
00F71465  movaps      xmm1,xmm0  

This code executes in 280 milliseconds, an 2.8x improvement over the original, non-SSE code. Interestingly, the vast majority of the time is spent within the integer division and subsequent floating point conversion. Removing the code shaves 200 milliseconds off the naiive implementation, to 600ms, and improves the sse implementation to 175ms, a 3.4x improvement. This improvement is significant, but only really affects vector computation heavy code, which should have been already instruction parallelized. Using an SSE version of a vector also requires special container considerations, which might require significant code rework. Considering the actual level of code that would benefit from this optimization, and the amount of code that would be required to change out all the containers and allocators, it may not even be a good cost/benefit ratio.