# vector.rb version 0.1, 2002/4/21
# simple vector caliculator
# by Tomokiyo Nomura
#
# Usage: this program is supposed to be used on the irb with --noinspect option.
# after irb prompt appeared, type 'require "vector.rb"'. To create a Matrix,
# type as follows.
#     irb> a = Matrix[[1,2],[3,4]]
# or to create a vector,
#     irb> b = Vector[[1,2]]
# in this case, try 'a + a', 'a * a', 'a * 2', 'b + b', 'b * 2', 'a * b'
# on the irb prompt. To calculate inner_product, use Vector.inner_product
# class method. ex,
#     irb> Vector.inner_product( b, b )
# To transpose matrix or vector, use transpose method.
#     irb> a.transpose
# To convert a vector to the matrix, use to_matrix method.
#     irb> b.to_matrix
#
# Each row or column of the matrix can be processed, to exchange rows,
#     irb> a.row.exchange(1,2)
# The following methods are usable for row or column calculation.
#     row_exchange!(i, j)                has 'ce!' shortcut 
#     row_add!( i, j, scalar )           has 'ca!' shortcut
#     row_subtract!( i, j, scalar )      has 'cs!' shortcut
#     row_multiply!( i, scalar )         has 'cm!' shortcut
#     row_divide!( i, scalar )           has 'cd!' shortcut
#     column_exchange!                   has 'ce!' shortcut 
#     column_add!( i, j, scalar )        has 'ca!' shortcut
#     column_subtract!( i, j, scalar )   has 'cs!' shortcut
#     column_multiply!( i, scalar )      has 'cm!' shortcut
#     column_divide!( i, scalar )        has 'cd!' shortcut

class Matrix
  def initialize( a )
    @a = a
  end
  attr_reader :a
  def to_s
    b = "["
    @a.each {|c| b += c.inspect + "\n"}
    b.chop + "]"
  end

  def Matrix.[]( *a )
    Matrix.new( a )
  end

  def Matrix.dim( row, column )
    (0...row).collect { Array.new(column).fill(0) }
  end
  def Matrix.add( x, y )
    row, column = x.size, x[0].size
    c = Matrix.dim(row, column)
    for i in 0...row
      for j in 0...column
        c[i][j] = x[i][j] + y[i][j]
      end
    end
    c
  end
  def Matrix.sub( x, y )
    row, column = x.size, x[0].size
    c = Matrix.dim(row, column)
    for i in 0...row
      for j in 0...column
        c[i][j] = x[i][j] - y[i][j]
      end
    end
    c
  end
  def Matrix.multiply( x, y )
    x_row, x_column = x.size, x[0].size
    y_column = y[0].size
    c = Matrix.dim(x_row, y_column)
    for i in 0...x_row
      for j in 0...y_column
        for k in 0...x_column
          c[i][j] += x[i][k] * y[k][j]
        end
      end
    end
    c
  end
  def Matrix.multiply_by_scalar( x, a )
    x_row, x_column = x.size, x[0].size
    c = Matrix.dim(x_row, x_column)
    for i in 0...x_row
      for j in 0...x_column
        c[i][j] = x[i][j] * a
      end
    end
    c
  end
  def Matrix.transpose( arry )
    row, column = arry.size, arry[0].size
    c = Matrix.dim( column, row )
    for i in 0...row
      for j in 0...column
        c[j][i] = arry[i][j]
      end
    end
    c
  end

  def +( other )
    Matrix.new( Matrix.add( self.a, other.a ) )
  end
  def -( other )
    Matrix.new( Matrix.sub( self.a, other.a ) )
  end
  def *( other )
    if other.kind_of?( Numeric )
      Matrix.new( Matrix.multiply_by_scalar( self.a, other ) )
    elsif other.type == Vector
      Vector.new( Matrix.multiply( self.a, other.a ) )
    else
      Matrix.new( Matrix.multiply( self.a, other.a ) )
    end
  end

  def row_exchange!( i, j )
    i -= 1; j -= 1
    @a[i], @a[j] = @a[j], @a[i]
    self
  end
  def row_add!( i, j, c = 1 )
    i -=1; j -= 1
    @a[i].each_index {|k| @a[i][k] += @a[j][k] * c}
    self
  end
  def row_subtract!( i, j, c = 1 )
    i -=1; j -= 1
    @a[i].each_index {|k| @a[i][k] -= @a[j][k] * c}
    self
  end
  def row_multiply!( i, c )
    i -=1; column = @a.size
    (0...column).each {|k| @a[i][k] *= c}
    self
  end
  def row_divide!( i, c )
    i -=1; column = @a[0].size
    (0...column).each {|k| @a[i][k] /= c.to_f}
    self
  end
  def transpose
    Matrix.new( Matrix.transpose( @a ) )
  end

  alias re! row_exchange!
  alias ra! row_add!
  alias rs! row_subtract!
  alias rm! row_multiply!
  alias rd! row_divide! 

  def column_exchange!( i, j )
    row = @a.size
    i -= 1; j -= 1
    for k in 0...row
      @a[k][i], @a[k][j] = @a[k][j], @a[k][i]
    end
    self
  end
  def column_add!( i, j, c = 1 )
    row = @a.size
    i -= 1; j -= 1
    for k in 0...row
      @a[k][i] += @a[k][j] * c
    end
    self
  end
  def column_subtract!( i, j, c = 1)
    row = @a.size
    i -= 1; j -= 1
    for k in 0...row
      @a[k][i] -= @a[k][j] * c
    end
    self
  end
  def column_multiply!( i, c )
    row = @a.size
    i -= 1
    for k in 0...row
      @a[k][i] *= c
    end
    self
  end
  def column_divide!( i, c )
    row = @a.size
    i -= 1
    for k in 0...row
      @a[k][i] /= c.to_f
    end
    self
  end

  alias ce! column_exchange!
  alias ca! column_add!
  alias cs! column_subtract!
  alias cm! column_multiply!
  alias cd! column_divide!

end

class Vector
  def initialize( a )
    @a = a
  end
  attr_reader :a
  def to_s
    @a.collect {|c| c[0]}.inspect
  end

  def Vector.[]( *a )
    Vector.new( a.collect{|c| c.to_a} )
  end

  def +( other )
    Vector.new( Matrix.add( self.a, other.a ) )
  end
  def -( other )
    Vector.new( Matrix.sub( self.a, other.a ) )
  end
  def *( other )
    if other.kind_of?( Numeric )
      Vector.new( Matrix.multiply_by_scalar( self.a, other ) )
    else
      Vector.new( Matrix.multiply( self.a, other.a ) )
    end
  end
  def to_matrix
    Matrix.new( @a )
  end

  def Vector.inner_product( x, y )
    p = 0
    x.a.each_index {|i| p += x.a[i][0] * y.a[i][0]}
    p
  end
end

# 'gaussjordan'
# this part of tha program describes methods to calculate inverse matrix
# and to solve simultaneous linear equation.
# Usage: to get the inverse matrix type a.inv or a.inverse.
# to solve the simultaneous linea equation type Matrix.solver( matrix, vector ).


class Matrix
  def Matrix.gaussjordan( arry )
    matrix = []
    arry.each{|c| matrix.push( c.dup )}
    n = matrix.size
    m = matrix[0].size
    for k in 0...n
      for i in k.succ...n
        if matrix[i][k].to_f.abs > matrix[k][k].to_f.abs
          matrix[i], matrix[k] = matrix[k], matrix[i]
        end
      end
      coeff = matrix[k][k].to_f
      for j in k.succ...m
        matrix[k][j] /= coeff
      end
      for i in 0...n
        if i == k
          next
        else
          coeff = matrix[i][k].to_f
          for j in k.succ...m
            matrix[i][j] -= matrix[k][j] * coeff
          end
        end
      end
    end
    matrix.collect{|c| c[n...m]}
  end
  def Matrix.unit( n )
    matrix = Matrix.dim( n, n )
    for i in 0...n
      matrix[i][i] = 1
    end
    matrix
  end

  def Matrix.solver( matrix, vector )
    arry = []
    matrix.a.each{|c| arry.push( c.dup )}
    v = vector.a
    v.each_index {|i| arry[i].push( v[i][0] )}
    Vector.new( Matrix.gaussjordan( arry ) )
  end
  def inverse
    n = @a.size
    unit = Matrix.unit( n )
    arry = @a.collect{|c| c + unit.shift}
    Matrix.new( Matrix.gaussjordan( arry ) )
  end
  alias inv inverse
end

# 'determinant'
# this part of the program describes the method to calculate the determinant
# of the matrix.
# Usage: Type 'a.det' or 'a.determinant' to get the determinant of the matrix
# 'a'.

class Matrix
  def determinant
    matrix = []
    @a.each{|c| matrix.push( c.dup )}
    n = matrix.size
    det = 1

    for k in 0...n
      if matrix[k][k] == 0
        for i in k.succ...n
          if matrix[i][k] != 0
            matrix[i], matrix[k] = matrix[k], matrix[i]
            det *= -1
            break
          end
        end
      end
      akk = matrix[k][k].to_f
      if akk == 0
        return 0
      else
        det *= akk
      end
      for i in k.succ...n
        coeff = matrix[i][k] / akk
        for j in k.succ...n
          matrix[i][j] -= matrix[k][j] * coeff
        end
      end
    end
    det
  end
  
  alias det determinant
end

# 'jacobi'
# This part of the program describes the method to get the eigenvalue and the
# eigenvector of the symmetric matrix.
# Usage: Type 'a.eva' or 'a.eigenvalue' to get the list of the eigenvalue of the
# symmetric matrix 'a'. Type 'a.evct' or 'a.eigenvector' to get the list of
# the eigenvectors of the symmetric matrix 'a'.

class Matrix
  Iterate_Max = 50
  Epsilon = 1e-15
  def Matrix.jacobi( arry )
    a = arry.collect{|c| c.dup }
    n = a.size
    w = (0...n).collect{ Array.new(n).fill(0) }
    w.each_index{|i| w[i][i] = 1}
    for iterate in 0..Iterate_Max
      amax = 0
      for i in 0...n
        for j in i.succ...n
          if amax < a[i][j].to_f.abs
            amax = a[i][j].to_f.abs
            imax = i
            jmax = j
          end
        end
      end
      if amax < Epsilon
        break
      end
      alpha = (a[imax][imax] - a[jmax][jmax])/2.0
      beta = Math.sqrt( alpha * alpha + amax * amax )
      sgn = (a[imax][imax] - a[jmax][jmax]) <=> 0
      cos = Math.sqrt( ( 1 + sgn * alpha / beta) / 2 )
      sin = sgn * amax / 2.0 / beta / cos
      for j in 0...n
        x, y = a[imax][j], a[jmax][j]
        a[imax][j] = x * cos + y * sin
        a[jmax][j] = -1 * x * sin + y * cos
      end
      for i in 0...n
        x, y = a[i][imax], a[i][jmax]
        a[i][imax] = x * cos + y * sin
        a[i][jmax] = -1 * x * sin + y * cos
      end
      for j in 0...n
        x, y = w[imax][j], w[jmax][j]
        w[imax][j] = x * cos + y * sin
        w[jmax][j] = -1 * x * sin + y * cos
      end
    end
    [ a, w ]
  end

  def eigenvalue
    a, w = Matrix.jacobi( @a )
    n = a.size
    (0...n).collect{|i| a[i][i] }
  end
  def eigenvector
    a, w = Matrix.jacobi( @a )
    n = w.size
    ret_val = []
    for i in 0...n
      v = (0...n).collect{ w[i].shift.to_a }
      ret_val.push( Vector.new( v ) )
    end
    ret_val
  end
  alias evl eigenvalue
  alias evct eigenvector
end