Atari Lynx Mode 7 table generator

gistfile1.txt

#include <cstdio>
#include <numbers>
#include <vector>

/*
	The prime factorisation of n can contain at most one prime factor that is greater than
	sqrt(n) as:
	
		(1) if it has any factors > sqrt(n) then either they're prime or they can
			decompose into other, smaller prime factors; and
		(2)	if more than one of those _prime_ factors were more than sqrt(n), then they
			would multiply together to make a number greater than n.

	Therefore an approach to factoring is to test all the primes up to
	sqrt(n) and take any residue as a single remaining prime.
*/
namespace {

/// One beyond the largest absolute value that either number output may have.
/// So I've opted to use the Lynx's **signed** multiplication (for reasons of
/// simplifying large values of the 32-bit total) here.
static constexpr uint32_t FactorLimit = 32768;

/// Uses the sieve of Eratosthenes to construct a list of primes that are strictly
/// less than the FactorLimit, and can subsequently vend them in descending order.
struct Primes {
	Primes() {
		std::vector<bool> not_prime(32768);

		int c = 1;
		while(++c < 32768) {
			if(not_prime[c]) continue;

			primes_.insert(primes_.begin(), c);
			int ic = c;
			while(ic < 32768) {
				not_prime[ic] = true;
				ic += c;
			}
		}
	}

	/// @returns A vector of all primes less than 65536, in descending order.
	const std::vector<int> primes() {
		return primes_;
	}

private:
	std::vector<int> primes_;
};
Primes primes;		// Construct only one of these. So it isn't worth using constexpr or traditional
					// templates to do this at compile time, as it will definitely be done only once
					// and isn't especially costly.

/// @returns The prime factorisation of @c value, in descending order.
std::vector<uint32_t> prime_factors(uint32_t value) {
	// Approach is as above: take out as many prime factors as possible as are less than 65536.
	// Then anything that's left is either '1' or a bigger prime.
	std::vector<uint32_t> result;
	for(auto prime: primes.primes()) {
		while(!(value %prime)) {
			result.push_back(prime);
			value /= prime;
		}
	}
	if(value > 1) result.insert(result. begin(), value);
	return result;
}
	
/// Given a list of prime factors, attempts to combine them into two 15-bit numbers that can be multiplied together
/// to produce the same product as the primes.
std::optional<std::pair<uint16_t, uint16_t>> factor_pair(const std::vector<uint32_t> &prime_factors) {
	uint32_t factors[2] = {1, 1};

	// Attempt 1: allocate to the current least.
	// Usually works.
	for(const auto factor: prime_factors) {
		const size_t target = factors[1] < factors[0];
		factors[target] *= factor;
		if(factors[target] >= FactorLimit) break;
	}

	if(factors[0] < FactorLimit && factors[1] < FactorLimit) {
		return std::make_pair<uint16_t>(factors[0], factors[1]);
	}

	// Failing that, try an exhaustive search.
	// Sometimes succeeds where the above failed.
	for(size_t c = 0; c < (1 << prime_factors.size()); c++) {
		factors[0] = factors[1] = 1;

		auto index = c;
		for(const auto factor: prime_factors) {
			factors[index & 1] *= factor;
			if(factors[index & 1] >= FactorLimit) break;
			index >>= 1;
		}

		if(factors[0] < FactorLimit && factors[1] < FactorLimit) {
			return std::make_pair<uint16_t>(factors[0], factors[1]);
		}
	}

	return std::nullopt;
}
}

int main(int argc, char *argv[]) {
	static constexpr int TicksPerCircle = 256;		// As per your desired number of units in a complete circle.
													// Doesn't need to be a power of two.

	// For bookkeeping only; the total error for a solution is the Manhattan distance between the interpolant
	// found and the one desired. The frequency of each value of error is counted for user information.
	int errors[256]{};

	for(int c = 0; c < TicksPerCircle; c++) {
		// Convert to an angle in radians and determine the ideal integer representation.
		const float angle = (float(c) / float(TicksPerCircle - 1)) * std::numbers::pi_v<float> * 2.0f;
		const uint16_t vector[2] = {
			uint16_t(int16_t(0.5f + sin(angle) * 256.0f)),	// i.e. round, don't truncate; use 8:8 fixed point.
			uint16_t(int16_t(0.5f + cos(angle) * 256.0f))
		};

		/// Searches for a soution that hits exactly @c target; if one is found then it will be output to
		/// the console. Its @c error is also recorded into the bookkeeping totals for final output.
		const auto attempt = [&](uint32_t target, const int error) -> bool {
			// Simple observation: larger numbers are harder to factorise both in general and within
			// the constraints of trying to boil down to two 16-bit numbers to multiply. So
			int multiplier = 1;
			if(int32_t(target) < 0) {
				target = -target;
				multiplier = -1;
			}

			const auto log = [&](const uint16_t first, const uint16_t second) {
				printf("0x%08x -> 0x%08x: ", (vector[0] << 16) | vector[1], int16_t(first * multiplier) * int16_t(second));
				printf("0x%04x * 0x%04x [%d]\n", uint16_t(first * multiplier), second, error);
				++errors[error];
			};

			if(target < FactorLimit) {
				// 0 and 1 aren't otherwise resolvable since neither has a prime decomposition
				// per my 1-isn't-prime logic. But then I might as well extend the test.
				log(1, target);
				return true;
			}
			const auto pair = factor_pair(prime_factors(target));
			if(!pair.has_value()) {
				return false;
			}

			log(pair->first, pair->second);
			return true;
		};


		// Prefer getting the exact answer.
		//
		// That won't always be possible; e.g. sometimes the true adder is itself just a large prime number.
		//
		int offset = 0;		// b0: major axis;
							// b1: direction along that axis;
							// b2: direction along orthogonal axis;
							// b3–: distance along orthogonal axis.
		int distance = 0;
		while(true) {
			// Compute point to test.
			const int axis = offset & 1;
			const int direction = offset & 2 ? -1 : 1;
			const int orthogonal_direction = offset & 4 ? -1 : 1;
			const int orthogonal_distance = offset >> 3;

			uint16_t test_vector[2] = {
				vector[0], vector[1]
			};

			// Walk directly along the major axis, per the direction.
			// From there walk along the orthogonal.

			test_vector[axis] += direction * distance;
			test_vector[axis ^ 1] += orthogonal_direction * orthogonal_distance;

			const uint32_t target = (test_vector[0] << 16) | test_vector[1];
			if(attempt(target, distance + orthogonal_distance)) break;

			++ offset;
			if((offset >> 3) > distance) {
				offset = 0;
				++distance;
			}
		}
	}


	for(int c = 0; c < 256; c++) {
		if(!errors[c]) continue;

		if(!c) printf("Exact results: %d\n", errors[c]);
		else printf("Off-by-%d/256: %d\n", c, errors[c]);
	}
	return 0;
}