subroutine exponential(x,n,nit1,nit2,c,M,k1,k2,tol,pi1,pi2,beta1,beta2,pi1pen,pi2pen,beta1pen,beta2pen,beta0,lrt,mlrt)
	implicit integer*4(i-k,n)
	implicit real*8(a-h,l-m,o-z)
	dimension x(1:n) 
	dimension tau1(1:n),tau2(1:n),f1(1:n),f2(1:n),f0(1:n)
	i = 0
	tmp1 = 0
	tmp2 = 0
	tmp3 = 0 
  102   continue
	i = i + 1
	oldpi1 = pi1
	oldpi2 = pi2
	oldbeta1 = beta1
	oldbeta2 = beta2
	do 103  k = 1, n 
	         tmp1 = x(k)*beta1
		 tmp2 = x(k)*beta2
		 tmp3 = pi1*beta1*exp(-tmp1) + pi2*beta2*exp(-tmp2)
		 tau1(k) = pi1*beta1*exp(-tmp1) / tmp3
		 tau2(k) = 1 - tau1(k)
  103   continue
	tmp4 = 0
	tmp5 = 0
	tmp6 = 0
	do 104  j = 1, n  
		 tmp4 = tmp4 + tau1(j)
		 tmp5 = tmp5 + tau1(j)*x(j)
		 tmp6 = tmp6 + tau2(j)*x(j)
  104   continue	  
	pi1 = tmp4 / n
	pi2 = 1 - pi1
	beta1 = tmp4 / tmp5
	beta2 = (n - tmp4) / tmp6
	if (beta1.ge.M) beta1 = M
	if (beta1.le.1/M) beta1 = 1/M
	if (beta2.ge.M) beta2 = M
	if (beta2.le.1/M) beta2 = 1/M
	dtmp1 = sqrt((beta1-oldbeta1)*(beta1-oldbeta1))/sqrt(beta1*beta1)
	dtmp2 = sqrt((beta2-oldbeta2)*(beta2-oldbeta2))/sqrt(beta2*beta2)
	dtmp3 = sqrt((pi1-oldpi1)*(pi1-oldpi1))/sqrt(pi1*pi1)
	dtmp4 = sqrt((pi2-oldpi2)*(pi2-oldpi2))/sqrt(pi2*pi2)   
	nit1 = i
	if (i.ge.k1) dtmp4 = -dtmp1 + -dtmp2 + -dtmp3
	if ((dtmp1+dtmp2+dtmp3+dtmp4).ge.tol) goto 102
	lrt = 0
        do 105  j = 1, n
        	tmp1 = x(j)*beta1
        	tmp2 = x(j)*beta2
        	tmp3 = x(j)*beta0
        	f1(j) = pi1*exp(-tmp1 )*beta1
        	f2(j) = pi2*exp(-tmp2 )*beta2
        	f0(j) = beta0*exp(-tmp3)
        	lrt = lrt + log(f1(j)+f2(j))  
        	lrt = lrt - log(f0(j))
  105   continue
  	lrt = 2*lrt
      i = 0
	pi1pen = pi1
	pi2pen = pi2
	beta1pen = beta1
	beta2pen = beta2
  106   continue
	i = i + 1
	oldpi1 = pi1pen
	oldpi2 = pi2pen
	oldbeta1 = beta1pen
	oldbeta2 = beta2pen
	do 107  k = 1, n 
	         tmp1 = x(k)*beta1pen
		 tmp2 = x(k)*beta2pen
		 tmp3 = pi1pen*beta1pen*exp(-tmp1) + pi2pen*beta2pen*exp(-tmp2)
		 tau1(k) = pi1pen*beta1pen*exp(-tmp1) / tmp3
		 tau2(k) = 1 - tau1(k)
  107   continue
	tmp4 = 0
	tmp5 = 0
	tmp6 = 0
	do 108  j = 1, n  
		 tmp4 = tmp4 + tau2(j)
		 tmp5 = tmp5 + tau1(j)*x(j)
		 tmp6 = tmp6 + tau2(j)*x(j)
  108   continue	 
	pi1pen = 1 - ((tmp4+c) / (n+2*c))
	pi2pen = 1 - pi1pen
	beta1pen = (n - tmp4) / tmp5
	beta2pen = tmp4 / tmp6
	if (beta1pen.ge.M) beta1pen = M
	if (beta1pen.le.1/M) beta1pen = 1/M
	if (beta2pen.ge.M) beta2pen = M
	if (beta2pen.le.1/M) beta2pen = 1/M
	dtmp1 = sqrt((beta1pen-oldbeta1)*(beta1pen-oldbeta1))/sqrt(beta1pen*beta1pen)
	dtmp2 = sqrt((beta2pen-oldbeta2)*(beta2pen-oldbeta2))/sqrt(beta2pen*beta2pen)
	dtmp3 = sqrt((pi1pen-oldpi1)*(pi1pen-oldpi1))/sqrt(pi1pen*pi1pen)
	dtmp4 = sqrt((pi2pen-oldpi2)*(pi2pen-oldpi2))/sqrt(pi2pen*pi2pen)   
	nit2 = i
	if (i.ge.k2) dtmp4 = -dtmp1 + -dtmp2 + -dtmp3
	if ((dtmp1+dtmp2+dtmp3+dtmp4).ge.tol) goto 106
	mlrt = 0
        do 109  j = 1, n
        	tmp1 = x(j)*beta1pen
        	tmp2 = x(j)*beta2pen
        	tmp3 = x(j)*beta0
        	f1(j) = pi1pen*exp(-tmp1 )*beta1pen
        	f2(j) = pi2pen*exp(-tmp2 )*beta2pen 
        	f0(j) = beta0*exp(-tmp3)
        	mlrt = mlrt + log(f1(j)+f2(j))  
        	mlrt = mlrt - log(f0(j))
  109   continue
        mlrt = mlrt + c*log(4*pi1pen*pi2pen)
  	mlrt = 2*mlrt
	return
	end