module Algebra::ElementaryDivisor

Public Class Methods

factorize(eds) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 163
def self.factorize(eds)
  prev = nil
  qeds = []
  eds.each do |f|
    qeds << if prev
              f / prev
            else
              f
            end
    prev = f
  end

  qeds_facts = qeds.collect(&:factorize)
  result = []
  prev = nil
  qeds_facts.each do |fac|
    prev = if prev
             prev * fac
           else
             fac
           end
    result << prev
  end
  result
end

Public Instance Methods

_coeff(sdd) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 102
def _coeff(sdd)
  if ground <= Algebra::Polynomial
    sdd.lc
  elsif ground.field?
    sdd
  elsif ground.euclidian?
    if ground <= Integer && sdd < 0
      - ground.unity
    elsif sdd.unit?
      sdd
    else
      ground.unity
    end
  else
    ground.unity
  end
end
e_diagonalize() click to toggle source
# File lib/algebra/elementary-divisor.rb, line 56
def e_diagonalize
  dup.e_diagonalize!
end
e_diagonalize!() click to toggle source
# File lib/algebra/elementary-divisor.rb, line 60
def e_diagonalize!
  size = rsize < csize ? rsize : csize
  0.upto size - 1 do |d|
    el_sweep!(d)
  end
  self
end
e_inverse() click to toggle source
# File lib/algebra/elementary-divisor.rb, line 40
def e_inverse
  raise 'not square matrix' unless rsize == csize
  size = rsize
  edm = to_triplet.e_diagonalize
  if edm.body[size - 1, size - 1].unity?
    edm.right * edm.left
  else
    puts 'elementary divisor matrix:'; edm.body.display
    raise 'not invertible'
  end
end
el_sweep!(d) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 68
def el_sweep!(d)
  loop do
    1 while horizontal_sweep!(d) || vertical_sweep!(d)
    if idx = sq_find_nondiv(d + 1, self[d, d])
      i, _j = idx
      mix_r!(d, i)
    else
      break
    end
  end
  unless self[d, d].zero?
    #       row!(d)[d] = self[d, d].monic
    divide_r!(d, _coeff(self[d, d]))
  end
  self
end
el_sweep_old!(d) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 85
def el_sweep_old!(d)
  loop do
    sr = vertical_sweep!(d)
    sc = horizontal_sweep!(d)
    break unless sr || sc
    next unless idx = sq_find_nondiv(d + 1, self[d, d])
    i, _j = idx
    mix_r!(d, i)
    redo
  end
  unless self[d, d].zero?
    #       row!(d)[d] = self[d, d].monic
    divide_r!(d, _coeff(self[d, d]))
  end
  self
end
elementary_divisor() click to toggle source
# File lib/algebra/elementary-divisor.rb, line 52
def elementary_divisor
  e_diagonalize.diag
end
horizontal_sweep!(d) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 153
def horizontal_sweep!(d)
  m = transpose
  if m.vertical_sweep!(d)
    replace m.transpose
    true
  else
    false
  end
end
sq_find_nondiv(d, x) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 120
def sq_find_nondiv(d, x)
  (d...rsize).each do |i|
    (d...csize).each do |j|
      return [i, j] unless (self[i, j] % x).zero?
    end
  end
  nil
end
vertical_sweep!(d) click to toggle source
# File lib/algebra/elementary-divisor.rb, line 129
def vertical_sweep!(d)
  sw = false
  if self[d, d].zero? && (i = ((d + 1)...rsize).find { |i1| !self[i1, d].zero? })
    sswap_r!(d, i)
    sw = true
  end
  unless self[d, d].zero?
    (d + 1...rsize).each do |i0|
      next if self[i0, d].zero?
      q = self[i0, d] / self[d, d]
      sw = true
      mix_r!(i0, d, -q)
      # high speed but lost genericity (wo'nt work on triplet)
      #         r = row!(i0)
      #         (d...csize).each do |j0|; r[j0] -=  q * self[d, j0]; end
      unless self[i0, d].zero?
        sswap_r!(d, i0)
        redo
      end
    end
  end
  sw
end