1 /**
2 Defines a squared Matrix with column-major ordering from 2x2 to 4x4 size
3 */
4 module zmath.matrix;
5 
6 import zmath.vector;
7 
8 alias Mat2f = Matrix!(float, 2); /// 2D squared matrix of floats
9 alias Mat3f = Matrix!(float, 3); /// 3D squared matrix of floats
10 alias Mat4f = Matrix!(float, 4); ///	4D squared matrix of floats
11 
12 alias Mat2d = Matrix!(double, 2); /// 2D squared matrix of doubles
13 alias Mat3d = Matrix!(double, 3); /// 3D squared matrix of doubles
14 alias Mat4d = Matrix!(double, 4); ///	4D squared matrix of doubles
15 
16 alias Mat2r = Matrix!(real, 2); /// 2D squared matrix of reals
17 alias Mat3r = Matrix!(real, 3); /// 3D squared matrix of reals
18 alias Mat4r = Matrix!(real, 4); ///	4D squared matrix of reals
19 
20 /**
21  * Defines a squared Matrix of n = 2, 3 or 4 size, like a linear array of numbers in a row-major order
22  */
23 public struct Matrix(T, size_t dim_)
24       if (__traits(isFloating, T) && dim_ >= 2 && dim_ <= 4)
25 {
26 nothrow:
27   static enum size_t dim = dim_; /// Matrix Dimension
28   static enum size_t cells = dim * dim; /// Matrix number of cells
29 
30   private alias VRow = Vector!(T, dim);
31   private alias VCol = Vector!(T, dim);
32 
33   union {
34     private T[cells] cell; /// Matrix like of a array of cells in row-major
35     private T[dim][dim] c; /// Matrix like of a 2d array of cells in row-major
36     private VRow[dim_] row; /// Matrix like of a array of rows vectors
37   }
38 
39   // Consts
40   static if (dim == 2) { // 2x2
41     // dfmt off
42     /** Zero R2 matrix */
43     public static enum Matrix!(T,2) ZERO = Matrix!(T,2)([
44         0, 0,
45         0, 0
46     ]);
47     /** Identity R2 matrix */
48     public static enum Matrix!(T,2) IDENTITY = Matrix!(T,2)([
49         1, 0,
50         0, 1
51     ]);
52     // dfmt on
53   } else static if (dim == 3) { // 3x3
54     // dfmt off
55     /** Zero R3 matrix */
56     public static enum Matrix!(T,3) ZERO = Matrix!(T,3)([
57         0, 0, 0,
58         0, 0, 0,
59         0, 0, 0
60     ]);
61     /** Identity R3 matrix */
62     public static enum Matrix!(T,3) IDENTITY = Matrix!(T,3)([
63         1, 0, 0,
64         0, 1, 0,
65         0, 0, 1
66     ]);
67     // dfmt on
68   } else static if (dim == 4) { // 4x4
69     // dfmt off
70     /** Zero R4 matrix */
71     public static enum Matrix!(T,4) ZERO = Matrix!(T,4)([
72         0, 0, 0, 0,
73         0, 0, 0, 0,
74         0, 0, 0, 0,
75         0, 0, 0, 0
76     ]);
77     /** Identity R4 matrix */
78     public static enum Matrix!(T,4) IDENTITY = Matrix!(T,4)([
79         1, 0, 0, 0,
80         0, 1, 0, 0,
81         0, 0, 1, 0,
82         0, 0, 0, 1
83     ]);
84     // dfmt on
85   }
86 
87   // Basic properties *****************************************************
88 
89   /**
90    * Returns i, j cell
91    */
92   @nogc T opIndex(size_t row, size_t col) pure const {
93     return this.c[row][col];
94   }
95 
96   /**
97    * Assigns a new cell value
98    */
99   @nogc void opIndexAssign(K)(K val, size_t row, size_t col) if (is(K : real)) {
100     this.c[row][col] = val;
101   }
102 
103   /**
104    * Returns i row vector
105    */
106   @nogc VRow opIndex(size_t i) pure const {
107     return this.row[i];
108   }
109 
110   /**
111    * Assigns a new row vector
112    */
113   @nogc void opIndexAssign(K)(K v, size_t i) if (isVector!(K)) {
114     static assert(v.dim <= dim, "Vector has a bigger dimension that this matrix.");
115     static if (!is(K == VRow)) {
116       this.row[i] = cast(VRow) v;
117     } else {
118       this.row[i] = v;
119     }
120   }
121 
122   // Operations ***************************************
123 
124   /**
125    * Define Equality
126    */
127   @nogc bool opEquals()(auto ref const Matrix rhs) pure const {
128     if (this.row[0] != rhs.row[0] || this.row[1] != rhs.row[1]) {
129       return false;
130     }
131     static if (dim >= 3) {
132       if (this.row[2] != rhs.row[2]) {
133         return false;
134       }
135     } else static if (dim >= 4) {
136       if (this.row[3] != rhs.row[3]) {
137         return false;
138       }
139     }
140     return true;
141   }
142 
143   /**
144   * Approximated equality
145   */
146   @nogc bool approxEqual()(auto ref const Matrix rhs, T maxRelDiff = 1e-2, T maxAbsDiff = 1e-05) pure const {
147     if (!this.row[0].approxEqual(rhs.row[0], maxRelDiff, maxAbsDiff)
148         || !this.row[1].approxEqual(rhs.row[1], maxRelDiff, maxAbsDiff)) {
149       return false;
150     }
151     static if (dim >= 3) {
152       if (!this.row[2].approxEqual(rhs.row[2], maxRelDiff, maxAbsDiff)) {
153         return false;
154       }
155     } else static if (dim >= 4) {
156       if (!this.row[3].approxEqual(rhs.row[3], maxRelDiff, maxAbsDiff)) {
157         return false;
158       }
159     }
160     return true;
161   }
162 
163   /**
164   * Define unary operators + and -
165   */
166   @nogc Matrix opUnary(string op)() pure const if (op == "+" || op == "-") {
167     import std.conv : to;
168 
169     string makeMix(size_t d, string op) {
170       string ret = "return Matrix( [";
171       foreach (i; 0 .. (d * d)) {
172         ret ~= op ~ "cell[" ~ to!string(i) ~ "] ";
173         if (i != (d * d - 1))
174           ret ~= ",";
175       }
176       ret ~= "] );";
177       return ret;
178     }
179 
180     static if (dim == 2) {
181       mixin(makeMix(2, op));
182     } else static if (dim == 3) {
183       mixin(makeMix(3, op));
184     } else static if (dim == 4) {
185       mixin(makeMix(4, op));
186     }
187   }
188 
189   /**
190   * Define binary operator + and -
191   */
192   @nogc Matrix opBinary(string op)(auto ref const Matrix rhs) pure const
193       if ((op == "+" || op == "-") && dim >= 2 && dim <= 4) {
194     import std.conv : to;
195 
196     string makeMixOp(size_t d, string op) {
197       string ret = "return Matrix( [";
198       foreach (i; 0 .. (d * d)) {
199         ret ~= "cell[" ~ to!string(i) ~ "] " ~ op ~ "rhs.cell[" ~ to!string(i) ~ "] ";
200         if (i != (d * d - 1))
201           ret ~= ",";
202       }
203       ret ~= "] );";
204       return ret;
205     }
206 
207     mixin(makeMixOp(dim, op));
208   }
209 
210   /**
211   * Define Scalar multiplication
212   */
213   @nogc Matrix opBinary(string op)(in real rhs) pure const if (op == "*" || op == "/") {
214     T[this.cells] ret = mixin("this.cell[]" ~ op ~ "cast(T) rhs");
215     return Matrix(ret);
216   }
217 
218   /**
219   * Define Scalar multiplication
220   */
221   @nogc Matrix opBinaryRight(string op)(in real rhs) pure const if (op == "*" || op == "/") {
222     T[this.cells] ret = mixin("this.cell[]" ~ op ~ "cast(T) rhs");
223     return Matrix(ret);
224   }
225 
226 
227   /**
228   * Define Matrix Product
229   */
230   @nogc Matrix opBinary(string op)(auto ref const Matrix rhs) pure const if (op == "*") {
231     Matrix mat;
232     foreach (size_t i; 0 .. dim) { // Runs along result rows
233       foreach (size_t j; 0 .. dim) { // Runs along result columns
234         static if (dim == 2) {
235           mat[i, j] = this[i, 0] * rhs[0, j] + this[i, 1] * rhs[1, j];
236         } else static if (dim == 3) {
237           mat[i, j] = this[i, 0] * rhs[0, j] + this[i, 1] * rhs[1, j] + this[i, 2] * rhs[2, j];
238         } else static if (dim == 4) {
239           mat[i, j] = this[i, 0] * rhs[0, j] + this[i, 1] * rhs[1, j] + this[i, 2] * rhs[2, j] + this[i, 3] * rhs[3, j];
240         }
241       }
242     }
243 
244     return mat;
245   }
246 
247   /**
248   * Define Matrix x Vector
249   */
250   @nogc VCol opBinary(string op)(auto ref const VRow rhs) pure const if (op == "*") {
251     static assert(rhs.dim == dim, "The vector dimension must be equal to Matrix dimmension");
252 
253     VCol ret;
254     foreach (size_t i; 0 .. dim) { // Runs along result coords
255       static if (dim == 2) {
256         ret[i] = this[i, 0] * rhs[0] + this[i, 1] * rhs[1];
257       } else static if (dim == 3) {
258         ret[i] = this[i, 0] * rhs[0] + this[i, 1] * rhs[1] + this[i, 2] * rhs[2];
259       } else static if (dim == 4) {
260         ret[i] = this[i, 0] * rhs[0] + this[i, 1] * rhs[1] + this[i, 2] * rhs[2] + this[i, 3] * rhs[3];
261       }
262     }
263     return ret;
264   }
265 
266   /**
267   * Define Vector x Matrix
268   */
269   @nogc VCol opBinaryRight(string op)(auto ref const VCol lhs) pure const if (op == "*") {
270     static assert(lhs.dim == dim, "The vector dimension must be equal to Matrix dimmension");
271     VCol ret;
272     foreach (size_t j; 0 .. dim) { // Runs along result coords
273       static if (dim == 2) {
274         ret[j] = this[0, j] * lhs[0] + this[1, j] * lhs[1];
275       } else static if (dim == 3) {
276         ret[j] = this[0, j] * lhs[0] + this[1, j] * lhs[1] + this[2, j] * lhs[2];
277       } else static if (dim == 4) {
278         ret[j] = this[0, j] * lhs[0] + this[1, j] * lhs[1] + this[2, j] * lhs[2] + this[3, j] * lhs[3];
279       }
280     }
281     return ret;
282   }
283 
284   /**
285   * Returns transposed matrix
286   *
287   * It can be used to "convert" Matrix to a row ordered matrix
288   */
289   @nogc Matrix transpose() pure const {
290     Matrix mat;
291     foreach (i; 0 .. dim) { // Runs along result rows
292       foreach (j; 0 .. dim) { // Runs along result columns
293         mat[j, i] = this[i, j];
294       }
295     }
296     return mat;
297   }
298 
299   /**
300   * Returns Determinant of this matrix
301   */
302   @nogc T determinant() pure const {
303     static if (dim == 2) {
304       return this.c[0][0] * this.c[1][1] - (this.c[1][0] * this.c[0][1]);
305     } else static if (dim == 3) {
306       T aei = this.c[0][0] * this.c[1][1] * this.c[2][2];
307       T bfg = this.c[0][1] * this.c[1][2] * this.c[2][0];
308       T cdh = this.c[0][2] * this.c[1][0] * this.c[2][1];
309       T afh = this.c[0][0] * this.c[1][2] * this.c[2][1];
310       T bdi = this.c[0][1] * this.c[1][0] * this.c[2][2];
311       T ceg = this.c[0][2] * this.c[1][1] * this.c[2][0];
312       return aei + bfg + cdh - afh - bdi - ceg;
313     } else {
314       // dfmt off
315       return (
316            this.c[0][0] * this.c[1][1] - this.c[0][1] * this.c[1][0])
317         * (this.c[2][2] * this.c[3][3] - this.c[2][3] * this.c[3][2])
318         - (this.c[0][0] * this.c[1][2] - this.c[0][2] * this.c[1][0])
319         * (this.c[2][1] * this.c[3][3] - this.c[2][3] * this.c[3][1])
320         + (this.c[0][0] * this.c[1][3] - this.c[0][3] * this.c[1][0])
321         * (this.c[2][1] * this.c[3][2] - this.c[2][2] * this.c[3][1])
322         + (this.c[0][1] * this.c[1][2] - this.c[0][2] * this.c[1][1])
323         * (this.c[2][0] * this.c[3][3] - this.c[2][3] * this.c[3][0])
324         - (this.c[0][1] * this.c[1][3] - this.c[0][3] * this.c[1][1])
325         * (this.c[2][0] * this.c[3][2] - this.c[2][2] * this.c[3][0])
326         + (this.c[0][2] * this.c[1][3] - this.c[0][3] * this.c[1][2])
327         * (this.c[2][0] * this.c[3][1] - this.c[2][1] * this.c[3][0]);
328       // dfmt on
329     }
330   }
331 
332   /**
333   * Calcs this matrix inverse
334   * Returns: A matrix full of NaNs if this matrix isn't inverible (!isOk =1)
335   */
336   @nogc Matrix inverse() pure const {
337     import std.math : approxEqual;
338 
339     Matrix mat;
340     auto det = this.determinant();
341     if (approxEqual(det, 0)) // Not invertible matrix
342       return mat; // At this point Mat it's full of NaNs and a not valid Matrix
343 
344     static if (dim == 2) {
345       mat[0, 0] = this[1, 1];
346       mat[1, 1] = this[0, 0];
347       mat[0, 1] = -this[0, 1];
348       mat[1, 0] = -this[1, 0];
349     } else if (dim == 3) {
350       mat[0, 0] = this[1, 1] * this[2, 2] - this[1, 2] * this[2, 1];
351       mat[0, 1] = this[0, 2] * this[2, 1] - this[0, 1] * this[2, 2];
352       mat[0, 2] = this[0, 1] * this[1, 2] - this[0, 2] * this[1, 1];
353 
354       mat[1, 0] = this[1, 2] * this[2, 0] - this[1, 0] * this[2, 2];
355       mat[1, 1] = this[0, 0] * this[2, 2] - this[0, 2] * this[2, 0];
356       mat[1, 2] = this[0, 2] * this[1, 0] - this[0, 0] * this[1, 2];
357 
358       mat[2, 0] = this[1, 0] * this[2, 1] - this[1, 1] * this[2, 2];
359       mat[2, 1] = this[0, 1] * this[2, 0] - this[0, 0] * this[2, 1];
360       mat[2, 2] = this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];
361     } else { // dim = 4
362       // dftm off
363       mat[0, 0] = this[1, 1] * (this[2, 2] * this[3, 3] - this[2, 3] * this[3, 2])
364         + this[1, 2] * (this[2, 3] * this[3, 1] - this[2, 1] * this[3, 3])
365         + this[1, 3] * (this[2, 1] * this[3, 2] - this[2, 2] * this[3, 1]);
366       mat[0, 1] = this[2, 1] * (this[0, 2] * this[3, 3] - this[0, 3] * this[3, 2])
367         + this[2, 2] * (this[0, 3] * this[3, 1] - this[0, 1] * this[3, 3])
368         + this[2, 3] * (this[0, 1] * this[3, 2] - this[0, 2] * this[3, 1]);
369       mat[0, 2] = this[3, 1] * (this[0, 2] * this[1, 3] - this[0, 3] * this[1, 2])
370         + this[3, 2] * (this[0, 3] * this[1, 1] - this[0, 1] * this[1, 3])
371         + this[3, 3] * (this[0, 1] * this[1, 2] - this[0, 2] * this[1, 1]);
372       mat[0, 3] = this[0, 1] * (this[1, 3] * this[2, 2] - this[1, 2] * this[2, 3])
373         + this[0, 2] * (this[1, 1] * this[2, 3] - this[1, 3] * this[2, 1])
374         + this[0, 3] * (this[1, 2] * this[2, 1] - this[1, 1] * this[2, 2]);
375 
376       mat[1, 0] = this[1, 2] * (this[2, 0] * this[3, 3] - this[2, 3] * this[3, 0])
377         + this[1, 3] * (this[2, 2] * this[3, 0] - this[2, 0] * this[3, 2])
378         + this[1, 0] * (this[2, 3] * this[3, 2] - this[2, 2] * this[3, 3]);
379       mat[1, 1] = this[2, 2] * (this[0, 0] * this[3, 3] - this[0, 3] * this[3, 0])
380         + this[2, 3] * (this[0, 2] * this[3, 0] - this[0, 0] * this[3, 2])
381         + this[2, 0] * (this[0, 3] * this[3, 2] - this[0, 2] * this[3, 3]);
382       mat[1, 2] = this[3, 2] * (this[0, 0] * this[1, 3] - this[0, 3] * this[1, 0])
383         + this[3, 3] * (this[0, 2] * this[1, 0] - this[0, 0] * this[1, 2])
384         + this[3, 0] * (this[0, 3] * this[1, 2] - this[0, 2] * this[1, 3]);
385       mat[1, 3] = this[0, 2] * (this[1, 3] * this[2, 0] - this[1, 0] * this[2, 3])
386         + this[0, 3] * (this[1, 0] * this[2, 2] - this[1, 2] * this[2, 0])
387         + this[0, 0] * (this[1, 2] * this[2, 3] - this[1, 3] * this[2, 2]);
388 
389       mat[2, 0] = this[1, 3] * (this[2, 0] * this[3, 1] - this[2, 1] * this[3, 0])
390         + this[1, 0] * (this[2, 1] * this[3, 3] - this[2, 3] * this[3, 1])
391         + this[1, 1] * (this[2, 3] * this[3, 0] - this[2, 0] * this[3, 3]);
392       mat[2, 1] = this[2, 3] * (this[0, 0] * this[3, 1] - this[0, 1] * this[3, 0])
393         + this[2, 0] * (this[0, 1] * this[3, 3] - this[0, 3] * this[3, 1])
394         + this[2, 1] * (this[0, 3] * this[3, 0] - this[0, 0] * this[3, 3]);
395       mat[2, 2] = this[3, 3] * (this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0])
396         + this[3, 0] * (this[0, 1] * this[1, 3] - this[0, 3] * this[1, 1])
397         + this[3, 1] * (this[0, 3] * this[1, 0] - this[0, 0] * this[1, 3]);
398       mat[2, 3] = this[0, 3] * (this[1, 1] * this[2, 0] - this[1, 0] * this[2, 1])
399         + this[0, 0] * (this[1, 3] * this[2, 1] - this[1, 1] * this[2, 3])
400         + this[0, 1] * (this[1, 0] * this[2, 3] - this[1, 3] * this[2, 0]);
401 
402       mat[3, 0] = this[1, 0] * (this[2, 2] * this[3, 1] - this[2, 1] * this[3, 2])
403         + this[1, 1] * (this[2, 0] * this[3, 2] - this[2, 2] * this[3, 0])
404         + this[1, 2] * (this[2, 1] * this[3, 0] - this[2, 0] * this[3, 1]);
405       mat[3, 1] = this[2, 0] * (this[0, 2] * this[3, 1] - this[0, 1] * this[3, 2])
406         + this[2, 1] * (this[0, 0] * this[3, 2] - this[0, 2] * this[3, 0])
407         + this[2, 2] * (this[0, 1] * this[3, 0] - this[0, 0] * this[3, 1]);
408       mat[3, 2] = this[3, 0] * (this[0, 2] * this[1, 1] - this[0, 1] * this[1, 2])
409         + this[3, 1] * (this[0, 0] * this[1, 2] - this[0, 2] * this[1, 0])
410         + this[3, 2] * (this[0, 1] * this[1, 0] - this[0, 0] * this[1, 1]);
411       mat[3, 3] = this[0, 0] * (this[1, 1] * this[2, 2] - this[1, 2] * this[2, 1])
412         + this[0, 1] * (this[1, 2] * this[2, 0] - this[1, 0] * this[2, 2])
413         + this[0, 2] * (this[1, 0] * this[2, 1] - this[1, 1] * this[2, 0]);
414       // dfmt on
415     }
416 
417     return mat / det;
418   }
419 
420   // Misc***********************************************************************
421 
422   /**
423   * Checks that the matrix not have a weird NaN value
424   */
425   @nogc bool isOk() pure const {
426     foreach (c; this.row) {
427       if (!c.isOk())
428         return false;
429     }
430     return true;
431   }
432 
433   /**
434   * Checks that the matrix have finite values
435   */
436   @nogc bool isFinite() pure const {
437     foreach (c; this.row) {
438       if (!c.isFinite())
439         return false;
440     }
441     return true;
442   }
443 
444   /**
445    * Return a pointer of the internal array
446    */
447   @property @nogc T* ptr() pure nothrow {
448     return this.cell.ptr;
449   }
450 
451   unittest {
452     auto m = Mat3f.IDENTITY;
453     auto pt = m.ptr;
454     pt[0] = 5;
455     assert(m[0, 0] == 5);
456     assert(pt[0] == 5);
457   }
458 
459   /**
460   * Casting operation that allow convert between Matrix types
461   */
462   @nogc Tout opCast(Tout)() pure const if (isMatrix!(Tout) && Tout.dim >= dim) {
463     static assert(isMatrix!(Tout), "This type not is a Matrix");
464     static assert(Tout.dim >= dim, "Original Matrix bigger that destiny" ~ " Vector");
465     Tout newMat;
466     auto i = 0;
467     static if (is(typeof(newMat[0, 0]) == typeof(this[0, 0]))) {
468       for (; i < dim; i++)
469         newMat[i] = this[i];
470     } else {
471       for (; i < dim; i++) //Vector Casting auto expands rows
472         newMat[i] = cast(newMat.VCol) this[i];
473     }
474 
475     // Expands a small matrix to a bigger dimension with a full 0's columns
476     for (; i < Tout.dim; i++)
477       newMat[i] = cast(newMat.VCol) newMat.VCol.ZERO;
478 
479     return newMat;
480   }
481 
482   /**
483   * Returns a visual representation of this matrix
484   */
485   string toString() const {
486     import std.conv : to;
487 
488     string ret;
489     try {
490       foreach (row; 0 .. dim) {
491         ret ~= "|" ~ this.row[row].toString();
492         /+foreach (col; 0 .. dim) {
493           ret ~= to!string(this.c[row][col]);
494           if (col < (dim - 1))
495             ret ~= ", ";
496         }
497         +/
498         ret ~= "|";
499         if (row < (dim - 1))
500           ret ~= "\n";
501       }
502     } catch (Exception ex) {
503       assert(false); // This never should happen
504     }
505     return ret;
506   }
507 }
508 
509 /**
510 * Say if a thing it's a Matrix
511 */
512 template isMatrix(T) {
513   //immutable bool isMatrix = is(T == Matrix);
514   immutable bool isMatrix = __traits(compiles, () {
515     T t;
516     static assert(T.dim >= 2 && T.dim <= 4);
517     auto cell = t.cell;
518     auto cell00 = t[0, 0];
519     auto col0 = t[0];
520     auto coln = t[T.dim];
521     t.isOk();
522 
523     // TODO : Should test for methods ?
524   });
525 }
526 
527 static assert(Mat2f.sizeof == 4 * 2*2);
528 static assert(Mat3d.sizeof == 8 * 3*3);
529 static assert(Mat4f.sizeof == 4 * 4*4);
530 
531 unittest {
532   // Fundamental Matrixes
533   auto z = Mat4d.ZERO;
534   auto ide = Mat4d.IDENTITY;
535 
536   foreach (ind; 0 .. z.cells) {
537     assert(z.cell[ind] == 0.0L);
538   }
539 
540   foreach (i; 0 .. ide.dim) {
541     foreach (j; 0 .. ide.dim) {
542       if (i == j) {
543         assert(ide[i, j] == 1.0L);
544       } else {
545         assert(ide[i, j] == 0.0L);
546       }
547     }
548   }
549 }
550 
551 unittest {
552   auto m = Mat4f.IDENTITY;
553   assert(m[0, 0] == 1);
554   assert(m[1, 1] == 1);
555   assert(m[2, 2] == 1);
556   assert(m[3, 3] == 1);
557   assert(m[3, 0] == 0);
558   m[3, 0] = 5;
559   assert(m[3, 0] == 5);
560 }
561 
562 unittest {
563   // Testing opIndexAssign of rows
564   auto m = Mat3r.IDENTITY;
565   const v = Vec2r(-1, -2);
566   const v2 = Vec3f(1, 2, 3);
567   m[1] = v; // Ha! internal vector casting
568   m[2] = v2;
569   assert(m[1, 0] == -1);
570   assert(m[1, 1] == -2);
571   assert(m[1, 2] == 0);
572   assert(m[2, 0] == 1);
573   assert(m[2, 1] == 2);
574   assert(m[2, 2] == 3);
575 }
576 
577 unittest {
578   // Check == and !=
579   auto z = Mat4r.ZERO;
580   auto ide = Mat4r.IDENTITY;
581   assert(z != ide);
582   assert(ide == ide);
583   assert(ide.approxEqual(ide));
584   assert(!ide.approxEqual(z));
585 }
586 
587 unittest {
588   // Check + - unary operators
589   auto n_ide = -Mat4r.IDENTITY;
590   foreach (i; 0 .. Mat4r.dim) {
591     foreach (j; 0 .. Mat4r.dim) {
592       if (i == j) {
593         assert(n_ide[i, j] == -1.0L);
594       } else {
595         assert(n_ide[i, j] == 0.0L);
596       }
597     }
598   }
599 }
600 
601 unittest {
602   // Check Addition and subtraction
603   // dfmt off
604   auto a = Mat4r([1, 1, 1, 1,
605                 1, 1, 1, 1,
606                 1, 1, 1, 1,
607                 1, 1, 1, 1]);
608   // dfmt on
609   auto b = a + Mat4r.IDENTITY;
610   auto c = Mat4r.ZERO - a;
611 
612   foreach (i; 0 .. Mat4r.dim) {
613     foreach (j; 0 .. Mat4r.dim) {
614       assert(c[i, j] == -1.0);
615       if (i == j) {
616         assert(b[i, j] == 2.0L);
617       } else {
618         assert(b[i, j] == 1.0L);
619       }
620     }
621   }
622 }
623 
624 unittest {
625   // Scalar multiplication
626   // dfmt off
627   auto a = Mat4f([2, 1, 1, 1,
628                 1, 2, 1, 1,
629                 1, 1, 2, 1,
630                 1, 1, 1, 2]);
631   // dfmt on
632   a = a * 2;
633   foreach (i; 0 .. Mat4f.dim) {
634     foreach (j; 0 .. Mat4f.dim) {
635       if (i == j) {
636         assert(a[i, j] == 4.0);
637       } else {
638         assert(a[i, j] == 2.0);
639       }
640     }
641   }
642 
643   // dfmt off
644   auto b = Mat4f([
645       2, 1, 1, 1,
646       1, 2, 1, 1,
647       1, 1, 2, 1,
648       1, 1, 1, 2]);
649   // dfmt on
650 
651   b = b * 2;
652   assert(a == b);
653 }
654 
655 unittest {
656   // Check Matrix multiplication - Neutral element
657   // dfmt off
658   const a = Mat4f([
659       2, 1, 1, 1,
660       1, 2, 1, 1,
661       1, 1, 2, 1,
662       1, 1, 1, 2]);
663   // dfmt on
664   const axi = a * Mat4f.IDENTITY;
665   const ixa = Mat4f.IDENTITY * a;
666   assert(a == axi, "Multiplication by identity isn't correct");
667   assert(axi == ixa);
668 
669   // dfmt off
670   const a3 = Mat3r([
671       1,4,7,
672       2,5,8,
673       3,6,9]);
674   const b3 = Mat3r([
675       9,6,3,
676       8,5,2,
677       7,4,1]);
678   const shouldBeC3 = Mat3r([
679        90, 54, 18,
680       114, 69, 24,
681       138, 84, 30]);
682   const shouldBeD3 = Mat3r([
683       30, 84, 138,
684       24, 69, 114,
685       18, 54,  90]);
686   // dfmt on
687   const c3 = a3 * b3;
688   const d3 = b3 * a3;
689   assert(c3 == shouldBeC3, "A * B should be \n" ~ shouldBeC3.toString() ~ " but was\n" ~ c3.toString());
690   assert(c3 != d3);
691   assert(d3 == shouldBeD3, "B * A should be \n" ~ shouldBeD3.toString() ~ " but was\n" ~ d3.toString());
692 }
693 
694 unittest {
695   // dfmt off
696   const mat = Mat2f([
697       0, -1,
698       1,  0]);
699   // dfmt on
700   const vec = Vec2f(1, 2);
701   const vecl = vec * mat;
702   const vecr = mat * vec;
703   assert(vec != vecl);
704   assert(vec != vecr);
705   assert(vecl != vecr);
706   assert(vecr == Vec2f(-2, 1));
707   assert(vecl == Vec2f(2, -1));
708 
709   // TODO : Do it with other type of matrix
710 }
711 
712 unittest { // Check Matrix transposition
713   auto ide = Mat4f.IDENTITY;
714   auto ide_t = ide.transpose();
715   assert(ide == ide_t);
716   auto a = Mat3r([1, 4, 7, 2, 5, 8, 3, 6, 9]);
717   auto a_t = a.transpose();
718   auto should_a_t = Mat3r([1, 2, 3, 4, 5, 6, 7, 8, 9]);
719   assert(a_t == should_a_t);
720 }
721 
722 unittest { // Check Matrix Determinant
723   auto ide = Mat4f.IDENTITY;
724   auto ide_t = ide.transpose;
725   assert(ide_t.determinant == 1);
726   assert(ide.determinant == 1);
727   // dfmt off
728   auto n = Mat4r([1, 20, 1, 0,
729                 1,  5, 1, 1,
730                 1, -1, 3, 1,
731                 1, 10, 1, 1]);
732   // dfmt on
733   auto d_n = n.determinant;
734   assert(d_n == -10);
735   auto b = Mat4r.IDENTITY;
736   b = b * 2;
737   b = n + b;
738   auto d_b = b.determinant;
739   auto prod = n * b;
740   auto d_prod = prod.determinant;
741   assert(d_n * d_b == d_prod);
742 
743   auto ide3 = Mat3f.IDENTITY;
744   assert(ide3.determinant == 1);
745   auto ide2 = Mat2f.IDENTITY;
746   assert(ide2.determinant == 1);
747 }
748 
749 unittest { // Check matrix inversion
750   // dfmt off
751   auto n = Mat4r([
752       1, 20, 1, 0,
753       1,  5, 1, 1,
754       1, -1, 3, 1,
755       1, 10, 1, 1
756   ]);
757   // dfmt on
758   auto n_inv = n.inverse;
759   // dfmt off
760   assert(n_inv == Mat4r([
761       1,  5.1, -0.5, -4.6,
762       0, -0.2,    0,  0.2,
763       0, -1.1,  0.5,  0.6,
764     -1,   -2,    0,    3
765   ]), "n_inv\n" ~ n_inv.toString());
766   // dfmt on
767   auto a = Mat3r([1, 4, 7, 2, 5, 8, 3, 6, 9]);
768   auto a_inv = a.inverse;
769   assert(!a_inv.isOk()); //not have inverse, so returns a matrix full of NaN
770   // dfmt off
771   auto m = Mat2r([
772       10, 5,
773       5, 2
774   ]);
775   // dfmt on
776   auto m_inv = m.inverse;
777   // dfmt off
778   assert (m_inv == Mat2r([
779         -.4, 1,
780         1, -2
781   ]));
782   // dfmt on
783 }
784 
785 unittest {
786   auto ide = Mat4d.IDENTITY;
787   assert(ide.isOk());
788   assert(ide.isFinite());
789   Mat3r full_nan;
790   assert(!full_nan.isOk());
791   assert(!full_nan.isFinite());
792 }
793 
794 unittest {
795   // Check casting
796   auto mat2f = Mat2f.IDENTITY;
797   mat2f[1, 0] = 5;
798   auto mat2r = cast(Mat2r) mat2f;
799   auto mat2d = cast(Mat2d) mat2f;
800   auto mat4d = cast(Mat4d) mat2r;
801   assert(is(typeof(mat2r) == Mat2r));
802   assert(is(typeof(mat2d) == Mat2d));
803   assert(is(typeof(mat4d) == Mat4d));
804 
805   assert(mat2f[0, 0] == mat4d[0, 0]);
806   assert(mat2f[0, 1] == mat4d[0, 1]);
807   assert(mat2f[1, 0] == mat4d[1, 0]);
808   assert(mat2f[1, 1] == mat4d[1, 1]);
809 }
810 
811 unittest {
812   auto ide = Mat4d.IDENTITY;
813   auto v = Vec3f.X_AXIS;
814   assert(isMatrix!(typeof(ide)));
815   assert(!isMatrix!(typeof(v)));
816 }