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:

PrintMatrix(transpose(transpose(PAM1[1..10])[1..10]),'% .5f ');