1 module gfm.math.matrix;
2 
3 import std.math,
4        std.typetuple,
5        std.traits,
6        std..string,
7        std.typecons,
8        std.conv;
9 
10 import gfm.math.vector,
11        gfm.math.shapes,
12        gfm.math.quaternion;
13 
14 /// Generic non-resizeable matrix with R rows and C columns.
15 /// Intended for 3D use (size 3x3 and 4x4).
16 /// Important: <b>Matrices here are in row-major order whereas OpenGL is column-major.</b>
17 /// Params:
18 ///   T = type of elements
19 ///   R = number of rows
20 ///   C = number of columns
21 align(1) struct Matrix(T, size_t R, size_t C)
22 {
23     align(1):
24     public
25     {
26         static assert(R >= 1u && C >= 1u);
27 
28         alias Vector!(T, C) row_t;
29         alias Vector!(T, R) column_t;
30 
31         enum bool isSquare = (R == C);
32 
33         // fields definition
34         union
35         {
36             T[C*R] v;        // all elements
37             row_t[R] rows;   // all rows
38             T[C][R] c;       // components
39         }
40 
41         this(U...)(U values) pure nothrow @nogc
42         {
43             static if ((U.length == C*R) && allSatisfy!(isTConvertible, U))
44             {
45                 // construct with components
46                 foreach(int i, x; values)
47                     v[i] = x;
48             }
49             else static if ((U.length == 1) && (isAssignable!(U[0])) && (!is(U[0] : Matrix)))
50             {
51                 // construct with assignment
52                 opAssign!(U[0])(values[0]);
53             }
54             else static assert(false, "cannot create a matrix from given arguments");
55         }
56 
57         /// Construct a matrix from columns.
58         static Matrix fromColumns(column_t[] columns) pure nothrow @nogc
59         {
60             assert(columns.length == C);
61             Matrix res;
62             for (size_t i = 0; i < R; ++i)
63                 for (size_t j = 0; j < C; ++j)
64                 {
65                    res.c[i][j] = columns[j][i];
66                 }
67             return res;
68         }
69 
70         /// Construct a matrix from rows.
71         static Matrix fromRows(row_t[] rows) pure nothrow @nogc
72         {
73             assert(rows.length == R);
74             Matrix res;
75             res.rows[] = rows[];
76             return res;
77         }
78 
79         /// Construct matrix with a scalar.
80         this(U)(T x) pure nothrow @nogc
81         {
82             for (size_t i = 0; i < _N; ++i)
83                 v[i] = x;
84         }
85 
86         /// Assign with a samey matrice.
87         ref Matrix opAssign(U : Matrix)(U x) pure nothrow @nogc
88         {
89             for (size_t i = 0; i < R * C; ++i)
90                 v[i] = x.v[i];
91             return this;
92         }
93 
94         /// Assign from other small matrices (same size, compatible type).
95         ref Matrix opAssign(U)(U x) pure nothrow @nogc
96             if (is(typeof(U._isMatrix))
97                 && is(U._T : _T)
98                 && (!is(U: Matrix))
99                 && (U._R == R) && (U._C == C))
100         {
101             for (size_t i = 0; i < R * C; ++i)
102                 v[i] = x.v[i];
103             return this;
104         }
105 
106         /// Assign with a static array of size R * C.
107         ref Matrix opAssign(U)(U x) pure nothrow @nogc
108             if ((isStaticArray!U)
109                 && is(typeof(x[0]) : T)
110                 && (U.length == R * C))
111         {
112             for (size_t i = 0; i < R * C; ++i)
113                 v[i] = x[i];
114             return this;
115         }
116 
117         /// Assign with a dynamic array of size R * C.
118         ref Matrix opAssign(U)(U x) pure nothrow @nogc
119             if ((isDynamicArray!U)
120                 && is(typeof(x[0]) : T))
121         {
122             assert(x.length == R * C);
123             for (size_t i = 0; i < R * C; ++i)
124                 v[i] = x[i];
125             return this;
126         }
127 
128         /// Return a pointer to content.
129         inout(T)* ptr() pure inout nothrow @property @nogc
130         {
131             return v.ptr;
132         }
133 
134         /// Returns: column j as a vector.
135         column_t column(size_t j) pure const nothrow @nogc
136         {
137             column_t res = void;
138             for (size_t i = 0; i < R; ++i)
139                 res.v[i] = c[i][j];
140             return res;
141         }
142 
143         /// Returns: row i as a vector.
144         row_t row(size_t i) pure const nothrow @nogc
145         {
146             return rows[i];
147         }
148 
149         /// Covnerts to pretty string.
150         string toString() const nothrow
151         {
152             try
153                 return format("%s", v);
154             catch (Exception e)
155                 assert(false); // should not happen since format is right
156         }
157 
158         /// Matrix * vector multiplication.
159         column_t opBinary(string op)(row_t x) pure const nothrow @nogc if (op == "*")
160         {
161             column_t res = void;
162             for (size_t i = 0; i < R; ++i)
163             {
164                 T sum = 0;
165                 for (size_t j = 0; j < C; ++j)
166                 {
167                     sum += c[i][j] * x.v[j];
168                 }
169                 res.v[i] = sum;
170             }
171             return res;
172         }
173 
174         /// Matrix * matrix multiplication.
175         auto opBinary(string op, U)(U x) pure const nothrow @nogc
176             if (is(typeof(U._isMatrix)) && (U._R == C) && (op == "*"))
177         {
178             Matrix!(T, R, U._C) result = void;
179 
180             for (size_t i = 0; i < R; ++i)
181             {
182                 for (size_t j = 0; j < U._C; ++j)
183                 {
184                     T sum = 0;
185                     for (size_t k = 0; k < C; ++k)
186                         sum += c[i][k] * x.c[k][j];
187                     result.c[i][j] = sum;
188                 }
189             }
190             return result;
191         }
192 
193         /// Matrix add and substraction.
194         Matrix opBinary(string op, U)(U other) pure const nothrow @nogc
195             if (is(U : Matrix) && (op == "+" || op == "-"))
196         {
197             Matrix result = void;
198 
199             for (size_t i = 0; i < R; ++i)
200             {
201                 for (size_t j = 0; j < C; ++j)
202                 {
203                     mixin("result.c[i][j] = c[i][j] " ~ op ~ " other.c[i][j];");
204                 }
205             }
206             return result;
207         }
208 
209         /// Assignment operator with another samey matrix.
210         ref Matrix opOpAssign(string op, U)(U operand) pure nothrow @nogc if (is(U : Matrix))
211         {
212             mixin("Matrix result = this " ~ op ~ " operand;");
213             return opAssign!Matrix(result);
214         }
215 
216         /// Assignment operator with another samey matrix.
217         ref Matrix opOpAssign(string op, U)(U operand) pure nothrow @nogc if (isConvertible!U)
218         {
219             Matrix conv = operand;
220             return opOpAssign!op(conv);
221         }
222 
223         /// Cast to other matrix types.
224         /// If the size are different, the resulting matrix is truncated
225         /// and/or filled with identity coefficients.
226         U opCast(U)() pure const nothrow @nogc if (is(typeof(U._isMatrix)))
227         {
228             U res = U.identity();
229             enum minR = R < U._R ? R : U._R;
230             enum minC = C < U._C ? C : U._C;
231             for (size_t i = 0; i < minR; ++i)
232                 for (size_t j = 0; j < minC; ++j)
233                 {
234                     res.c[i][j] = cast(U._T)(c[i][j]);
235                 }
236             return res;
237         }
238 
239         bool opEquals(U)(U other) pure const nothrow @nogc if (is(U : Matrix))
240         {
241             for (size_t i = 0; i < R * C; ++i)
242                 if (v[i] != other.v[i])
243                     return false;
244             return true;
245         }
246 
247         bool opEquals(U)(U other) pure const nothrow @nogc
248             if ((isAssignable!U) && (!is(U: Matrix)))
249         {
250             Matrix conv = other;
251             return opEquals(conv);
252         }
253 
254         // +matrix, -matrix, ~matrix, !matrix
255         Matrix opUnary(string op)() pure const nothrow @nogc if (op == "+" || op == "-" || op == "~" || op == "!")
256         {
257             Matrix res = void;
258             for (size_t i = 0; i < N; ++i)
259                 mixin("res.v[i] = " ~ op ~ "v[i];");
260             return res;
261         }
262 
263         /// Convert 3x3 rotation matrix to quaternion.
264         /// See_also: 3D Math Primer for Graphics and Game Development.
265         U opCast(U)() pure const nothrow @nogc if (is(typeof(U._isQuaternion))
266                                                    && is(U._T : _T)
267                                                    && (_R == 3) && (_C == 3))
268         {
269             T fourXSquaredMinus1 = c[0][0] - c[1][1] - c[2][2];
270             T fourYSquaredMinus1 = c[1][1] - c[0][0] - c[2][2];
271             T fourZSquaredMinus1 = c[2][2] - c[0][0] - c[1][1];
272             T fourWSquaredMinus1 = c[0][0] + c[1][1] + c[2][2];
273 
274             int biggestIndex = 0;
275             T fourBiggestSquaredMinus1 = fourWSquaredMinus1;
276 
277             if(fourXSquaredMinus1 > fourBiggestSquaredMinus1)
278             {
279                 fourBiggestSquaredMinus1 = fourXSquaredMinus1;
280                 biggestIndex = 1;
281             }
282 
283             if(fourYSquaredMinus1 > fourBiggestSquaredMinus1)
284             {
285                 fourBiggestSquaredMinus1 = fourYSquaredMinus1;
286                 biggestIndex = 2;
287             }
288 
289             if(fourZSquaredMinus1 > fourBiggestSquaredMinus1)
290             {
291                 fourBiggestSquaredMinus1 = fourZSquaredMinus1;
292                 biggestIndex = 3;
293             }
294 
295             T biggestVal = sqrt(fourBiggestSquaredMinus1 + 1) / 2;
296             T mult = 1 / (biggestVal * 4);
297 
298             U quat;
299             switch(biggestIndex)
300             {
301                 case 1:
302                     quat.w = (c[1][2] - c[2][1]) * mult;
303                     quat.x = biggestVal;
304                     quat.y = (c[0][1] + c[1][0]) * mult;
305                     quat.z = (c[2][0] + c[0][2]) * mult;
306                     break;
307 
308                 case 2:
309                     quat.w = (c[2][0] - c[0][2]) * mult;
310                     quat.x = (c[0][1] + c[1][0]) * mult;
311                     quat.y = biggestVal;
312                     quat.z = (c[1][2] + c[2][1]) * mult;
313                     break;
314 
315                 case 3:
316                     quat.w = (c[0][1] - c[1][0]) * mult;
317                     quat.x = (c[2][0] + c[0][2]) * mult;
318                     quat.y = (c[1][2] + c[2][1]) * mult;
319                     quat.z = biggestVal;
320                     break;
321 
322                 default: // biggestIndex == 0
323                     quat.w = biggestVal;
324                     quat.x = (c[1][2] - c[2][1]) * mult;
325                     quat.y = (c[2][0] - c[0][2]) * mult;
326                     quat.z = (c[0][1] - c[1][0]) * mult;
327                     break;
328             }
329 
330             return quat;
331         }
332 
333         /// Converts a 4x4 rotation matrix to quaternion.
334         U opCast(U)() pure const nothrow @nogc if (is(typeof(U._isQuaternion))
335                                                   && is(U._T : _T)
336                                                   && (_R == 4) && (_C == 4))
337         {
338             auto m3 = cast(mat3!T)(this);
339             return cast(U)(m3);
340         }
341 
342         /// Matrix inversion is provided for 2x2, 3x3 and 4x4 floating point matrices.
343 
344         static if (isSquare && isFloatingPoint!T && R == 2)
345         {
346             /// Returns: inverse of matrix.
347             Matrix inverse() pure const nothrow @nogc
348             {
349                 T invDet = 1 / (c[0][0] * c[1][1] - c[0][1] * c[1][0]);
350                 return Matrix( c[1][1] * invDet, -c[0][1] * invDet,
351                                    -c[1][0] * invDet,  c[0][0] * invDet);
352             }
353         }
354 
355         static if (isSquare && isFloatingPoint!T && R == 3)
356         {
357             /// Returns: inverse of matrix.
358             Matrix inverse() pure const nothrow @nogc
359             {
360                 T det = c[0][0] * (c[1][1] * c[2][2] - c[2][1] * c[1][2])
361                       - c[0][1] * (c[1][0] * c[2][2] - c[1][2] * c[2][0])
362                       + c[0][2] * (c[1][0] * c[2][1] - c[1][1] * c[2][0]);
363                 T invDet = 1 / det;
364 
365                 Matrix res = void;
366                 res.c[0][0] =  (c[1][1] * c[2][2] - c[2][1] * c[1][2]) * invDet;
367                 res.c[0][1] = -(c[0][1] * c[2][2] - c[0][2] * c[2][1]) * invDet;
368                 res.c[0][2] =  (c[0][1] * c[1][2] - c[0][2] * c[1][1]) * invDet;
369                 res.c[1][0] = -(c[1][0] * c[2][2] - c[1][2] * c[2][0]) * invDet;
370                 res.c[1][1] =  (c[0][0] * c[2][2] - c[0][2] * c[2][0]) * invDet;
371                 res.c[1][2] = -(c[0][0] * c[1][2] - c[1][0] * c[0][2]) * invDet;
372                 res.c[2][0] =  (c[1][0] * c[2][1] - c[2][0] * c[1][1]) * invDet;
373                 res.c[2][1] = -(c[0][0] * c[2][1] - c[2][0] * c[0][1]) * invDet;
374                 res.c[2][2] =  (c[0][0] * c[1][1] - c[1][0] * c[0][1]) * invDet;
375                 return res;
376             }
377         }
378 
379         static if (isSquare && isFloatingPoint!T && R == 4)
380         {
381             /// Returns: inverse of matrix.
382             Matrix inverse() pure const nothrow @nogc
383             {
384                 T det2_01_01 = c[0][0] * c[1][1] - c[0][1] * c[1][0];
385                 T det2_01_02 = c[0][0] * c[1][2] - c[0][2] * c[1][0];
386                 T det2_01_03 = c[0][0] * c[1][3] - c[0][3] * c[1][0];
387                 T det2_01_12 = c[0][1] * c[1][2] - c[0][2] * c[1][1];
388                 T det2_01_13 = c[0][1] * c[1][3] - c[0][3] * c[1][1];
389                 T det2_01_23 = c[0][2] * c[1][3] - c[0][3] * c[1][2];
390 
391                 T det3_201_012 = c[2][0] * det2_01_12 - c[2][1] * det2_01_02 + c[2][2] * det2_01_01;
392                 T det3_201_013 = c[2][0] * det2_01_13 - c[2][1] * det2_01_03 + c[2][3] * det2_01_01;
393                 T det3_201_023 = c[2][0] * det2_01_23 - c[2][2] * det2_01_03 + c[2][3] * det2_01_02;
394                 T det3_201_123 = c[2][1] * det2_01_23 - c[2][2] * det2_01_13 + c[2][3] * det2_01_12;
395 
396                 T det = - det3_201_123 * c[3][0] + det3_201_023 * c[3][1] - det3_201_013 * c[3][2] + det3_201_012 * c[3][3];
397                 T invDet = 1 / det;
398 
399                 T det2_03_01 = c[0][0] * c[3][1] - c[0][1] * c[3][0];
400                 T det2_03_02 = c[0][0] * c[3][2] - c[0][2] * c[3][0];
401                 T det2_03_03 = c[0][0] * c[3][3] - c[0][3] * c[3][0];
402                 T det2_03_12 = c[0][1] * c[3][2] - c[0][2] * c[3][1];
403                 T det2_03_13 = c[0][1] * c[3][3] - c[0][3] * c[3][1];
404                 T det2_03_23 = c[0][2] * c[3][3] - c[0][3] * c[3][2];
405                 T det2_13_01 = c[1][0] * c[3][1] - c[1][1] * c[3][0];
406                 T det2_13_02 = c[1][0] * c[3][2] - c[1][2] * c[3][0];
407                 T det2_13_03 = c[1][0] * c[3][3] - c[1][3] * c[3][0];
408                 T det2_13_12 = c[1][1] * c[3][2] - c[1][2] * c[3][1];
409                 T det2_13_13 = c[1][1] * c[3][3] - c[1][3] * c[3][1];
410                 T det2_13_23 = c[1][2] * c[3][3] - c[1][3] * c[3][2];
411 
412                 T det3_203_012 = c[2][0] * det2_03_12 - c[2][1] * det2_03_02 + c[2][2] * det2_03_01;
413                 T det3_203_013 = c[2][0] * det2_03_13 - c[2][1] * det2_03_03 + c[2][3] * det2_03_01;
414                 T det3_203_023 = c[2][0] * det2_03_23 - c[2][2] * det2_03_03 + c[2][3] * det2_03_02;
415                 T det3_203_123 = c[2][1] * det2_03_23 - c[2][2] * det2_03_13 + c[2][3] * det2_03_12;
416 
417                 T det3_213_012 = c[2][0] * det2_13_12 - c[2][1] * det2_13_02 + c[2][2] * det2_13_01;
418                 T det3_213_013 = c[2][0] * det2_13_13 - c[2][1] * det2_13_03 + c[2][3] * det2_13_01;
419                 T det3_213_023 = c[2][0] * det2_13_23 - c[2][2] * det2_13_03 + c[2][3] * det2_13_02;
420                 T det3_213_123 = c[2][1] * det2_13_23 - c[2][2] * det2_13_13 + c[2][3] * det2_13_12;
421 
422                 T det3_301_012 = c[3][0] * det2_01_12 - c[3][1] * det2_01_02 + c[3][2] * det2_01_01;
423                 T det3_301_013 = c[3][0] * det2_01_13 - c[3][1] * det2_01_03 + c[3][3] * det2_01_01;
424                 T det3_301_023 = c[3][0] * det2_01_23 - c[3][2] * det2_01_03 + c[3][3] * det2_01_02;
425                 T det3_301_123 = c[3][1] * det2_01_23 - c[3][2] * det2_01_13 + c[3][3] * det2_01_12;
426 
427                 Matrix res = void;
428                 res.c[0][0] = - det3_213_123 * invDet;
429                 res.c[1][0] = + det3_213_023 * invDet;
430                 res.c[2][0] = - det3_213_013 * invDet;
431                 res.c[3][0] = + det3_213_012 * invDet;
432 
433                 res.c[0][1] = + det3_203_123 * invDet;
434                 res.c[1][1] = - det3_203_023 * invDet;
435                 res.c[2][1] = + det3_203_013 * invDet;
436                 res.c[3][1] = - det3_203_012 * invDet;
437 
438                 res.c[0][2] = + det3_301_123 * invDet;
439                 res.c[1][2] = - det3_301_023 * invDet;
440                 res.c[2][2] = + det3_301_013 * invDet;
441                 res.c[3][2] = - det3_301_012 * invDet;
442 
443                 res.c[0][3] = - det3_201_123 * invDet;
444                 res.c[1][3] = + det3_201_023 * invDet;
445                 res.c[2][3] = - det3_201_013 * invDet;
446                 res.c[3][3] = + det3_201_012 * invDet;
447                 return res;
448             }
449         }
450 
451         /// Returns: transposed matrice.
452         Matrix!(T, C, R) transposed() pure const nothrow @nogc
453         {
454             Matrix!(T, C, R) res;
455             for (size_t i = 0; i < C; ++i)
456                 for (size_t j = 0; j < R; ++j)
457                     res.c[i][j] = c[j][i];
458             return res;
459         }
460 
461         static if (isSquare && R > 1)
462         {
463             /// In-place translate by (v, 1)
464             void translate(Vector!(T, R-1) v) pure nothrow @nogc
465             {
466                 for (size_t i = 0; i < R; ++i)
467                 {
468                     T dot = 0;
469                     for (size_t j = 0; j + 1 < C; ++j)
470                         dot += v.v[j] * c[i][j];
471 
472                     c[i][C-1] += dot;
473                 }
474             }
475 
476             /// Make a translation matrix.
477             static Matrix translation(Vector!(T, R-1) v) pure nothrow @nogc
478             {
479                 Matrix res = identity();
480                 for (size_t i = 0; i + 1 < R; ++i)
481                     res.c[i][C-1] += v.v[i];
482                 return res;
483             }
484 
485             /// In-place matrix scaling.
486             void scale(Vector!(T, R-1) v) pure nothrow
487             {
488                 for (size_t i = 0; i < R; ++i)
489                     for (size_t j = 0; j + 1 < C; ++j)
490                         c[i][j] *= v.v[j];
491             }
492 
493             /// Make a scaling matrix.
494             static Matrix scaling(Vector!(T, R-1) v) pure nothrow @nogc
495             {
496                 Matrix res = identity();
497                 for (size_t i = 0; i + 1 < R; ++i)
498                     res.c[i][i] = v.v[i];
499                 return res;
500             }
501         }
502 
503         // rotations are implemented for 3x3 and 4x4 matrices.
504         static if (isSquare && (R == 3 || R == 4) && isFloatingPoint!T)
505         {
506             public static Matrix rotateAxis(size_t i, size_t j)(T angle) pure nothrow @nogc
507             {
508                 Matrix res = identity();
509                 const T cosa = cos(angle);
510                 const T sina = sin(angle);
511                 res.c[i][i] = cosa;
512                 res.c[i][j] = -sina;
513                 res.c[j][i] = sina;
514                 res.c[j][j] = cosa;
515                 return res;
516             }
517 
518             /// Returns: rotation matrix along axis X
519             alias rotateAxis!(1, 2) rotateX;
520 
521             /// Returns: rotation matrix along axis Y
522             alias rotateAxis!(2, 0) rotateY;
523 
524             /// Returns: rotation matrix along axis Z
525             alias rotateAxis!(0, 1) rotateZ;
526 
527             /// Similar to the glRotate matrix, however the angle is expressed in radians
528             /// See_also: $(LINK http://www.cs.rutgers.edu/~decarlo/428/gl_man/rotate.html)
529             static Matrix rotation(T angle, vec3!T axis) pure nothrow @nogc
530             {
531                 Matrix res = identity();
532                 const T c = cos(angle);
533                 const oneMinusC = 1 - c;
534                 const T s = sin(angle);
535                 axis = axis.normalized();
536                 T x = axis.x,
537                   y = axis.y,
538                   z = axis.z;
539                 T xy = x * y,
540                   yz = y * z,
541                   xz = x * z;
542 
543                 res.c[0][0] = x * x * oneMinusC + c;
544                 res.c[0][1] = x * y * oneMinusC - z * s;
545                 res.c[0][2] = x * z * oneMinusC + y * s;
546                 res.c[1][0] = y * x * oneMinusC + z * s;
547                 res.c[1][1] = y * y * oneMinusC + c;
548                 res.c[1][2] = y * z * oneMinusC - x * s;
549                 res.c[2][0] = z * x * oneMinusC - y * s;
550                 res.c[2][1] = z * y * oneMinusC + x * s;
551                 res.c[2][2] = z * z * oneMinusC + c;
552                 return res;
553             }
554         }
555 
556         // 4x4 specific transformations for 3D usage
557         static if (isSquare && R == 4 && isFloatingPoint!T)
558         {
559             /// Returns: orthographic projection.
560             static Matrix orthographic(T left, T right, T bottom, T top, T near, T far) pure nothrow @nogc
561             {
562                 T dx = right - left,
563                   dy = top - bottom,
564                   dz = far - near;
565 
566                 T tx = -(right + left) / dx;
567                 T ty = -(top + bottom) / dy;
568                 T tz = -(far + near)   / dz;
569 
570                 return Matrix(2 / dx,   0,      0,    tx,
571                                 0,    2 / dy,   0,    ty,
572                                 0,      0,    2 / dz, tz,
573                                 0,      0,      0,     1);
574             }
575 
576             /// Returns: perspective projection.
577             static Matrix perspective(T FOVInRadians, T aspect, T zNear, T zFar) pure nothrow @nogc
578             {
579                 T f = 1 / tan(FOVInRadians / 2);
580                 T d = 1 / (zNear - zFar);
581 
582                 return Matrix(f / aspect, 0,                  0,                    0,
583                                        0, f,                  0,                    0,
584                                        0, 0, (zFar + zNear) * d, 2 * d * zFar * zNear,
585                                        0, 0,                 -1,                    0);
586             }
587 
588             /// Returns: "lookAt" projection.
589             /// Thanks to vuaru for corrections.
590             static Matrix lookAt(vec3!T eye, vec3!T target, vec3!T up) pure nothrow @nogc
591             {
592                 vec3!T Z = (eye - target).normalized();
593                 vec3!T X = cross(-up, Z).normalized();
594                 vec3!T Y = cross(Z, -X);
595 
596                 return Matrix(-X.x,        -X.y,        -X.z,      dot(X, eye),
597                                Y.x,         Y.y,         Y.z,     -dot(Y, eye),
598                                Z.x,         Z.y,         Z.z,     -dot(Z, eye),
599                                0,           0,           0,        1);
600             }
601 
602             /// Extract frustum from a 4x4 matrice.
603             Frustum!T frustum() pure const nothrow @nogc
604             {
605                 auto left   = Plane!T(row(3) + row(0));
606                 auto right  = Plane!T(row(3) - row(0));
607                 auto top    = Plane!T(row(3) - row(1));
608                 auto bottom = Plane!T(row(3) + row(1));
609                 auto near   = Plane!T(row(3) + row(2));
610                 auto far    = Plane!T(row(3) - row(2));
611                 return Frustum!T(left, right, top, bottom, near, far);
612             }
613 
614         }
615     }
616 
617     private
618     {
619         alias T _T;
620         enum _R = R;
621         enum _C = C;
622         enum bool _isMatrix = true;
623 
624         template isAssignable(T)
625         {
626             enum bool isAssignable =
627                 is(typeof(
628                 {
629                     T x;
630                     Matrix m = x;
631                 }()));
632         }
633 
634         template isConvertible(T)
635         {
636             enum bool isConvertible = (!is(T : Matrix)) && isAssignable!T;
637         }
638 
639         template isTConvertible(U)
640         {
641             enum bool isTConvertible = is(U : T);
642         }
643 
644         template isRowConvertible(U)
645         {
646             enum bool isRowConvertible = is(U : row_t);
647         }
648 
649         template isColumnConvertible(U)
650         {
651             enum bool isColumnConvertible = is(U : column_t);
652         }
653     }
654 
655     public
656     {
657         /// Returns: an identity matrice.
658         /// Note: the identity matrix, while only meaningful for square matrices,
659         /// is also defined for non-square ones.
660         static Matrix identity() pure nothrow @nogc
661         {
662             Matrix res = void;
663             for (size_t i = 0; i < R; ++i)
664                 for (size_t j = 0; j < C; ++j)
665                     res.c[i][j] = (i == j) ? 1 : 0;
666             return res;
667         }
668 
669         /// Returns: a constant matrice.
670         static Matrix constant(U)(U x) pure nothrow @nogc
671         {
672             Matrix res = void;
673 
674             for (size_t i = 0; i < R * C; ++i)
675                 res.v[i] = cast(T)x;
676             return res;
677         }
678     }
679 }
680 
681 // GLSL is a big inspiration here
682 // we defines types with more or less the same names
683 template mat2x2(T) { alias Matrix!(T, 2u, 2u) mat2x2; }
684 template mat3x3(T) { alias Matrix!(T, 3u, 3u) mat3x3; }
685 template mat4x4(T) { alias Matrix!(T, 4u, 4u) mat4x4; }
686 
687 // WARNING: in GLSL, first number is _columns_, second is rows
688 // It is the opposite here: first number is rows, second is columns
689 // With this convention mat2x3 * mat3x4 -> mat2x4.
690 template mat2x3(T) { alias Matrix!(T, 2u, 3u) mat2x3; }
691 template mat2x4(T) { alias Matrix!(T, 2u, 4u) mat2x4; }
692 template mat3x2(T) { alias Matrix!(T, 3u, 2u) mat3x2; }
693 template mat3x4(T) { alias Matrix!(T, 3u, 4u) mat3x4; }
694 template mat4x2(T) { alias Matrix!(T, 4u, 2u) mat4x2; }
695 template mat4x3(T) { alias Matrix!(T, 4u, 3u) mat4x3; }
696 
697 alias mat2x2 mat2;
698 alias mat3x3 mat3;  // shorter names for most common matrices
699 alias mat4x4 mat4;
700 
701 private string definePostfixAliases(string type)
702 {
703     return "alias " ~ type ~ "!byte "   ~ type ~ "b;\n"
704          ~ "alias " ~ type ~ "!ubyte "  ~ type ~ "ub;\n"
705          ~ "alias " ~ type ~ "!short "  ~ type ~ "s;\n"
706          ~ "alias " ~ type ~ "!ushort " ~ type ~ "us;\n"
707          ~ "alias " ~ type ~ "!int "    ~ type ~ "i;\n"
708          ~ "alias " ~ type ~ "!uint "   ~ type ~ "ui;\n"
709          ~ "alias " ~ type ~ "!long "   ~ type ~ "l;\n"
710          ~ "alias " ~ type ~ "!ulong "  ~ type ~ "ul;\n"
711          ~ "alias " ~ type ~ "!float "  ~ type ~ "f;\n"
712          ~ "alias " ~ type ~ "!double " ~ type ~ "d;\n"
713          ~ "alias " ~ type ~ "!real "   ~ type ~ "L;\n";
714 }
715 
716 // define a lot of type names
717 mixin(definePostfixAliases("mat2"));
718 mixin(definePostfixAliases("mat3"));
719 mixin(definePostfixAliases("mat4"));
720 mixin(definePostfixAliases("mat2x2"));
721 mixin(definePostfixAliases("mat2x3"));
722 mixin(definePostfixAliases("mat2x4"));
723 mixin(definePostfixAliases("mat3x2"));
724 mixin(definePostfixAliases("mat3x3"));
725 mixin(definePostfixAliases("mat3x4"));
726 mixin(definePostfixAliases("mat4x2"));
727 mixin(definePostfixAliases("mat4x3"));
728 mixin(definePostfixAliases("mat4x4"));
729 
730 unittest
731 {
732     mat2i x = mat2i(0, 1,
733                     2, 3);
734     assert(x.c[0][0] == 0 && x.c[0][1] == 1 && x.c[1][0] == 2 && x.c[1][1] == 3);
735 
736     vec2i[2] cols = [vec2i(0, 2), vec2i(1, 3)];
737     mat2i y = mat2i.fromColumns(cols[]);
738     assert(y.c[0][0] == 0 && y.c[0][1] == 1 && y.c[1][0] == 2 && y.c[1][1] == 3);
739     y = mat2i.fromRows(cols[]);
740     assert(y.c[0][0] == 0 && y.c[1][0] == 1 && y.c[0][1] == 2 && y.c[1][1] == 3);
741     y = y.transposed();
742 
743     assert(x == y);
744     x = [0, 1, 2, 3];
745     assert(x == y);
746 
747 
748     mat2i z = x * y;
749     assert(z == mat2i([2, 3, 6, 11]));
750     vec2i vz = z * vec2i(2, -1);
751     assert(vz == vec2i(1, 1));
752 
753     mat2f a = z;
754     mat2d ad = a;
755     ad += a;
756     mat2f w = [4, 5, 6, 7];
757     z = cast(mat2i)w;
758     assert(w == z);
759 
760     {
761         mat2x3f A;
762         mat3x4f B;
763         mat2x4f C = A * B;
764     }
765 }