class Algebra::SquareMatrix

Constants

EigenSystem
Matrices

Public Class Methods

const(x) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 715
def self.const(x)
  matrix{|i, j| i == j ? ground.regulate(x) : ground.zero }
end
create(ground, n, m = nil) click to toggle source
Calls superclass method Algebra::MatrixAlgebra::create
# File lib/algebra/matrix-algebra.rb, line 731
def self.create(ground, n, m = nil)
  if m && n != m
    raise "type error to create SquareMatrix"
  end
  klass = super(ground, n, n)
  klass.sysvar(:size, n)
  klass
end
det(aa)
Alias for: determinant
determinant(aa) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 752
def self.determinant(aa)
  r = create(Integer, aa.size)
  m = r.new(aa)
  m.determinant
end
Also aliased as: det
matrices() click to toggle source
# File lib/algebra/matrix-algebra.rb, line 723
def self.matrices
  Matrices
end
regulate(x) click to toggle source
Calls superclass method Algebra::MatrixAlgebra::regulate
# File lib/algebra/matrix-algebra.rb, line 742
def self.regulate(x)
  if x.is_a? superclass #SquareMatrix
    x
  elsif y = ground.regulate(x)
    const(y)
  else
    super
  end
end
symmetric(*a) click to toggle source
# File lib/algebra/linear-algebra.rb, line 49
def self.symmetric(*a)
  k = size * (size + 1) / 2
  raise "the size of #{a.inspect} must be <= #{k}" if a.size > k
  matrix do |i, j|
    d = (i - j).abs
    a[(i < j ? i : j) + (2 * size + 1 - d) * d / 2] || ground.zero
  end
end
unity() click to toggle source
# File lib/algebra/matrix-algebra.rb, line 711
def self.unity
  matrix{|i, j| i == j ? ground.unity : ground.zero }
end

Public Instance Methods

/(other) click to toggle source
Calls superclass method Algebra::MatrixAlgebra#/
# File lib/algebra/matrix-algebra.rb, line 762
  def /(other)
    #It might be better to use regulate ...
    case other
#    when ground, Numeric
#      matrix{|i, j| self[i, j] / other}
    when self.class.superclass #SquareMatrix
      self * other.inverse
    else
      super
    end
  end
_char_matrix(mop) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 839
def _char_matrix(mop)
  # assume mop.ground.ground == self.ground
  mop.unity * mop.ground.var - self
end
char_matrix(polyring) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 844
def char_matrix(polyring)
  matrix_over_poly = Algebra.SquareMatrix(polyring, size)
  _char_matrix(matrix_over_poly)
end
char_polynomial(polyring) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 849
def char_polynomial(polyring)
  matrix_over_poly = Algebra.SquareMatrix(polyring, size)
  _char_matrix(matrix_over_poly).determinant
end
char_polynomial0(obj = "x") click to toggle source
# File lib/algebra/matrix-algebra.rb, line 831
def char_polynomial0(obj = "x")
  require "algebra/polynomial"
  a = Algebra.Polynomial(ground, obj)
  k = Algebra.SquareMatrix(a, size)
  (k.unity * a.var - self).determinant
  #    (k.unity * a.var - k.matrix{|i, j| self[i,j]}).determinant
end
const(x) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 719
def const(x)
  self.class.const(x)
end
det()
Alias for: determinant
determinant() click to toggle source
# File lib/algebra/gaussian-elimination.rb, line 313
def determinant
  if ground.field?
    determinant_by_elimination
  elsif ground.euclidian?
    determinant_by_elimination_euclidian
  else
    determinant!
  end
end
Also aliased as: determinant!, det
determinant!()
Alias for: determinant
determinant_by_elimination() click to toggle source

ground ring must be a field

# File lib/algebra/gaussian-elimination.rb, line 260
def determinant_by_elimination
  m = dup
  det = ground.unity
  each_j do |d|
    if i = (d...size).find{|i| !m[i, d].zero?}
      if i != d
        m.sswap_r!(d, i)
        det = -det
      end
      c = m[d, d]
      det *= c
      (d+1...size).each do |i0|
        r = m.row!(i0)
        q = ground.unity * m[i0, d] / c #this lets the entries be in ground
        (d+1...size).each do |j0|
          r[j0] -= q * m[d, j0]
        end
      end
    else
      return ground.zero
    end
  end
  det
end
determinant_by_elimination_euclidian() click to toggle source

ground ring must be an euclidian domain

# File lib/algebra/gaussian-elimination.rb, line 286
def determinant_by_elimination_euclidian#_new
  m = dup
  det = ground.unity
  each_j do |d|
    if i = (d...size).find{|i| !m[i, d].zero?}
      if i != d
        m.sswap_r!(d, i)
        det = -det
      end
      (d+1...size).each do |i0|
        r = m.row!(i0)
        q = m[i0, d] / m[d, d]
        (d...size).each do |j0|; r[j0] -=  q * m[d, j0]; end
        unless m[i0, d].zero?
          m.sswap_r!(d, i0)
          det = -det
          redo
        end
      end
      det *= m[d, d]
    else
      return ground.zero
    end
  end
  det
end
determinant_by_elimination_euclidian_old() click to toggle source
# File lib/algebra/gaussian-elimination.rb, line 249
def determinant_by_elimination_euclidian_old
  m = dup
  inv, k = m.left_eliminate_euclidian!
  s = ground.unity
  (0...size).each do |i|
    s *= m[i, i]
  end
  s / k
end
determinant_by_elimination_old() click to toggle source

def inverse_euclidian; left_inverse_euclidian; end

# File lib/algebra/gaussian-elimination.rb, line 239
def determinant_by_elimination_old
  m = dup
  inv, k = m.left_eliminate!
  s = ground.unity
  (0...size).each do |i|
    s *= m[i, i]
  end
  s / k
end
diagonalize() click to toggle source
# File lib/algebra/linear-algebra.rb, line 69
def diagonalize
  chp = char_polynomial(Algebra.Polynomial(ground, 't'))

  facts = chp.factorize
  mdf, mods, fas, roots, proots = chp.decompose(facts)
  nu = Algebra.SquareMatrix(mdf, size)

  espaces = {}
  evalues = []
  evectors = []
  roots0 = []
  k = 0
  facts.each do |f, _n|
    dim = f.deg
    if dim <= 0 # 0 might be occurred
      next
    elsif dim == 1
      evs = [-f[0] / f[1]]
    else
      evs = roots[k, dim] # .reverse
      roots0.push [f, evs]
    end
    k += dim

    evs.each do |evalue|
      ks = nu.const(evalue) - self
      kb = ks.kernel_basis
      espaces[evalue] = kb

      kb.each do |v|
        evalues.push evalue
        evectors.push v
      end
    end
  end

  raise "can't Diagonalize" if evalues.size < size
  ttype = Algebra.SquareMatrix(mdf, size)
  r = ttype.collect_column { |i| evectors[i] }
  #    :extfield, :roots, :tmatrix, :evalues, :addelms, :evectors,
  #                   :espaces, :chpoly, :facts
  e = EigenSystem.new(mdf, roots0, r, evalues, proots, evectors,
                      espaces, chp, facts)
  def e.to_ary
    to_a
  end
  e
end
initialzie(array, conv = false) click to toggle source
Calls superclass method
# File lib/algebra/matrix-algebra.rb, line 727
def initialzie(array, conv = false)
  super
end
inverse() click to toggle source
# File lib/algebra/matrix-algebra.rb, line 789
def inverse
  #gaussian-elimination.rb
  if ground.field? &&
      (!ground.respond_to?(:reducible) || ground.reducible)
    inverse_field
  else
    a = adjoint
    d = determinant
    unless d.unit?
      warn("#{a}/#{d} may not be done.")
    end
    a / d
  end
end
inverse_field() click to toggle source
# File lib/algebra/gaussian-elimination.rb, line 230
def inverse_field
  n, k, pi = dup.left_eliminate!
  if pi != size
    raise "Error: invert non invertible SquareMatrix"
  end
  n
end
perm(n = size, stack = []) { |y;| ... } click to toggle source
# File lib/algebra/matrix-algebra.rb, line 817
def perm(n = size, stack = [])
  (0...n).each do |x|
    unless stack.include? x
      stack.push x
      if stack.size < n
        perm(n, stack) do |y|; yield y; end
      else
        yield stack
      end
      stack.pop
    end
  end
end
sign(a) click to toggle source
# File lib/algebra/matrix-algebra.rb, line 804
def sign(a)
  n = a.size
  a = a.dup
  b = (0...n).collect{|i| a.index(i)}
  s = 1
  (0...(n-1)).each do |i|
    if (j = a[i]) != i
      a[b[i]], b[j], s = j, b[i], -s
    end
  end
  s
end
size() click to toggle source
# File lib/algebra/matrix-algebra.rb, line 740
def size; self.class.size; end
solve_eigen_value_problem() click to toggle source
# File lib/algebra/linear-algebra.rb, line 118
def solve_eigen_value_problem
  puts 'A = '; display; puts

  e = diagonalize
  puts "Charactoristic Poly.: #{e.chpoly} => #{e.facts}"; puts
  puts 'Algebraic Numbers:'
  e.roots.each do |po, rs|
    puts "#{rs.join(', ')} : roots of #{po} == 0"
  end; puts

  puts 'EigenSpaces: '
  e.evalues.uniq.each do |ev|
    puts "W_{#{ev}} = <#{e.espaces[ev].join(', ')}>"
  end; puts

  puts 'Trans. Matrix:'
  puts 'P ='
  e.tmatrix.display; puts

  puts 'P^-1 * A * P = '; (e.tmatrix.inverse * self * e.tmatrix).display; puts
end