Compiler Explorer

Source code

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

// uncommend to toggle implementation
#define TOGGLE_AOS     (1<<0)
#define TOGGLE_SOA     (1<<1)
#define TOGGLE_float3  (1<<2)
#define TOGGLE_float4  (1<<3)
#define TOGGLE_vec4    (1<<4)
#define TOGGLE_floats3 (1<<5)

#ifdef __FAST_MATH__
// otherwise it won't get inlined on RV
#define sqrtf __builtin_sqrtf
#endif

#if (N) & (1<<0)
struct planet {
	float x, y, z;
	float vx, vy, vz;
	float mass;
};

void
advance_AOS(struct planet *bodies, float dt, size_t nbodies)
{
	float rdt = 1/dt;
	for (size_t i = 0; i < nbodies; ++i) {
		struct planet *b = &bodies[i];
		float x1 = b->x, y1 = b->y, z1 = b->z;
		float vx1 = b->vx, vy1 = b->vy, vz1 = b->vz;
		float mass1 = b->mass;
		for (size_t j = i + 1; j < nbodies; ++j) {
			struct planet *b2 = &bodies[j];
			float dx = x1 - b2->x;
			float dy = y1 - b2->y;
			float dz = z1 - b2->z;
			float normsq = dx * dx + dy * dy + dz * dz;
			float dist = sqrtf(normsq);
			float mag = rdt * normsq * dist;
			float m2 = b2->mass * mag;
			float m1 = mass1 * mag;
			vx1 -= dx * m2;
			vy1 -= dy * m2;
			vz1 -= dz * m2;
			b2->vx += dx * m1;
			b2->vy += dy * m1;
			b2->vz += dz * m1;
		}
		b->vx = vx1;
		b->vy = vy1;
		b->vz = vz1;
	}
	for (size_t i = 0; i < nbodies; ++i) {
		struct planet * b = &bodies[i];
		b->x += dt * b->vx;
		b->y += dt * b->vy;
		b->z += dt * b->vz;
	}
}
#endif

#if (N) & (1<<1)
struct planets {
	float *x, *y, *z;
	float *vx, *vy, *vz;
	float *mass;
};

void
advance_SOA(struct planets bodies, float dt, size_t nbodies)
{
	float *x = bodies.x, *y = bodies.y, *z = bodies.z;
	float *vx = bodies.vx, *vy = bodies.vy, *vz = bodies.vz;
	float *mass = bodies.mass;
	float rdt = 1/dt;

for (size_t i = 0; i < nbodies; ++i) {
		float xi = x[i], yi = y[i], zi = z[i];
		float vxi = vx[i], vyi = vy[i], vzi = vz[i];
		float massi = mass[i];
		for (size_t j = i + 1; j < nbodies; ++j) {
			float dx = xi - x[j];
			float dy = yi - y[j];
			float dz = zi - z[j];
			float normsq = dx * dx + dy * dy + dz * dz;
			float dist = sqrtf(normsq); 
			float mag = rdt * normsq * dist;
			float mi = massi * mag;
			float mj = mass[j] * mag;
			vxi -= dx * mj;
			vyi -= dy * mj;
			vzi -= dz * mj;
			vx[j] += dx * mi;
			vy[j] += dy * mi;
			vz[j] += dz * mi;
		}
		vx[i] = vxi;
		vy[i] = vyi;
		vz[i] = vzi;
	}
	for (size_t i = 0; i < nbodies; ++i) {
		x[i] += dt * vx[i];
		y[i] += dt * vy[i];
		z[i] += dt * vz[i];
	}
}
#endif

#if (N & (1<<2)) || (N & (1<<5))
typedef struct {
	float x, y, z;
} float3;

static float3  f3_add(float3 a, float3 b) { return (float3){a.x+b.x, a.y+b.y, a.z+b.z}; }
static float3  f3_sub(float3 a, float3 b) { return (float3){a.x-b.x, a.y-b.y, a.z-b.z}; }
static float3  f3_mul(float3 a, float3 b) { return (float3){a.x*b.x, a.y*b.y, a.z*b.z}; }
static float3  f3_muls(float3 a, float b) { return (float3){a.x*b, a.y*b, a.z*b}; }
static float f3_sum(float3 a) { return a.x + a.y + a.z; }
#endif
#if (N) & (1<<2)

struct f3planet {
	float3 pos, vel;
	float mass;
};

void
advance_float3(struct f3planet *bodies, float dt, size_t nbodies)
{
	float rdt = 1/dt;
	for (size_t i = 0; i < nbodies; ++i) {
		struct f3planet *b = &bodies[i];
		float3 pos1 = b->pos;
		float3 vel1 = b->vel;
		float mass1 = b->mass;
		for (size_t j = i + 1; j < nbodies; ++j) {
			struct f3planet *b2 = &bodies[j];
			float3 diff = f3_sub(pos1, b2->pos);
			float normsq = f3_sum(f3_mul(diff, diff));
			float dist = sqrtf(normsq);
			float mag = rdt * normsq * dist;
			vel1 = f3_sub(vel1, f3_muls(diff, mag * b2->mass));
			b2->vel = f3_add(b2->vel, f3_muls(diff, mag * mass1));
		}
		b->vel = vel1;
	}
	for (size_t i = 0; i < nbodies; ++i) {
		struct f3planet * b = &bodies[i];
		b->pos = f3_add(b->pos, f3_muls(b->vel, dt));
	}
}
#endif

#if (N) & (1<<3)

typedef struct {
	float x, y, z, w;
} float4;

static float4  f4_add(float4 a, float4 b) { return (float4){a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w}; }
static float4  f4_sub(float4 a, float4 b) { return (float4){a.x-b.x, a.y-b.y, a.z-b.z, a.w-b.w}; }
static float4  f4_mul(float4 a, float4 b) { return (float4){a.x*b.x, a.y*b.y, a.z*b.z, a.w*b.w}; }
static float4  f4_muls(float4 a, float b) { return (float4){a.x*b, a.y*b, a.z*b, a.w*b}; }
static float f4_sum3(float4 a) { return a.x + a.y + a.z; }

struct f4planet {
	float4 pos, vel;
	float mass;
};

void
advance_float4(struct f4planet *bodies, float dt, size_t nbodies)
{
	float rdt = 1/dt;
	for (size_t i = 0; i < nbodies; ++i) {
		struct f4planet *b = &bodies[i];
		float4 pos1 = b->pos;
		float4 vel1 = b->vel;
		float mass1 = b->mass;
		for (size_t j = i + 1; j < nbodies; ++j) {
			struct f4planet *b2 = &bodies[j];
			float4 diff = f4_sub(pos1, b2->pos);
			float normsq = f4_sum3(f4_mul(diff, diff));
			float dist = sqrtf(normsq);
			float mag = rdt * normsq * dist;
			vel1 = f4_sub(vel1, f4_muls(diff, mag * b2->mass));
			b2->vel = f4_add(b2->vel, f4_muls(diff, mag * mass1));
		}
		b->vel = vel1;
	}
	for (size_t i = 0; i < nbodies; ++i) {
		struct f4planet * b = &bodies[i];
		b->pos = f4_add(b->pos, f4_muls(b->vel, dt));
	}
}
#endif

#if (N) & (1<<4)

typedef float vec4 __attribute__ ((vector_size(16)));

static float v4_sum(vec4 v) { return v[0] + v[1] + v[2]; }

struct v4planet {
	vec4 pos, vel;
	float mass;
};

void
advance_vec4(struct v4planet *bodies, float dt, size_t nbodies)
{
	float rdt = 1/dt;
	for (size_t i = 0; i < nbodies; ++i) {
		struct v4planet *b = &bodies[i];
		vec4 pos1 = b->pos;
		vec4 vel1 = b->vel;
		float mass1 = b->mass;
		for (size_t j = i + 1; j < nbodies; ++j) {
			struct v4planet *b2 = &bodies[j];
			vec4 diff = pos1 - b2->pos;
			vec4 sq = diff*diff;
			float normsq = v4_sum(sq);
			float dist = sqrtf(normsq);
			float mag = rdt * normsq * dist;
			vel1    -= diff * (mag * b2->mass);
			b2->vel += diff * (mag * mass1);
		}
		b->vel = vel1;
	}
	for (size_t i = 0; i < nbodies; ++i) {
		struct v4planet * b = &bodies[i];
		b->pos += b->vel * dt;
	}
}
#endif

#if (N) & (1<<5)

typedef struct {
	float *x, *y, *z;
} floats3;

struct fs3planet {
	floats3 pos, vel;
	float *mass;
};
static float3 fs3(floats3 a, size_t i) { return (float3){a.x[i], a.y[i], a.z[i]}; }
static void   fs3_store(floats3 a, size_t i, float3 b) { a.x[i] = b.x; a.y[i] = b.y; a.z[i] = b.z; }

void
advance_floats3(struct fs3planet b, float dt, size_t nbodies)
{
	float rdt = 1/dt;
	for (size_t i = 0; i < nbodies; ++i) {
		float3 pos1 = fs3(b.pos, i);
		float3 vel1 = fs3(b.vel, i);
		float mass1 = b.mass[i];
		for (size_t j = i + 1; j < nbodies; ++j) {
			float3 pos = fs3(b.pos, j);
			float3 diff = f3_sub(pos1, pos);
			float normsq = f3_sum(f3_mul(diff, diff));
			float dist = sqrtf(normsq);
			float mag = rdt * normsq * dist;
			vel1 = f3_sub(vel1, f3_muls(diff, mag * b.mass[j]));
			fs3_store(b.vel, j, f3_add(fs3(b.vel, j), f3_muls(diff, mag * mass1)));
		}
		fs3_store(b.vel, i, vel1);
	}
	for (size_t i = 0; i < nbodies; ++i) {
		fs3_store(b.pos, i, f3_add(fs3(b.pos, i), f3_muls(fs3(b.vel, i), dt)));
	}
}
#endif