Calculate Residuals and extra info for Arcosphere Balancing

arcospheres.c

#include <stdbool.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#define DIM 8
struct vec {
	int32_t x[DIM];
};

struct vec a = { .x = {0, -1, 0, 1, 0, 1, 0, -1} };
struct vec b = { .x = {1, -1, 0, 0, 1, -1, 0, 0} };
struct vec c = { .x = {0, 1, -1, -1, 1, 0, 0, 0} };
struct vec d = { .x = {0, 0, 0, -1, -1, 1, 1, 0} };
struct vec e = { .x = {-1, 0, 0, 0, 0, -1, 1, 1} };
struct vec f = { .x = {1, 0, 1, 0, -1, 0, -1, 0} };
struct vec g = { .x = {0, 0, -1, 1, 0, 0, -1, 1} };
struct vec h = { .x = {-1, 1, 1, 0, 0, 0, 0, -1} };
struct vec i = { .x = {-1, -1, 1, -1, 1, 1, -1, 1} };

struct vec u[7] = { 0 };
struct vec zero = {0};
int32_t dots[7] = {0};

static void vec_add(struct vec *d, int32_t c1, struct vec *s1, int32_t c2,
	struct vec *s2)
{
	for (size_t i = 0; i < DIM; i++)
		d->x[i] = c1*s1->x[i] + c2*s2->x[i];
}

static int32_t l1(struct vec *s) {
	return abs(s->x[0]) + abs(s->x[1]) + abs(s->x[2]) + abs(s->x[3])
		+ abs(s->x[4]) + abs(s->x[5]) + abs(s->x[6]) + abs(s->x[7]);
}

static int32_t dot(struct vec *a, struct vec *b) {
	int32_t res = 0;

	for (size_t i = 0; i < DIM; ++i)
		res += a->x[i]*b->x[i];

	return res;
}

static void residual_floor(struct vec *r, struct vec *src) {
	vec_add(r, 0, r, 1, src);

	for (size_t k = 0; k < 7; ++k) {
		int32_t T = dot(src, &u[k]);
		vec_add(r, 1, r, -T/dots[k], &u[k]);
	}
}

static void residual_round(struct vec *r, struct vec *src) {
	vec_add(r, 0, r, 1, src);

	for (size_t k = 0; k < 7; ++k) {
		int32_t T = dot(src, &u[k]);
		if ((T < 0) == (dots[k] < 0)) {
			T = T + dots[k]/2;
		}
		else {
			T = T - dots[k]/2;
		}
		vec_add(r, 1, r, -T/dots[k], &u[k]);
	}
}

static void print_vec(struct vec *v) {
	printf("[%d, %d, %d, %d, %d, %d, %d, %d]", v->x[0], v->x[1], v->x[2], v->x[3],
		v->x[4], v->x[5], v->x[6], v->x[7]);
}

static bool in_residual_region(struct vec *vec) {
	for (size_t k = 0; k < 7; ++k) {
		int32_t d = dot(&u[k], vec);
		if (abs(d) >= abs(dots[k])) {
			return false;
		}
	}

	return true;
}

static void (*search_test)(struct vec *s);

/*
 * Change this between residual_floor and residual_round to compare the difference
 * between the two rounding methods on the max error L1 and max error per arcosphere
 */
static void (*get_residual)(struct vec *r, struct vec *src) = residual_round;

size_t tried = 0;
void iterate_volumes(size_t coord) {
	static struct vec s = {0};

	if (coord == 7) {
		int32_t sum = s.x[0] + s.x[1] + s.x[2] + s.x[3] + s.x[4] + s.x[5] + s.x[6];
		s.x[7] = -sum;
		tried += 1;
		search_test(&s);
	}
	else {
		for (int32_t k = -4; k < 5; ++k) {
			s.x[coord] = k;
			iterate_volumes(coord+1);
		}
	}
}

static int32_t max_l1 = 0;
static struct vec max_vec = {0};
static size_t outside = 0;
void max_l1_test(struct vec *s) {
	struct vec r;

	if (!in_residual_region(s)) {
		outside += 1;
		return;
	}

	get_residual(&r, s);

	int32_t sum = l1(&r);

	if (sum > max_l1) {
		max_l1 = sum;
		memcpy(&max_vec, &r, sizeof(r));
	}
}

static size_t has_max_n = 0;
void equals_max_l1(struct vec *s) {
	struct vec r;

	if (!in_residual_region(s))
		return;

	get_residual(&r, s);
	int32_t sum = l1(&r);

	if (sum == max_l1) {
		print_vec(&r);
		printf("\n");
		has_max_n += 1;
	}
}

static int32_t max_coord_val = 0;
static size_t max_coord_coord = 0;
static struct vec max_coord_vec = {0};
void max_coord_test(struct vec *s) {
	struct vec r;

	if (!in_residual_region(s))
		return;

	get_residual(&r, s);

	if (abs(r.x[max_coord_coord]) > max_coord_val) {
		max_coord_val = abs(r.x[max_coord_coord]);
		memcpy(&max_coord_vec, &r, sizeof(r));
	}
}

int main(void) {
	size_t k;

	vec_add(&u[0], 1, &zero, 1, &c);
	vec_add(&u[1], 1, &zero, 1, &e);

	vec_add(&u[2], 1, &c, 1, &e);
	vec_add(&u[2], 1, &u[2], 1, &a);
	vec_add(&u[2], 1, &u[2], 1, &a);

	vec_add(&u[3], 1, &zero, 1, &b);
	vec_add(&u[4], 1, &zero, 1, &g);
	vec_add(&u[5], 1, &b, 1, &g);
	vec_add(&u[5], 1, &u[5], 1, &d);
	vec_add(&u[5], 1, &u[5], 1, &d);

	vec_add(&u[6], 1, &zero, 1, &i);

	for (k = 0; k < 7; ++k) {
		dots[k] = dot(&u[k], &u[k]);
	}

	for (k = 0; k < 7; ++k) {
		printf("%zu: ", k);
		print_vec(&u[k]);
		printf(", len = %d\n", dots[k]);
	}

	search_test = max_l1_test;
	iterate_volumes(0);
	printf("Examined %zu vectors\n", tried);
	printf("%zu vectors were outside the residual region\n", outside);
	printf("-> %zu total points inside the residual region\n", tried - outside);
	printf("Max L1: %d at ", max_l1);
	print_vec(&max_vec);
	printf("\n");

	search_test = equals_max_l1;
	iterate_volumes(0);

	printf("%zu vectors achieve the max L1\n", has_max_n);

	search_test = max_coord_test;
	for (k = 0; k < DIM; ++k) {
		max_coord_val = 0;
		max_coord_coord = k;
		iterate_volumes(0);
		printf("Max imbalance of %d at coord %zu, ex ", max_coord_val, max_coord_coord);
		print_vec(&max_coord_vec);
		printf("\n");
	}

	return 0;
}