//Two sample Kolmogorov Smirnov test
//Written by Raymond N. Greenwell, March 2000
#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <fstream.h>
#include <iomanip.h>
int n1,n2;//number of data per test
double x[1001],y[1001],//data values,
alpha=.01,
correct, correct1;
double KS();
double max(double, double);
double min(double, double);
double pvalue(double);
int lcd(int,int);

int main()
{
int ii,//test number
jj,//data number
kk,n2count,
nu,//degrees of freedom for t distribution
nuval,
reject,//number of times p<alpha
count;//number of times p=pmin
double T,//test statistic
rejected,//total number of times H0 rejected
val,
crit,//critical KS value
pmin,//minimum p>=alpha
pmax,//maximum p<alpha
pmax1,pmin1, //with continuity correction
pmax2,pmin2, //using Kim correction
badness;
crit=sqrt(-.5*log(alpha/2));
cout << "n1 " << "n2" << " alphaL " << " alphaR " << " alphaL' " << " alphaR' "<< endl;
//for (n1=40; n1<=160; n1=n1+40)
//cout << "correct " << " pmax " << "pmin" << endl;
for (n1=40; n1<=160; n1=n1+40)
{
for (n2count=1; n2count<=5; n2count++)
{
switch(n2count)
{case 1: n2=n1/4; break;
case 2: n2=2*n1/5; break;
case 3: n2=n1/2; break;
case 4: n2=3*n1/4; break;
case 5: n2=n1; break;
default: n2=n1;
}
val=lcd(n1,n2)*sqrt((n1+n2)/(1.0*n1*n2));
if (n1%n2==0) correct=n1*0.032/(pow(n2,1.5));
else correct=0.25/(sqrt(n1));
pmax1=pvalue(ceil(val*(crit-correct))/val+correct);
pmin1=pvalue(floor(val*(crit-correct))/val+correct);
correct1=correct;
if (n1>2*n2) correct=1/(2*sqrt(n1));
else if (n1%n2==0) correct=2/(3*sqrt(n1));
else correct=2/(5*sqrt(n1));
pmax2=pvalue(ceil(val*(crit-correct))/val+correct);
pmin2=pvalue(floor(val*(crit-correct))/val+correct);
//cout << correct << " " << pmax2 << " " << pmin2 << endl;
pmax=pvalue(ceil(val*(crit))/val);
pmin=pvalue(floor(val*(crit))/val);
//cout << n1 << " & "<< n2 << " & "<< pmax << " & "<< pmax1 <<" & "<< pmax2 <<" \\ "<< endl;
cout << n1 << " & "<< n2 << " & "<< pmin<<" & "<< pmin1 <<" & "<< pmin2<<" \\\\ "<< endl;
//cout << "val = " << val << " crit = " << crit <<" Our correction = " << correct1 << " Kim's correction = " << correct << endl;
}//end for n2
}//end for n2
}

double KS()
//Find maximum value of difference between the two
//sample distribution functions.
{
double fmax=0;//largest value found so far
int count1=1,//data number
count2=1;
while ((count1<=n1) && (count2<=n2))
{
if (x[count1]<y[count2])
{fmax=max(fmax, fabs((count1-1.0)/n1-(count2-1.0)/n2));
fmax=max(fmax, fabs(count1/(1.0*n1)-(count2-1.0)/n2));
count1++;
}
else{if (x[count1]>y[count2])
{fmax=max(fmax, fabs((count1-1.0)/n1-(count2-1.0)/n2));
fmax=max(fmax, fabs((count1-1.0)/n1-count2/(1.0*n2)));
count2++;
}
else {fmax=max(fmax, fabs((count1-1.0)/n1-(count2-1.0)/n2));
fmax=max(fmax, fabs(count1/(1.0*n1)-count2/(1.0*n2)));
fmax=max(fmax, fabs((count1-1.0)/n1-(count2)/(1.0*n2)));
fmax=max(fmax, fabs(count1/(1.0*n1)-(count2-1.0)/(1.0*n2)));
count1++;count2++;
}//end second else
}//end first else
}//end while
return fmax;
}


double max(double aa, double bb)
{double mm;
mm=aa;
if (bb>aa) mm=bb;
return mm;}

double min(double aa, double bb)
{double mm;
mm=bb;
if (bb>aa) mm=aa;
return mm;}


double pvalue(double x)
//Calculate pvalue based on KS coefficient
//See Fisz, p. 396
{
return 2*exp(-2*x*x)-2*exp(-8*x*x);
}

int lcd(int n1, int n2)
{int r,//gcd
m,n;
m=n1; n=n2;
do
{r=m%n;
m=n;
n=r;
}
while (r>0);
return n1*n2/m;
}//end lcd