Skip to content

Instantly share code, notes, and snippets.

@dadeba
Last active Jun 27, 2019
Embed
What would you like to do?
g5: bare metal C code
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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment