kly=function(time,cens)
{
################################
# based on Kochar, Lam and Yip #
# Lifetime Data Analysis,2002  #
# time=time to first event     #
# cens=censor variable         #
#      1=one type of event     #
#      2=another type of event #
#      0=the rest of events and#
#        censored records      #
################################       
dd0=table(time,cens)
if (sum(cens==0)==0) {dd=cbind(rep(0,nrow(dd0)),dd0)}
else {dd=dd0}
if (sum(cens==1)==0 | sum(cens==2)==0 | sum(cens>2)>0) 
{stop ("Codes for events should be 1 or 2")}
n=length(time)
sigma=sqrt(n*sum(cens>0))
nt=apply(dd,1,sum)
ntt=cumsum(nt)
atrisk=n-ntt
lnti=(dd[,3]-dd[,2])
lnt=cumsum(lnti)
k=nrow(dd)
maxlnt=vector("numeric",k)

for (i in 1:k)
{
maxlnt[i]=max(abs(lnt-lnt[i]))
}

cnw=sqrt(n)*max(maxlnt)/sigma
## for C, (without absolute value)
pcnw=function(cc)
{
k=c(0:15)
parti=(-1)^k*exp(-pi^2*(2*k+1)^2/(8*cc^2))/(2*k+1)
p=sum(parti)*4/pi
return(p)
}

## for Cstar (with absolute value)
ph=function(z)
{
f1=function(x1,x2,k,xi)
{
ki=c(1:k) 
x=c((x1*xi):(x2*xi))/xi 
n=length(x)
h=vector("numeric",n)

for (i in 1:n) 
{ 
h[i]=sum(8*(-1)^(ki-1)*ki^2*dnorm(ki*x[i]))
} 

parth=sum(h)/xi
return(parth)}

if (z<=0.1) 
ph=f1(0.01,z,1000,100)
else if (z<=0.5)
ph=f1(0.01,0.1,1000,100)+f1(0.11,z,100,100)
else ph=f1(0.01,0.1,1000,100)+f1(0.11,0.5,100,100)+
f1(0.51,z,30,1000)
return(ph)
}

pvalue=1-ph(cnw)
if (pvalue>=0.00001) r=data.frame(cnw,pvalue)
if (pvalue<0.00001) r=data.frame(cnw,pvalue="p-value<0.00001")
row.names(r)=""
return(r)
}