require 'test/unit'
require 'algebra'
require 'algebra/jordan-form'

class TestJordanForm < Test::Unit::TestCase
  # "Jordan Hyouzyunkei I & Isa, Tokudai Syuppan Kai # pp.96-97
  def setup
    @A = [
      [ # 0
        [4, 1, 1, 1],
        [1, 5, 3, -1],
        [-1, -2, 1, -1],
        [-1, -1, -1, 2]
      ],
      [ # 1
        [-3, -3, 4, 2],
        [5, 1, -4, -6],
        [4, -1, -2, -5],
        [-5, -3, 4, 4]
      ],
      [ # 2
        [3, -1, 3],
        [-1, 3, 3],
        [1, -1, 5]
      ],
      [ # 3
        [-1, 1, 2, -1],
        [-5, 3, 4, -2],
        [3, -1, 0, 1],
        [5, -2, -2, 3]
      ],
      [ # 4
        [1, 2, 0, 0],
        [0, 1, 2, 0],
        [0, 0, 1, 2],
        [0, 0, 0, 1]
      ],
      [ # 5
        [3, 0, 2, 1],
        [0, 3, 0, 2],
        [0, 0, 3, 0],
        [0, 0, 0, 3]
      ],
      [ # 6
        [1, 0, 0, 0],
        [2, 1, 0, 0],
        [0, 2, 1, 0],
        [0, 0, 2, 1]
      ],
      [ # 7
        [2, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 2, 0],
        [5, 0, 0, 2]
      ],
      [ # 8
        [0, 0, 8],
        [1, 0, 4],
        [0, 1, -2]
      ],

      #### Original
      [ # 9
        [1, 1, 0],
        [1, 1, 1],
        [0, 1, 1]
      ],
      [ # 10
        [1, 0, 8],
        [1, 1, 4],
        [0, 1, -2]
      ],
      [ # 11 288sec.
        [0, 0, 1, 0],
        [1, 1, 1, 1],
        [1, 1, 1, 1],
        [0, 1, 1, 1]
      ]
    ]
    @Ai = (0...@A.size).to_a
    @PR = Algebra.Polynomial(Rational, 'x')
  end

  def test_jordan_form0
    # puts
    @Ai.each do |i|
      m = @A[i]
      #    Ai[0..9].each do |i|; m = A[i]
      # puts "##{i}"
      _test_jordan_form(m)
    end
  end

  def _test_jordan_form(a)
    m = Algebra.SquareMatrix(Rational, a.size)[*a]
    # m.display; puts
    asr = Algebra.SquareMatrix(@PR, a.size)

    me = m._char_matrix(asr).to_quint.e_diagonalize
    elem_divs = me.body.diag

    # p elem_divs; puts
    # me.body.display;puts

    facts = Algebra::ElementaryDivisor.factorize(elem_divs)
    dl = facts.last
    field, modulus, facts, roots, elms = dl.pi.decompose(dl)
    # p [facts, roots, elms]

    pfield = Algebra.Polynomial(field, 'x')
    pfm_alg = Algebra.SquareMatrix(pfield, a.size)

    jm = Algebra::JordanForm.construct(elem_divs, facts, field, pfield)
    jme = jm._char_matrix(pfm_alg).to_quint.e_diagonalize

    mebody = pfm_alg.convert_from(me.body)
    meright = pfm_alg.convert_from(me.right)
    meleft = pfm_alg.convert_from(me.left)

    assert_equal(mebody, jme.body)

    sR = (meright * jme.righti).i2o.evaluateR(jm)
    tL = (jme.lefti * meleft).i2o.evaluateL(jm)

    m = m.convert_to(jm.class)

    k = tL * m * sR

    # k.display;puts
    assert_equal(k, jm)
  end

  def test_jordan_form
    puts
    @Ai[0..0].each do |i|
      a = @A[i]
      # puts "##{i}"
      m = Algebra.SquareMatrix(Rational, a.size)[*a]
      # m.display; puts

      jf, p, q, field, modulus = Algebra::JordanForm.decompose(m)
      # print "extension = "
      # p modulus
      # jf.display; puts
      #
      # print "P =\n"; p.display; puts
      # print "P^-1 =\n"; q.display; puts

      m = m.convert_to(jf.class)
      assert_equal(jf, p * m * q)
    end
  end
end
