Optimize Projection::determinant()

[raylibfan]

Describe the project you are working on

Research in Godot / Raylib / blender 3d matrix math performance differences.

Describe the problem or limitation you are having in your project

Projection determinant calculation can be found in module https://github.com/godotengine/godot/blob/master/core/math/projection.cpp LINE 40 and looks like:

real_t Projection::determinant() const {
	return columns[0][3] * columns[1][2] * columns[2][1] * columns[3][0] - columns[0][2] * columns[1][3] * columns[2][1] * columns[3][0] -
	columns[0][3] * columns[1][1] * columns[2][2] * columns[3][0] + columns[0][1] * columns[1][3] * columns[2][2] * columns[3][0] +
	columns[0][2] * columns[1][1] * columns[2][3] * columns[3][0] - columns[0][1] * columns[1][2] * columns[2][3] * columns[3][0] -
	columns[0][3] * columns[1][2] * columns[2][0] * columns[3][1] + columns[0][2] * columns[1][3] * columns[2][0] * columns[3][1] +
	columns[0][3] * columns[1][0] * columns[2][2] * columns[3][1] - columns[0][0] * columns[1][3] * columns[2][2] * columns[3][1] -
	columns[0][2] * columns[1][0] * columns[2][3] * columns[3][1] + columns[0][0] * columns[1][2] * columns[2][3] * columns[3][1] +
	columns[0][3] * columns[1][1] * columns[2][0] * columns[3][2] - columns[0][1] * columns[1][3] * columns[2][0] * columns[3][2] -
	columns[0][3] * columns[1][0] * columns[2][1] * columns[3][2] + columns[0][0] * columns[1][3] * columns[2][1] * columns[3][2] +
	columns[0][1] * columns[1][0] * columns[2][3] * columns[3][2] - columns[0][0] * columns[1][1] * columns[2][3] * columns[3][2] -
	columns[0][2] * columns[1][1] * columns[2][0] * columns[3][3] + columns[0][1] * columns[1][2] * columns[2][0] * columns[3][3] +
	columns[0][2] * columns[1][0] * columns[2][1] * columns[3][3] - columns[0][0] * columns[1][2] * columns[2][1] * columns[3][3] -
	columns[0][1] * columns[1][0] * columns[2][2] * columns[3][3] + columns[0][0] * columns[1][1] * columns[2][2] * columns[3][3];
}

It takes 72 multiplication to calculate 4x4 matrix determinant as one can see.

Meanwhile, accordinly to this Wikipedia's article https://en.wikipedia.org/wiki/Laplace_expansion, this calculation can be reduced to (4 + 4 * (3 + 3 * 2)) = 40 multiplication by decreasing matrix size from 4x4 to 2x2 using minors.

Such algorithm is used in blender 3d application source code as for example: https://github.com/blender/blender/blob/main/source/blender/blenlib/intern/math_matrix_c.cc from LINE 1837

So we can replace those existing code fragment from above by something like this (slower variant):

real_t Projection::determinant const {
									
	return (columns[0][0] * ((columns[1][1] * (columns[2][2] * columns[3][3] - columns[2][3] * columns[3][2]) - columns[2][1] * (columns[1][2] * columns[3][3] - columns[1][3] * columns[3][2]) + columns[3][1] * (columns[1][2] * columns[2][3] - columns[1][3] * columns[2][2]))) -
			columns[1][0] * ((columns[0][1] * (columns[2][2] * columns[3][3] - columns[2][3] * columns[3][2]) - columns[2][1] * (columns[0][2] * columns[3][3] - columns[0][3] * columns[3][2]) + columns[3][1] * (columns[0][2] * columns[2][3] - columns[0][3] * columns[2][2]))) +
			columns[2][0] * ((columns[0][1] * (columns[1][2] * columns[3][3] - columns[1][3] * columns[3][2]) - columns[1][1] * (columns[0][2] * columns[3][3] - columns[0][3] * columns[3][2]) + columns[3][1] * (columns[0][2] * columns[1][3] - columns[0][3] * columns[1][2]))) -
			columns[3][0] * ((columns[0][1] * (columns[1][2] * columns[2][3] - columns[1][3] * columns[2][2]) - columns[1][1] * (columns[0][2] * columns[2][3] - columns[0][3] * columns[2][2]) + columns[2][1] * (columns[0][2] * columns[1][3] - columns[0][3] * columns[1][2]))) );
	
}

or this (faster variant):

real_t Projection::determinant() const {
									
	real_t m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15;

	m0 = columns[0][0]; m4 = columns[1][0]; m8 =  columns[2][0]; m12 = columns[3][0];
	m1 = columns[0][1]; m5 = columns[1][1]; m9 =  columns[2][1]; m13 = columns[3][1];
	m2 = columns[0][2]; m6 = columns[1][2]; m10 = columns[2][2]; m14 = columns[3][2];
	m3 = columns[0][3]; m7 = columns[1][3]; m11 = columns[2][3]; m15 = columns[3][3];

	return (m0 *  ((m5 * (m10 * m15 - m11 * m14) - m9 * (m6 * m15 - m7 * m14) + m13 * (m6 * m11 - m7 * m10))) -
		m4 *  ((m1 * (m10 * m15 - m11 * m14) - m9 * (m2 * m15 - m3 * m14) + m13 * (m2 * m11 - m3 * m10))) +
		m8 *  ((m1 * (m6  * m15 - m7  * m14) - m5 * (m2 * m15 - m3 * m14) + m13 * (m2 * m7  - m3 * m6))) -
		m12 * ((m1 * (m6  * m11 - m7  * m10) - m5 * (m2 * m11 - m3 * m10) + m9  * (m2 * m7  - m3 * m6))));
	
}

Describe the feature / enhancement and how it helps to overcome the problem or limitation

We can achieve 1.7x to 1.8x faster determinant calculation.

I wrote some simple test application to check perfomance boost. Look on green and dark green results on screenshot:

Describe how your proposal will work, with code, pseudo-code, mock-ups, and/or diagrams

One should just replace old Projection::determinant() implementation with faster one from above.

If this enhancement will not be used often, can it be worked around with a few lines of script?

No.

Is there a reason why this should be core and not an add-on in the asset library?

This enhancement will be used often because its in the core math module.

[sockeye-d]

How often is this method actually called? I don't think that it shouldn't be implemented but would users actually see a real performance gain?

[Calinou]

How often is this method actually called? I don't think that it shouldn't be implemented but would users actually see a real performance gain?

I doubt it'll benefit performance much in most projects, but this kind of optimization generally saves on CPU cycles and therefore reduce power consumption on mobile devices.

[raylibfan]

One can think of the Projection class in more wide viewing angle. It's just a regular 4x4 matrix of float (or double) with handy set of built-in functions by its nature. So, from mathematical point of view, its usage limited only by developer fantasy.

So, for example, just by 4 calls of determinant() method and Projection class constructor user can find:

1. Solution of a system of linear equations by Cramer's rule (We did this when I was in high school). Not every Godot project should be a game. It can be some kind of equation solver for lazy student with Godot gui. Maybe exported on a smartphone. In those scenario perfomance is not the key factor;
2. An intersection point of three geometric planes (or triangles or quads) in 3d space;
3. An intersection point of one geometric plane (or triangle or quad) and vector (like in raycasting).
   In examples 2 and 3 perfomance is the key.

In 3d game it can be an archer summoning cloud of hundreds of arrow from the sky (like in Path of Exile) or hundreds of archers shoots with one arrow each (like in Total War series). To calculate a position (Vector3) array where all of the arrows hits the 3d-mesh (landscape, or boss, strucrures, e.t.c.) determinant() function can be used hundreds of times per second per mesh triangle (it's normal if to be clear).