時間tのときの、六面体のどの点をすでに通ったかという状態の確率を求めるのに、漸化式を導出してそれを解くために固有値を求めて、それに対応する固有ベクトルを求めようとしているのだった。
今回は、そのうち、有理数の固有値に対応する固有ベクトルを求める。
まず、行列を縦ベクトルに分解して、そのうち一次独立のものを選ぶ。これにはGramの行列式を使う。一次従属のベクトルがあるはずなので、その代わりに残りのベクトルとすべて直交するベクトルを選ぶ。
理屈はすぐに思いついたが、オーバーフローしたりしてなかなかうまくいかなかった。最後にPerlで組んで、結果こうなった。
λ=0
0 0 13 -2 -11 0 4 13 -2 -11 0 0 0 0 0 0 0 4 0 13 -2 0 0 0 0
0 0 0 1 -1 0 -2 0 1 -1 0 0 0 0 0 0 0 -2 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 -1 2 -1 2 -1 -1 0 -1 0 0 2 0 0 0λ=1
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 -1 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0λ=2
2 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0λ=3
61 61 87 48 48 139 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48
0 0 2 -1 -1 6 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
数があわないが。
use strict;
use bigint;my @A = (
[1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,1,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,1,0,2,0,0,0,0,1,2,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,2,0,0,0,0,2,0,1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,2,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,2,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,1,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3,2,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3]
);main();
sub main() {
my $tA = mat_transpose(\@A);
for my $e(0..3) {
my $B = subtract_unit_mat($tA, $e);
my $C = extract_independent_vectors($B);
my $D = orthogonal_vectors($C);
mat_print($D);
check_orthogonal($D, $B);
}
}sub mat_transpose {
my $B = [ ];
my $a = $_[0];
for(my $i = 0; $i <= $#{$a}; $i++) {
for(my $j = 0; $j <= $#{$a->[$i]}; $j++) {
$B->[$j] = [ ] if $i == 0;
$B->[$j]->[$i] = $a->[$i]->[$j];
}
}
return $B;
}sub subtract_unit_mat {
my $B = [ ];
my ($A, $e) = @_;
for(my $i = 0; $i <= $#{$A}; $i++) {
$B->[$i] = [ ];
for(my $j = 0; $j <= $#{$A->[$i]}; $j++) {
$B->[$i]->[$j] = $A->[$i]->[$j];
if($i == $j) {
$B->[$i]->[$j] -= $e;
}
}
}
return $B;
}sub extract_independent_vectors {
my $B = [ ];
my ($A) = @_;
push @{$B}, $A->[0];
for(my $i = 1; $i <= $#{$A}; $i++) {
my $G = Gram($B, $A->[$i]);
if(is_independent($G)) {
push @{$B}, $A->[$i];
}
}
return $B;
}sub orthogonal_vectors {
my $B = [ ];
my ($A) = @_;
my $n = $#{$A->[0]} + 1;
my $m = $n - $#{$A} - 1;
for(my $i = 0; $i < $n; $i++) {
my $e = [ map { $_ == $i ? 1 : 0 } 0..$n-1 ];
my $v = extract_orthogonal_vector($A, $e);
if(!is_zero_vector($v)) {
push @{$B}, $v;
push @{$A}, $v;
if($#{$A} == $n - 1) {
return $B;
}
}
}
return $B;
}sub check_orthogonal {
my ($D, $B) = @_;
for(my $i = 0; $i <= $#{$D}; $i++) {
my $ok = 1;
for(my $j = 0; $j < $#{$B}; $j++) {
if(inner_product($D->[$i], $B->[$j]) != 0) {
$ok = 0;
last;
}
}
print $ok ? "ok!" : "NG", "\n";
}
}sub extract_orthogonal_vector {
my ($A, $e) = @_;
my $B = [ ];
my $m = $#{$A} + 1;
my $n = $#{$e} + 1;
for(my $i = 0; $i < $m; $i++) {
$B->[$i] = [ ];
for(my $j = 0; $j < $m; $j++) {
$B->[$i]->[$j] = inner_product($A->[$i], $A->[$j]);
}
$B->[$i]->[$m] = inner_product($A->[$i], $e);
}
my $t = solve_linear_equation($B);
my $v = [ map { $_ * $t->[0]->[1] } @{$e} ];
for(my $i = 0; $i < $m; $i++) {
for(my $j = 0; $j < $n; $j++) {
$v->[$j] -= $t->[$i]->[0] * $A->[$i]->[$j];
}
}
vec_normalize($v);
return $v;
}sub solve_linear_equation {
my ($A) = @_;
my $n = $#{$A} + 1;
my $B = mat_copy($A);
for(my $i = 0; $i < $n - 1; $i++) {
if($B->[$i]->[$i] == 0) {
my $j = 0;
for($j = $i + 1; $j <= $n; $j++) {
last if $B->[$j]->[$i] != 0;
}
for(my $k = $i; $k <= $n; $k++) {
$B->[$i]->[$k] += $B->[$j]->[$k];
}
}
for(my $j = $i + 1; $j < $n; $j++) {
if($B->[$j]->[$i] != 0) {
my ($f1, $f2);
if($B->[$j]->[$i] % $B->[$i]->[$i] == 0) {
$f1 = 1;
$f2 = $B->[$j]->[$i] / $B->[$i]->[$i];
}
else {
my $d = GCD($B->[$j]->[$i], $B->[$i]->[$i]);
$f1 = $B->[$i]->[$i] / $d;
$f2 = $B->[$j]->[$i] / $d;
}
for(my $k = $i + 1; $k <= $n; $k++) {
$B->[$j]->[$k] = $f1 * $B->[$j]->[$k]
- $f2 * $B->[$i]->[$k];
}
}
}
}
for(my $i = $n - 1; $i > 0; $i--) {
for(my $j = 0; $j < $i; $j++) {
if($B->[$j]->[$i] != 0) {
my ($f1, $f2);
if($B->[$j]->[$i] % $B->[$i]->[$i] == 0) {
$f1 = 1;
$f2 = $B->[$j]->[$i] / $B->[$i]->[$i];
}
else {
my $d = GCD($B->[$j]->[$i], $B->[$i]->[$i]);
$f1 = $B->[$i]->[$i] / $d;
$f2 = $B->[$j]->[$i] / $d;
}
$B->[$j]->[$n] = $f1 * $B->[$j]->[$n]
- $f2 * $B->[$i]->[$n];
for(my $k = $j; $k < $i; $k++) {
$B->[$j]->[$k] *= $f1;
}
}
}
}
for(my $i = 0; $i < $n; $i++) {
my $d = GCD($B->[$i]->[$i], $B->[$i]->[$n]);
if($d > 1) {
$B->[$i]->[$i] /= $d;
$B->[$i]->[$n] /= $d;
}
}
my $m = $B->[0]->[0];
for(my $i = 1; $i < $n; $i++) {
my $d = GCD($m, $B->[$i]->[$i]);
$m *= $B->[$i]->[$i] / $d;
}
my $t = [ ];
for(my $i = 0; $i < $n; $i++) {
$t->[$i] = [ $B->[$i]->[$n] * ($m / $B->[$i]->[$i]), $m ];
}
return $t;
}sub Gram {
my ($A, $v) = @_;
my $n = $#{$A} + 2;
my $G = [ ];
for(my $i = 0; $i < $n - 1; $i++) {
$G->[$i] = [ ];
for(my $j = 0; $j < $n - 1; $j++) {
$G->[$i]->[$j] = inner_product($A->[$i], $A->[$j]);
}
$G->[$i]->[$n-1] = inner_product($A->[$i], $v);
}
$G->[$n-1] = [ ];
for(my $j = 0; $j < $n - 1; $j++) {
$G->[$n-1]->[$j] = inner_product($v, $A->[$j]);
}
$G->[$n-1]->[$n-1] = inner_product($v, $v);
return $G;
}sub is_independent {
my ($A) = @_;
my $n = $#{$A} + 1;
my $B = mat_copy($A);
for(my $i = 0; $i < $n; $i++) {
if($B->[$i]->[$i] == 0) {
my $j = 0;
for($j = $i + 1; $j < $n; $j++) {
last if $B->[$j]->[$i] != 0;
}
return 0 if $j == $n;
for(my $k = $i; $k < $n; $k++) {
$B->[$i]->[$k] += $B->[$j]->[$k];
}
}
for(my $j = $i + 1; $j < $n; $j++) {
if($B->[$j]->[$i] != 0) {
my ($f1, $f2);
if($B->[$j]->[$i] % $B->[$i]->[$i] == 0) {
$f1 = 1;
$f2 = $B->[$j]->[$i] / $B->[$i]->[$i];
}
else {
my $d = GCD($B->[$j]->[$i], $B->[$i]->[$i]);
$f1 = $B->[$i]->[$i] / $d;
$f2 = $B->[$j]->[$i] / $d;
}
for(my $k = $i + 1; $k < $n; $k++) {
$B->[$j]->[$k] = $f1 * $B->[$j]->[$k]
- $f2 * $B->[$i]->[$k];
}
}
}
}
return 1;
}sub GCD {
my ($a, $b) = @_;
if($a == 0) {
return $b;
}
if($b == 0) {
return $a;
}
if($a < 0) {
$a = -$a;
}
if($b < 0) {
$b = -$b;
}
my ($a1, $b1);
while(1) {
if(($a1 = $a % $b) == 0) {
return $b;
}
if(($b1 = $b % $a1) == 0) {
return $a1;
}
$a = $a1;
$b = $b1;
}
}
sub is_zero_vector {
my $v = $_[0];
for my $e(@{$v}) {
return 0 if $e != 0;
}
return 1;
}sub inner_product {
my $sum = 0;
my ($a, $b) = @_;
for(my $i = 0; $i <= $#{$a}; $i++) {
$sum += $a->[$i] * $b->[$i];
}
return $sum;
}sub vec_normalize {
my $v = $_[0];
my $d = 0;
for my $e(@{$v}) {
$d = GCD($e, $d);
}
if($d > 1) {
for(my $i = 0; $i <= $#{$v}; $i++) {
$v->[$i] /= $d;
}
}
}sub mat_copy {
my $A = $_[0];
my $B = [ ];
for(my $i = 0; $i <= $#{$A}; $i++) {
$B->[$i] = [ ];
@{$B->[$i]} = @{$A->[$i]};
}
return $B;
}sub rat_print {
my $r = $_[0];
if($r->[1] == 1) {
print $r->[0];
}
else {
print $r->[0], "/", $r->[1];
}
}sub vec_print {
my $v = $_[0];
print join(' ', @{$v}), "\n";
}sub mat_print {
my $a = $_[0];
for(my $i = 0; $i <= $#{$a}; $i++) {
vec_print $a->[$i];
}
}
