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