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