const std = @import("std");
const math = std.math;
const testing = std.testing;
const Vector = @import("vector.zig").Vector;

pub fn Matrix(comptime N: usize) type {
  return packed struct {
    const Self = @This();
    pub const Scalar = f32;

    values: [N][N]Scalar,

    /// Initialies a matrix.
    pub fn init(values: [N][N]Scalar) Self {
      return .{ .values = values };
    }

    /// Creates a matrix filled with the given value.
    pub fn filled(n: Scalar) Self {
      var result: Self = undefined;

      comptime var i = 0;
      inline while (i < N) : (i += 1) {
        comptime var j = 0;
        inline while (j < N) : (j += 1) {
          result.values[i][j] = n;
        }
      }

      return result;
    }

    /// A zero initialized matrix.
    pub const ZERO: Self = comptime Self.filled(0);

    /// An identity initialized matrix.
    pub const IDENTITY: Self = comptime brk: {
      var result: Self = undefined;

      comptime var i = 0;
      while (i < N) : (i += 1) {
        comptime var j = 0;
        while (j < N) : (j += 1) {
          result.values[i][j] = if (i == j) 1 else 0;
        }
      }

      break :brk result;
    };

    /// Transposes this matrix
    pub fn transpose(self: Self) Self {
      var result: Self = undefined;

      comptime var i = 0;
      inline while (i < N) : (i += 1) {
        comptime var j = 0;
        inline while (j < N) : (j += 1) {
          result.values[j][i] = self.values[i][j];
        }
      }

      return result;
    }

    pub fn invert(src: Self) ?Self {
        // https://github.com/stackgl/gl-mat4/blob/master/invert.js
        const a = @bitCast([16]f32, src.values);

        const a00 = a[0];
        const a01 = a[1];
        const a02 = a[2];
        const a03 = a[3];
        const a10 = a[4];
        const a11 = a[5];
        const a12 = a[6];
        const a13 = a[7];
        const a20 = a[8];
        const a21 = a[9];
        const a22 = a[10];
        const a23 = a[11];
        const a30 = a[12];
        const a31 = a[13];
        const a32 = a[14];
        const a33 = a[15];

        const b00 = a00 * a11 - a01 * a10;
        const b01 = a00 * a12 - a02 * a10;
        const b02 = a00 * a13 - a03 * a10;
        const b03 = a01 * a12 - a02 * a11;
        const b04 = a01 * a13 - a03 * a11;
        const b05 = a02 * a13 - a03 * a12;
        const b06 = a20 * a31 - a21 * a30;
        const b07 = a20 * a32 - a22 * a30;
        const b08 = a20 * a33 - a23 * a30;
        const b09 = a21 * a32 - a22 * a31;
        const b10 = a21 * a33 - a23 * a31;
        const b11 = a22 * a33 - a23 * a32;

        // Calculate the determinant
        var det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;

        if (std.math.approxEqAbs(f32, det, 0, 1e-8)) {
            return null;
        }
        det = 1.0 / det;

        const out = [16]f32{
            (a11 * b11 - a12 * b10 + a13 * b09) * det, // 0
            (a02 * b10 - a01 * b11 - a03 * b09) * det, // 1
            (a31 * b05 - a32 * b04 + a33 * b03) * det, // 2
            (a22 * b04 - a21 * b05 - a23 * b03) * det, // 3
            (a12 * b08 - a10 * b11 - a13 * b07) * det, // 4
            (a00 * b11 - a02 * b08 + a03 * b07) * det, // 5
            (a32 * b02 - a30 * b05 - a33 * b01) * det, // 6
            (a20 * b05 - a22 * b02 + a23 * b01) * det, // 7
            (a10 * b10 - a11 * b08 + a13 * b06) * det, // 8
            (a01 * b08 - a00 * b10 - a03 * b06) * det, // 9
            (a30 * b04 - a31 * b02 + a33 * b00) * det, // 10
            (a21 * b02 - a20 * b04 - a23 * b00) * det, // 11
            (a11 * b07 - a10 * b09 - a12 * b06) * det, // 12
            (a00 * b09 - a01 * b07 + a02 * b06) * det, // 13
            (a31 * b01 - a30 * b03 - a32 * b00) * det, // 14
            (a20 * b03 - a21 * b01 + a22 * b00) * det, // 15
        };
        return Self{
            .values = @bitCast([4][4]f32, out),
        };
    }

    /// Multiplies 2 matrices together.
    pub fn mul(self: Self, other: Self) Self {
      var result: Self = undefined;

      const a = self.transpose();
      const b = other;

      comptime var i = 0;
      inline while (i < N) : (i += 1) {
        comptime var j = 0;
        inline while (j < N) : (j += 1) {
          const row: @Vector(N, Scalar) = a.values[j];
          const column: @Vector(N, Scalar) = b.values[i];
          const products: [N]Scalar = row * column;

          var sum = @floatCast(f32, 0);
          comptime var k = 0;
          inline while (k < N) : (k += 1) {
            sum += products[k];
          }

          result.values[i][j] = sum;
        }
      }

      return result;
    }

    /// Multiplies 2 matrices together.
    pub fn mulAssign(self: *Self, other: Self) void {
      const a = self.transpose();
      const b = other;

      comptime var i = 0;
      inline while (i < N) : (i += 1) {
        comptime var j = 0;
        inline while (j < N) : (j += 1) {
          const row: @Vector(N, Scalar) = a.values[j];
          const column: @Vector(N, Scalar) = b.values[i];
          const products: [N]Scalar = row * column;

          var sum = @floatCast(f32, 0);
          comptime var k = 0;
          inline while (k < N) : (k += 1) {
            sum += products[k];
          }

          self.values[i][j] = sum;
        }
      }
    }

    const VecN = Vector(N);
    pub fn apply(self: Self, vec: VecN) VecN {
      var result: VecN = undefined;

      comptime var i = 0;
      inline while (i < N) : (i += 1) {
        result.values[i] = 0;

        comptime var j = 0;
        inline while (j < N) : (j += 1) {
          result.values[i] += self.values[i][j] * vec.values[j];
        }
      }

      return result;
    }
  };
}

fn expectMatrixEqual(comptime N: usize, expected: [N][N]f32, actual: Matrix(N)) void {
  comptime var i = 0;
  inline while (i < N) : (i += 1) {
    comptime var j = 0;
    inline while (j < N) : (j += 1) {
      testing.expectEqual(expected[i][j], actual.values[i][j]);
    }
  }
}

test "matrix initialization" {
  const A = Matrix(2).init(.{ .{ 1, 2 }, .{ 3, 4 } });
  expectMatrixEqual(2, [2][2]f32{ .{ 1, 2, }, .{ 3, 4 } }, A);
}

test "matrix zero" {
  expectMatrixEqual(3, Matrix(3).filled(0).values, Matrix(3).ZERO);
}

test "matrix identity" {
  expectMatrixEqual(3, [3][3]f32{ .{ 1, 0, 0 }, .{ 0, 1 , 0 }, .{ 0, 0, 1 } }, Matrix(3).IDENTITY);
}

test "matrix transpose" {
  const actual = Matrix(3).init(.{ .{ 1, 2, 3 }, .{ 4, 5, 6 }, .{ 7, 8, 9 } });
  const expected = [3][3]f32{ .{ 1, 4, 7 }, .{  2, 5, 8 }, .{ 3, 6, 9 } };
  expectMatrixEqual(3, expected, actual.transpose());
}

test "matrix multiplication" {
  const A = Matrix(2).init(.{ .{ 1, 2 }, .{ 3, 4 } });
  const B = Matrix(2).init(.{ .{ 6, 7 }, .{ 8, 9 } });
  const expected = [2][2]f32{ .{ 27, 40 }, .{ 35, 52 } };
  expectMatrixEqual(2, expected, A.mul(B));
}