#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NMAX 1024

#pragma goose func
double rsqrt(double r2)
{
  return 1.0/sqrt(r2);
}

#pragma goose func
double g_p3m(double re)
{
  double func, cppfrc;

  if(re<1.0){    
    func=re*(224.+re*re*(-224.+re*(70.+re*(48.-re*21.))))/(35.*4.0);
    cppfrc = 1.0-re*re*func;
  }else if(re<2.0){
    func=(-224.+re*(896.+re*(-840.+re*(224.+re*(70.+re*(-48.+re*7.))))))/(35.*4.0);    
    cppfrc = 1.0-12./(35.*4.0)-re*re*func;
  }else{
    func = 0.0;
    cppfrc = 0.0;
  }

  return cppfrc;
}

main(int argc, char **argv)
{
  static double x[NMAX][3];
  static double a[NMAX][3];
  static double m[NMAX];
  static double eps2i[NMAX];
  
  double r2, mrinv, mr3inv, rinv, eps2;
  int i, j, k,n;
  double dx[3];
  double gfunc,etainv,eta,re; 

  if(argc == 2){
    sscanf(argv[1],"%d",&n);
  }else{
    n = 16;
  }
  printf("n %d\n",n);

  eps2 = (1.0/128.0)*(1.0/128.0);
  for(i=0;i<n;i++){
    x[i][0] = drand48()-0.5;
    x[i][1] = drand48()-0.5;
    x[i][2] = drand48()-0.5;
    m[i] = 1.0/((double)n);
    if(i<(n/2)){
      eps2i[i] = eps2;
    }else{
      eps2i[i] = 2.0*eps2;
    }
  }
  eta = 0.3;
  etainv = 1.0/eta;

#pragma omp parallel for private(j,dx,r2,rinv,mrinv,mr3inv,k,re,gfunc)
#pragma goose parallel for precision ("double") loopcounter(i, j) result (a[i][0..2])
  for(i=0;i<n;i++) {
    for(k=0;k<3;k++) a[i][k] = 0.0;
    for (j=0;j<n;j++) {
        for(k=0;k<3;k++){
	  dx[k] = x[j][k] - x[i][k];
	  if(dx[k] > 0.5) dx[k] = dx[k] - 1.0;
	  if(dx[k] < -0.5) dx[k] = dx[k] + 1.0;
	}
	r2 = dx[0]*dx[0]+dx[1]*dx[1]+dx[2]*dx[2]+0.5*(eps2i[i]+eps2i[j]);
        rinv = rsqrt(r2);
	re = r2 * rinv * etainv ;
	gfunc = g_p3m(re);
        mrinv = m[j]*rinv;
        mr3inv = gfunc*mrinv*rinv*rinv;	
	a[i][0] += mr3inv * dx[0];
	a[i][1] += mr3inv * dx[1];
	a[i][2] += mr3inv * dx[2];      
    }
  }

  for(i=0;i<n;i++){
    printf("i %d a %22.14e %22.14e %22.14e\n",i,a[i][0],a[i][1],a[i][2]);
  }
}