require 'matrix' require 'mathn' A=Matrix[ [1,2,3], [3,1,2], [2,3,1] ] n=A.row_size abort '!error: not full ranked!' if A.det.zero? #see (2.6) and Remark 2.1 of [Pan-Reif:1985] norm=(A*A.transpose).map(&:abs).row_vectors.map{|r| r.to_a.inject(:+)}.max x=Matrix.scalar(n,1.0/norm)*A.transpose p x=0xFF.times.inject(x){|x,_| (Matrix.scalar(n,2.0)-x*A)*x} p Matrix.unit(n)-A*x