#! perl -slw
use strict;
use List::Util qw[ sum ];

my @matrix =(
[0.54, 0.03, 0.03, 0.03, 0.08, 0.04, 0.05, 0.05, 0.02, 0.02, 0.03, 0.04, 0.03, 0.01, 0.06, 0.11, 0.06, 0.01, 0.01,0.06],
[0.02, 0.57, 0.03, 0.01, 0.02, 0.06, 0.03, 0.01, 0.04, 0.01, 0.01, 0.11, 0.02, 0.01, 0.01, 0.02, 0.02, 0.01, 0.01, 0.01],
[0.01, 0.02, 0.48, 0.06, 0.01, 0.03, 0.03, 0.03, 0.05, 0.01, 0.01, 0.03, 0.01, 0.00, 0.01, 0.04, 0.03, 0.00, 0.02, 0.01],
[0.02, 0.01, 0.09, 0.57, 0.01, 0.04, 0.11, 0.02, 0.03, 0.00, 0.00, 0.03, 0.01, 0.00, 0.02, 0.04, 0.02, 0.00, 0.01, 0.00],
[0.01, 0.00, 0.00, 0.00, 0.57, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.01, 0.00, 0.01, 0.01, 0.01, 0.01, 0.01],
[0.02, 0.04, 0.03, 0.02, 0.01, 0.42, 0.05, 0.01, 0.04, 0.01, 0.01, 0.05, 0.02, 0.00, 0.01, 0.02, 0.02, 0.01, 0.01, 0.01],
[0.04, 0.04, 0.05, 0.13, 0.01, 0.11, 0.52, 0.02, 0.03, 0.01, 0.01, 0.07, 0.01, 0.00, 0.02, 0.04, 0.04, 0.01, 0.01, 0.01],
[0.05, 0.02, 0.05, 0.03, 0.02, 0.02, 0.02, 0.77, 0.02, 0.00, 0.00, 0.02, 0.01, 0.00, 0.01, 0.04, 0.01, 0.01, 0.00, 0.01],
[0.01, 0.01, 0.02, 0.01, 0.01, 0.02, 0.01, 0.00, 0.54, 0.00, 0.00, 0.01, 0.01, 0.02, 0.01, 0.01, 0.01, 0.02, 0.04, 0.00],
[0.02, 0.01, 0.01, 0.00, 0.03, 0.01, 0.01, 0.00, 0.01, 0.50, 0.10, 0.01, 0.09, 0.04, 0.01, 0.01, 0.03, 0.02, 0.01, 0.18],
[0.03, 0.02, 0.01, 0.01, 0.03, 0.03, 0.01, 0.00, 0.02, 0.13, 0.60, 0.02, 0.19, 0.09, 0.02, 0.02, 0.02, 0.02, 0.03, 0.08],
[0.03, 0.14, 0.06, 0.03, 0.01, 0.09, 0.06, 0.02, 0.05, 0.01, 0.01, 0.50, 0.02, 0.01, 0.02, 0.04, 0.03, 0.01, 0.01, 0.01],
[0.01, 0.01, 0.00, 0.00, 0.01, 0.02, 0.00, 0.00, 0.01, 0.03, 0.05, 0.01, 0.42, 0.02, 0.00, 0.01, 0.01, 0.01, 0.01, 0.02],
[0.01, 0.01, 0.00, 0.00, 0.02, 0.01, 0.00, 0.00, 0.03, 0.02, 0.04, 0.00, 0.04, 0.61, 0.00, 0.01, 0.01, 0.07, 0.14, 0.01],
[0.03, 0.01, 0.01, 0.01, 0.00, 0.01, 0.01, 0.01, 0.02, 0.00, 0.01, 0.01, 0.01, 0.00, 0.72, 0.03, 0.01, 0.00, 0.00, 0.01],
[0.07, 0.03, 0.06, 0.04, 0.05, 0.04, 0.03, 0.03, 0.02, 0.01, 0.01, 0.04, 0.02, 0.01, 0.04, 0.43, 0.09, 0.02, 0.01, 0.01],
[0.04, 0.02, 0.04, 0.02, 0.04, 0.03, 0.03, 0.01, 0.02, 0.02, 0.02, 0.03, 0.03, 0.01, 0.02, 0.09, 0.52, 0.01, 0.01, 0.04],
[0.00, 0.00, 0.00, 0.00, 0.01, 0.00, 0.00, 0.00, 0.01, 0.00, 0.00, 0.00, 0.00, 0.02, 0.00, 0.00, 0.00, 0.71, 0.02, 0.00],
[0.00, 0.01, 0.01, 0.00, 0.01, 0.01, 0.00, 0.00, 0.05, 0.01, 0.01, 0.01, 0.01, 0.10, 0.00, 0.01, 0.01, 0.06, 0.61, 0.01],
[0.05, 0.01, 0.01, 0.01, 0.07, 0.02, 0.02, 0.01, 0.01, 0.20, 0.07, 0.01, 0.06, 0.03, 0.01, 0.02, 0.06, 0.02, 0.02, 0.51]
);

my @freq = (
    0.08477395, 0.05103334, 0.03837665, 0.05740129, 0.01256165,
    0.03471746, 0.06883297, 0.07907659, 0.02077239, 0.06813598,
    0.08906677, 0.06501537, 0.02318017, 0.03823936, 0.04036729,
    0.05668318, 0.05721648, 0.00871765, 0.02785317, 0.07797831,
);

my @lnPAM = map[ map log( $_ || 1e-308 ), @$_ ], @matrix;
my $k = 1;
my $mutRate = sum map{ $freq[ $_ ] * ( 1 - $matrix[ $_ ][ $_ ] ) } 0..19;
my $iter = 0;

# The following loop goes until the difference reach the computer's
# precision.
while( abs( $mutRate - 0.01) > 1e-16 ) {
    ++$iter;
    $k *= $mutRate * 100;
    map $_ /= ( $mutRate * 100), @$_ for @lnPAM;
    @matrix = map[ map exp( $_ ), @$_ ], @lnPAM;
    $mutRate = sum map{ $freq[ $_ ] * ( 1 - $matrix[ $_ ][ $_ ] ) } 0..19;
    printf "After %d iteration(s), k=%f, RateMutation=%f\n", $iter, $k, $mutRate;
}

__END__
C:\test>junk
After 1 iteration(s), k=44.513577, RateMutation=0.013418
After 2 iteration(s), k=59.727332, RateMutation=0.010018
After 3 iteration(s), k=59.836899, RateMutation=0.010000
After 4 iteration(s), k=59.837487, RateMutation=0.010000
After 5 iteration(s), k=59.837490, RateMutation=0.010000
After 6 iteration(s), k=59.837490, RateMutation=0.010000
After 7 iteration(s), k=59.837490, RateMutation=0.010000