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 }