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

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

main(int argc, char **argv)
{
  static double x[NMAX][3];
  static double v[NMAX][3];
  static double a[NMAX][3];
  static double adot[NMAX][3];
  static double a0[NMAX][3];
  static double adot0[NMAX][3];
  static double m[NMAX];
  static double pot[NMAX];
  static double ti[NMAX];  

  double r2, eps2;
  double xdotv,rinv,r3inv,r5inv,xdotvr5inv;
  int i, j, k,n;
  double dx[3],dv[3],x1j[3],v1j[3];
  double dt,dt2half,dt3over6,time;
      
  if(argc == 2){
    sscanf(argv[1],"%d",&n);
  }else{
    n = 17;
  }
  printf("n %d\n",n);

  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);
  }
  for(i=0;i<n;i++){
    v[i][0] = drand48()-0.5;
    v[i][1] = drand48()-0.5;
    v[i][2] = drand48()-0.5;
    a0[i][0] = drand48()-0.5;
    a0[i][1] = drand48()-0.5;
    a0[i][2] = drand48()-0.5;
    adot0[i][0] = drand48()-0.5;
    adot0[i][1] = drand48()-0.5;
    adot0[i][2] = drand48()-0.5;
    ti[i] = 0.25;
  }
  eps2 = 1.0/1024.0;
  eps2 *= eps2;
  time = 0.5;

#pragma omp parallel for private(j,k,dx,dv,xdotv,r2,rinv,r3inv,r5inv,xdotvr5inv)
#pragma goose parallel for precision ("double") loopcounter(i, j) \
  result (a[i][0..2],adot[i][0..2],pot[i])

  for(i=0;i<n;i++) {
    for(k=0;k<3;k++) a[i][k] = 0.0;
    pot[i] = 0.0;
    for (j=0;j<n;j++) {
      r2 = eps2;
      xdotv = 0.0;
      for(k=0;k<3;k++){
	dx[k] = x[j][k] - x[i][k];
	dv[k] = v[j][k] - v[i][k];
	r2 = r2 + dx[k] * dx[k];
	xdotv = xdotv + dx[k]*dv[k];
      }
      rinv = rsqrt(r2);
      //      if(i==j) rinv = 0.0;
      r3inv = rinv*rinv*rinv;
      r5inv = r3inv*rinv*rinv;
      xdotvr5inv = 3.0*xdotv*r5inv;
      if(i!=j){
	for(k=0;k<3;k++){
	  a[i][k] += m[j] * r3inv * dx[k];
          adot[i][k] += m[j] * (r3inv * dv[k] - xdotvr5inv * dx[k]); 
	}
	pot[i] -= m[j]*rinv;
      }
    }
  }

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