karlin_function(pss, pff,N,d,k) 
{
   # Output: asymptotic value of P(Lk>d) using karlin's formula# 
   # Input: N: total length of sequence,#
   #pss=P(S|S), pff=P(f|f), transition matrix of Markov Chain#
# begin of the subroutine #
karlin3_function(pss, pff,N,d,n1,n2,n3) 
 {
   # Output: asymptotic value of P(Lk>d) using karlin's formula# 
   #n1 : number singleton mismatch within the subseq.#
   #n2 : number doublet mismatch within the subseq.#
   #n3 : number triplet mismatch within the subseq.#
       psf_1-pff
       pfs_1-pss
       temp_c(pss, psf ,pfs ,pff)
       A_matrix(temp,2,2)
       n_2
       cc_matrix(0,n,n)
       B_matrix(c(1,1,0,0),nrow=2)
       {
         for(i in 1:n)   
         for(j in 1:n) 
        {
          cc[i,j]_A[i,j]*B[i,j] 
        }
         lamda_max(eigen(cc)$values)
         phi0_eigen(cc)$vectors[,1]
         kesi0_t(eigen(t(cc))$vectors[,1])
         b_kesi0%*%phi0
         phi_phi0/b
         paiA_matrix(c(pfs/(psf+pfs),psf/(psf+pfs)),nrow=1) 
         paiB_matrix(c(1,0),nrow=1) 
         a_matrix(0,1,n)
          {
            for(j in 1:n)
              {
                a[1,j]_paiA[1,j]*paiB[1,j]*phi[j]
              }
          }
         r_(sum(a))*(sum(kesi0))/(lamda)
         cat ("lamda=",lamda, "\n")
         cat ("r=",r, "\n")
       }
   n0_n1+n2+n3
   E_(log(N)+(n0)*(log(log(N))-log(-log(lamda))))/(-log(lamda))
   s_d-E
   t_n1+2*n2+3*n3
   1-exp(((lamda)-1)*((lamda)^(s))*(((1-(lamda))/(lamda))^t)*(r)/(factorial(n1)*factorial(n2)*factorial(n3)))
 }
# end of the subroutine #
 if (k==0) karlin3(pss, pff,N,d,0,0,0) else
if (k==1) karlin3(pss, pff,N,d,1,0,0) else
if (k==2) karlin3(pss,pff,N,d,2,0,0)+ karlin3(pss,pff,N,d,0,1,0)  else
if (k==3) karlin3(pss,pff,N,d,3,0,0)+karlin3(pss,pff,N,d,1,1,0)+ karlin3(pss,pff,N,d,0,0,1) else
if (k>3) cat("This program allows only at most 3 mismatches","\n")
}