Scala Matrix Inversion

Uh, yeah, I'd really need a quick input from someone without creator's eyes. Something's wrong in here, according to my scalacheck tests... but I don't really know enough about it to know where it's wrong.

case class Matrix(_1: (Float, Float, Float, Float), _2: (Float, Float, Float, Float),
                  _3: (Float, Float, Float, Float), _4: (Float, Float, Float, Float)) extends Immutable {
  def invert = {
    val _11 = _2._2 * _3._3 * _4._4 - _2._2 * _3._4 * _4._3 - _3._2 * _2._3 * _4._4
      +_3._2 * _2._4 * _4._3 + _4._2 * _2._3 * _3._4 - _4._2 * _2._4 * _3._3
    val _21 = -_2._1 * _3._3 * _4._4 + _2._1 * _3._4 * _4._3 + _3._1 * _2._3 * _4._4
      -_3._1 * _2._4 * _4._3 - _4._1 * _2._3 * _3._4 + _4._1 * _2._4 * _3._3
    val _31 = _2._1 * _3._2 * _4._4 - _2._1 * _3._4 * _4._2 - _3._1 * _2._2 * _4._4
      +_3._1 * _2._4 * _4._2 + _4._1 * _2._2 * _3._4 - _4._1 * _2._4 * _3._2
    val _41 = -_2._1 * _3._2 * _4._3 + _2._1 * _3._3 * _4._2 + _3._1 * _2._2 * _4._3
      -_3._1 * _2._3 * _4._2 - _4._1 * _2._2 * _3._3 + _4._1 * _2._3 * _3._2
    val _12 = -_1._2 * _3._3 * _4._4 + _1._2 * _3._4 * _4._3 + _3._2 * _1._3 * _4._4
      -_3._2 * _1._4 * _4._3 - _4._2 * _1._3 * _3._4 + _4._2 * _1._4 * _3._3
    val _22 = _1._1 * _3._3 * _4._4 - _1._1 * _3._4 * _4._3 - _3._1 * _1._3 * _4._4
      +_3._1 * _1._4 * _4._3 + _4._1 * _1._3 * _3._4 - _4._1 * _1._4 * _3._3
    val _32 = -_1._1 * _3._2 * _4._4 + _1._1 * _3._4 * _4._2 + _3._1 * _1._2 * _4._4
      -_3._1 * _1._4 * _4._2 - _4._1 * _1._2 * _3._4 + _4._1 * _1._4 * _3._2
    val _42 = _1._1 * _3._2 * _4._3 - 开发者_JS百科_1._1 * _3._3 * _4._2 - _3._1 * _1._2 * _4._3
      +_3._1 * _1._3 * _4._2 + _4._1 * _1._2 * _3._3 - _4._1 * _1._3 * _3._2
    val _13 = _1._2 * _2._3 * _4._4 - _1._2 * _2._4 * _4._3 - _2._2 * _1._3 * _4._4
      +_2._2 * _1._4 * _4._3 + _4._2 * _1._3 * _2._4 - _4._2 * _1._4 * _2._3
    val _23 = -_1._1 * _2._3 * _4._4 + _1._1 * _2._4 * _4._3 + _2._1 * _1._3 * _4._4
      -_2._1 * _1._4 * _4._3 - _4._1 * _1._3 * _2._4 + _4._1 * _1._4 * _2._3
    val _33 = _1._1 * _2._2 * _4._4 - _1._1 * _2._4 * _4._2 - _2._1 * _1._2 * _4._4
      +_2._1 * _1._4 * _4._2 + _4._1 * _1._2 * _2._4 - _4._1 * _1._4 * _2._2
    val _43 = -_1._1 * _2._2 * _4._3 + _1._1 * _2._3 * _4._2 + _2._1 * _1._2 * _4._3
      -_2._1 * _1._3 * _4._2 - _4._1 * _1._2 * _2._3 + _4._1 * _1._3 * _2._2
    val _14 = -_1._2 * _2._3 * _3._4 + _1._2 * _2._4 * _3._3 + _2._2 * _1._3 * _3._4
      -_2._2 * _1._4 * _3._3 - _3._2 * _1._3 * _2._4 + _3._2 * _1._4 * _2._3
    val _24 = _1._1 * _2._3 * _3._4 - _1._1 * _2._4 * _3._3 - _2._1 * _1._3 * _3._4
      +_2._1 * _1._4 * _3._3 + _3._1 * _1._3 * _2._4 - _3._1 * _1._4 * _2._3
    val _34 = -_1._1 * _2._2 * _3._4 + _1._1 * _2._4 * _3._2 + _2._1 * _1._2 * _3._4
      -_2._1 * _1._4 * _3._2 - _3._1 * _1._2 * _2._4 + _3._1 * _1._4 * _2._2
    val _44 = _1._1 * _2._2 * _3._3 - _1._1 * _2._3 * _3._2 - _2._1 * _1._2 * _3._3
      +_2._1 * _1._3 * _3._2 + _3._1 * _1._2 * _2._3 - _3._1 * _1._3 * _2._2

    val det = _1._1 * _11 + _1._2 * _21 + _1._3 * _31 + _1._4 * _41
    if (det == 0) this
    else Matrix(
      (_11, _12, _13, _14),
      (_21, _22, _23, _24),
      (_31, _32, _33, _34),
      (_41, _42, _43, _44)
    ) * (1 / det)
  }

  def *(f: Float) = Matrix(
    (_1._1 * f, _1._2 * f, _1._3 * f, _1._4 * f),
    (_2._1 * f, _2._2 * f, _2._3 * f, _2._4 * f),
    (_3._1 * f, _3._2 * f, _3._3 * f, _3._4 * f),
    (_4._1 * f, _4._2 * f, _4._3 * f, _4._4 * f)
  )
}

Also, can I load this Matrix into OpenGL or do I have to transpose it first. I really always get confused about this maths.

---

Inverting a matrix is usually a bad idea, because the calculations can be ill-conditioned.

If you want to solve a system of equations it's a better idea to do it using something like LU decomposition and forward-backward substitution, especially if you can reuse the decomposition to solve for several right hand side vectors.

This link shows a Java example for Gaussian elimination with pivoting.

Here's another thought: maybe you can just use Java libraries like Apache Commons Math, the successor to JAMA, in your application?

If you have a particular case in mind, I'd recommend entering it into Wolfram Alpha so you can see what the answer should be before you start coding.

---

I'm pretty sure Simplex3D implements this calculation (and very likely it's done correctly there).

---

If you want to play with numerics - by all means go ahead and do it yourself - you have some good suggestions from Jesper and duffymo (inverting matrices is not useful in practice - look into LU decomposition).

If however if you just want to Get Stuff Done^TM look into Scalala and Scalalab.

Either way you will need Linear Algebra background knowledge - which is incredibly useful math for lots of fields.

---

Look at Invertible matrix: Analytic solution on Wikipedia. The whole bunch of calculations at the top compute the adjugate of the matrix, from which the determinant is calculated, and the inverse is then 1 / det times the adjugate matrix.

The whole calculation is written out explicitly for a 4 x 4 matrix in your code, so if there's a bug in it will take some effort to check the whole thing. The Wikipedia articles explain how it's supposed to work.