gt:         m[0] = m[1] = m[3] = 0;
context: /*
   *   __               .__       .__  ._____.           
   * _/  |_  _______  __|__| ____ |  | |__\_ |__   ______
   * \   __\/  _ \  \/  /  |/ ___\|  | |  || __ \ /  ___/
   *  |  | (  <_> >    <|  \  \___|  |_|  || \_\ \\___ \ 
   *  |__|  \____/__/\_ \__|\___  >____/__||___  /____  >
   *                   \/       \/             \/     \/ 
   *
   * Copyright (c) 2006-2011 Karsten Schmidt
   * 
   * This library is free software; you can redistribute it and/or
   * modify it under the terms of the GNU Lesser General Public
   * License as published by the Free Software Foundation; either
   * version 2.1 of the License, or (at your option) any later version.
   * 
   * http://creativecommons.org/licenses/LGPL/2.1/
   * 
   * This library is distributed in the hope that it will be useful,
   * but WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   * Lesser General Public License for more details.
   * 
   * You should have received a copy of the GNU Lesser General Public
   * License along with this library; if not, write to the Free Software
   * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
   */
  
  package toxi.geom;
  
  import toxi.math.MathUtils;
  
  /**
   * Implements a simple row-major 4x4 matrix class, all matrix operations are
   * applied to new instances. Use {@link #transpose()} to convert from
   * column-major formats...
   */
  public class Matrix4x4 {
  
      private static final Matrix4x4 TEMP = new Matrix4x4();
  
      /**
       * Given a MxM array "matrix0", this function replaces it with the LU
       * decomposition of a row-wise permutation of itself. The input parameters
       * are "matrix0" and "dimen". The array "matrix0" is also an output
       * parameter. The vector "row_perm[4]" is an output parameter that contains
       * the row permutations resulting from partial pivoting. The output
       * parameter "even_row_xchg" is 1 when the number of row exchanges is even,
       * or -1 otherwise. Assumes data type is always double.
       * 
       * This function is similar to luDecomposition, except that it is tuned
       * specifically for 4x4 matrices.
       * 
       * Reference: Press, Flannery, Teukolsky, Vetterling,
       * _Numerical_Recipes_in_C_, Cambridge University Press, 1988, pp 40-45.
       * 
       * @param matrix0
       * @param row_perm
       * @param width
       * @return true if the matrix is nonsingular, or false otherwise.
       */
      static boolean LUDecomposition(double[] matrix0, int[] row_perm, int width) {
          double row_scale[] = new double[width];
          // Determine implicit scaling information by looping over rows
          {
              int i, j;
              int ptr, rs;
              double big, temp;
  
              ptr = 0;
              rs = 0;
  
              // For each row ...
              i = width;
              while (i-- != 0) {
                  big = 0.0;
  
                  // For each column, find the largest element in the row
                  j = width;
                  while (j-- != 0) {
                      temp = matrix0[ptr++];
                      temp = Math.abs(temp);
                      if (temp > big) {
                          big = temp;
                      }
                  }
                  // Is the matrix singular?
                  if (big == 0.0) {
                      return false;
                  }
                  row_scale[rs++] = 1.0 / big;
              }
          }
  
          {
              int j;
              int mtx = 0;
  
              // For all columns, execute Crout's method
              for (j = 0; j < width; j++) {
                  int i, imax, k;
                  int target, p1, p2;
                  double sum, big, temp;
  
                  // Determine elements of upper diagonal matrix U
                  for (i = 0; i < j; i++) {
                      target = mtx + (width * i) + j;
                      sum = matrix0[target];
                      k = i;
                      p1 = mtx + (width * i);
                      p2 = mtx + j;
                      while (k-- != 0) {
                          sum -= matrix0[p1] * matrix0[p2];
                          p1++;
                          p2 += width;
                      }
                      matrix0[target] = sum;
                  }
  
                  // Search for largest pivot element and calculate
                  // intermediate elements of lower diagonal matrix L.
                  big = 0.0;
                  imax = -1;
                  for (i = j; i < width; i++) {
                      target = mtx + (width * i) + j;
                      sum = matrix0[target];
                      k = j;
                      p1 = mtx + (width * i);
                      p2 = mtx + j;
                      while (k-- != 0) {
                          sum -= matrix0[p1] * matrix0[p2];
                          p1++;
                          p2 += width;
                      }
                      matrix0[target] = sum;
  
                      // Is this the best pivot so far?
                      if ((temp = row_scale[i] * Math.abs(sum)) >= big) {
                          big = temp;
                          imax = i;
                      }
                  }
  
                  if (imax < 0) {
                      throw new RuntimeException();
                  }
  
                  // Is a row exchange necessary?
                  if (j != imax) {
                      // Yes: exchange rows
                      k = width;
                      p1 = mtx + (width * imax);
                      p2 = mtx + (width * j);
                      while (k-- != 0) {
                          temp = matrix0[p1];
                          matrix0[p1++] = matrix0[p2];
                          matrix0[p2++] = temp;
                      }
  
                      // Record change in scale factor
                      row_scale[imax] = row_scale[j];
                  }
  
                  // Record row permutation
                  row_perm[j] = imax;
  
                  // Is the matrix singular
                  if (matrix0[(mtx + (width * j) + j)] == 0.0) {
                      return false;
                  }
  
                  // Divide elements of lower diagonal matrix L by pivot
                  if (j != width - 1) {
                      temp = 1.0 / (matrix0[(mtx + (width * j) + j)]);
                      target = mtx + (width * (j + 1)) + j;
                      i = (width - 1) - j;
                      while (i-- != 0) {
                          matrix0[target] *= temp;
                          target += width;
                      }
                  }
              }
          }
          return true;
      }
  
      public double[][] matrix;
      protected double[] temp = new double[4];
  
      public Matrix4x4() {
          init();
          matrix[0][0] = 1;
          matrix[1][1] = 1;
          matrix[2][2] = 1;
          matrix[3][3] = 1;
      }
  
      public Matrix4x4(double v11, double v12, double v13, double v14,
              double v21, double v22, double v23, double v24, double v31,
              double v32, double v33, double v34, double v41, double v42,
              double v43, double v44) {
          init();
          double[] m = matrix[0];
          m[0] = v11;
          m[1] = v12;
          m[2] = v13;
          m[3] = v14;
          m = matrix[1];
          m[0] = v21;
          m[1] = v22;
          m[2] = v23;
          m[3] = v24;
          m = matrix[2];
          m[0] = v31;
          m[1] = v32;
          m[2] = v33;
          m[3] = v34;
          m = matrix[3];
          m[0] = v41;
          m[1] = v42;
          m[2] = v43;
          m[3] = v44;
      }
  
      /**
       * Initialising constructor from a 1d array. Assumes row-major ordering
       * (column index increases faster).
       * 
       * @param array
       */
      public Matrix4x4(double[] array) {
          if (array.length != 9 && array.length != 16) {
              throw new RuntimeException("Array.length must == 9 or 16");
          }
          init();
          if (array.length == 16) {
              matrix[0][0] = array[0];
              matrix[0][1] = array[1];
              matrix[0][2] = array[2];
              matrix[0][3] = array[3];
  
              matrix[1][0] = array[4];
              matrix[1][1] = array[5];
              matrix[1][2] = array[6];
              matrix[1][3] = array[7];
  
              matrix[2][0] = array[8];
              matrix[2][1] = array[9];
              matrix[2][2] = array[10];
              matrix[2][3] = array[11];
  
              matrix[3][0] = array[12];
              matrix[3][1] = array[13];
              matrix[3][2] = array[14];
              matrix[3][3] = array[15];
          } else {
              matrix[0][0] = array[0];
              matrix[0][1] = array[1];
              matrix[0][2] = array[2];
  
              matrix[1][0] = array[3];
              matrix[1][1] = array[4];
              matrix[1][2] = array[5];
  
              matrix[2][0] = array[6];
              matrix[2][1] = array[7];
              matrix[2][2] = array[8];
  
              matrix[3][0] = array[9];
              matrix[3][1] = array[10];
              matrix[3][2] = array[11];
              matrix[3][3] = 1;
          }
      }
  
      public Matrix4x4(Matrix4x4 m) {
          init();
          for (int i = 0; i < 4; i++) {
              double[] mi = matrix[i];
              double[] mmi = m.matrix[i];
              mi[0] = mmi[0];
              mi[1] = mmi[1];
              mi[2] = mmi[2];
              mi[3] = mmi[3];
          }
      }
  
      public Matrix4x4 add(Matrix4x4 rhs) {
          Matrix4x4 result = new Matrix4x4(this);
          return result.addSelf(rhs);
      }
  
      public Matrix4x4 addSelf(Matrix4x4 m) {
          for (int i = 0; i < 4; i++) {
              double[] mi = matrix[i];
              double[] rhsm = m.matrix[i];
              mi[0] += rhsm[0];
              mi[1] += rhsm[1];
              mi[2] += rhsm[2];
              mi[3] += rhsm[3];
          }
          return this;
      }
  
      /**
       * Creates a copy of the given vector, transformed by this matrix.
       * 
       * @param v
       * @return transformed vector
       */
      public Vec3D applyTo(ReadonlyVec3D v) {
          return applyToSelf(new Vec3D(v));
      }
  
      public Vec3D applyToSelf(Vec3D v) {
          for (int i = 0; i < 4; i++) {
              double[] m = matrix[i];
              temp[i] = v.x * m[0] + v.y * m[1] + v.z * m[2] + m[3];
          }
          v.set((float) temp[0], (float) temp[1], (float) temp[2]).scaleSelf(
                  (float) (1.0 / temp[3]));
          return v;
      }
  
      public Matrix4x4 copy() {
          return new Matrix4x4(this);
      }
  
      public Matrix4x4 getInverted() {
          return new Matrix4x4(this).invert();
      }
  
      public Matrix4x4 getRotatedAroundAxis(ReadonlyVec3D axis, double theta) {
          return new Matrix4x4(this).rotateAroundAxis(axis, theta);
      }
  
      public Matrix4x4 getRotatedX(double theta) {
          return new Matrix4x4(this).rotateX(theta);
      }
  
      public Matrix4x4 getRotatedY(double theta) {
          return new Matrix4x4(this).rotateY(theta);
      }
  
      public Matrix4x4 getRotatedZ(double theta) {
          return new Matrix4x4(this).rotateZ(theta);
      }
  
      public Matrix4x4 getTransposed() {
          return new Matrix4x4(this).transpose();
      }
  
      public Matrix4x4 identity() {
          double[] m;
          m = matrix[0];
          m[1] = m[2] = m[3] = 0;
          m = matrix[1];
          m[0] = m[2] = m[3] = 0;
          m = matrix[2];
  