# Matrix Algebra Triplet
#
#   by Shin-ichiro Hara
#
# Version 1.0 (2001.11.01)

require 'algebra/gaussian-elimination'
require 'algebra/elementary-divisor'

module Algebra
  class MatrixAlgebra
    def to_triplet
      Algebra::MatrixAlgebraTriplet.new(self)
    end

    def to_quint
      Algebra::MatrixAlgebraQuint.new(self)
    end
  end

  class MatrixAlgebraTriplet
    include GaussianElimination
    include ElementaryDivisor

    attr_reader :left, :body, :right, :ground
    def initialize(matrix, left = nil, right = nil)
      @body = matrix
      @type = @body.class
      @ground = @type.ground
      @left_type = Algebra.SquareMatrix(@type.ground, @type.rsize)
      @right_type = Algebra.SquareMatrix(@type.ground, @type.csize)
      @left = left ? left : @left_type.unity
      @right = right ? right : @right_type.unity
    end

    def to_ary
      [body, left, right]
    end

    alias to_a to_ary

    def dup
      self.class.new(body.dup, left.dup, right.dup)
    end

    def transpose
      self.class.new(body.transpose, right.transpose, left.transpose)
    end

    def replace(other)
      initialize(other.body, other.left, other.right)
      self
    end

    def display
      puts '============= begin'
      @left.display
      puts
      @body.display
      puts
      @right.display
      puts '============= end'
    end

    def [](i, j)
      @body[i, j]
    end

    def rsize
      @body.rsize
    end

    def csize
      @body.csize
    end

    def each_i(&block)
      @body.each_i(&block)
    end

    def each_j(&block)
      @body.each_j(&block)
    end

    def row!(i)
      @body.row!(i)
    end

    # ElementaryOpeartion
    def sswap_r!(i, j)
      @left.sswap_r!(i, j)
      @body.sswap_r!(i, j)
      self
    end

    def swap_r!(i, j)
      @left.swap_r!(i, j)
      @body.swap_r!(i, j)
      self
    end

    #  def swap_r(i, j)
    #    dup.swap_r!(i, j)
    #  end

    def swap_c!(i, j)
      @body.swap_c!(i, j)
      @right.swap_c!(i, j)
      self
    end

    #  def swap_c(i, j)
    #    dup.swap_c!(i, j)
    #  end

    def multiply_r!(i, c)
      @left.multiply_r!(i, c)
      @body.multiply_r!(i, c)
      self
    end

    #  def multiply_r(i, c)
    #    dup.multiply_r!(i, c)
    #  end

    def multiply_c!(j, c)
      @body.multiply_c!(j, c)
      @right.multiply_c!(j, c)
      self
    end

    #  def multiply_c(j, c)
    #    dup.multiply_c!(j, c)
    #  end

    def divide_r!(i, c)
      @left.divide_r!(i, c)
      @body.divide_r!(i, c)
    end

    #  def divide_r(i, c)
    #    dup.divide_r!(i, c)
    #  end

    def divide_c!(j, c)
      @body.divide_c!(j, c)
      @right.divide_c!(j, c)
    end

    #  def divide_c(j, c)
    #    dup.divide_c!(j, c)
    #  end

    def mix_r!(i, j, c = nil)
      @left.mix_r!(i, j, c)
      @body.mix_r!(i, j, c)
      self
    end

    #  def mix_r(i, j, c = nil)
    #    dup.mix_r!(i, j, c)
    #  end

    def mix_c!(i, j, c = nil)
      @body.mix_c!(i, j, c)
      @right.mix_c!(i, j, c)
      self
    end

    #  def mix_c(i, j, c = nil)
    #    dup.mix_c!(i, j, c)
    #  end

    def left_eliminate!
      #    inv = Algebra.SquareMatrix(ground, rsize).unity
      k = ground.unity
      pi = 0
      each_j do |j|
        next unless i = (pi...rsize).find { |i1| !self[i1, j].zero? }
        if i != pi
          swap_r!(pi, i) # ; inv.swap_r!(pi, i)
          k = -k
        end
        c = ground.unity / self[pi, j] # this lets the entries be in ground
        multiply_r!(pi, c) # ; inv.multiply_r!(pi, c)
        k *= c
        each_i do |i0|
          next if i0 == pi
          d = self[i0, j] # / self[pi, j]
          mix_r!(i0, pi, -d) # ; inv.mix_r!(i0, pi, -d)
        end
        pi += 1
      end
      [left, k]
    end
  end

  class MatrixAlgebraQuint < MatrixAlgebraTriplet
    attr_reader :lefti, :righti
    def initialize(matrix, left = nil, right = nil, lefti = nil, righti = nil)
      super(matrix, left, right)
      @lefti = lefti ? lefti : @left_type.unity
      @righti = righti ? righti : @right_type.unity
    end

    def to_ary
      [body, left, right, lefti, righti]
    end

    alias to_a to_ary

    def dup
      self.class.new(body.dup, left.dup, right.dup, lefti.dup, righti.dup)
    end

    def transpose
      self.class.new(body.transpose, right.transpose, left.transpose,
                     righti.transpose, lefti.transpose)
    end

    def replace(other)
      initialize(other.body, other.left, other.right, other.lefti, other.righti)
      self
    end

    # ElementaryOpeartion
    def sswap_r!(i, j)
      @lefti.swap_c!(i, j) # sswap_c doesn't exist
      super
    end

    def swap_r!(i, j)
      @lefti.swap_c!(i, j)
      super
    end

    def swap_c!(i, j)
      @righti.swap_r!(i, j)
      super
    end

    def multiply_r!(i, c)
      @lefti.divide_c!(i, c)
      super
    end

    def multiply_c!(j, c)
      @righti.divide_r!(j, c)
      super
    end

    def divide_r!(i, c)
      @lefti.multiply_c!(i, c)
      super
    end

    def divide_c!(j, c)
      @righti.multiply_r!(j, c)
      super
    end

    def mix_r!(i, j, c = nil)
      @lefti.mix_c!(j, i, c ? -c : -1)
      super
    end

    def mix_c!(i, j, c = nil)
      @righti.mix_r!(j, i, c ? -c : -1)
      super
    end
  end
end

if $PROGRAM_NAME == __FILE__
  require 'algebra/m-polynomial'
  require 'algebra/polynomial'
  require 'algebra/rational'
  require 'algebra/matrix-algebra'
  #  class Rational# < Numeric
  #    def inspect; to_s; end
  #  end
  #  include Algebra

  G = Algebra.MatrixAlgebra(Rational, 2, 2)
  a = G[[1, 1], [1, -1]]
  a0 = Algebra::MatrixAlgebraTriplet.new(a)
  a0.display
  a0.left_eliminate!
  puts '-------'
  a0.display
  #  p a.kernel_basis
end
