1 /**
2 dlsl.matrix
3 
4 
5 Authors: Peter Particle ( based on gl3n by David Herberth )
6 License: MIT
7 
8 Note: All methods marked with pure are weakly pure since, they all access an instance member.
9 All static methods are strongly pure.
10 */
11 
12 module dlsl.matrix;
13 
14 import dlsl.vector;
15 
16 import std.conv	: to;
17 import std.math	: sqrt, sin, cos, tan, PI;
18 import std.array : join;
19 import std.traits : isNumeric, isFloatingPoint, isArray;
20 import std.string : format, rightJustify;
21 import std.algorithm : max, min, reduce;
22 
23 
24 /// predefined matrix types, the first number represents the number of cols 
25 /// and the second the number of columns, if there's just one it's a nxn matrix.
26 /// All of these matrices are floating - point matrices.
27 alias Matrix!( float, 2, 2 ) mat2;
28 alias Matrix!( float, 2, 3 ) mat23;
29 alias Matrix!( float, 2, 4 ) mat24;
30 
31 alias Matrix!( float, 3, 2 ) mat32;
32 alias Matrix!( float, 3, 3 ) mat3;
33 alias Matrix!( float, 3, 4 ) mat34;
34 
35 alias Matrix!( float, 4, 2 ) mat42;
36 alias Matrix!( float, 4, 3 ) mat43;
37 alias Matrix!( float, 4, 4 ) mat4;
38 
39 
40 enum float	deg2rad	= 0.0174532925199432957692369076849f;
41 enum float	rad2deg	= 57.295779513082320876798154814105f;
42 
43 version(NoReciprocalMul) {
44 	private enum rmul = false;
45 }	else {
46 	private enum rmul = true;
47 }
48 
49 /// If T is a matrix, this evaluates to true, otherwise false
50 template isMatrix( T ) {  enum isMatrix = is( typeof( isMatrixImpl( T.init )));  }
51 private void isMatrixImpl( T, int cols, int rows )( Matrix!( T, cols, rows ) mat ) {}
52 
53 
54 /// Base template for all matrix types. See predefined matrix types
55 /// Params:
56 /// type = the value type of each matrix element
57 /// cols = count of columns of the matrix
58 /// rows = count of rows of the matrix
59 /// - - - 
60 struct Matrix( type, int numCols, int numRows ) if(( numCols > 1 ) && ( numRows > 1 )) {
61 	alias type valueType;						/// Holds the internal primitive type of the matrix;
62 	alias Vector!( type, rows ) vectorType;		/// Holds the internal vectorType of the matrix;
63 
64 	/// Holds the matrix $( RED column - major ) in memory.
65 	vectorType[ cols ] data;					/// Each Column is Vector of length rows
66 	alias data this;							/// Treat this matrix as an nxm array
67 	alias cols = numCols;						/// Holds the number of cols;
68 	alias rows = numRows;						/// Holds the number of columns;
69 
70 
71 	unittest {		// Construction through aliased this static array of Vectors
72 
73 		mat2 m2 = mat2( 0.0f, 1.0f, 2.0f, 3.0f );
74 		m2 = mat2( [ 0.0f, 1.0f ], [ 2.0f, 3.0f ] );
75 		m2 = mat2( vec2( 0.0f, 1.0f ), vec2( 2.0f, 3.0f ));
76 		assert( m2[ 0 ][ 0 ] == 0.0f );
77 		assert( m2[ 0 ][ 1 ] == 1.0f );
78 		assert( m2[ 1 ][ 0 ] == 2.0f );
79 		assert( m2[ 1 ][ 1 ] == 3.0f );
80 		
81 		m2[ 0 ] = [ 2.0, 2.0 ];
82 		m2[ 0 ] = vec2( 2.0, 2.0 );
83 		//m2[ 0 .. 1 ] = [ [ 2.0f, 2.0f ] ];
84 		m2[ 0 .. 1 ] = [ vec2( 2.0, 2.0 ) ];
85 
86 		assert( m2 == [ [ 2.0f, 2.0f ], [ 2.0f, 3.0f ] ] );
87 
88 		mat3 m3 = mat3( 0.0f, 0.1f, 0.2f, 1.0f, 1.1f, 1.2f, 2.0f, 2.1f, 2.2f );
89 		assert( m3[ 0 ][ 1 ] == 0.1f );
90 		assert( m3[ 2 ][ 0 ] == 2.0f );
91 		assert( m3[ 1 ][ 2 ] == 1.2f );
92 
93 		m3[ 0 ][ 0 .. $ ] = 0.0f;
94 		assert(  m3 == [[ 0.0f, 0.0f, 0.0f ],
95 						[ 1.0f, 1.1f, 1.2f ], 
96 						[ 2.0f, 2.1f, 2.2f ] ] );
97 		//m3[ 1 .. 3 ]  = [ [ 1, 1, 1 ], [ 2, 2, 2 ] ];
98 		m3[ 1 .. 3 ] =  [ vec3( 1, 1, 1 ), vec3( 2, 2, 2 ) ];
99 		assert(  m3 == [[ 0.0f, 0.0f, 0.0f ],
100 						[ 1.0f, 1.0f, 1.0f ], 
101 						[ 2.0f, 2.0f, 2.0f ]] );
102 
103 		mat4 m4 = mat4( 0.0f, 0.1f, 0.2f, 0.3f,
104 							 1.0f, 1.1f, 1.2f, 1.3f,
105 							 2.0f, 2.1f, 2.2f, 2.3f,
106 							 3.0f, 3.1f, 3.2f, 3.3f );
107 		assert( m4[ 0 ][ 3 ] == 0.3f );
108 		assert( m4[ 1 ][ 1 ] == 1.1f );
109 		assert( m4[ 2 ][ 0 ] == 2.0f );
110 		assert( m4[ 3 ][ 2 ] == 3.2f );
111 
112 		m4[ 2 ][ 1 .. 3 ] = [ 1.0f, 2.0f ];
113 		assert(  m4 == [[ 0.0f, 0.1f, 0.2f, 0.3f ],
114 						[ 1.0f, 1.1f, 1.2f, 1.3f ],
115 						[ 2.0f, 1.0f, 2.0f, 2.3f ],
116 						[ 3.0f, 3.1f, 3.2f, 3.3f ]] );
117 
118 	}
119 
120 
121 	/// Returns the pointer to the stored values as OpenGL requires it.
122 	/// Note this will return a pointer to a $( RED column - major ) matrix, 
123 	/// $( RED this is the OpneGL convention and expected in programs via Uniforms or UBOs ).
124 	@property auto ptr() { return data[ 0 ].ptr; }
125 
126 	/// Returns the current matrix formatted as flat string.
127 	@property string asString() {  return format( "%s", data );  }
128 	alias asString toString;
129 
130 	/// Returns the current matrix as pretty formatted string.
131 	/// TODO : Check This
132 	@property string asPrettyString() {
133 		string fmtr = "%s";
134 
135 		size_t rjust = max( format( fmtr, reduce!( max )( data[] )).length,
136 								  format( fmtr, reduce!( min )( data[] )).length ) - 1;
137 
138 		string[] outer_parts;
139 		foreach( valueType[] col; data ) {
140 			string[] inner_parts;
141 			foreach( valueType row; col ) {
142 				inner_parts ~= rightJustify( format( fmtr, row ), rjust );
143 			}
144 			outer_parts ~= " [ " ~ join( inner_parts, ", " ) ~ " ]";
145 		}
146 
147 		return "[ " ~ join( outer_parts, "\n" )[ 1 .. $ ] ~ " ]";
148 	}
149 	alias asPrettyString toPrettyString;
150 
151 	@safe pure nothrow :
152 	template isCompatibleMatrix( T ) {  enum isCompatibleMatrix = is( typeof( isCompatibleMatrixImpl( T.init )));  }
153 	static void isCompatibleMatrixImpl( int col, int row )( Matrix!( valueType, col, row ) mat ) {}
154 
155 	template isCompatibleVector( T ) {  enum isCompatibleVector = is( typeof( isCompatibleVectorImpl( T.init )));  }
156 	static void isCompatibleVectorImpl( int dim )( Vector!( valueType, dim ) vec ) {}
157 
158 	
159 	/// TODO: Fix Construction with combination of numeric, static and dynamic array, vector and matrix
160 	private void construct( int i, T, Tail... )( T head, Tail tail ) {
161 		static if( i >= cols * rows ) {
162 			static assert( false, "constructor has too many arguments" );
163 		}	else static if( is( T : valueType )) {
164 			data[ i / rows ][ i % rows ] = head;
165 			construct!( i + 1 )( tail );
166 		}	else static if( is( T == Vector!( valueType, rows ))) {
167 			static if( i % rows == 0 ) {
168 				data[ i / rows ] = head;
169 				construct!( i + T.dimension )( tail );
170 			}	else {
171 				static assert( false, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit" );
172 			}
173 		}	/*else {
174 			static assert( false, "Matrix constructor argument must be of type " ~ valueType.stringof ~ " or Vector, not " ~ T.stringof );
175 		}*/
176 	}
177 
178 	private void construct( int i )() {}	// terminate
179 
180 
181 	/// Constructs the matrix:
182 	/// If a single value is passed, the matrix will be cleared with this value ( each column in each col will contain this value ).
183 	/// If a matrix with more cols and columns is passed, the matrix will be the upper left nxm matrix.
184 	/// If a matrix with less cols and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
185 	/// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
186 	/// Examples:
187 	/// - - - 
188 	/// mat2 m2 = mat2( 0.0f ); // mat2 m2 = mat2( 0.0f, 0.0f, 0.0f, 0.0f );
189 	/// mat3 m3 = mat3( m2 ); // mat3 m3 = mat3( 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f );
190 	/// mat3 m3_2 = mat3( vec3( 1.0f, 2.0f, 3.0f ), 4.0f, 5.0f, 6.0f, vec3( 7.0f, 8.0f, 9.0f ));
191 	/// mat4 m4 = mat4.identity; // just an identity matrix
192 	/// mat3 m3_3 = mat3( m4 ); // mat3 m3_3 = mat3.identity
193 	/// - - - 
194 	this( Args... )( Args args ) {  construct!( 0 )( args );  }
195 
196 	/// Construct a Matrix from another Matrix, equal sized or bigger
197 	this( T )( T mat ) if( isMatrix!T && ( T.rows >= rows ) && ( T.cols >= cols )) {
198 		for( int c = 0; c < cols; ++c ) {
199 			for( int r = 0; r < rows; ++r ) {
200 				data[ c ][ r ] = mat.data[ c ][ r ];
201 			}
202 		}
203 	}
204 
205 
206 	/// Construct a Matrix from another Matrix
207 	this( T )( T mat ) if( isMatrix!T ) {
208 		this( 1 );	// Could be optimized with not initializing and setting missing values within static if range checks
209 		//alias Cols = cols < T.cols ? cols : T.cols;
210 		//alias Rows = rows < T.rows ? rows : T.rows;
211 		foreach( c; 0 .. cols < T.cols ? cols : T.cols ) {
212 			foreach( r; 0 .. rows < T.rows ? rows : T.rows ) {
213 				data[ c ][ r ] = mat.data[ c ][ r ];
214 			}
215 		}
216 	}
217 
218 	/// Construct a Matrix from a single value, for non- and square matrices
219 	/// see GLSL 4.5 Spec, section 5.4.2 Vector and Matrix Constructors 
220 	this()( valueType value ) {
221 		clear( 0 );
222 		foreach( rc; 0 .. cols < rows ? cols : rows ) {
223 			data[ rc ][ rc ] = value;
224 		}
225 	}
226 
227 	/// Sets all values of the matrix to value ( each column in each col will contain this value ).
228 	void clear( valueType value ) {
229 		foreach( ref vec; data )
230 			vec[] = value;
231 	}
232 
233 	unittest {		/// Matrix.isvalid, Constructing
234 
235 		mat2 m2 = mat2( 1.0f, 1.0f, vec2( 2.0f, 2.0f ));
236 		assert( m2.data == [[ 1.0f, 1.0f ], [ 2.0f, 2.0f ]] );
237 		m2.clear( 3.0f );
238 		assert( m2.data == [[ 3.0f, 3.0f ], [ 3.0f, 3.0f ]] );
239 		assert( m2.isvalid );
240 		m2.clear( float.nan );
241 		assert( !m2.isvalid );
242 		m2.clear( float.infinity );
243 		assert( !m2.isvalid );
244 		m2.clear( 0.0f );
245 		assert( m2.isvalid );
246 
247 		mat3 m3 = mat3( 1.0f );
248 		assert( m3.data == [[ 1.0f, 0.0f, 0.0f ],
249 							[ 0.0f, 1.0f, 0.0f ],
250 							[ 0.0f, 0.0f, 1.0f ]] );
251 
252 		mat4 m4 = mat4( vec4( 	1.0f, 1.0f, 1.0f, 1.0f ),
253 								2.0f, 2.0f, 2.0f, 2.0f,
254 								3.0f, 3.0f, 3.0f, 3.0f,
255 							 vec4( 4.0f, 4.0f, 4.0f, 4.0f ));
256 		assert( m4.data == [[ 1.0f, 1.0f, 1.0f, 1.0f ],
257 							[ 2.0f, 2.0f, 2.0f, 2.0f ],
258 							[ 3.0f, 3.0f, 3.0f, 3.0f ],
259 							[ 4.0f, 4.0f, 4.0f, 4.0f ]] );
260 		assert( mat3( m4 ).data == [[ 1.0f, 1.0f, 1.0f ],
261 									[ 2.0f, 2.0f, 2.0f ],
262 									[ 3.0f, 3.0f, 3.0f ]] );
263 		assert( mat2( mat3( m4 )).data == [[ 1.0f, 1.0f ],
264 														[ 2.0f, 2.0f ] ] );
265 		assert( mat2( m4 ).data == mat2( mat3( m4 )).data );
266 		assert( mat4( mat3( m4 )).data ==  [[ 1.0f, 1.0f, 1.0f, 0.0f ],
267 											[ 2.0f, 2.0f, 2.0f, 0.0f ],
268 											[ 3.0f, 3.0f, 3.0f, 0.0f ],
269 											[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
270 
271 		Matrix!( float, 2, 3 ) mt1 = Matrix!( float, 2, 3 )( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f );
272 		Matrix!( float, 3, 2 ) mt2 = Matrix!( float, 3, 2 )( 6.0f, - 1.0f, 3.0f, 2.0f, 0.0f, - 3.0f );
273 
274 		assert( mt1.data == [[ 1.0f, 2.0f, 3.0f ], [ 4.0f, 5.0f, 6.0f ]] );
275 		assert( mt2.data == [[ 6.0f, - 1.0f ], [ 3.0f, 2.0f ], [ 0.0f, - 3.0f ]] );
276 	}
277 
278 
279 	//void opAssign( T )( T array )  if( is( T array == U[][], U )) { //if( isArray( A )) {
280 	//	pragma( msg, "Instantiat opAssign" );
281 	//	data[] = array[];
282 	//}
283 
284 	static if( cols == rows ) {
285 		/// static create identity matrix
286 		static @property Matrix identity() {
287 			return Matrix( 1 );
288 		}
289 
290 		/// Returns a identity matrix.
291 		@property bool isIdentity() {
292 			return data == Matrix( 1 ).data;
293 		}
294 
295 		/// Unittest identity, isIdentity
296 		unittest {
297 			mat2 m2 = mat2( 1.0f );
298 			assert( m2.isIdentity );
299 			assert( m2.data == mat2( 1 ).data );
300 
301 			m2 = mat2.identity;
302 			assert( m2.data == [[ 1.0f, 0.0f ],
303 								[ 0.0f, 1.0f ]] );
304 
305 			mat3 m3 = mat3.identity;
306 			assert( m3.isIdentity );
307 			assert( m3.data == mat3( 1.0f ).data );
308 
309 			mat4 m4 = mat4( 2.0f );
310 			assert( m4.data == mat4( 2.0f ).data );
311 
312 			m4 = mat4.identity;
313 			assert( m4.isIdentity );
314 			assert( m4.data == [[ 1.0f, 0.0f, 0.0f, 0.0f ],
315 								[ 0.0f, 1.0f, 0.0f, 0.0f ],
316 								[ 0.0f, 0.0f, 1.0f, 0.0f ],
317 								[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
318 			assert( m4.data == m4.identity.data );
319 		}
320 
321 
322 
323 
324 		/////////////////
325 		// Translation //
326 		/////////////////
327 		static if( cols == 2 ) {			// homogeneous 1D translations
328 			/// static construction of a translation matrix
329 			static Matrix translation( valueType x ) {
330 				Matrix result = Matrix.identity;
331 				result.data[ 1 ][ 0 ] = x;
332 				return result;
333 			}
334 
335 			/// translate an existing matrix
336 			Matrix translate( valueType x ) {
337 				this = Matrix.translation( x ) * this;
338 				return this;
339 			}
340 		}
341 
342 		else static if( cols == 3 ) {		// homogeneous 2D translations
343 			/// static construction of a translation matrix from an array/vector
344 			static Matrix translation( valueType[ 2 ] vec )  {
345 				Matrix result = Matrix.identity;
346 				result.data[ 2 ][ 0..2 ] = vec;
347 				return result;
348 			}
349 
350 			/// static construction of a translation matrix from two coordinates
351 			static Matrix translation( valueType x, valueType y )  {
352 				return translation( [ x, y ] );
353 			}
354 
355 			/// translate an existing matrix with an array/vector
356 			Matrix translate( valueType[ 2 ] vec ) {
357 				this = Matrix.translation( vec ) * this;
358 				return this;
359 			}
360 
361 			/// translate an existing matrix with two scalars
362 			Matrix translate( valueType x, valueType y ) {
363 				this = Matrix.translation( [ x, y ] ) * this;
364 				return this;
365 			}
366 		}
367 
368 		else static if( cols == 4 ) {		// homogeneous 3D translations
369 			/// static construction of a translation matrix from an array/vector
370 			static Matrix translation( valueType[ 3 ] vec )  {
371 				Matrix result = Matrix.identity;
372 				result.data[ 3 ][ 0..3 ] = vec;
373 				return result;
374 			}
375 
376 			/// static construction of a translation matrix from three coordinates
377 			static Matrix translation( valueType x, valueType y, valueType z )  {
378 				return translation( [ x, y, z ] );
379 			}
380 
381 			/// translate an existing matrix with a vector
382 			Matrix translate( valueType[ 3 ] vec ) {
383 				this = Matrix.translation( vec ) * this;
384 				return this;
385 			}
386 
387 			/// translate an existing matrix with three scalars
388 			Matrix translate( valueType x, valueType y, valueType z ) {
389 				this = Matrix.translation( [ x, y, z ] ) * this;
390 				return this;
391 			}
392 
393 			unittest {
394 				mat4 m4 = mat4( 1.0f );
395 				assert( m4.translation( 1.0f, 2.0f, 3.0f ).data == mat4.translation( 1.0f, 2.0f, 3.0f ).data );
396 				assert( mat4.translation( 1.0f, 2.0f, 3.0f ).data  ==  [[ 1.0f, 0.0f, 0.0f, 0.0f ],
397 																		[ 0.0f, 1.0f, 0.0f, 0.0f ],
398 																		[ 0.0f, 0.0f, 1.0f, 0.0f ],
399 																		[ 1.0f, 2.0f, 3.0f, 1.0f ] ] );
400 				assert( mat4.identity.translate( 0.0f, 1.0f, 2.0f ).data == mat4.translation( 0.0f, 1.0f, 2.0f ).data );
401 			}
402 		}
403 
404 		/// extract translation part of the matrix in a copy 
405 		Matrix translation() {
406 			Matrix result = Matrix.identity;
407 			result[ $-1 ][ 0 .. $-1 ] = data[ $-1 ][ 0 .. $-1 ];
408 			return result;
409 		}
410 
411 		unittest  {
412 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
413 							3.0f, 4.0f, 5.0f,
414 							6.0f, 7.0f, 1.0f );
415 			assert( m3.translation.data == [ [ 1.0f, 0.0f, 2.0f ], [ 0.0f, 1.0f, 5.0f ], [ 0.0f, 0.0f, 1.0f ] ] );
416 
417 
418 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
419 							 4.0f,  5.0f,  6.0f,  7.0f,
420 							 8.0f,  9.0f, 10.0f, 11.0f,
421 							12.0f, 13.0f, 14.0f,  1.0f );
422 			assert( m4.translation.data == [ [ 1.0f, 0.0f, 0.0f,  3.0f ],
423 											 [ 0.0f, 1.0f, 0.0f,  7.0f ],
424 											 [ 0.0f, 0.0f, 1.0f, 11.0f ],
425 											 [ 0.0f, 0.0f, 0.0f,  1.0f ] ] );
426 		}
427 
428 
429 
430 		//////////////
431 		// Rotation // - integer rotations might not be useful, but should they be disabled?
432 		//////////////
433 		static if( cols == 2 || cols == 3 ) {		// non- and homogeneous 2D rotations		
434 			/// static construction of a rotation matrix
435 			static Matrix rotation( real angle ) {
436 				Matrix result = Matrix.identity;
437 				valueType cosAngle = to!valueType( cos( angle ));
438 				valueType sinAngle = to!valueType( sin( angle ));
439 				result.data[ 0 ][ 0 ] = cosAngle;
440 				result.data[ 1 ][ 0 ] = - sinAngle;
441 				result.data[ 0 ][ 1 ] = sinAngle;
442 				result.data[ 1 ][ 1 ] = cosAngle;
443 				return result;
444 			}
445 
446 			/// rotate an existing matrix
447 			Matrix rotate( real angle ) {
448 				this = rotation( angle ) * this;
449 				return this;
450 			}
451 		}
452 
453 		static if( cols >= 3 ) {						// non- and homogeneous 3D rotations
454 			/// static construction of a rotation matrix angle and axis
455 			static Matrix rotation( real angle, Vector!( valueType, 3 ) axis ) {
456 				Matrix result = Matrix.identity;
457 
458 				auto length = axis.length;
459 				if( length != 1 ) {
460 					version( NoReciprocalMul ) {
461 						axis /= length;
462 					}	else	{
463 						auto invLength = 1.0 / length;
464 						axis *= invLength;
465 					}
466 				}
467 
468 				const real cosAngle = cos( angle );
469 				const real sinAngle = sin( angle );
470 
471 				Vector!( valueType, 3 ) temp = ( 1 - cosAngle ) * axis;
472 
473 				result.data[ 0 ][ 0 ] = to!valueType( cosAngle	+	temp.x * axis.x );
474 				result.data[ 0 ][ 1 ] = to!valueType(				temp.x * axis.y + sinAngle * axis.z );
475 				result.data[ 0 ][ 2 ] = to!valueType(				temp.x * axis.z - sinAngle * axis.y );
476 				result.data[ 1 ][ 0 ] = to!valueType(				temp.y * axis.x - sinAngle * axis.z );
477 				result.data[ 1 ][ 1 ] = to!valueType( cosAngle	+	temp.y * axis.y );
478 				result.data[ 1 ][ 2 ] = to!valueType(				temp.y * axis.z + sinAngle * axis.x );
479 				result.data[ 2 ][ 0 ] = to!valueType(				temp.z * axis.x + sinAngle * axis.y );
480 				result.data[ 2 ][ 1 ] = to!valueType(				temp.z * axis.y - sinAngle * axis.x );
481 				result.data[ 2 ][ 2 ] = to!valueType( cosAngle	+	temp.z * axis.z );
482 
483 				return result;
484 			}
485 
486 			/// static construction of a rotation matrix angle and axis coordinates
487 			static Matrix rotation( real angle, valueType x, valueType y, valueType z ) {
488 				return Matrix.rotation( angle, Vector!( valueType, 3 )( x, y, z ));
489 			}
490 
491 			/// rotate an existing matrix with an angle around an axis
492 			Matrix rotate( real angle, Vector!( valueType, 3 ) axis ) {
493 				this = rotation( angle, axis ) * this;
494 				return this;
495 			}
496 
497 			/// rotate an existing matrix with an angle around axis coordinates
498 			Matrix rotate( real angle, valueType x, valueType y, valueType z  ) {
499 				this = rotation( angle, x, y, z ) * this;
500 				return this;
501 			}
502 
503 			/// static construction of a rotation matrix with an angle around a canonical axis
504 			private static Matrix rotationA( ubyte a, ubyte b )( real angle ) {
505 				Matrix result = Matrix.identity;
506 
507 				valueType cosAngle = to!valueType( cos( angle ));
508 				valueType sinAngle = to!valueType( sin( angle ));
509 
510 				result.data[ a ][ a ] = cosAngle;
511 				result.data[ b ][ a ] = - sinAngle;
512 				result.data[ a ][ b ] = sinAngle;
513 				result.data[ b ][ b ] = cosAngle;
514 
515 				return result;
516 			}
517 
518 			/// static construction of a rotation matrix with an angle around X axis 
519 			static Matrix rotationX( real angle ) { return rotationA!( 1, 2 )( angle );	} /// A-canonical = X
520 
521 			/// static construction of a rotation matrix with an angle around Y axis
522 			static Matrix rotationY( real angle ) { return rotationA!( 2, 0 )( angle );	} /// A-canonical = Y
523 
524 			/// static construction of a rotation matrix with an angle around Z axis
525 			static Matrix rotationZ( real angle ) { return rotationA!( 0, 1 )( angle );	} /// A-canonical = Z
526 
527 			/// rotate an existing matrix with an angle around X axis
528 			Matrix rotateX( real angle ) {
529 				this = rotationX( angle ) * this;
530 				return this;
531 			}
532 
533 			/// rotate an existing matrix with an angle around Y axis
534 			Matrix rotateY( real angle ) {
535 				this = rotationY( angle ) * this;
536 				return this;
537 			}
538 
539 			/// rotate an existing matrix with an angle around Z axis
540 			Matrix rotateZ( real angle ) {
541 				this = rotationZ( angle ) * this;
542 				return this;
543 			}
544 
545 			unittest {
546 				assert( mat4.rotationX( 0 ).data == [[ 1.0f,   0.0f,   0.0f, 0.0f ],
547 													 [ 0.0f,   1.0f, - 0.0f, 0.0f ],
548 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
549 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
550 				assert( mat4.rotationY( 0 ).data == [[ 1.0f,   0.0f,   0.0f, 0.0f ],
551 													 [ 0.0f,   1.0f,   0.0f, 0.0f ],
552 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
553 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
554 				assert( mat4.rotationZ( 0 ).data == [[ 1.0f, - 0.0f,   0.0f, 0.0f ],
555 													 [ 0.0f,   1.0f,   0.0f, 0.0f ],
556 													 [ 0.0f,   0.0f,   1.0f, 0.0f ],
557 													 [ 0.0f,   0.0f,   0.0f, 1.0f ]] );
558 
559 				//mat4 rotX = mat4.identity;
560 				//rotX.rotateX( 1 );
561 				//assert( mat4.rotationX( 1 ).data == rotX.data );
562 				//assert( rotX.data == mat4.identity.rotateX( 1 ).data );
563 				//assert( rotX.data == mat4.rotation( 1, vec3( 1.0f, 0.0f, 0.0f )).data );
564 
565 				//mat4 rotY = mat4.identity;
566 				//rotY.rotateY( 2 );
567 				//assert( mat4.rotationY( 2 ).data == rotY.data );
568 				//assert( rotY.data == mat4.identity.rotateY( 2 ).data );
569 				//assert( rotY.data == mat4.rotation( 2, vec3( 0.0f, 1.0f, 0.0f )).data );
570 
571 				//mat4 rotZ = mat4.identity;
572 				//rotZ.rotateZ( 3 );
573 				//assert( mat4.rotationZ( 3 ).data == rotZ.data );
574 				//assert( rotZ.data == mat4.identity.rotateZ( 3 ).data );
575 				//assert( rotZ.data == mat4.rotation( 3, vec3( 0.0f, 0.0f, 1.0f )).data );
576 			}
577 		}
578 
579 		/// extract rotation part of the matrix in a copy 
580 		// TODO(pp): fix proper extraction, currently works only right if determinant = 1
581 		Matrix rotation() {
582 			Matrix result = data;
583 			result[ $-1 ][ 0..$-1 ] = 0;
584 			result[ 0..$-1 ][ $-1 ] = 0;
585 			result[ $-1 ][ $-1 ] = 1;
586 			return result;
587 		}
588 
589 		unittest {
590 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
591 							3.0f, 4.0f, 5.0f,
592 							6.0f, 7.0f, 1.0f );
593 			assert( m3.rotation.data == [[ 0.0f, 1.0f, 0.0f ],
594 										 [ 3.0f, 4.0f, 0.0f ],
595 										 [ 0.0f, 0.0f, 1.0f ]] );
596 
597 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
598 							 4.0f,  5.0f,  6.0f,  7.0f,
599 							 8.0f,  9.0f, 10.0f, 11.0f,
600 							12.0f, 13.0f, 14.0f,  1.0f );
601 
602 			assert( m4.rotation.data == [[  0.0f,  1.0f,  2.0f, 0.0f ],
603 										 [  4.0f,  5.0f,  6.0f, 0.0f ],
604 										 [  8.0f,  9.0f, 10.0f, 0.0f ],
605 										 [  0.0f,  0.0f,  0.0f, 1.0f ]] );
606 		}
607 
608 
609 
610 		/////////////
611 		// Scaling //
612 		/////////////
613 		static if( cols == 2 ) {						// homogeneous 1D scaling
614 			/// static construction of a scaling matrix
615 			static Matrix scaling( valueType x ) {
616 				Matrix result = Matrix.identity;
617 				result.data[ 0 ][ 0 ] = x;
618 				return result;
619 			}
620 
621 			/// scale an existing matrix 
622 			Matrix scale( valueType x ) {
623 				this = Matrix.scaling( x ) * this;
624 				return this;
625 			}
626 		}
627 
628 		static if( cols == 2 || cols == 3 ) {		// non- and homogeneous 2D scaling
629 			/// static construction of a scaling matrix from an array/vector
630 			static Matrix scaling( valueType[ 2 ] vec ) {
631 				Matrix result = Matrix.identity;
632 				result.data[ 0 ][ 0 ] = vec[ 0 ];
633 				result.data[ 1 ][ 1 ] = vec[ 1 ];
634 				return result;
635 			}
636 
637 			/// static construction of a scaling matrix from two scalars
638 			static Matrix scaling( valueType x, valueType y ) {
639 				return scaling( [ x, y ] );
640 			}
641 
642 			/// scale an existing matrix with an array/vector 
643 			Matrix scale( valueType[ 2 ] vec ) {
644 				this = Matrix.scaling( vec ) * this;
645 				return this;
646 			}
647 
648 			/// scale an existing matrix with two scalars
649 			Matrix scale( valueType x, valueType y ) {
650 				this = Matrix.scaling( [ x, y ] ) * this;
651 				return this;
652 			}
653 		}
654 
655 		else static if( cols >= 3 ) {
656 			/// static construction of a scaling matrix from an array/vector
657 			static Matrix scaling( valueType[ 3 ] vec ) {
658 				Matrix result = Matrix.identity;
659 				result.data[ 0 ][ 0 ] = vec[ 0 ];
660 				result.data[ 1 ][ 1 ] = vec[ 1 ];
661 				result.data[ 2 ][ 2 ] = vec[ 2 ];
662 				return result;
663 			}
664 
665 			/// static construction of a scaling matrix from three scalars
666 			static Matrix scaling( valueType x, valueType y, valueType z ) {
667 				return scaling( [ x, y, z ] );
668 			}
669 
670 			/// scale an existing matrix with an array/vector
671 			Matrix scale( valueType[ 3 ] vec ) {
672 				this = Matrix.scaling( vec ) * this;
673 				return this;
674 			}
675 
676 			/// scale an existing matrix with three scalars
677 			Matrix scale( valueType x, valueType y, valueType z ) {
678 				this = Matrix.scaling( [ x, y, z ] ) * this;
679 				return this;
680 			}
681 
682 			unittest {
683 				mat4 m4 = mat4( 1.0f );
684 				assert( m4.scaling( 0.0f, 1.0f, 2.0f ).data == mat4.scaling( 0.0f, 1.0f, 2.0f ).data );
685 				assert( mat4.scaling( 0.0f, 1.0f, 2.0f ).data  ==  [[ 0.0f, 0.0f, 0.0f, 0.0f ],
686 																	[ 0.0f, 1.0f, 0.0f, 0.0f ],
687 																	[ 0.0f, 0.0f, 2.0f, 0.0f ],
688 																	[ 0.0f, 0.0f, 0.0f, 1.0f ] ] );
689 				assert( mat4.identity.scale( 0.0f, 1.0f, 2.0f ).data == mat4.scaling( 0.0f, 1.0f, 2.0f ).data );
690 			}
691 		}
692 
693 		/// extract scaling part of the matrix in a copy
694 		Matrix scale() { 
695 			Matrix result = Matrix.identity;
696 			foreach( rc; 0 .. cols - 1 )
697 				result[ rc ][ rc ] = data[ rc ][ rc ];
698 			return result;
699 		}
700 
701 		unittest {
702 			mat3 m3 = mat3( 0.0f, 1.0f, 2.0f,
703 							3.0f, 4.0f, 5.0f,
704 							6.0f, 7.0f, 1.0f );
705 
706 			assert( m3.scale.data  ==  [[ 0.0f, 0.0f, 0.0f ], 
707 										[ 0.0f, 4.0f, 0.0f ],
708 										[ 0.0f, 0.0f, 1.0f ] ] );
709 
710 
711 			mat4 m4 = mat4(  0.0f,  1.0f,  2.0f,  3.0f,
712 							 4.0f,  5.0f,  6.0f,  7.0f,
713 							 8.0f,  9.0f, 10.0f, 11.0f,
714 							12.0f, 13.0f, 14.0f,  1.0f );
715 
716 			assert( m4.scale.data  ==  [[ 0.0f, 0.0f,  0.0f, 0.0f ],
717 										[ 0.0f, 5.0f,  0.0f, 0.0f ],
718 										[ 0.0f, 0.0f, 10.0f, 0.0f ],
719 										[ 0.0f, 0.0f,  0.0f, 1.0f ]] );
720 		}
721 	}	// end static if( cols == rows )
722 
723 
724 	/// Component-wise binary matrix-scalar operation: addition, subtraction, multiplication, division
725 	auto opBinary( string op, T )( T s ) const if( isNumeric!T && (( op == "+" ) || ( op == "-" ) || ( op == "*" ) || ( op == "/" ))) {
726 		Matrix result;
727 		mixin( "result[0] = data[0] " ~ op ~ " s;" );
728 		mixin( "result[1] = data[1] " ~ op ~ " s;" );
729 		static if( cols >= 3 )  mixin( "result[2] = data[2] " ~ op ~ " s;" );
730 		static if( cols == 4 )  mixin( "result[3] = data[3] " ~ op ~ " s;" );
731 		return result;
732 	}
733 
734 	/// Component-wise binary scalar-matrix operation: addition, subtraction, multiplication, division
735 	auto opBinaryRight( string op, T )( T s ) const if( isNumeric!T && (( op == "+" ) || ( op == "-" ) || ( op == "*" ) || ( op == "/" ))) {
736 		Matrix result;
737 		mixin( "result[0] = s " ~ op ~ "data[0];" );
738 		mixin( "result[1] = s " ~ op ~ "data[1];" );
739 		static if( cols >= 3 )  mixin( "result[2] = s " ~ op ~ "data[2];" );
740 		static if( cols == 4 )  mixin( "result[3] = s " ~ op ~ "data[3];" );
741 		return result;
742 	}
743 
744 	unittest {		/// Matrix-Scalar Operations
745 		mat2 m2 = mat2( 1.0f, 2.0f, 3.0f, 4.0f );
746 		assert(( m2 * 2.0f ).data == [ [ 2.0f, 4.0f ], [ 6.0f, 8.0f ] ] );
747 		assert(( 2 * m2 ).data == ( m2 * 2 ).data );
748 
749 		mat3 m3 = mat3( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f );
750 		assert(( m3 * 2 ).data == [ [ 2.0f, 4.0f, 6.0f ], [ 8.0f, 10.0f, 12.0f ], [ 14.0f, 16.0f, 18.0f ] ] );
751 		assert(( 2 * m3 ).data == ( m3 * 2 ).data );
752 
753 		//TODO: tests for mat4, mat34
754 	}
755 
756 
757 	/// Matrix-vector multiplication
758 	auto opBinary( string op : "*", V )( V vec ) const if( isCompatibleVector!V && V.dimension == cols ) {
759 		vectorType result;
760 		result.clear( 0 );
761 		foreach( c; 0 .. cols )
762 			foreach( r; 0 .. rows )
763 				result[ r ] += data[ c ][ r ] * vec[ c ];
764 		return result;
765 	}
766 
767 	/// Vector-Matrix multiplication
768 	auto opBinaryRight( string op : "*", V : vectorType )( V vec ) const if( isVector!V && V.dimension == rows ) {
769 		auto transposedMatrix = transposed();
770 		return transposedMatrix * vec;
771 	}
772 
773 	unittest {		/// matrix-vector, vector-matrix multiplication
774 		mat2 m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
775 		vec2 v2 = vec2( 2.0f, 2.0f );
776 		assert(( m2 * v2 ) == [ 16.0f, 24.0f ] );
777 		assert(( v2 * m2 ) == [ 12.0f, 28.0f ] );
778 
779 		mat3 m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
780 		vec3 v3 = vec3( 2.0f, 2.0f, 2.0f );
781 		assert(( m3 * v3 ) == [ 48.0f, 60.0f, 72.0f ] );
782 		assert(( v3 * m3 ) == [ 24.0f, 60.0f, 96.0f ] );
783 
784 		mat4 m4 = mat4( 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32 );
785 		vec4 v4 = vec4( 2, 2, 2, 2 );
786 		assert(( m4 * v4 ) == [ 112, 128, 144, 160 ] );
787 		assert(( v4 * m4 ) == [  40, 104, 168, 232 ] );
788 
789 		//TODO: tests for mat4, mat34, mat43, mat24, mat 42 
790 	}
791 
792 
793 	/// matrix-matrix component-wise operations addition, subtraction and division, using vector-vector operations of all colums
794 	Matrix opBinary( string op )( Matrix mat ) const if( op == "+" || op == "-" || op == "/" ) {
795 		Matrix result;
796 		mixin( "result[0] = data[0] " ~ op ~ "ma t[0];" );
797 		mixin( "result[1] = data[1] " ~ op ~ "ma t[1];" );
798 		static if( cols >= 3 )  mixin( "result[2] = data[2] " ~ op ~ "ma t[2];" );
799 		static if( cols == 4 )  mixin( "result[3] = data[3] " ~ op ~ "ma t[3];" );
800 		return result;
801 	}
802 
803 	/// matrix-matrix multiplication, using matrix-vector multiplication for each column of mat
804 	Matrix opBinary( string op : "*", M )( M mat ) const if( isCompatibleMatrix!M && cols == M.rows ) {
805 		Matrix!( valueType, cols, M.cols ) result;
806 		result[0] = this * mat[0];
807 		result[1] = this * mat[1];
808 		static if( M.cols >= 3 )  result[2] = this * mat[2];
809 		static if( M.cols == 4 )  result[3] = this * mat[3];
810 		return result;
811 	}
812 
813 
814 	unittest {		/// Matrix-Matrix Operations
815 		mat2 m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
816 		assert(( m2 * m2 ) == [ [ 28.0f, 40.0f ], [ 60.0f, 88.0f ] ] );
817 		assert(( m2 - m2 ) == [ [ 0.0f, 0.0f ], [ 0.0f, 0.0f ] ] );
818 		mat2 m2_2 = m2 + m2;
819 		assert(  m2_2  == [ [ 4.0f, 8.0f ], [ 12.0f, 16.0f ] ] );
820 		assert(( m2_2 / m2 ) == [ vec2( 2.0f, 2.0f ), vec2( 2.0f, 2.0f ) ] );
821 
822 		mat3 m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
823 		assert(( m3 * m3 ) == [ [ 120.0f, 144.0f, 168.0f ], [ 264.0f, 324.0f, 384.0f ], [ 408.0f, 504.0f, 600.0f ] ] );
824 		assert(( m3 - m3 ) == [ [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ] ] );
825 		assert(( m3 + m3 ) == [ [ 4.0f, 8.0f, 12.0f ], [ 16.0f, 20.0f, 24.0f ], [ 28.0f, 32.0f, 36.0f ] ] );
826 		assert(( m3 / m3 ) == [ [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ], vec3( 1.0f, 1.0f, 1.0f ) ] );
827 
828 		//TODO: tests for mat4, mat34
829 	}
830 
831 
832 	void opOpAssign( string op, T )( T val ) {
833 		mixin( "this = this " ~ op ~  "val;" );
834 	}
835 
836 	unittest {		/// Matrix Unary Operations
837 		mat2  m2 = mat2( 2.0f, 4.0f, 6.0f, 8.0f );
838 		m2 += m2; assert( m2 == [ [ 4.0f, 8.0f ], [ 12.0f, 16.0f ] ] );
839 		m2 /= m2; assert( m2 == [ [ 1.0f, 1.0f ], [ 1.0f, 1.0f ] ] );
840 		m2 -= m2; assert( m2 == [ [ 0.0f, 0.0f ], [ 0.0f, 0.0f ] ] );
841 
842 		mat3  m3 = mat3( 2.0f, 4.0f, 6.0f, 8.0f, 10.0f, 12.0f, 14.0f, 16.0f, 18.0f );
843 		m3 += m3; assert( m3 == [ [ 4.0f, 8.0f, 12.0f ], [ 16.0f, 20.0f, 24.0f ], [ 28.0f, 32.0f, 36.0f ] ] );
844 		m3 /= m3; assert( m3 == [ [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ], [ 1.0f, 1.0f, 1.0f ] ] );
845 		m3 -= m3; assert( m3 == [ [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ], [ 0.0f, 0.0f, 0.0f ] ] );
846 
847 		//TODO: tests for mat4, mat34
848 	}
849 
850 
851 	bool opCast( T : bool )() const {
852 		return isvalid;
853 	}
854 
855 	unittest {
856 		assert( mat2( 1.0f, 2.0f, 1.0f, 1.0f ) == mat2( 1.0f, 2.0f, 1.0f, 1.0f ));
857 		assert( mat2( 1.0f, 2.0f, 1.0f, 1.0f ) != mat2( 1.0f, 1.0f, 1.0f, 1.0f ));
858 
859 		assert( mat3( 1.0f ) == mat3( 1.0f ));
860 		assert( mat3( 1.0f ) != mat3( 2.0f ));
861 
862 		assert( mat4( 1.0f ) == mat4( 1.0f ));
863 		assert( mat4( 1.0f ) != mat4( 2.0f ));
864 
865 		assert( !( mat4( float.nan )));
866 		if( mat4( 1.0f )) { }
867 		else { assert( false ); }
868 	}
869 
870 }
871 
872 
873 /////////////////
874 /// Transpose ///
875 /////////////////
876 auto transpose( Matrix )( Matrix mat ) if( isMatrix!Matrix ) {
877 	Matrix!( Matrix.valueType, Matrix.rows, Matrix.cols ) result;
878 
879 	foreach( c; 0 .. Matrix.cols ) {
880 		foreach( r; 0 .. Matrix.rows ) {
881 			result[ r ][ c ] = mat[ c ][ r ];
882 		}
883 	}
884 
885 	return result;
886 }
887 
888 /// Unittest transpose 
889 unittest {
890 	mat2 m2 = mat2( 1.0f );
891 	assert( m2.transpose.data == mat2( 1.0f ).data );
892 
893 	assert( m2.transpose.data  ==  [[ 1.0f, 0.0f ],
894 									[ 0.0f, 1.0f ]] );
895 	assert( m2.data == mat2.identity.data );
896 
897 	mat3 m3 = mat3( 1.1f, 1.2f, 1.3f,
898 					2.1f, 2.2f, 2.3f,
899 					3.1f, 3.2f, 3.3f );
900 	m3 = m3.transpose;
901 	assert( m3.data == [[ 1.1f, 2.1f, 3.1f ],
902 						[ 1.2f, 2.2f, 3.2f ],
903 						[ 1.3f, 2.3f, 3.3f ]] );
904 
905 	mat4 m4 = mat4( 2.0f );
906 	assert( m4.transpose.data == mat4( 2.0f ).data );
907 }
908 
909 
910 
911 
912 /////////////////
913 // Determinant //
914 /////////////////
915 auto determinant( Matrix )( Matrix mat ) if( isMatrix!Matrix && Matrix.cols == Matrix.rows ) {
916 	static if( mat.cols == 2 ) {
917 		return ( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] );
918 	}
919 
920 	else static if( mat.cols == 3 ) {
921 		return (  mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ]
922 				+ mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ]
923 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ]
924 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ]
925 				- mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ]
926 				- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] );
927 	}
928 
929 	else {
930 		return (  mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
931 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ]
932 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ]
933 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
934 				+ mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ]
935 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ]
936 				+ mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ]
937 				- mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ]
938 				+ mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ]
939 				- mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ]
940 				+ mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ]
941 				- mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] );
942 	}
943 }
944 
945 ////////////
946 // Invert //
947 ////////////
948 auto invert( Matrix )( Matrix mat ) if( isMatrix!Matrix && Matrix.cols == Matrix.rows ) {
949 
950 	static if( isFloatingPoint!( Matrix.valueType ) && rmul ) {
951 		Matrix.valueType d = 1 / mat.determinant;
952 		enum op = "*";
953 	}	else {
954 		Matrix.valueType d = mat.determinant;
955 		enum op = "/";
956 	}
957 
958 	Matrix result;
959 	alias vectorType = Matrix.vectorType;
960 
961 	static if( Matrix.cols == 2 )
962 		mixin( "
963 			result = [
964 				vectorType(   mat[ 1 ][ 1 ], - mat[ 0 ][ 1 ] ) " ~ op ~ " d,
965 				vectorType( - mat[ 1 ][ 0 ],   mat[ 0 ][ 0 ] ) " ~ op ~ " d ];" );
966 	
967 	else static if( Matrix.cols == 3 )
968 		mixin( "
969 			result = [
970 				vectorType((  mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d,
971 							( mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d,
972 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] ) " ~ op ~ " d ),
973 				vectorType((  mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] - mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d,
974 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d,
975 							( mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] - mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] ) " ~ op ~ " d ),
976 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d,
977 							( mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] - mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d,
978 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] ) " ~ op ~ " d ) ];" );
979 
980 	else
981 		mixin( "
982 			result = [
983 				vectorType((  mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ]
984 							- mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
985 							( mat[ 0 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ]
986 							- mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
987 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 2 ]
988 							- mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
989 							( mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ]
990 							- mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] ) " ~ op ~ " d ),
991 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ]
992 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
993 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ]
994 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
995 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 0 ]
996 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 2 ] ) " ~ op ~ " d,
997 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ]
998 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d ),
999 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
1000 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1001 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
1002 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 3 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1003 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 3 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 1 ]
1004 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 3 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1005 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 3 ] + mat[ 0 ][ 3 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ]
1006 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 3 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 3 ] * mat[ 2 ][ 0 ] - mat[ 0 ][ 3 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] ) " ~ op ~ " d ),
1007 				vectorType((  mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] + mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] + mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ]
1008 							- mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] - mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] - mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1009 							( mat[ 0 ][ 0 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 1 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 0 ] + mat[ 0 ][ 2 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 1 ]
1010 							- mat[ 0 ][ 0 ] * mat[ 2 ][ 2 ] * mat[ 3 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 2 ][ 0 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 2 ][ 1 ] * mat[ 3 ][ 0 ] ) " ~ op ~ " d,
1011 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 1 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 2 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 0 ]
1012 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 3 ][ 2 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 3 ][ 0 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 3 ][ 1 ] ) " ~ op ~ " d,
1013 							( mat[ 0 ][ 0 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 2 ] + mat[ 0 ][ 1 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 0 ] + mat[ 0 ][ 2 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 1 ]
1014 							- mat[ 0 ][ 0 ] * mat[ 1 ][ 2 ] * mat[ 2 ][ 1 ] - mat[ 0 ][ 1 ] * mat[ 1 ][ 0 ] * mat[ 2 ][ 2 ] - mat[ 0 ][ 2 ] * mat[ 1 ][ 1 ] * mat[ 2 ][ 0 ] ) " ~ op ~ " d ) ];" );	
1015 
1016 	return result;
1017 }
1018 
1019 
1020 unittest {
1021 	mat2 m2 = mat2( 1.0f, 2.0f, vec2( 3.0f, 4.0f ));
1022 	assert( m2.determinant == - 2.0f );
1023 	assert( m2.invert.data ==  [[ - 2.0f,   1.0f ],
1024 								[   1.5f, - 0.5f ]  );
1025 
1026 	mat3 m3 = mat3( 1.0f, - 2.0f,   3.0f,
1027 					7.0f, - 1.0f,   0.0f,
1028 					3.0f,   2.0f, - 4.0f );
1029 	assert( m3.determinant == - 1.0f );
1030 	assert( m3.invert.data ==  [[ -  4.0f,  2.0f, -  3.0f ],
1031 								[ - 28.0f, 13.0f, - 21.0f ],
1032 								[ - 17.0f,  8.0f, - 13.0f ]] );
1033 
1034 	mat4 m4 = mat4(	  1.0f,   2.0f,   3.0f,   4.0f,
1035 					- 2.0f,   1.0f,   5.0f, - 2.0f,
1036 					  2.0f, - 1.0f,   7.0f,   1.0f,
1037 					  3.0f, - 3.0f,   2.0f,   0.0f );
1038 	assert( m4.determinant == - 8.0f );
1039 	assert( m4.invert.data ==  [[   6.875f,   7.875f, - 11.75f,  11.125f ],
1040 								[   6.625f,   7.625f, - 11.25f,  10.375f ],
1041 								[ - 0.375f, - 0.375f,    0.75f, - 0.625f ],
1042 								[ - 4.5f,   - 5.5f,      8.0f,  - 7.5f ]] );
1043 }
1044 
1045 
1046 /// query if any entry is nan
1047 bool isnan( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1048 	foreach( const ref vec; mat )
1049 		if( vec.isnan )
1050 			return true;
1051 	return false;
1052 }
1053 
1054 
1055 /// query if any entry is inf
1056 bool isinf( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1057 	foreach( const ref vec; mat )
1058 		if( vec.isinf )
1059 			return true;
1060 	return false;
1061 }
1062 
1063 
1064 /// query if all entries are not nan and not inf
1065 bool isvalid( genType )( const ref genType mat ) if( isMatrix!genType && isFloatingPoint!( genType.valueType )) {
1066 	return !( mat.isinf || mat.isnan );
1067 }
1068 
1069 
1070 unittest {		/// Matrix.isvalid
1071 
1072 	mat2 m2 = mat2( 1.0f, 1.0f, vec2( 2.0f, 2.0f ));
1073 	assert( m2.data == [[ 1.0f, 1.0f ], [ 2.0f, 2.0f ]] );
1074 	m2.clear( 3.0f );
1075 	assert( m2.data == [[ 3.0f, 3.0f ], [ 3.0f, 3.0f ]] );
1076 	assert( m2.isvalid );
1077 	m2.clear( float.nan );
1078 	assert( !m2.isvalid );
1079 	m2.clear( float.infinity );
1080 	assert( !m2.isvalid );
1081 	m2.clear( 0.0f );
1082 	assert( m2.isvalid );
1083 
1084 	mat3 m3 = mat3( 1.0f );
1085 	assert( m3.data == [[ 1.0f, 0.0f, 0.0f ],
1086 						[ 0.0f, 1.0f, 0.0f ],
1087 						[ 0.0f, 0.0f, 1.0f ]] );
1088 
1089 	mat4 m4 = mat4( vec4( 	1.0f, 1.0f, 1.0f, 1.0f ),
1090 							2.0f, 2.0f, 2.0f, 2.0f,
1091 							3.0f, 3.0f, 3.0f, 3.0f,
1092 						 vec4( 4.0f, 4.0f, 4.0f, 4.0f ));
1093 	assert( m4.data == [[ 1.0f, 1.0f, 1.0f, 1.0f ],
1094 						[ 2.0f, 2.0f, 2.0f, 2.0f ],
1095 						[ 3.0f, 3.0f, 3.0f, 3.0f ],
1096 						[ 4.0f, 4.0f, 4.0f, 4.0f ]] );
1097 	assert( mat3( m4 ).data == [[ 1.0f, 1.0f, 1.0f ],
1098 								[ 2.0f, 2.0f, 2.0f ],
1099 								[ 3.0f, 3.0f, 3.0f ]] );
1100 	assert( mat2( mat3( m4 )).data == [[ 1.0f, 1.0f ],
1101 													[ 2.0f, 2.0f ] ] );
1102 	assert( mat2( m4 ).data == mat2( mat3( m4 )).data );
1103 	assert( mat4( mat3( m4 )).data ==  [[ 1.0f, 1.0f, 1.0f, 0.0f ],
1104 										[ 2.0f, 2.0f, 2.0f, 0.0f ],
1105 										[ 3.0f, 3.0f, 3.0f, 0.0f ],
1106 										[ 0.0f, 0.0f, 0.0f, 1.0f ]] );
1107 
1108 	Matrix!( float, 2, 3 ) mt1 = Matrix!( float, 2, 3 )( 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f );
1109 	Matrix!( float, 3, 2 ) mt2 = Matrix!( float, 3, 2 )( 6.0f, - 1.0f, 3.0f, 2.0f, 0.0f, - 3.0f );
1110 
1111 	assert( mt1.data == [[ 1.0f, 2.0f, 3.0f ], [ 4.0f, 5.0f, 6.0f ]] );
1112 	assert( mt2.data == [[ 6.0f, - 1.0f ], [ 3.0f, 2.0f ], [ 0.0f, - 3.0f ]] );
1113 }	
1114 
1115 
1116 
1117