#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <assert.h>
#include <mpi.h>
/* #include "problem2.c" */


/* 
 * SuperCon2001 Honsen Program
 * by KaLa
 *
 */



#define EOC 0x7fffffff
#define NEWTYPE 0
#define OLDTYPE 1

#define PX(a) (pxy[(a)*2])
#define PY(a) (pxy[(a)*2+1])
#define PVX(a) (pvxy[(a)*2])
#define PVY(a) (pvxy[(a)*2+1])


/* #define SC2001OUT */
#ifdef DEBUG
#define TRACE(a,b) (my_proc==3?printf(a,b): 0 )
#else
#define TRACE(a,b) 
#endif
typedef struct
{
  int x;
  int n_box;
  int box;
} ROW;

typedef struct
{
  int y;
  int n_star;
  int star;
  int next;
} BOX;

double DT;
int n_steps;
int n_particles;
void supercon2k1_get_config(int problem_num, int * size_p, 
			   int * n_steps_p, double * dt_p);
void supercon2k1_get_data(int problem_num, double * a);
void supercon2k1_verify(int problem_num, double * a);
double round(double x);
void do_step();
void create_array();
int find_row(int *row,int rx);
int find_box(int nbox,BOX **pbox,BOX **pbeforebox,int ry,int size);
void calc_gravity(int i, int j, double *ax, double *ay);


double *pmass;
double *pxy;
double *pvxy;
double *pay;
double *pax;

ROW *prow;
BOX *pbox;
int *next_star;
int n_row;
int n_col;

int my_proc;
int n_procs;

int start_index[33];

int main(int argc, char ** argv)
{
  int problem_num, i;
  double *a,*ap;

  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD,&n_procs);
  MPI_Comm_rank(MPI_COMM_WORLD,&my_proc);


  if(argc < 2){
    fprintf(stderr,"bad args...\n");
    exit(1);
  }
  problem_num = atoi(argv[1]);

  supercon2k1_get_config(problem_num,&n_particles,&n_steps,&DT);

/*
  if(my_proc == 0){
    printf("problem_num=%d: n_particle=%d, n_steps=%d, DT=%g\n",
	 problem_num, n_particles, n_steps, DT);
  }
*/
  a = (double *)malloc(sizeof(double) * 5 * n_particles);

  pmass = (double *)malloc(sizeof(double) * n_particles);
  pxy = (double *)malloc(sizeof(double) * n_particles * 2);
  pvxy = (double *)malloc(sizeof(double) * n_particles *2);
  pax = (double *)malloc(sizeof(double) * n_particles);
  pay = (double *)malloc(sizeof(double) * n_particles);

  next_star = (int *)malloc(sizeof(int) * n_particles);
  
  prow = (ROW *)malloc(sizeof(ROW) * n_particles);
  pbox = (BOX *)malloc(sizeof(BOX) * n_particles *2); /* ! */

  if(a == NULL || pxy==NULL || pmass==NULL 
     || pvxy==NULL || pax==NULL || pay==NULL || next_star==NULL 
     || prow==NULL || pbox == NULL){
    fprintf(stderr,"no memory...\n");
    exit(1);
  }

  supercon2k1_get_data(problem_num,a);

  ap = a;
  for(i = 0; i < n_particles; i++) {
    pmass[i] = *ap++;
    PX(i) = *ap++;
    PY(i) = *ap++;
    PVX(i) = *ap++;
    PVY(i) = *ap++;
    
  }


  for(i = 0; i< n_procs; i++){
    start_index[i] = (n_particles/n_procs)*i;
  }
  start_index[n_procs] = n_particles;

  for(i = 0; i < n_steps; i++){
    do_step();
#ifdef SC2001OUT
    if(my_proc==0)
      {
	int xx;
	for (xx=0; xx<n_particles; xx++)   /* don't worry the order of */
	  printf("%d %d %5f %5f\n", my_proc ,i, PX(xx), PY(xx));  /* the output */
      }
#endif
  }
  
  if(my_proc == 0){
    ap = a;
    for(i = 0; i < n_particles; i++) {
      *ap++ = pmass[i];
      *ap++ = PX(i);
      *ap++ = PY(i);
      *ap++ = PVX(i);
      *ap++ = PVY(i);
    }
    supercon2k1_verify(problem_num, (double *)a);
  }

  MPI_Finalize();
  return 0;
}
int Nbox;
int newbox(void)
{
  if(Nbox<n_particles || Nbox>=n_particles*4)
  {
    fprintf(stderr,"FATAL ERROR IN NEWBOX\n");
    exit(-1);
  }
  return Nbox++;
}
void calc_star_box(int si,BOX *sbox,double *ax,double *ay)
{
  double dx,dy;
  double norm2,cos,sin,acc;
  int ji;
  
  for(ji=sbox->star; ji!=EOC; ji=next_star[ji]){
    if(si==ji) continue;
    TRACE("ji:%d\n",ji);
    dx = PX(ji) - PX(si);
    dy = PY(ji) - PY(si);
    
    if(dx > 0) dx = ceil(dx); else dx = floor(dx); 
    if(dy > 0) dy = ceil(dy); else dy = floor(dy); 

    norm2 = dx * dx + dy * dy;
    assert(norm2);
    if(norm2 < 1024.0) {
      cos = dx / sqrt(norm2);
      sin = dy / sqrt(norm2);
      acc = pmass[ji] / norm2;
      *ax += round(acc * cos);
      *ay += round(acc * sin);
       }/* end of ji loop */
  }
}
#if 0
#define TIME() cl2=clock();if(my_proc==3)printf("%d\n",(int)(cl2-cl));cl=cl2
#else
#define TIME()
#endif

void do_step()
{
  int si;
  clock_t cl,cl2;
  int rx,ry;
  int row,box;
  BOX *pb,*pbbef;
  int rank;
  int i3;
  int tmprow,tmpbox;
  cl=clock();
/*  printf("CLOCK START :%d\n",my_proc);*/

/*  fprintf(stderr,"PROC%d\n",my_proc);*/

  TRACE("loop:%d\n",0);
  create_array();
  TRACE("array created:%d\n",0);

/*  fprintf(stderr,"PROC(2)%d\n",my_proc);*/

  TIME();
/*  printf("CA OK!\n");*/
  for(si=start_index[my_proc];si<start_index[my_proc+1];si++){
    TRACE("si:%d\n",si);
/*    if(my_proc==0)printf("%d\n",si);*/
    rx = (int)floor(PX(si)/32.0);
    ry = (int)floor(PY(si)/32.0);
    pax[si] = 0.0;
    pay[si] = 0.0;
    
    for(tmprow=-1;tmprow<=1;tmprow++){
      if(find_row(&row,rx+tmprow)==NEWTYPE) continue;
      if(find_box(row,&pb,&pbbef,ry,prow[row].n_box)==NEWTYPE)
	{
	  if(find_box(row,&pb,&pbbef,ry-1,prow[row].n_box)==OLDTYPE)
	    calc_star_box(si,pb,&pax[si],&pay[si]);
	  if(find_box(row,&pb,&pbbef,ry+1,prow[row].n_box)==OLDTYPE)
	    calc_star_box(si,pb,&pax[si],&pay[si]);
	}
      else
	{
	  calc_star_box(si,pb,&pax[si],&pay[si]);
	  if(pbbef->y==ry-1 && pbbef->n_star>0)
	    calc_star_box(si,pbbef,&pax[si],&pay[si]);
	  if(pb->next>=0 && pbox[pb->next].y==ry+1)
	    calc_star_box(si,&pbox[pb->next],&pax[si],&pay[si]);
	}
    }
  }
  TRACE("ax calced:%d\n",0);
  TIME();
  for(si=start_index[my_proc];si<start_index[my_proc+1];si++){
    PX(si) += PVX(si) * DT;
    PY(si) += PVY(si) * DT;
    PVX(si) += pax[si] * DT;
    PVY(si) += pay[si] * DT;
  }
  TIME();
  
  for(rank=0;rank<n_procs;rank++)
  {
    MPI_Bcast(pxy+(start_index[rank]*2),(start_index[rank+1]-start_index[rank])*2,
        MPI_DOUBLE,rank,MPI_COMM_WORLD);
    MPI_Bcast(pvxy+(start_index[rank]*2),(start_index[rank+1]-start_index[rank])*2,
        MPI_DOUBLE,rank,MPI_COMM_WORLD);
  }

  TIME();
}

int find_row(int *row,int rx)
{
  int left_row,right_row,middle_row;

  if( n_row == 0 ){
    *row = 0;
    return NEWTYPE;
  }

  right_row = n_row-1;
  left_row = 0;
  while(1){
    if(right_row < left_row){
      *row = left_row;
      return NEWTYPE;
    }
    if(right_row == left_row){
      if(prow[left_row].x==rx){
	*row=left_row;
	return OLDTYPE;
      }
      else if(prow[left_row].x>rx){
	*row = left_row;
	return NEWTYPE;
      }
      else{ /* prow[left_row].x < rx */
	*row = left_row+1;
	return NEWTYPE;
      }
    }

    middle_row = (left_row+right_row)/2;
    if(prow[middle_row].x < rx){
      left_row = middle_row + 1;
    }
    else if(prow[middle_row].x > rx){
      right_row = middle_row -1;
    }
    else{
      *row = middle_row;
      return OLDTYPE;
    }
  }
}

int find_box(int nline,BOX **pboxo,BOX **pbeforebox,int ry,int size)
{
  int nbox;
  BOX *p,*pbefore;
  int count=size;
  int nbef;
  pbefore=&pbox[prow[nline].box];
  nbef=prow[nline].box;
  for(nbox=pbefore->next;nbox>=0;nbox=p->next,count--){
    if(nbox>=n_particles*4 || count<=0)
      {
/*	fprintf(stderr,"UG\n");*/
	break;
      }
    p=&pbox[nbox];
    if(ry<p->y)
      {
	*pboxo=pbefore;
	return NEWTYPE;
      }
    else if(ry==p->y)
      {
	*pboxo=p;
	*pbeforebox=pbefore;
	return OLDTYPE;
      }
    pbefore=p;
    nbef=nbox;
    }
  *pboxo=pbefore;
  return NEWTYPE;
}

void create_array()
{
 
  int si;
  int nbox;
  int rx,ry;
  int i;
  int insrow;
  BOX *pinsbox,*pbefbox,*pinsbox2;
  int T=0;
  int i5;
  n_row = 0;
  n_col = 0;
  Nbox=n_particles;
  for(i=0;i<n_particles;i++)
    {
      pbox[i].next=-1;
      pbox[i].n_star=0;
    }
  for(si = 0 ;si < n_particles; si++){
    TRACE("si:%d\n",si);
/*    fprintf(stderr,"si:%d\n",si);*/
/*    MPI_Bcast(&i5,1,MPI_INT,0,MPI_COMM_WORLD);*/
    
/*    if(my_proc==0)printf("%d\n",si);*/
    rx = (int)floor(PX(si)/32.0);
    ry = (int)floor(PY(si)/32.0);

    i5=find_row(&insrow,rx);
#define   TT()

    if(i5 == OLDTYPE){
      TT();
      if(find_box(insrow,&pinsbox,&pbefbox,ry,prow[insrow].n_box) == OLDTYPE){
	TT();
	next_star[si] = pinsbox->star;
	TT();
	pinsbox->star = si;
	TT();
	pinsbox->n_star++;
	TT();
      }
      else{ /* LINE=old BOX=new */
	int nbox;
	TT();
	nbox=newbox();
	pinsbox2=&pbox[nbox];
	pinsbox2->next=pinsbox->next;
	pinsbox2->y=ry;
	pinsbox2->n_star=1;
	pinsbox2->star=si;
	pinsbox->next=nbox;

	next_star[si]=EOC;
	prow[insrow].n_box++;
      }
    }
    else{ /* LINE=box=new */
      TT();
      memmove(&prow[insrow+1],&prow[insrow],
	      sizeof(ROW)*(n_row-insrow));
      prow[insrow].x=rx;
      TT();
      prow[insrow].n_box=1;
      prow[insrow].box=n_row;
      TT();
/*      if(n_row>=n_particles)
	{
	  printf("FATAL ERROR!\n");
	}*/
/*	printf("CREATE:%d\n",n_row);*/
      
	nbox=newbox();
/*      if(pbox[n_row].next>=0)
	printf("ERROR! %d pboxhead=%d\n",n_row,pbox[n_row].next);*/
	pbox[n_row].next=nbox;
	pinsbox=&pbox[nbox];
	pinsbox->y=ry;
	pinsbox->n_star=1;
	pinsbox->star=si;
	pinsbox->next=-1;

      TT();
      next_star[si]=EOC;
      n_row++;
    }
  } /* end of si loop */
}


double round(double x){ 
  if (x < 0.0) {
    return - floor (-x + 0.5);
  } else {
    return floor (x + 0.5);
  }
}