#include <string.h>#include <stdio.h>#include <stdlib.h>#include <math.h>void make_box1(double *x, int N, double scal, int bs, 	       int *box, int *lis, int *mxi) {  int i, ix;  for (i=0;i<=bs;i++) box[i]=-1;  for (i=0;i<=bs;i++) mxi[i]=0;  for(i=0;i<N;i++) {ix=(int)(x[i]*scal);lis[i]=box[ix]; box[ix]=i; mxi[ix]++;}  for (i=1;i<=bs;i++) mxi[i]+=mxi[i-1];}void make_box2(double **x, int dim, int N, int comp1, int comp2, int bs, int inveps, 	       int **box, int *lis) {  int d,i,ix,iy,ixy;  int ib=bs-1;  double **xx;  xx=(double **)calloc(dim, sizeof(double*));   for (d=0;d<dim;d++) {    xx[d]=(double *)calloc(N, sizeof(double));     memcpy(xx[d],x[d],N*sizeof(double));  }  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) box[ix][iy] = -1;    for (i=0;i<N;i++) {    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;    lis[i]=box[ix][iy];box[ix][iy]=i;  }  i=-1;  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) {    ixy=box[ix][iy];    while(ixy>=0) {      i++;      for (d=0;d<dim;d++) x[d][i]=xx[d][ixy];      ixy=lis[ixy];    }    box[ix][iy]=-1;  }  for (i=0;i<N;i++) {    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;    lis[i]=box[ix][iy];box[ix][iy]=i;  }  for (d=0;d<dim;d++) free(xx[d]); free(xx);}void make_box2ind(double **x, int dim, int N, int comp1, int comp2, int bs, int inveps, 		  int *ind,  int **box, int *lis) {  int d,i,ix,iy,ixy;  int ib=bs-1;  double **xx;  xx=(double **)calloc(dim, sizeof(double*));   for (d=0;d<dim;d++) {    xx[d]=(double *)calloc(N, sizeof(double));     memcpy(xx[d],x[d],N*sizeof(double));  }  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) box[ix][iy] = -1;    for (i=0;i<N;i++) {    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;    lis[i]=box[ix][iy];box[ix][iy]=i;  }  i=-1;  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) {    ixy=box[ix][iy];    while(ixy>=0) {      i++;      for (d=0;d<dim;d++) x[d][i]=xx[d][ixy];      ind[ixy]=i;ixy=lis[ixy];    }    box[ix][iy]=-1;  }  for (i=0;i<N;i++) {    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;    lis[i]=box[ix][iy];box[ix][iy]=i;  }  for (d=0;d<dim;d++) free(xx[d]); free(xx);}int neiE1(double *x, int i, double scal, int bs, double eps,  int *box, int *lis, int *mxi) {  double dd;  int mi,mp,mm,nx=0;  double xc=x[i];  mp=int((xc+eps)*scal); if(mp>bs) mp=bs;  mm=int((xc-eps)*scal); if(mm<0)  mm=0;  mi=box[mp];  while(mi>=0) {    dd=x[mi]-xc;if(fabs(dd)<=eps) nx++;    mi=lis[mi];  }  if(mm>=mp) return nx-1;  mi=box[mm];  while(mi>=0) {    dd=xc-x[mi];if(fabs(dd)<=eps) nx++;    mi=lis[mi];  }  nx+=mxi[mp-1]-mxi[mm];  return nx-1;}int neiE(double **x, int i, int comp1, int comp2, int dim, int bs, double epsgrid, double eps, int **box, int *lis) {  int ix,iy,ix1,iy1,ix2,jj,step,d;  int el,nx,ib=bs-1;  double dd,dy;  double *xx;  xx=(double *)calloc(dim,sizeof(double));  for (d=0;d<dim;d++) xx[d]=x[d][i];  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;  jj=0; nx=0;  while (eps>epsgrid*(jj-1)) {    step = (jj) ? 2*jj : 1;    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {      ix2=ix1&ib;      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {	el=box[ix2][iy1&ib];	while (el != -1) { 	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;	  if (dd<=eps) nx++;	  el=lis[el];	}      }    }    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {      ix2=ix1&ib;      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {	el=box[ix2][iy1&ib];	while (el != -1) { 	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;	  if (dd<=eps) nx++;	  el=lis[el];	}      }    }    jj++;    if (jj==(bs/2)) break;  }  if ( jj==(bs/2) ) { //half of the layer    for (ix1=ix-jj;ix1<ix+jj;ix1++) {      ix2=ix1&ib; iy1=iy-jj;      el=box[ix2][iy1&ib];      while (el != -1) { 	dd=fabs(xx[0]-x[0][el]); 	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;	if (dd<=eps) nx++;	el=lis[el];      }    }    ix1=ix-jj; ix2=ix1&ib;    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {      el=box[ix2][iy1&ib];      while (el != -1) { 	dd=fabs(xx[0]-x[0][el]); 	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;	if (dd<=eps) nx++;	el=lis[el];      }    }  }  free(xx);  return nx-1;}void neiEK(double **x, int i, int comp1, int comp2, int dim, int K,	   int bs, double epsgrid, double *eps, int **box, int *lis,	   int *nx) {  int ix,iy,ix1,iy1,ix2,jj,step,d,ik;  int el,ib=bs-1;  double dd,dy;  double *xx;  xx=(double *)calloc(dim,sizeof(double));  for (d=0;d<dim;d++) xx[d]=x[d][i];  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;  jj=0;   for (ik=0;ik<K;ik++) nx[ik]=0;  while (eps[K-1]>epsgrid*(jj-1)) {    step = (jj) ? 2*jj : 1;    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {      ix2=ix1&ib;      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {	el=box[ix2][iy1&ib];	while (el != -1) { 	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;; 	  for (ik=0;ik<K;ik++) {	    if (dd<=eps[ik]) nx[ik]++;	  }	  el=lis[el];	}      }    }    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {      ix2=ix1&ib;      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {	el=box[ix2][iy1&ib];	while (el != -1) { 	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;	  for (ik=0;ik<K;ik++) {	    if (dd<=eps[ik]) nx[ik]++;	  }	  el=lis[el];	}      }    }    jj++;    if (jj==(bs/2)) break;  }  if ( jj==(bs/2) ) { //half of the layer    for (ix1=ix-jj;ix1<ix+jj;ix1++) {      ix2=ix1&ib; iy1=iy-jj;      el=box[ix2][iy1&ib];      while (el != -1) { 	dd=fabs(xx[0]-x[0][el]); 	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;	for (ik=0;ik<K;ik++) {	  if (dd<=eps[ik]) nx[ik]++;	}	el=lis[el];      }    }    ix1=ix-jj; ix2=ix1&ib;    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {      el=box[ix2][iy1&ib];      while (el != -1) { 	dd=fabs(xx[0]-x[0][el]); 	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;	for (ik=0;ik<K;ik++) {	  if (dd<=eps[ik]) nx[ik]++;	}	el=lis[el];      }    }  }  free(xx);  for (ik=0;ik<K;ik++) nx[ik]-=1;}void neiK(double **x, int dim, int comp1, int comp2, int i, 	  int bs, double epsgrid, int K, int **box, int *lis,	  int *nn) {  double *dn,*xx;  int k,ix,iy,ix1,iy1,ix2,jj,step,ib=bs-1;  int el;  double dd,dy;  int d;  dn=(double*)calloc(K+1,sizeof(double));  xx=(double*)calloc(dim,sizeof(double));  for(k=0;k<dim;k++) xx[k]=x[k][i];  dn[0]=0;  for(k=1;k<=K;k++) dn[k]=1e30;  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;  jj=0;   while (dn[K]>epsgrid*(jj-1)) {    step = (jj) ? 2*jj : 1;    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {      ix2=ix1&ib;      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {	el=box[ix2][iy1&ib];	while (el != -1) { if (el!=i) {	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;	  if(dd<dn[K]) {	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }	    dn[k+1]=dd; nn[k+1]=el;	  }	} el=lis[el];	}      }    }    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {      ix2=ix1&ib;      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {	el=box[ix2][iy1&ib];	while (el != -1) { if (el!=i) {	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;	  if ( (dy=fabs(xx[1]-x[1][el]))>dd ) dd=dy;	  if(dd<dn[K]) {	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }	    dn[k+1]=dd; nn[k+1]=el;	  } 	} el=lis[el];	}      }    }    jj++;    if (jj==bs/2) break;  }  if ( jj==(bs/2) ) { //half of the layer    for (ix1=ix-jj;ix1<ix+jj;ix1++) {      ix2=ix1&ib; iy1=iy-jj;      el=box[ix2][iy1&ib];      while (el != -1) { 	if (el!=i) {	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;	  if(dd<dn[K]) {	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }	    dn[k+1]=dd; nn[k+1]=el;	  }	}	el=lis[el];      }    }    ix1=ix-jj; ix2=ix1&ib;    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {      el=box[ix2][iy1&ib];      while (el != -1) { 	if (el!=i) {	  dd=fabs(xx[0]-x[0][el]); 	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;	  if(dd<dn[K]) {	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }	    dn[k+1]=dd; nn[k+1]=el;	  }	} 	el=lis[el];      }    }  }  free(dn);free(xx);}/*********/void mi2(double **x, int N, int K,	 double *psi, 	 double *scal,	 double *mic, double *mir) {    int i,k,*n1,*n2;  double *xc,dx;  double *eps,Eps;  int maxdim;  double dxy1,dxy2;  int *nn;  int d;    int BOX,BOX1;  int **box,*lis; // two dimensional boxes  int **box1; // onedimensional boxes  int **lis1; // lists for one dimensions  int **mxi; //accumulative lists of points in oned boxes  double epsilon;  int inveps;  int dim=2;  double *phi;  phi=(double*)calloc(K+1,sizeof(double));  for (i=1;i<=K;i++) phi[i]=psi[i]-(dim-1)/(double(i));   //   nn=(int*)calloc(K+1,sizeof(int));  xc=(double*)calloc(dim+1,sizeof(double));  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;  epsilon=4.0/BOX;  inveps=BOX/4;  BOX1=N-5;  box1=(int**)calloc(dim,sizeof(int*));   lis1=(int**)calloc(dim,sizeof(int*));   mxi=(int**)calloc(dim,sizeof(int*));   for (d=0;d<dim;d++) {    box1[d]=(int*)calloc(BOX1+1,sizeof(int));     lis1[d]=(int*)calloc(N,sizeof(int));     mxi[d]=(int*)calloc(BOX1+1,sizeof(int));  }   box=(int**)calloc(BOX,sizeof(int*));  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));  lis=(int*)calloc(N,sizeof(int));  eps=(double*)calloc(dim,sizeof(double));  n1=(int*)calloc(dim,sizeof(int));  n2=(int*)calloc(dim,sizeof(int));  make_box2(x,dim,N,0,1,BOX,inveps,box,lis); //for searching neighbours in prodict space  for (d=0;d<dim;d++) make_box1(x[d],N,scal[d],BOX1,box1[d],lis1[d],mxi[d]);          dxy1=dxy2=0.0;  for (i=0;i<N;i++) {    for (d=0;d<dim;d++) xc[d]=x[d][i];    neiK(x,dim,0,dim-1,i,BOX,epsilon,K,box,lis,nn);        Eps=0;maxdim=-1;    for (d=0;d<dim;d++) {      eps[d]=0;      for(k=1;k<=K;k++) {if( (dx=fabs(xc[d]-x[d][nn[k]]))>eps[d] ) eps[d]=dx; }      if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}    }	    for (d=0;d<dim;d++) {      n2[d]=neiE1(x[d],i,scal[d],BOX1,eps[d],box1[d],lis1[d],mxi[d]);      if (d==maxdim) { n1[d]=n2[d]; dxy1+=psi[n1[d]]; }      else { 	n1[d]=neiE1(x[d],i,scal[d],BOX1,Eps,box1[d],lis1[d],mxi[d]); 	dxy1+=psi[n1[d]+1];      }      dxy2+=psi[n2[d]];    }    //    fprintf(stdout,"%f %f %d %d %d %d %d %f %d\n",xc[0],xc[1],(eps[0]<=eps[1]),n1[0],n1[1],n2[0],n2[1],Eps,nn[K]);  }  dxy1/=N;*mic=psi[N]+psi[K]-dxy1;  dxy2/=N;*mir=psi[N]+phi[K]-dxy2;    free(xc);free(nn);  for (i=0;i<BOX;i++) free(box[i]); free(box);  free(lis);  for (d=0;d<dim;d++) {    free(box1[d]);free(lis1[d]);free(mxi[d]);  }  free(box1);free(lis1);free(mxi);  free(eps);free(n1);free(n2);  free(phi);}void mi2c(double **x, int N, int K,	  double *psi, 	  double *scal,	  double *mic) {    int i,k,*n1;  double *xc,dx;  double *eps,Eps;  int maxdim;  double dxy1;  int *nn;  int d;    int BOX,BOX1;  int **box,*lis; // two dimensional boxes  int **box1; // onedimensional boxes  int **lis1; // lists for one dimensions  int **mxi; //accumulative lists of points in oned boxes  double epsilon;  int inveps;  int dim=2;  nn=(int*)calloc(K+1,sizeof(int));  xc=(double*)calloc(dim+1,sizeof(double));  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;  epsilon=4.0/BOX;  inveps=BOX/4;  BOX1=N-5;  box1=(int**)calloc(dim,sizeof(int*));   lis1=(int**)calloc(dim,sizeof(int*));   mxi=(int**)calloc(dim,sizeof(int*));   for (d=0;d<dim;d++) {    box1[d]=(int*)calloc(BOX1+1,sizeof(int));     lis1[d]=(int*)calloc(N,sizeof(int));     mxi[d]=(int*)calloc(BOX1+1,sizeof(int));  }   box=(int**)calloc(BOX,sizeof(int*));  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));  lis=(int*)calloc(N,sizeof(int));  eps=(double*)calloc(dim,sizeof(double));  n1=(int*)calloc(dim,sizeof(int));  make_box2(x,dim,N,0,1,BOX,inveps,box,lis); //for searching neighbours in prodict space  for (d=0;d<dim;d++) make_box1(x[d],N,scal[d],BOX1,box1[d],lis1[d],mxi[d]);          dxy1=0.0;  for (i=0;i<N;i++) {    for (d=0;d<dim;d++) xc[d]=x[d][i];    neiK(x,dim,0,dim-1,i,BOX,epsilon,K,box,lis,nn);        Eps=0;maxdim=-1;    for (d=0;d<dim;d++) {      eps[d]=0;      for(k=1;k<=K;k++) {if( (dx=fabs(xc[d]-x[d][nn[k]]))>eps[d] ) eps[d]=dx; }      if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}    }	    for (d=0;d<dim;d++) {      n1[d]=neiE1(x[d],i,scal[d],BOX1,Eps,box1[d],lis1[d],mxi[d]);       if (d==maxdim) { 	dxy1+=psi[n1[d]];      } else { 	dxy1+=psi[n1[d]+1];      }    }    //    fprintf(stdout,"%f %f %d %d %d %f %d\n",xc[0],xc[1],(eps[0]<=eps[1]),n1[0],n1[1],Eps,nn[K]);  }  dxy1/=N;*mic=psi[N]+psi[K]-dxy1;    free(xc);free(nn);  for (i=0;i<BOX;i++) free(box[i]); free(box);  free(lis);  for (d=0;d<dim;d++) {    free(box1[d]);free(lis1[d]);free(mxi[d]);  }  free(box1);free(lis1);free(mxi);  free(eps);free(n1);}void mi2r(double **x, int N, int K,	  double *psi, 	  double *scal,	  double *mir) {    int i,k,*n2;  double *xc,dx;  double *eps;  double dxy2;  int *nn;  int d;    int BOX,BOX1;  int **box,*lis; // two dimensional boxes  int **box1; // onedimensional boxes