SIMD implementation of dot-product for real numbers

realdotproduct.d

F dotProduct(F)(in F[] a, in F[] b) @trusted
if(is(F == double) || is(F == float))
in
{
	assert(a.length == b.length);
}
out(result)
{
	assert(!result.isNaN);
} 
body
{
	enum N = MaxVectorSizeof / F.sizeof;
	alias VP = const(Vector!(F[N]))*;

	auto ap = a.ptr, bp = b.ptr; 
	const n = a.length;
	const end = ap + n;
	const ste = ap + (n & -N);
	const sbe = ap + (n & -N*4);
	
	Vector!(F[N]) s0 = 0, s1 = 0, s2 = 0, s3 = 0;	

	for(; ap < sbe; ap+=4*N, bp+=4*N)
	{
		s0 += load(cast(VP)ap+0) * load(cast(VP)bp+0);
		s1 += load(cast(VP)ap+1) * load(cast(VP)bp+1);
		s2 += load(cast(VP)ap+2) * load(cast(VP)bp+2);
		s3 += load(cast(VP)ap+3) * load(cast(VP)bp+3);
	}

	s0 = (s0+s1)+(s2+s3);
		
	for(; ap < ste; ap+=N, bp+=N)
		s0 += load(cast(VP)ap+0) * load(cast(VP)bp+0);

	//reduce to scalar
	F s = toScalarSum(s0);

	for(; ap < end; ap++, bp++)
		s += ap[0] * bp[0];
	
	return s;
}