PAM1 := MutationMatrix:
lnPAM1 := ln(MutationMatrix):
k := 1:
RateMutation := sum(Freq[i] * (1 - PAM1[i,i]), i= 1..20):
iter := 0:

# The following loop goes until the difference reach the computer's
# precision.
while abs(RateMutation - 0.01) > DBL_EPSILON do
    iter := iter + 1;
    k:= k * RateMutation * 100:
    lnPAM1 := lnPAM1 / (RateMutation * 100):
    PAM1 := exp(lnPAM1):
    RateMutation := sum(Freq[i] * (1 - PAM1[i,i]), i= 1..20);
    printf('After %d iteration(s), k=%f, RateMutation=%f
', iter, k, RateMutation);
od:

##</code><code>##

After 1 iteration(s), k=44.530219, RateMutation=0.014324
After 2 iteration(s), k=63.785069, RateMutation=0.010027
After 3 iteration(s), k=63.959878, RateMutation=0.010000
After 4 iteration(s), k=63.960985, RateMutation=0.010000
After 5 iteration(s), k=63.960992, RateMutation=0.010000
After 6 iteration(s), k=63.960992, RateMutation=0.010000
After 7 iteration(s), k=63.960992, RateMutation=0.010000
After 8 iteration(s), k=63.960992, RateMutation=0.010000


##</code><code>##

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 @frequency = ([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 $mutrate=0;
my $matrix2={};
my $e;
$e=1;$e/=2 while $e/2+1>1;

my @PAM1= @matrix;

my @lnPAM1=();
LINE: for (my $i=0; $i<20; $i++){
    for (my $j=0; $j<20; $j++){
        if ($matrix[$i][$j] > 0){
            $lnPAM1[$i][$j] = log $matrix[$i][$j];
        }else{ next LINE;
        }
    }
}


for (my $i=0; $i<20; $i++){
    $mutrate += $frequency[0][$i]*(1-$matrix[$i][$i]);
}

my $k=1;
my $iter = 0;

while (($mutrate-0.01) > $e){
    $iter++;
    $k= $k * $mutrate * 100;
    for (my $i=0; $i<20; $i++){
        for (my $j=0; $j<20; $j++){
            $lnPAM1[$i][$j] = $lnPAM1[$i][$j]/($mutrate * 100);
            $PAM1[$i][$j]=exp($lnPAM1[$i][$j]);
            $mutrate += $frequency[0][$i]*(1-$PAM1[$i][$i]);
            
        }
    }
    
    print " iteration: $iter, k: $k, MutationRate:$mutrate\n";
    
}