[ruby-dev:27731] Re: Matrix class is broken without mathn

原です。

In [ruby-dev:27700]
>けいじゅ＠いしつかです.

>あと, Matrix#detだけでなく, Matirx#rankも同じ問題があるのですが, 
>いくらなんでもこちらをFloatで計算しろとは言えないのではないかと?

その主張は同意します。だからrationalをrequireして、Rational
で(quoで)計算しろ、ということなんです。

Float行列にしろInteger行列にしろrankを計算したい人がそもそも
非常に少ないのでまっいっかというのもあります(^^;。元々、
Float行列のrankの値なんか信用できないので、Integer行列のrank
を知りたい人は、rational.rbをrequireする手間は惜しくないので
はないかと。逆に、quo式はrational.rbをrequireするだけで済む
だけ、更に各成分Rational()しなくてはならない現行よりアドバン
テージがある。

さて、以下にInteger行列のrankを求める互除法版を載せます。
(やはりRational との混交可です。)

ここで提案ですが、rankもquo式にし、このrank_e（名前は改めて
考えるとして）も取り込んで、マニュアルに

  【整数行列を扱う場合の注意】
    整数行列の階数を計算させるときは、Matrix#rank_eを使いましょ
    う。rational.rbをrequireしている状態でMatrix#rankを使って
    も正しい値が得られます。

としたらどうでしょうか。

#!/usr/local/bin/ruby
require "./matrix"
require "rational"

class Matrix
  def rank_e
    a = to_a
    a_column_size = column_size
    a_row_size = row_size
    pi = 0
    (0 ... a_column_size).each do |j|
      if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
        if i != pi
          a[pi], a[i] = a[i], a[pi]
        end
        (pi + 1 ... a_row_size).each do |k|
          q = a[k][j] / a[pi][j]
          (pi ... a_column_size).each do |j0|
            a[k][j0] -= q * a[pi][j0]
          end
          if k > pi && !a[k][j].zero?
            a[k], a[pi] = a[pi], a[k]
            redo
          end
        end
        pi += 1
      end
    end
    pi
  end
end

class Test
  attr_reader :value
  def initialize(name)
    @name = name
    @value = nil
    @time = 0.0
  end

  def stime(v = nil)
    t = Time.now
    @value = yield
    @time += Time.now - t
    if v && v != @value
      puts "Test Failed: #{@value} != #{v}"
      exit
    end
  end

  def to_s
    sprintf("%-30s: %3.3fsec.", @name, @time)
  end
end

def typical(row, col, rank, range = 100)
  s = 0
 (0...row).map{|i|
    (0...col).map{|j| (i == j && s < rank) ? (s += 1; 1) : 0}
  }
end

def test(rsize, csize, r, num, range = 100)
  tr = Test.new("rank(rational)")
  te = Test.new("rank_e(int)")
  t1 = Test.new("rank_e(int & rational)")
  tf = Test.new("rank(float)")

  num.times do
    left = (0...rsize).map{(0...rsize).map{rand(2*range+1)-range}}
    middle = typical(rsize, csize, r)
    right = (0...csize).map{(0...csize).map{rand(2*range+1)-range}}

    m = Matrix[*left]*Matrix[*middle]*Matrix[*right]
    # rank m == r (almost)
    a = m.to_a
    m1 = Matrix[*a.map{|b| b.map{|c|
          rand(2).zero? ? Rational(c, 1) : c}}]
    mr = Matrix[*a.map{|b| b.map{|c| Rational(c, 1)}}]
    mf = Matrix[*a.map{|b| b.map{|c| Float(c)}}]

    tr.stime{Matrix[*mr].rank}
    te.stime(tr.value){Matrix[*m].rank_e}
    t1.stime(tr.value){Matrix[*m1].rank_e}
    tf.stime(){Matrix[*mf].rank}
  end

  puts "--- #{rsize}x#{csize}, rank(#{r}), #{num}times ---"
  puts tr, te, t1, tf
  puts
end

test(4, 3, 3, 100)
test(4, 5, 3, 100)
test(5, 5, 4, 100)
test(10, 10, 7, 10)
test(9, 11, 7, 10)
test(11, 9, 7, 10)


結果：
--- 4x3, rank(3), 100times ---
rank(rational)                : 0.470sec.
rank(rational)                : 0.044sec.(with raional.so)
rank_e(int)                   : 0.188sec.
rank_e(int & rational)        : 0.597sec.
rank(float)                   : 0.029sec.

--- 4x5, rank(3), 100times ---
rank(rational)                : 1.003sec.
rank(rational)                : 0.210sec.(with raional.so)
rank_e(int)                   : 0.227sec.
rank_e(int & rational)        : 1.187sec.
rank(float)                   : 0.073sec.

--- 5x5, rank(4), 100times ---
rank(rational)                : 1.995sec.
rank(rational)                : 0.386sec.(with raional.so)
rank_e(int)                   : 0.422sec.
rank_e(int & rational)        : 2.349sec.
rank(float)                   : 0.032sec.

--- 10x10, rank(7), 10times ---
rank(rational)                : 3.016sec.
rank(rational)                : 1.003sec.(with raional.so)
rank_e(int)                   : 0.284sec.
rank_e(int & rational)        : 3.671sec.
rank(float)                   : 0.018sec.

--- 9x11, rank(7), 10times ---
rank(rational)                : 3.030sec.
rank(rational)                : 1.067sec.(with raional.so)
rank_e(int)                   : 0.283sec.
rank_e(int & rational)        : 3.365sec.
rank(float)                   : 0.022sec.

--- 11x9, rank(7), 10times ---
rank(rational)                : 3.008sec.
rank(rational)                : 1.010sec.(with raional.so)
rank_e(int)                   : 0.218sec.
rank_e(int & rational)        : 3.639sec.
rank(float)                   : 0.018sec.



あ、そういえば matrix.rbのtypo見つけました。

--- matrix.rb~	2005-11-17 01:10:30.000000000 +0900
+++ matrix.rb	2005-11-17 10:57:56.000000000 +0900
@@ -740,7 +740,7 @@ class Matrix
           i = k
           exists = true
           begin
-            if (i += 1) > a_row_size - 1
+            if (i += 1) > a_column_size - 1
               exists = false
               break
             end