longd2_function(ps, pss, pff, n,d)
{
#this program compute the probability P(L2>=d)#
#L2 is the 2-interrupted run within a sequence of length n# 
# ps is the initial prob. for "s", and pss is the transition prob of 1 given 1#
 NN_(2*d^3+10*d)/12+1
pf_1-ps
 psf_1-pss
 pfs_1-pff
c0_matrix(0,1,NN-1)
p0_cbind(1,c0)
c1_matrix(1,NN-1,1)
U_rbind(c1,0)

 M_matrix(0,NN,NN)
 for(ir in 0:(d-1)){
   for(jr in 0:ir){
     for(kr in 0:ir){
for(ic in 0:(d-1)){
          for(jc in 0:ic){
            for(kc in 0:ic){
           rid_((2*ir^3+10*ir)/12+(jr-1)*(2*ir-jr+2)/2+kr-jr+2)*(kr>jr)*(jr>0)+
              ((2*ir^3+10*ir)/12+((jr-1)*(2*ir-jr+2)/2+2)*(jr>0))*(kr==0)+
                1*(kr==0)*(jr==0)
           cid_((2*ic^3+10*ic)/12+(jc-1)*(2*ic-jc+2)/2+kc-jc+2)*(kc>jc)*(jc>0)+
((2*ic^3+10*ic)/12+((jc-1)*(2*ic-jc+2)/2+2)*(jc>0))*(kc==0)+
                1*(kc==0)*(jc==0)
if(ir>0){
if(jr>1 && kr>jr &&ic==ir+1 && jc==jr+1 && kc==kr+1){M[rid,cid]_pss}
if(jr>1 && kr==0 && ic==ir+1 && jc==jr+1 && kc==0){M[rid,cid]_pss}
if(jr==0 && kr==0 && ic==ir+1 && kc==0 && jc==0){M[rid,cid]_pss}
if(jr==0 && kr==0 && ic==ir+1 && jc==1 && kc==0){M[rid,cid]_psf}
if(jr>1&& kr==0&& ic==ir+1 && jc==1 && kc==jr+1){M[rid,cid]_psf}
if(jr>1 && kr>jr && ic==kr && jc==1 && kc==jr+1){M[rid,cid]_psf}
if(jr==1 && kr==0 && ic==ir+1 && jc==1 && kc==2){M[rid,cid]_pff}
if(jr==1 && kr>jr+1 && ic==kr && jc==1 && kc==2){M[rid,cid]_pff}
if(jr==1 && kr==2 && ic==2 && jc==1 && kc==2){M[rid,cid]_pff}
if(jr==1 && kr==0 && ic==ir+1 && jc==2 && kc==0){M[rid,cid]_pfs}
if(jr==1 && kr>jr && ic==ir+1 && jc==jr+1 && kc==kr+1){M[rid,cid]_pfs}}
if(ir==0 && jr==0 && kr==0 && ic==1 && jc==0 && kc==0){M[rid,cid]_ps}
if(ir==0 && jr==0 && kr==0 && ic==1 && jc==1 && kc==0){M[rid,cid]_pf}
}}}}}}
for(i in 1:NN) { M[i,NN]_1-sum(M[i,]) }
Mn_M%*%M
for (i in 3:n) {Mn_Mn%*%M }
1-p0%*%Mn%*%U
}