class Matrix

def *(m) # m is matrix or vector or number
  case(m)
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e * m }
    }
    return new_matrix(rows, column_size)
  when Vector
    m = self.class.column_vector(m)
    r = self * m
    return r.column(0)
  when Matrix
    Matrix.Raise ErrDimensionMismatch if column_size != m.row_size
    rows = Array.new(row_size) {|i|
      Array.new(m.column_size) {|j|
        (0 ... column_size).inject(0) do |vij, k|
          vij + self[i, k] * m[k, j]
        end
      }
    }
    return new_matrix(rows, m.column_size)
  else
    return apply_through_coercion(m, __method__)
  end
end

def **(other)
  case other
  when Integer
    x = self
    if other <= 0
      x = self.inverse
      return self.class.identity(self.column_size) if other == 0
      other = -other
    end
    z = nil
    loop do
      z = z ? z * x : x if other[0] == 1
      return z if (other >>= 1).zero?
      x *= x
    end
  when Numeric
    v, d, v_inv = eigensystem
    v * self.class.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
  else
    Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class
  end
end

def +(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class
  when Vector
    m = self.class.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end
  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
  rows = Array.new(row_size) {|i|
    Array.new(column_size) {|j|
      self[i, j] + m[i, j]
    }
  }
  new_matrix rows, column_size
end

def -(m)
  case m
  when Numeric
    Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class
  when Vector
    m = self.class.column_vector(m)
  when Matrix
  else
    return apply_through_coercion(m, __method__)
  end
  Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
  rows = Array.new(row_size) {|i|
    Array.new(column_size) {|j|
      self[i, j] - m[i, j]
    }
  }
  new_matrix rows, column_size
end

def /(other)
  case other
  when Numeric
    rows = @rows.collect {|row|
      row.collect {|e| e / other }
    }
    return new_matrix(rows, column_size)
  when Matrix
    return self * other.inverse
  else
    return apply_through_coercion(other, __method__)
  end
end

def ==(other)
  return false unless Matrix === other &&
                      column_size == other.column_size # necessary for empty matrices
  rows == other.rows
end

def [](i, j)
  @rows.fetch(i){return nil}[j]
end

def []=(i, j, v)
  @rows[i][j] = v
end

def clone
  new_matrix @rows.map(&:dup), column_size
end

def coerce(other)
  case other
  when Numeric
    return Scalar.new(other), self
  else
    raise TypeError, "#{self.class} can't be coerced into #{other.class}"
  end
end

def collect(&block) # :yield: e
  return to_enum(:collect) unless block_given?
  rows = @rows.collect{|row| row.collect(&block)}
  new_matrix rows, column_size
end

def column(j) # :yield: e
  if block_given?
    return self if j >= column_size || j < -column_size
    row_size.times do |i|
      yield @rows[i][j]
    end
    self
  else
    return nil if j >= column_size || j < -column_size
    col = Array.new(row_size) {|i|
      @rows[i][j]
    }
    Vector.elements(col, false)
  end
end

def column_vectors
  Array.new(column_size) {|i|
    column(i)
  }
end

def conjugate
  collect(&:conjugate)
end

def determinant
  Matrix.Raise ErrDimensionMismatch unless square?
  m = @rows
  case row_size
    # Up to 4x4, give result using Laplacian expansion by minors.
    # This will typically be faster, as well as giving good results
    # in case of Floats
  when 0
    +1
  when 1
    + m[0][0]
  when 2
    + m[0][0] * m[1][1] - m[0][1] * m[1][0]
  when 3
    m0, m1, m2 = m
    + m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
    - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
    + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
  when 4
    m0, m1, m2, m3 = m
    + m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
    - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
    + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
    - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
    + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
    - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
    + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
    - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
    + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
    - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
    + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
    - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
  else
    # For bigger matrices, use an efficient and general algorithm.
    # Currently, we use the Gauss-Bareiss algorithm
    determinant_bareiss
  end
end

def determinant_bareiss
  size = row_size
  last = size - 1
  a = to_a
  no_pivot = Proc.new{ return 0 }
  sign = +1
  pivot = 1
  size.times do |k|
    previous_pivot = pivot
    if (pivot = a[k][k]) == 0
      switch = (k+1 ... size).find(no_pivot) {|row|
        a[row][k] != 0
      }
      a[switch], a[k] = a[k], a[switch]
      pivot = a[k][k]
      sign = -sign
    end
    (k+1).upto(last) do |i|
      ai = a[i]
      (k+1).upto(last) do |j|
        ai[j] =  (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
      end
    end
  end
  sign * pivot
end

def determinant_e
  warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
  rank
end

def diagonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  each(:off_diagonal).all?(&:zero?)
end

def each(which = :all) # :yield: e
  return to_enum(:each, which) unless block_given?
  last = column_size - 1
  case which
  when :all
    block = Proc.new
    @rows.each do |row|
      row.each(&block)
    end
  when :diagonal
    @rows.each_with_index do |row, row_index|
      yield row.fetch(row_index){return self}
    end
  when :off_diagonal
    @rows.each_with_index do |row, row_index|
      column_size.times do |col_index|
        yield row[col_index] unless row_index == col_index
      end
    end
  when :lower
    @rows.each_with_index do |row, row_index|
      0.upto([row_index, last].min) do |col_index|
        yield row[col_index]
      end
    end
  when :strict_lower
    @rows.each_with_index do |row, row_index|
      [row_index, column_size].min.times do |col_index|
        yield row[col_index]
      end
    end
  when :strict_upper
    @rows.each_with_index do |row, row_index|
      (row_index+1).upto(last) do |col_index|
        yield row[col_index]
      end
    end
  when :upper
    @rows.each_with_index do |row, row_index|
      row_index.upto(last) do |col_index|
        yield row[col_index]
      end
    end
  else
    Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end

def each_with_index(which = :all) # :yield: e, row, column
  return to_enum(:each_with_index, which) unless block_given?
  last = column_size - 1
  case which
  when :all
    @rows.each_with_index do |row, row_index|
      row.each_with_index do |e, col_index|
        yield e, row_index, col_index
      end
    end
  when :diagonal
    @rows.each_with_index do |row, row_index|
      yield row.fetch(row_index){return self}, row_index, row_index
    end
  when :off_diagonal
    @rows.each_with_index do |row, row_index|
      column_size.times do |col_index|
        yield row[col_index], row_index, col_index unless row_index == col_index
      end
    end
  when :lower
    @rows.each_with_index do |row, row_index|
      0.upto([row_index, last].min) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :strict_lower
    @rows.each_with_index do |row, row_index|
      [row_index, column_size].min.times do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :strict_upper
    @rows.each_with_index do |row, row_index|
      (row_index+1).upto(last) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  when :upper
    @rows.each_with_index do |row, row_index|
      row_index.upto(last) do |col_index|
        yield row[col_index], row_index, col_index
      end
    end
  else
    Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
  end
  self
end

def eigensystem
  EigenvalueDecomposition.new(self)
end

def elements_to_f
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_f is deprecated, use map(&:to_f)"
  map(&:to_f)
end

def elements_to_i
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_i is deprecated, use map(&:to_i)"
  map(&:to_i)
end

def elements_to_r
  warn "#{caller(1)[0]}: warning: Matrix#elements_to_r is deprecated, use map(&:to_r)"
  map(&:to_r)
end

def empty?
  column_size == 0 || row_size == 0
end

def eql?(other)
  return false unless Matrix === other &&
                      column_size == other.column_size # necessary for empty matrices
  rows.eql? other.rows
end

def hash
  @rows.hash
end

def hermitian?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper).all? do |e, row, col|
    e == rows[col][row].conj
  end
end if 42.respond_to?(:conj)

def imaginary
  collect(&:imaginary)
end

def index(*args)
  raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2
  which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all
  return to_enum(:find_index, which, *args) unless block_given? || args.size == 1
  if args.size == 1
    value = args.first
    each_with_index(which) do |e, row_index, col_index|
      return row_index, col_index if e == value
    end
  else
    each_with_index(which) do |e, row_index, col_index|
      return row_index, col_index if yield e
    end
  end
  nil
end

def initialize(rows, column_size = rows[0].size)
  # No checking is done at this point. rows must be an Array of Arrays.
  # column_size must be the size of the first row, if there is one,
  # otherwise it *must* be specified and can be any integer >= 0
  @rows = rows
  @column_size = column_size
end

def inspect
  if empty?
    "#{self.class}.empty(#{row_size}, #{column_size})"
  else
    "#{self.class}#{@rows.inspect}"
  end
end

def inverse
  Matrix.Raise ErrDimensionMismatch unless square?
  self.class.I(row_size).send(:inverse_from, self)
end

def inverse_from(src) # :nodoc:
  last = row_size - 1
  a = src.to_a
  0.upto(last) do |k|
    i = k
    akk = a[k][k].abs
    (k+1).upto(last) do |j|
      v = a[j][k].abs
      if v > akk
        i = j
        akk = v
      end
    end
    Matrix.Raise ErrNotRegular if akk == 0
    if i != k
      a[i], a[k] = a[k], a[i]
      @rows[i], @rows[k] = @rows[k], @rows[i]
    end
    akk = a[k][k]
    0.upto(last) do |ii|
      next if ii == k
      q = a[ii][k].quo(akk)
      a[ii][k] = 0
      (k + 1).upto(last) do |j|
        a[ii][j] -= a[k][j] * q
      end
      0.upto(last) do |j|
        @rows[ii][j] -= @rows[k][j] * q
      end
    end
    (k+1).upto(last) do |j|
      a[k][j] = a[k][j].quo(akk)
    end
    0.upto(last) do |j|
      @rows[k][j] = @rows[k][j].quo(akk)
    end
  end
  self
end

def lower_triangular?
  each(:strict_upper).all?(&:zero?)
end

def lup
  LUPDecomposition.new(self)
end

def minor(*param)
  case param.size
  when 2
    row_range, col_range = param
    from_row = row_range.first
    from_row += row_size if from_row < 0
    to_row = row_range.end
    to_row += row_size if to_row < 0
    to_row += 1 unless row_range.exclude_end?
    size_row = to_row - from_row
    from_col = col_range.first
    from_col += column_size if from_col < 0
    to_col = col_range.end
    to_col += column_size if to_col < 0
    to_col += 1 unless col_range.exclude_end?
    size_col = to_col - from_col
  when 4
    from_row, size_row, from_col, size_col = param
    return nil if size_row < 0 || size_col < 0
    from_row += row_size if from_row < 0
    from_col += column_size if from_col < 0
  else
    Matrix.Raise ArgumentError, param.inspect
  end
  return nil if from_row > row_size || from_col > column_size || from_row < 0 || from_col < 0
  rows = @rows[from_row, size_row].collect{|row|
    row[from_col, size_col]
  }
  new_matrix rows, [column_size - from_col, size_col].min
end

def new_matrix(rows, column_size = rows[0].size) # :nodoc:
  self.class.send(:new, rows, column_size) # bypass privacy of Matrix.new
end

def normal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row_i, i|
    rows.each_with_index do |row_j, j|
      s = 0
      rows.each_with_index do |row_k, k|
        s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j]
      end
      return false unless s == 0
    end
  end
  true
end if 42.respond_to?(:conj)

def orthogonal?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_size.times do |j|
      s = 0
      row_size.times do |k|
        s += row[k] * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end

def permutation?
  Matrix.Raise ErrDimensionMismatch unless square?
  cols = Array.new(column_size)
  rows.each_with_index do |row, i|
    found = false
    row.each_with_index do |e, j|
      if e == 1
        return false if found || cols[j]
        found = cols[j] = true
      elsif e != 0
        return false
      end
    end
    return false unless found
  end
  true
end

def rank
  # We currently use Bareiss' multistep integer-preserving gaussian elimination
  # (see comments on determinant)
  a = to_a
  last_column = column_size - 1
  last_row = row_size - 1
  pivot_row = 0
  previous_pivot = 1
  0.upto(last_column) do |k|
    switch_row = (pivot_row .. last_row).find {|row|
      a[row][k] != 0
    }
    if switch_row
      a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
      pivot = a[pivot_row][k]
      (pivot_row+1).upto(last_row) do |i|
         ai = a[i]
         (k+1).upto(last_column) do |j|
           ai[j] =  (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
         end
       end
      pivot_row += 1
      previous_pivot = pivot
    end
  end
  pivot_row
end

def rank_e
  warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
  rank
end

def real
  collect(&:real)
end

def real?
  all?(&:real?)
end

def rect
  [real, imag]
end

def regular?
  not singular?
end

def round(ndigits=0)
  map{|e| e.round(ndigits)}
end

def row(i, &block) # :yield: e
  if block_given?
    @rows.fetch(i){return self}.each(&block)
    self
  else
    Vector.elements(@rows.fetch(i){return nil})
  end
end

def row_size
  @rows.size
end

def row_vectors
  Array.new(row_size) {|i|
    row(i)
  }
end

def singular?
  determinant == 0
end

def square?
  column_size == row_size
end

def symmetric?
  Matrix.Raise ErrDimensionMismatch unless square?
  each_with_index(:strict_upper) do |e, row, col|
    return false if e != rows[col][row]
  end
  true
end

def to_a
  @rows.collect(&:dup)
end

def to_s
  if empty?
    "#{self.class}.empty(#{row_size}, #{column_size})"
  else
    "#{self.class}[" + @rows.collect{|row|
      "[" + row.collect{|e| e.to_s}.join(", ") + "]"
    }.join(", ")+"]"
  end
end

def trace
  Matrix.Raise ErrDimensionMismatch unless square?
  (0...column_size).inject(0) do |tr, i|
    tr + @rows[i][i]
  end
end

def transpose
  return self.class.empty(column_size, 0) if row_size.zero?
  new_matrix @rows.transpose, row_size
end

def unitary?
  Matrix.Raise ErrDimensionMismatch unless square?
  rows.each_with_index do |row, i|
    column_size.times do |j|
      s = 0
      row_size.times do |k|
        s += row[k].conj * rows[k][j]
      end
      return false unless s == (i == j ? 1 : 0)
    end
  end
  true
end if 42.respond_to?(:conj)

def upper_triangular?
  each(:strict_lower).all?(&:zero?)
end

def zero?
  all?(&:zero?)
end