g5: bare metal C code

main_grav2.c

extern float base0;

float sqrtf(float);

void grav_kernel(float *x, float *y, float *z, float *m,
		 float *ax, float *ay, float *az, float *pt, float *res,
		 int n);

int main(int argc, char * argv[])
{
  float *base;
  float *x, *y, *z, *m, *ax, *ay, *az, *pt, *res;

  base = (float *)&base0;
  const int nword = 128;
  x   = base; base += nword;
  y   = base; base += nword;
  z   = base; base += nword;
  m   = base; base += nword;
  ax  = base; base += nword;
  ay  = base; base += nword;
  az  = base; base += nword;
  pt  = base; base += nword;
  res = base;

  grav_kernel(x, y, z, m, ax, ay, az, pt, res, nword);
  
  return 0;
}

void grav_kernel(float *x, float *y, float *z, float *m,
		 float *ax, float *ay, float *az, float *pt, float *res,
		 int n)
{
  float eps = 0.01f;
  float e2 = eps*eps;
  
  for(int i = 0; i < n; i++) {
    float xi = x[i];
    float yi = y[i];
    float zi = z[i];

    float a_x = 0.0f;
    float a_y = 0.0f;
    float a_z = 0.0f;
    float a_p = 0.0f;
    for(int j = 0; j < n; j++) {
      float xj = x[j];
      float yj = y[j];
      float zj = z[j];
      float mj = m[j];

      float dx, dy, dz;
      dx = xj - xi; 
      dy = yj - yi;
      dz = zj - zi;

      float r2 = dx*dx + dy*dy + dz*dz + e2;
      float r1i = 1.0f/sqrtf(r2);
      float r2i = r1i*r1i;
      float r1im = mj*r1i;
      float r3im = r1im*r2i;

      a_x += dx*r3im;
      a_y += dy*r3im;
      a_z += dz*r3im;
      a_p += -r1im;
    }

    ax[i] = a_x;
    ay[i] = a_y;
    az[i] = a_z;
    pt[i] = a_p;
  }

  float e_pot = 0.0;
  float eps_i = 1.0f/eps;
  for(int i = 0; i < n; i++){
    e_pot += m[i]*(pt[i] + m[i]*eps_i); 
  }  
  res[0] = e_pot*0.5f;
}

float sqrtf(float x)
{
  float t;
  asm("fsqrt.s %0, %1":"=rf"(t):"rf"(x));
  return t;
}