dvec.matrix source code

1 /** 
2  * This module contains the `Mat` templated struct representing matrices and
3  * their operations.
4  */
5 module dvec.matrix;
6 
7 import std.traits : isNumeric, isFloatingPoint;
8 import dvec.vector;
9 import dvec.vector_types;
10 
11 /** 
12  * Generic struct that represents a matrix with `rowCount` rows and `colCount`
13  * columns, holding elements of type `T`. A matrix must be at least 1x1.
14  */
15 struct Mat(T, size_t rowCount, size_t colCount) if (isNumeric!T && rowCount > 0 && colCount > 0) {
16     public static const size_t WIDTH = colCount;
17     public static const size_t HEIGHT = rowCount;
18 
19     /** 
20      * The internal static array storing the elements in row-major form.
21      */
22     public T[rowCount * colCount] data;
23 
24     /** 
25      * Constructs a matrix from the given elements.
26      * Params:
27      *   elements = The elements to place in the matrix. Should be in row-major
28      *              form.
29      */
30     public this(T[rowCount * colCount] elements) {
31         setData(elements);
32     }
33 
34     /** 
35      * Constructs a matrix from the given elements.
36      * Params:
37      *   elements = The elements to place in the matrix. Should be in row-major
38      *              form.
39      */
40     public this(T[] elements...) {
41         data[0 .. $] = elements[0 .. data.length];
42     }
43 
44     /** 
45      * Constructs a matrix from the given 2-dimensional array.
46      * Params:
47      *   elements = The elements to place in the matrix.
48      */
49     public this(T[colCount][rowCount] elements) {
50         static foreach (i; 0 .. rowCount) {
51             static foreach (j; 0 .. colCount) {
52                 this[i, j] = elements[i][j];
53             }
54         }
55     }
56 
57     /** 
58      * Constructs a matrix as a copy of another.
59      * Params:
60      *   other = The matrix to copy.
61      */
62     public this(Mat!(T, rowCount, colCount) other) {
63         setData(other.data);
64     }
65 
66     /** 
67      * Constructs a matrix with all elements initialized to the given value.
68      * Params:
69      *   value = The value to set all elements to.
70      */
71     public this(T value) {
72         static foreach (i; 0 .. data.length) data[i] = value;
73     }
74 
75     /** 
76      * Helper method to convert a 2-dimensional `[i, j]` index into a
77      * 1-dimensional `[i]` index.
78      * Params:
79      *   i = The row number.
80      *   j = The column number.
81      * Returns: A 1-dimensional index.
82      */
83     private static size_t convertToIndex(size_t i, size_t j) {
84         return colCount * i + j;
85     }
86 
87     /** 
88      * Gets the value at a specified location in the matrix. For example,
89      * with a `Mat!(float, 2, 2) m`, we can say `float v = m[0, 1];`
90      * Params:
91      *   i = The row number, starting from 0.
92      *   j = The column number, starting from 0.
93      * Returns: The value at the specified location.
94      */
95     public T opIndex(size_t i, size_t j) const {
96         return data[convertToIndex(i, j)];
97     }
98 
99     /** 
100      * Sets the value at a specified location in the matrix.
101      * Params:
102      *   value = The value to assign.
103      *   i = The row number, starting from 0.
104      *   j = The column number, starting from 0.
105      */
106     public void opIndexAssign(T value, size_t i, size_t j) {
107         data[convertToIndex(i, j)] = value;
108     }
109 
110     /** 
111      * Gets a specified row as a vector.
112      * Params:
113      *   row = The row number, starting from 0.
114      * Returns: The row.
115      */
116     public Vec!(T, colCount) getRow(size_t row) {
117         size_t idx = convertToIndex(row, 0);
118         return Vec!(T, colCount)(data[idx .. idx + colCount]);
119     }
120 
121     /** 
122      * Sets a specified row to the given vector of elements.
123      * Params:
124      *   row = The row number, starting from 0.
125      *   vector = The elements to set in the row.
126      * Returns: A reference to this matrix, for method chaining.
127      */
128     public ref Mat!(T, rowCount, colCount) setRow(size_t row, Vec!(T, colCount) vector) {
129         size_t idx = convertToIndex(row, 0);
130         data[idx .. idx + colCount] = vector.data;
131         return this;
132     }
133 
134     /** 
135      * Gets a specified column as a vector.
136      * Params:
137      *   col = The column number, starting from 0.
138      * Returns: The column.
139      */
140     public Vec!(T, rowCount) getCol(size_t col) const {
141         Vec!(T, rowCount) v;
142         static foreach (i; 0 .. rowCount) {
143             v[i] = this[i, col];
144         }
145         return v;
146     }
147 
148     /** 
149      * Sets a specified column to the given vector of elements.
150      * Params:
151      *   col = The column number, starting from 0.
152      *   vector = The elements to set in the column.
153      * Returns: A reference to this matrix, for method chaining.
154      */
155     public ref Mat!(T, rowCount, colCount) setCol(size_t col, Vec!(T, rowCount) vector) {
156         static foreach (i; 0 .. rowCount) {
157             this[i, col] = vector[i];
158         }
159         return this;
160     }
161 
162     /** 
163      * Sets all elements of this matrix using the given elements.
164      * Params:
165      *   elements = The elements to set.
166      * Returns: A reference to this matrix, for method chaining.
167      */
168     public ref Mat!(T, rowCount, colCount) setData(T[rowCount * colCount] elements) {
169         data[0 .. $] = elements[0 .. $];
170         return this;
171     }
172 
173     /** 
174      * Gets a copy of this matrix.
175      * Returns: A copy of this matrix.
176      */
177     public Mat!(T, rowCount, colCount) copy() const {
178         return Mat!(T, rowCount, colCount)(this);
179     }
180 
181     /** 
182      * Adds a given matrix to this one.
183      * Params:
184      *   other = The other matrix to add to this one.
185      * Returns: A reference to this matrix, for method chaining.
186      */
187     public ref Mat!(T, rowCount, colCount) add(Mat!(T, rowCount, colCount) other) {
188         static foreach (i; 0 .. data.length) data[i] += other.data[i];
189         return this;
190     }
191 
192     /** 
193      * Subtracts a given matrix from this one.
194      * Params:
195      *   other = The other matrix to subtract from this one.
196      * Returns: A reference to this matrix, for method chaining.
197      */
198     public ref Mat!(T, rowCount, colCount) sub(Mat!(T, rowCount, colCount) other) {
199         static foreach (i; 0 .. data.length) data[i] -= other.data[i];
200         return this;
201     }
202 
203     /** 
204      * Multiplies this matrix by a factor.
205      * Params:
206      *   factor = The factor to muliply by.
207      * Returns: A reference to this matrix, for method chaining.
208      */
209     public ref Mat!(T, rowCount, colCount) mul(T factor) {
210         static foreach (i; 0 .. data.length) data[i] *= factor;
211         return this;
212     }
213 
214     /** 
215      * Divides this matrix by a factor.
216      * Params:
217      *   factor = The factor to divide by.
218      * Returns: A reference to this matrix, for method chaining.
219      */
220     public ref Mat!(T, rowCount, colCount) div(T factor) {
221         static foreach (i; 0 .. data.length) data[i] /= factor;
222         return this;
223     }
224 
225     /** 
226      * Gets the [transpose](https://en.wikipedia.org/wiki/Transpose) of this
227      * matrix, and stores it in this matrix.
228      * Returns: A reference to this matrix, for method chaining.
229      */
230     public ref Mat!(T, rowCount, colCount) transpose() {
231         Mat!(T, colCount, rowCount) m;
232         static foreach (i; 0 .. rowCount) {
233             static foreach (j; 0 .. colCount) {
234                 m[j, i] = this[i, j];
235             }
236         }
237         return setData(m.data);
238     }
239 
240     /** 
241      * Computes the [matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication)
242      * of `this * other`.
243      * Params:
244      *   other = The matrix to multiply with this one.
245      * Returns: A new matrix of size `this.rowCount x other.colCount`.
246      */
247     public Mat!(T, rowCount, otherColCount) mul(T, size_t otherRowCount, size_t otherColCount)
248         (Mat!(T, otherRowCount, otherColCount) other) const {
249         Mat!(T, rowCount, otherColCount) m;
250         T sum;
251         static foreach (i; 0 .. rowCount) {
252             static foreach (j; 0 .. otherColCount) {
253                 sum = 0;
254                 static foreach (k; 0 .. colCount) {
255                     sum += this[i, k] * other[k, j];
256                 }
257                 m[i, j] = sum;
258             }
259         }
260         return m;
261     }
262 
263     /** 
264      * Multiplies a vector against this matrix.
265      * Params:
266      *   vector = The vector to multiply.
267      * Returns: The resultant transformed vector.
268      */
269     public Vec!(T, rowCount) mul(Vec!(T, colCount) vector) const {
270         Vec!(T, rowCount) result;
271         T sum;
272         static foreach (i; 0 .. rowCount) {
273             sum = 0;
274             static foreach (j; 0 .. colCount) {
275                 sum += this[i, j] * vector[j];
276             }
277             result[i] = sum;
278         }
279         return result;
280     }
281 
282     /** 
283      * Switches two rows in this matrix.
284      * Params:
285      *   rowI = A row to swap, starting from 0.
286      *   rowJ = A row to swap, starting from 0.
287      * Returns: A reference to this matrix, for method chaining.
288      */
289     public ref Mat!(T, rowCount, colCount) rowSwitch(size_t rowI, size_t rowJ) {
290         auto r = getRow(rowI);
291         setRow(rowI, getRow(rowJ));
292         setRow(rowJ, r);
293         return this;
294     }
295 
296     /** 
297      * Multiplies a row by a factor.
298      * Params:
299      *   row = The row number, starting from 0.
300      *   factor = The factor to multiply by.
301      * Returns: A reference to this matrix, for method chaining.
302      */
303     public ref Mat!(T, rowCount, colCount) rowMultiply(size_t row, T factor) {
304         size_t idx = convertToIndex(row, 0);
305         static foreach (i; 0 .. colCount) {
306             data[idx + i] *= factor;
307         }
308         return this;
309     }
310 
311     /** 
312      * Adds a row, multiplied by a factor, to another row.
313      * Params:
314      *   rowI = The row number to add to.
315      *   factor = The factor to multiply `rowJ` by.
316      *   rowJ = The row to add to `rowI`.
317      * Returns: A reference to this matrix, for method chaining.
318      */
319     public ref Mat!(T, rowCount, colCount) rowAdd(size_t rowI, T factor, size_t rowJ) {
320         auto row = getRow(rowJ);
321         row.mul(factor);
322         setRow(rowI, row);
323         return this;
324     }
325 
326     /** 
327      * Gets a submatrix of this matrix, with the given rows and columns removed.
328      * Params:
329      *   rows = The set of rows to remove.
330      *   cols = The set of columns to remove.
331      * Returns: A matrix with the given rows and columns removed.
332      */
333     public Mat!(T, rowCount - n, colCount - m) subMatrix(size_t n, size_t m)(size_t[n] rows, size_t[m] cols) const
334         if (rowCount - n > 0 && colCount - m > 0) {
335         // TODO: Improve efficiency with static stuff.
336         Mat!(T, rowCount - n, colCount - m) sub;
337         size_t subIdx = 0;
338         foreach (idx; 0 .. data.length) {
339             size_t row = idx / colCount;
340             size_t col = idx % colCount;
341             bool skip = false;
342             foreach (r; rows) {
343                 if (r == row) {
344                     skip = true;
345                     break;
346                 }
347             }
348             if (!skip) {
349                 foreach (c; cols) {
350                     if (c == col) {
351                         skip = true;
352                         break;
353                     }
354                 }
355             }
356             if (!skip) {
357                 sub.data[subIdx++] = this[row, col];
358             }
359         }
360         return sub;
361     }
362 
363     /** 
364      * Converts this matrix to a string.
365      * Returns: A string representation of the matrix.
366      */
367     public string toString() const {
368         import std.conv : to;
369         import std.algorithm : max;
370         import std.array : appender, replicate;
371         string[data.length] values;
372         size_t[colCount] columnWidths;
373         foreach (i; 0 .. data.length) {
374             values[i] = data[i].to!string;
375             size_t colIdx = i % colCount;
376             columnWidths[colIdx] = max(columnWidths[colIdx], values[i].length + 1);
377         }
378         auto s = appender!string;
379         foreach (r; 0 .. rowCount) {
380             s ~= "| ";
381             foreach (c; 0 .. colCount) {
382                 size_t idx = convertToIndex(r, c);
383                 size_t padAmount = columnWidths[c] - values[idx].length;
384                 string padding = replicate(" ", padAmount);
385                 s ~= padding ~ values[idx];
386                 if (c < colCount - 1) s ~= ", ";
387             }
388             s ~= " |";
389             if (r < rowCount - 1) s ~= "\n";
390         }
391         return s[];
392     }
393 
394     // Special methods for square matrices.
395     static if (rowCount == colCount) {
396         alias N = rowCount;
397 
398         /** 
399          * Gets an [identity matrix](https://en.wikipedia.org/wiki/Identity_matrix).
400          * Returns: An identity matrix.
401          */
402         public static Mat!(T, N, N) identity() {
403             Mat!(T, N, N) m;
404             static foreach (i; 0 .. N) {
405                 static foreach (j; 0 .. N) {
406                     m[i, j] = i == j ? 1 : 0;
407                 }
408             }
409             return m;
410         }
411 
412         /** 
413          * Gets the [determinant](https://en.wikipedia.org/wiki/Determinant)
414          * of this matrix.
415          * Returns: The determinant of this matrix.
416          */
417         public T det() const {
418             static if (N == 1) {
419                 return data[0];
420             } else static if (N == 2) {
421                 return data[0] * data[3] - data[1] * data[2];
422             } else {
423                 // Laplace expansion, taking i = 0.
424                 T sum = 0;
425                 static foreach (j; 0 .. N) {
426                     sum += (j % 2 == 0 ? 1 : -1) * this[0, j] * this.subMatrix([0], [j]).det();
427                 }
428                 return sum;
429             }
430         }
431 
432         /** 
433          * Determines if this matrix is invertible, which simply means a
434          * non-zero determinant.
435          * Returns: True if this matrix is invertible.
436          */
437         public bool invertible() const {
438             return det() != 0;
439         }
440 
441         /** 
442          * Gets a [cofactor matrix](https://en.wikipedia.org/wiki/Minor_(linear_algebra)#Cofactor_expansion_of_the_determinant).
443          * Returns: A reference to this matrix, for method chaining.
444          */
445         public ref Mat!(T, N, N) cofactor() {
446             static if (N > 1) {
447                 Mat!(T, N, N) c;
448                 static foreach (i; 0 .. N) {
449                     static foreach (j; 0 .. N) {
450                         c[i, j] = ((i + j) % 2 == 0 ? 1 : -1) * this.subMatrix([i], [j]).det();
451                     }
452                 }
453                 setData(c.data);
454             }
455             return this;
456         }
457 
458         /** 
459          * Gets an [adjugate matrix](https://en.wikipedia.org/wiki/Adjugate_matrix).
460          * Returns: A reference to this matrix, for method chaining.
461          */
462         public ref Mat!(T, N, N) adjugate() {
463             cofactor();
464             transpose();
465             return this;
466         }
467 
468         /** 
469          * Gets the inverse of this matrix.
470          * Returns: A reference to this matrix, for method chaining.
471          */
472         public ref Mat!(T, N, N) inv() {
473             T d = det();
474             adjugate();
475             div(d);
476             return this;
477         }
478 
479         // Special case for 3x3 floating-point matrices: linear transformations
480         static if (N == 3 && isFloatingPoint!T) {
481             /** 
482              * Applies a 2D rotation to this matrix.
483              * Params:
484              *   theta = The angle in radians.
485              * Returns: A reference to this matrix, for method chaining.
486              */
487             public ref Mat!(T, N, N) rotate(T theta) {
488                 import std.math : cos, sin;
489                 this = mul(Mat!(T, N, N)(
490                     cos(theta), -sin(theta), 0,
491                     sin(theta), cos(theta),  0,
492                     0,          0,           1
493                 ));
494                 return this;
495             }
496 
497             /** 
498              * Applies a 2D translation to this matrix.
499              * Params:
500              *   dx = The translation on the x-axis.
501              *   dy = The translation on the y-axis.
502              * Returns: A reference to this matrix, for method chaining.
503              */
504             public ref Mat!(T, N, N) translate(T dx, T dy) {
505                 this = mul(Mat!(T, N, N)(
506                     1, 0, dx,
507                     0, 1, dy,
508                     0, 0, 1
509                 ));
510                 return this;
511             }
512 
513             /** 
514              * Applies a 2D scaling to this matrix.
515              * Params:
516              *   sx = The scale factor on the x-axis.
517              *   sy = The scale factor on the y-axis.
518              * Returns: A reference to this matrix, for method chaining.
519              */
520             public ref Mat!(T, N, N) scale(T sx, T sy) {
521                 this = mul(Mat!(T, N, N)(
522                     sx, 0,  0,
523                     0,  sy, 0,
524                     0,  0,  1
525                 ));
526                 return this;
527             }
528 
529             /** 
530              * Applies a uniform 2D scaling to this matrix on all axes.
531              * Params:
532              *   s = The scale factor to apply.
533              * Returns: A reference to this matrix, for method chaining.
534              */
535             public ref Mat!(T, N, N) scale(T s) {
536                 return scale(s, s);
537             }
538 
539             /** 
540              * Applies a 2D shear to this matrix.
541              * Params:
542              *   sx = The shear factor on the x-axis.
543              *   sy = The shear factor on the y-axis.
544              * Returns: A reference to this matrix, for method chaining.
545              */
546             public ref Mat!(T, N, N) shear(T sx, T sy) {
547                 this = mul(Mat!(T, N, N)(
548                     1,  sx, 0,
549                     sy, 1,  0,
550                     0,  0,  1
551                 ));
552                 return this;
553             }
554 
555             /** 
556              * Maps the given 2D vector using this matrix.
557              * Params:
558              *   v = The vector to map.
559              * Returns: The resultant vector.
560              */
561             public Vec!(T, 2) map(Vec!(T, 2) v) {
562                 Vec!(T, 3) v1 = this.mul(Vec!(T, 3)(v[0], v[1], 1));
563                 return Vec!(T, 2)(v1[0], v1[1]);
564             }
565         }
566 
567         // Special case for 4x4 floating-point matrices: linear transformations
568         static if (N == 4 && isFloatingPoint!T) {
569             /** 
570              * Applies a 3D translation to this matrix.
571              * Params:
572              *   dx = The translation on the x-axis.
573              *   dy = The translation on the y-axis.
574              *   dz = The translation on the z-axis.
575              * Returns: A reference to this matrix, for method chaining.
576              */
577             public ref Mat!(T, N, N) translate(T dx, T dy, T dz) {
578                 this = mul(Mat!(T, N, N)(
579                     1, 0, 0, dx,
580                     0, 1, 0, dy,
581                     0, 0, 1, dz,
582                     0, 0, 0, 1
583                 ));
584                 return this;
585             }
586 
587             /** 
588              * Applies a 3D scaling to this matrix.
589              * Params:
590              *   sx = The scale factor on the x-axis.
591              *   sy = The scale factor on the y-axis.
592              *   sz = The scale factor on the z-axis.
593              * Returns: A reference to this matrix, for method chaining.
594              */
595             public ref Mat!(T, N, N) scale(T sx, T sy, T sz) {
596                 this = mul(Mat!(T, N, N)(
597                     sx, 0,  0,  0,
598                     0,  sy, 0,  0,
599                     0,  0,  sz, 0,
600                     0,  0,  0,  1
601                 ));
602                 return this;
603             }
604 
605             /** 
606              * Applies a uniform 3D scaling to this matrix on all axes.
607              * Params:
608              *   s = The scale factor to apply.
609              * Returns: A reference to this matrix, for method chaining.
610              */
611             public ref Mat!(T, N, N) scale(T s) {
612                 return scale(s, s, s);
613             }
614 
615             /** 
616              * Applies a rotation to this matrix about the x-axis.
617              * Params:
618              *   theta = The angle in radians.
619              * Returns: A reference to this matrix, for method chaining.
620              */
621             public ref Mat!(T, N, N) rotateX(T theta) {
622                 import std.math : cos, sin;
623                 this = mul(Mat!(T, N, N)(
624                     1, 0,          0,           0,
625                     0, cos(theta), -sin(theta), 0,
626                     0, sin(theta), cos(theta),  0,
627                     0, 0,          0,           1
628                 ));
629                 return this;
630             }
631 
632             /** 
633              * Applies a rotation to this matrix about the y-axis.
634              * Params:
635              *   theta = The angle in radians.
636              * Returns: A reference to this matrix, for method chaining.
637              */
638             public ref Mat!(T, N, N) rotateY(T theta) {
639                 import std.math : cos, sin;
640                 this = mul(Mat!(T, N, N)(
641                     cos(theta),  0, sin(theta), 0,
642                     0,           1, 0,          0,
643                     -sin(theta), 0, cos(theta), 0,
644                     0,           0, 0,          1
645                 ));
646                 return this;
647             }
648 
649             /** 
650              * Applies a rotation to this matrix about the z-axis.
651              * Params:
652              *   theta = The angle in radians.
653              * Returns: A reference to this matrix, for method chaining.
654              */
655             public ref Mat!(T, N, N) rotateZ(T theta) {
656                 import std.math : cos, sin;
657                 this = mul(Mat!(T, N, N)(
658                     cos(theta), -sin(theta), 0, 0,
659                     sin(theta), cos(theta),  0, 0,
660                     0,          0,           1, 0,
661                     0,          0,           0, 1
662                 ));
663                 return this;
664             }
665 
666             /** 
667              * Applies a 3D rotation to this matrix about the x, y, and then z
668              * axes.
669              * Params:
670              *   x = The angle to rotate about the x-axis.
671              *   y = The angle to rotate about the y-axis.
672              *   z = The angle to rotate about the z-axis.
673              * Returns: A reference to this matrix, for method chaining.
674              */
675             public ref Mat!(T, N, N) rotate(T x, T y, T z) {
676                 rotateX(x);
677                 rotateY(y);
678                 rotateZ(z);
679                 return this;
680             }
681         }
682     }
683 }
684 
685 // Aliases for common matrix types.
686 alias Mat2f = Mat!(float, 2, 2);
687 alias Mat3f = Mat!(float, 3, 3);
688 alias Mat4f = Mat!(float, 4, 4);
689 
690 alias Mat2d = Mat!(double, 2, 2);
691 alias Mat3d = Mat!(double, 3, 3);
692 alias Mat4d = Mat!(double, 4, 4);
693 
694 unittest {
695     import std.stdio;
696     import dvec.vector;
697 
698     auto m1 = Mat3d();
699     assert(m1.data.length == 9);
700 
701     auto m2 = Mat!(double, 2, 3)([1, 2, 3, 0, -6, 7]);
702     assert(m2.transpose.data == [1, 0, 2, -6, 3, 7]);
703 
704     auto m4 = Mat2d([1, 2, 3, 4]);
705     assert(m4.getRow(0).data == [1, 2]);
706     assert(m4.getRow(1).data == [3, 4]);
707     assert(m4.getCol(0).data == [1, 3]);
708     assert(m4.getCol(1).data == [2, 4]);
709     auto m5 = m4.mul(Mat2d([0, 1, 0, 0]));
710     assert(m5.data == [0, 1, 0, 3]);
711 
712     assert(Mat2d([[1, 2], [3, 4]]).data == [1, 2, 3, 4]);
713     assert(Mat!(double, 2, 3)([[1, 2, 3], [-1, -2, -3]]).data == [1, 2, 3, -1, -2, -3]);
714 
715     assert(Mat2d.identity.data == [1, 0, 0, 1]);
716     assert(Mat3d.identity.data == [1, 0, 0, 0, 1, 0, 0, 0, 1]);
717 
718     auto m8 = Mat!(double, 2, 3)([1, -1, 2, 0, -3, 1]);
719     Vec3d v1 = Vec3d(2, 1, 0);
720     assert(m8.mul(v1).data == [1, -3]);
721 
722     Vec3f p = Vec3f(0, 0, 1);
723     Mat3f tx = Mat3f([
724         [1, 0, 42],
725         [0, 1, 64],
726         [0, 0, 1]
727     ]);
728     auto transformed = tx.mul(p);
729     assert(transformed.data == [42, 64, 1]);
730 
731     auto p2 = Vec2f(0, 0);
732     auto tx2 = Mat3f.identity();
733     tx2.translate(42, 64);
734     auto transformed2 = tx2.map(p2);
735     assert(transformed2.data == [42, 64]);
736 
737     auto m9 = Mat!(double, 3, 4)([
738         1, 2, 3, 4,
739         5, 6, 7, 8,
740         9, 10, 11, 12
741     ]);
742     auto m10 = m9.subMatrix([2], [1]);
743     assert(m10.data == [
744         1, 3, 4,
745         5, 7, 8
746     ]);
747     assert(m10.subMatrix(cast(size_t[])[], [0]).data == [
748         3, 4,
749         7, 8
750     ]);
751 
752     assert(Mat2d([3, 7, 1, -4]).det == -19);
753     assert(Mat2d([1, 2, 3, 4]).det == -2);
754     assert(Mat3d([1, 2, 3, 4, 5, 6, 7, 8, 9]).det == 0);
755 
756     assert(Mat2d.identity.inv() == Mat2d.identity());
757     assert(Mat3f.identity.inv() == Mat3f.identity());
758     assert(Mat2d(4, 7, 2, 6).inv() == Mat2d(0.6, -0.7, -0.2, 0.4));
759     assert(Mat2d(-3, 1, 5, -2).inv() == Mat2d(-2, -1, -5, -3));
760     assert(Mat3d(1, 3, 3, 1, 4, 3, 1, 3, 4).inv() == Mat3d(7, -3, -3, -1, 1, 0, -1, 0, 1));
761 }