# 4×4 cofactor in SSE

The cofactor of a 4×4 matrix can be used to convert a "regular geometry" matrix into the matrix that transforms the normals. It's an alternative to the common inverse-transpose pattern. In this post I use row-major matrixes.

Here is reference scalar code from https://github.com/graphitemaster/normals_revisited

float minor(const float m[16], int r0, int r1, int r2, int c0, int c1, int c2) {
return m[4*r0+c0] * (m[4*r1+c1] * m[4*r2+c2] - m[4*r2+c1] * m[4*r1+c2]) -
m[4*r0+c1] * (m[4*r1+c0] * m[4*r2+c2] - m[4*r2+c0] * m[4*r1+c2]) +
m[4*r0+c2] * (m[4*r1+c0] * m[4*r2+c1] - m[4*r2+c0] * m[4*r1+c1]);
}

void cofactor(const float src[16], float dst[16]) {
dst[ 0] =  minor(src, 1, 2, 3, 1, 2, 3);
dst[ 1] = -minor(src, 1, 2, 3, 0, 2, 3);
dst[ 2] =  minor(src, 1, 2, 3, 0, 1, 3);
dst[ 3] = -minor(src, 1, 2, 3, 0, 1, 2);
dst[ 4] = -minor(src, 0, 2, 3, 1, 2, 3);
dst[ 5] =  minor(src, 0, 2, 3, 0, 2, 3);
dst[ 6] = -minor(src, 0, 2, 3, 0, 1, 3);
dst[ 7] =  minor(src, 0, 2, 3, 0, 1, 2);
dst[ 8] =  minor(src, 0, 1, 3, 1, 2, 3);
dst[ 9] = -minor(src, 0, 1, 3, 0, 2, 3);
dst[10] =  minor(src, 0, 1, 3, 0, 1, 3);
dst[11] = -minor(src, 0, 1, 3, 0, 1, 2);
dst[12] = -minor(src, 0, 1, 2, 1, 2, 3);
dst[13] =  minor(src, 0, 1, 2, 0, 2, 3);
dst[14] = -minor(src, 0, 1, 2, 0, 1, 3);
dst[15] =  minor(src, 0, 1, 2, 0, 1, 2);
}


With SSE it could be implemented as:

void cofactorSSE(const float src[16], float dst[16]) {

__m128 r0_0001 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(0, 0, 0, 1));
__m128 r1_0001 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(0, 0, 0, 1));
__m128 r2_0001 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(0, 0, 0, 1));
__m128 r3_0001 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(0, 0, 0, 1));
__m128 r0_1122 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(1, 1, 2, 2));
__m128 r1_1122 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(1, 1, 2, 2));
__m128 r2_1122 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(1, 1, 2, 2));
__m128 r3_1122 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(1, 1, 2, 2));
__m128 r0_2333 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(2, 3, 3, 3));
__m128 r1_2333 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(2, 3, 3, 3));
__m128 r2_2333 = _mm_shuffle_ps(r2, r2, _MM_SHUFFLE(2, 3, 3, 3));
__m128 r3_2333 = _mm_shuffle_ps(r3, r3, _MM_SHUFFLE(2, 3, 3, 3));
__m128 odd = _mm_set_ps(0.0, -0.0, 0.0, -0.0);
__m128 even = _mm_set_ps(-0.0, 0.0, -0.0, 0.0);

__m128 res0 = _mm_mul_ps(_mm_mul_ps(r1_0001, r2_1122), r3_2333);
__m128 res1 = _mm_mul_ps(_mm_mul_ps(r0_0001, r2_1122), r3_2333);
__m128 res2 = _mm_mul_ps(_mm_mul_ps(r0_0001, r1_1122), r3_2333);
__m128 res3 = _mm_mul_ps(_mm_mul_ps(r0_0001, r1_1122), r2_2333);
res0 = _mm_add_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_1122, r2_2333), r3_0001));
res1 = _mm_add_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_1122, r2_2333), r3_0001));
res2 = _mm_add_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_2333), r3_0001));
res3 = _mm_add_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_2333), r2_0001));
res0 = _mm_add_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_2333, r2_0001), r3_1122));
res1 = _mm_add_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_2333, r2_0001), r3_1122));
res2 = _mm_add_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_0001), r3_1122));
res3 = _mm_add_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_0001), r2_1122));
res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_2333, r2_1122), r3_0001));
res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_2333, r2_1122), r3_0001));
res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_1122), r3_0001));
res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_2333, r1_1122), r2_0001));
res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_1122, r2_0001), r3_2333));
res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_1122, r2_0001), r3_2333));
res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_0001), r3_2333));
res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_1122, r1_0001), r2_2333));
res0 = _mm_sub_ps(res0, _mm_mul_ps(_mm_mul_ps(r1_0001, r2_2333), r3_1122));
res1 = _mm_sub_ps(res1, _mm_mul_ps(_mm_mul_ps(r0_0001, r2_2333), r3_1122));
res2 = _mm_sub_ps(res2, _mm_mul_ps(_mm_mul_ps(r0_0001, r1_2333), r3_1122));
res3 = _mm_sub_ps(res3, _mm_mul_ps(_mm_mul_ps(r0_0001, r1_2333), r2_1122));

_mm_store_ps(&dst[0], _mm_xor_ps(res0, even));
_mm_store_ps(&dst[4], _mm_xor_ps(res1, odd));
_mm_store_ps(&dst[8], _mm_xor_ps(res2, even));
_mm_store_ps(&dst[12], _mm_xor_ps(res3, odd));
}


But it's not that fast compared to scalar code, only about 2.6 times faster (compiled with MSVC 2010, tested on Haswell). Of course, there are many shuffles. Also a couple of spills, which don't help, but there are not as many as it naively looks like there might be (MSVC changes the order and interleaves the shuffles and their uses somewhat).

Can this be done more efficiently? For example, a trick to save some shuffles? SSE through SSE3 can definitely be used, AVX(2) answers are also interesting (but I will have to dual-wield two versions of the function). AVX512 as a bonus but I doubt it will see much use.

• You could use restrict (C11) in your parameters: void cofactorSSE(const float [static restrict 16], float [static restrict 16]). Maybe it would help optimizations. I would be surprised if MSVC supported it, though. Commented Aug 6, 2019 at 16:40
• @CacahueteFrito isn't that a C thing technically? sorry for the C tag by the way someone else added it Commented Aug 6, 2019 at 16:43
• Ah sorry then. I saw the C tag :). GCC has it also in C++ as an extension (as __restrict__, but you can easily do the macro). Maybe other compilers also provide it as an extension. Commented Aug 6, 2019 at 16:43
• Also, if you care about performance, why not program the function in C (or Fortran or assembly) and just export it to C++ with an extern "C"? Commented Aug 6, 2019 at 16:48
• @CacahueteFrito I've tried __restrict now but it didn't do anything. I don't have an assembly version that's actually better, if you have one I can consider that Commented Aug 6, 2019 at 17:02