Advertisement

How to use unordered_map in C++

Started by February 07, 2025 12:40 AM
12 comments, last by RomanGenkhel 3 days, 23 hours ago

I'm trying to compile the following code related to unordered_map, but I get the error:

Error C2280 'std::_Uhash_compare<_Kty,_Hasher,_Keyeq>::_Uhash_compare(const std::_Uhash_compare<_Kty,_Hasher,_Keyeq> &)': attempting to reference a deleted function
#include <Eigen/Dense>
using namespace Eigen;

#include <iostream>
#include <vector>
#include <numeric>
#include <string>
#include <sstream>
#include <algorithm>
#include <array>
#include <map>
#include <unordered_map>

using namespace std;


// https://claude.ai/chat/3caf4077-28b5-497f-b704-1b0c336a104d
// https://codereview.stackexchange.com/questions/295232/calculating-the-determinant-of-a-matrix-using-a-purely-analytical-method-that-in


template<class T, size_t N>
class Vector_nD
{
public:
	array<T, N> components;

	T operator[](size_t index) const
	{
		return components[index];
	}

	Vector_nD<T, N>& operator=(const Vector_nD<T, N>& other)
	{
		for(size_t i = 0; i < N; i++)
		components[i] = other.components[i];
	}

	Vector_nD<T, N>(const array<T, N>& comps) : components(comps)
	{

	}

	Vector_nD<T, N>(const Vector_nD<T, N>& other) : components(other.components)
	{

	}

	Vector_nD<T, N>(void)
	{
		components.fill(0.0);
	}

	static signed char permutation_sign(const array<int, (N - 1)> perm)
	{
		static unordered_map<array<int, (N - 1)>, size_t> cache;

		if (cache.end() != cache.find(perm))
		{
			size_t inversion_count = cache[perm];
			
			if (inversion_count % 2 == 0)
				return 1;
			else
				return -1;
		}

		bool sign = true;
		size_t swap_count = 0;

		for (size_t i = 0; i < (N - 2); i++)
		{
			for (size_t j = i + 1; j < (N - 1); j++)
			{
				if (perm[i] > perm[j])
				{
					swap_count++;
					// a swap would be needed here if sorting
					sign = !sign;
				}
			}
		}

		cache[perm] = swap_count;

		if (swap_count % 2 == 0)
			return 1;
		else
			return -1;
	}



	// Hodge star operator, where k = n - 1
	static Vector_nD star(const vector<Vector_nD<T, N>>& vectors)
	{
		if (vectors.size() != (N - 1))
		{
			cout << "nD operation requires (n - 1) input vectors" << endl;
			return Vector_nD<T, N>();
		}

		array<T, N> result;

		for (size_t i = 0; i < N; i++)
			result[i] = 0.0;

		// These are the indices we'll use for each component calculation
		array<int, (N - 1)> base_indices;

		for (int i = 0; i < (N - 1); i++)
			base_indices[i] = i;

		// Skip k in our calculations - this is equivalent to removing the k-th column
		// For each permutation of the remaining (N - 1) indices
		for (int k = 0; k < N; k++)
		{
			do
			{
				signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);
				//signed char sign = lehmerPermutationSign(base_indices);

				// Calculate the product for this permutation
				T product = 1.0;
				ostringstream product_oss;

				for (int i = 0; i < (N - 1); i++)
				{
					const int col = base_indices[i];

					// Adjust column index if it's past k
					int actual_col = 0;

					if (col < k)
						actual_col = col;
					else
						actual_col = col + 1;

					product_oss << "v_{" << i << actual_col << "} ";

					product *= vectors[i][actual_col];

				//	sign2 = permutation_sign(k, i);
				}


			//	cout << "Signs: " << (int)sign << " " << (int)sign2 << endl;

				//if (sign == 1)
				//	cout << "x_{" << k << "} += " << product_oss.str() << endl;
				//else
				//	cout << "x_{" << k << "} -= " << product_oss.str() << endl;

				result[k] += sign * product;

			} while(next_permutation(
					base_indices.begin(), 
					base_indices.end()));

			// (N - 1)!
			//cout << x << endl;
		}

		// Flip handedness
		for (size_t k = 0; k < N; k++)
			if (k % 2 == 1)
				result[k] = -result[k];

		cout << endl;

		for (int k = 0; k < N; k++)
			cout << "result[" << k << "] = " << result[k] << endl;

		cout << endl;

		if (N == 3)
		{
			// Demonstrate the traditional cross product too
			double x = vectors[0][0];
			double y = vectors[0][1];
			double z = vectors[0][2];

			double rhs_x = vectors[1][0];
			double rhs_y = vectors[1][1];
			double rhs_z = vectors[1][2];

			double cross_x = y * rhs_z - rhs_y * z;
			double cross_y = z * rhs_x - rhs_z * x;
			double cross_z = x * rhs_y - rhs_x * y;

			cout << cross_x << " " << cross_y << " " << cross_z << endl << endl;
		}

		return Vector_nD(result);
	}

	static T dot_product(const Vector_nD<T, N>& a, const Vector_nD<T, N>& b)
	{
		return inner_product(
			a.components.begin(), 
			a.components.end(), 
			b.components.begin(), 0.0);
	}
};


template <class T, typename size_t N>
T determinant_nxn(const MatrixX<T>& m)
{
	if (m.cols() != m.rows())
	{
		cout << "Matrix must be square" << endl;
		return 0;
	}

	// We will use this N-vector later, in the dot product operation
	Vector_nD<T, N> a_vector;

	for (size_t i = 0; i < N; i++)
		a_vector.components[i] = m(0, i);

	// We will use these (N - 1) N-vectors later, in the Hodge star operation
	vector<Vector_nD<T, N>> input_vectors;

	for (size_t i = 1; i < N; i++)
	{
		Vector_nD<T, N> b_vector;

		for (size_t j = 0; j < N; j++)
			b_vector.components[j] = m(i, j);

		input_vectors.push_back(b_vector);
	}

	// Compute the Hodge star operation using (N - 1) N-vectors
	Vector_nD<T, N> result = Vector_nD<T, N>::star(input_vectors);

	// Compute the dot product
	T det = Vector_nD<T, N>::dot_product(a_vector, result);

	// These numbers should match
	cout << "Determinant:       " << det << endl;
	cout << "Eigen Determinant: " << m.determinant() << endl << endl;

	return det;
}


template<class T, size_t N>
class std::hash<Vector_nD<T, N>>
{
public:
	size_t operator()(const Vector_nD<T, N> &key) const
	{
		hash<T> hash_fn;

		size_t r = hash_fn(key.components[0]);

		for (size_t i = 1; i < N; i++)
			r ^= hash_fn(key.components[i]);

		return r;
	}
};


int main(int argc, char** argv)
{
	srand(static_cast<unsigned int>(time(0)));

	const size_t N = 9; // Anything larger than 12 takes eons to solve for

	MatrixX<double> m(N, N);

	for (size_t i = 0; i < N; i++)
	{
		for (size_t j = 0; j < N; j++)
		{
			m(i, j) = rand() / static_cast<double>(RAND_MAX);

			if (rand() % 2 == 0)
				m(i, j) = -m(i, j);
		}
	}

	determinant_nxn<double, N>(m);

	return 0;
}

for std::unordered_map<T, P>

If i not mistaken:
‘T’ used to calculate hash (some comparable key of data, for example hash of a person is his name, yes several people can have same name, but it`s help to distinguish people)
so
std::hash<int> a; // OK, ‘int’ could be hashed
std::hash<std::array<int,5>> b; // ERROR no hash function for array, compiler don`t know hot to calculate hash of an array

None

Advertisement

Seems to me your Vector_nD is missing a operator==. unordered_map requires both an hash, but also the elements to be equality comparable, to properly handle hash-collisions.

template<class T, size_t N>
class Vector_nD
{
public:
	array<T, N> components;
	
	bool operator==(const Vector_nD& other) const = default; // implement yourself if you don't have access to C++20
}

Thank you all. I finally got it to compile, and now that I have it compiling, it is slower than not using a cache lol

#include <Eigen/Dense>
using namespace Eigen;

#include <iostream>
#include <vector>
#include <numeric>
#include <string>
#include <sstream>
#include <algorithm>
#include <array>
#include <map>
#include <unordered_map>
#include <chrono>

using namespace std;


// https://claude.ai/chat/3caf4077-28b5-497f-b704-1b0c336a104d
// https://codereview.stackexchange.com/questions/295232/calculating-the-determinant-of-a-matrix-using-a-purely-analytical-method-that-in


template<class T, size_t N>
class Vector_nD;

template<class T, size_t N>
class hash_
{
public:
	size_t operator()(const array<int, (N - 1)>& key) const
	{
		hash<T> hash_fn;

		size_t r = hash_fn(key[0]);

		for (size_t i = 1; i < N - 1; i++)
			r ^= hash_fn(key[i]);

		return r;
	} 
};

template<class T, size_t N>
class Vector_nD
{
public:
	array<T, N> components;

	T operator[](size_t index) const
	{
		return components[index];
	}

	Vector_nD<T, N>& operator=(const Vector_nD<T, N>& other)
	{
		for(size_t i = 0; i < N; i++)
		components[i] = other.components[i];

		return *this;
	}

	bool operator==(const Vector_nD<T, N>& other) const
	{
		for (size_t i = 0; i < N; i++)
		{
			if (components[i] != other.components[i])
				return false;
		}

		return true;
	}


	Vector_nD<T, N>(const array<T, N>& comps) : components(comps)
	{

	}

	Vector_nD<T, N>(const Vector_nD<T, N>& other) : components(other.components)
	{

	}

	Vector_nD<T, N>(void)
	{
		components.fill(0.0);
	}

	static signed char permutation_sign(const array<int, (N - 1)> perm)
	{
		static map<array<int, (N - 1)>, size_t> cache;



		if (0)//cache.end() != cache.find(perm))
		{
			size_t inversion_count = cache[perm];
			
			if (inversion_count % 2 == 0)
				return 1;
			else
				return -1;
		}

		bool sign = true;
		size_t swap_count = 0;

		for (size_t i = 0; i < (N - 2); i++)
		{
			for (size_t j = i + 1; j < (N - 1); j++)
			{
				if (perm[i] > perm[j])
				{
					swap_count++;
					// a swap would be needed here if sorting
					sign = !sign;
				}
			}
		}

		//cache[perm] = swap_count;

		if (swap_count % 2 == 0)
			return 1;
		else
			return -1;
	}



	// Hodge star operator, where k = n - 1
	static Vector_nD star(const vector<Vector_nD<T, N>>& vectors)
	{
		if (vectors.size() != (N - 1))
		{
			cout << "nD operation requires (n - 1) input vectors" << endl;
			return Vector_nD<T, N>();
		}

		array<T, N> result;

		for (size_t i = 0; i < N; i++)
			result[i] = 0.0;

		// These are the indices we'll use for each component calculation
		array<int, (N - 1)> base_indices;

		for (int i = 0; i < (N - 1); i++)
			base_indices[i] = i;

		// Skip k in our calculations - this is equivalent to removing the k-th column
		// For each permutation of the remaining (N - 1) indices
		for (int k = 0; k < N; k++)
		{
			do
			{
				signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);
				//signed char sign = lehmerPermutationSign(base_indices);

				// Calculate the product for this permutation
				T product = 1.0;
				ostringstream product_oss;

				for (int i = 0; i < (N - 1); i++)
				{
					const int col = base_indices[i];

					// Adjust column index if it's past k
					int actual_col = 0;

					if (col < k)
						actual_col = col;
					else
						actual_col = col + 1;

					product_oss << "v_{" << i << actual_col << "} ";

					product *= vectors[i][actual_col];

				//	sign2 = permutation_sign(k, i);
				}


			//	cout << "Signs: " << (int)sign << " " << (int)sign2 << endl;

				//if (sign == 1)
				//	cout << "x_{" << k << "} += " << product_oss.str() << endl;
				//else
				//	cout << "x_{" << k << "} -= " << product_oss.str() << endl;

				result[k] += sign * product;

			} while(next_permutation(
					base_indices.begin(), 
					base_indices.end()));

			// (N - 1)!
			//cout << x << endl;
		}

		// Flip handedness
		for (size_t k = 0; k < N; k++)
			if (k % 2 == 1)
				result[k] = -result[k];

		cout << endl;

		for (int k = 0; k < N; k++)
			cout << "result[" << k << "] = " << result[k] << endl;

		cout << endl;

		if (N == 3)
		{
			// Demonstrate the traditional cross product too
			double x = vectors[0][0];
			double y = vectors[0][1];
			double z = vectors[0][2];

			double rhs_x = vectors[1][0];
			double rhs_y = vectors[1][1];
			double rhs_z = vectors[1][2];

			double cross_x = y * rhs_z - rhs_y * z;
			double cross_y = z * rhs_x - rhs_z * x;
			double cross_z = x * rhs_y - rhs_x * y;

			cout << cross_x << " " << cross_y << " " << cross_z << endl << endl;
		}

		return Vector_nD(result);
	}

	static T dot_product(const Vector_nD<T, N>& a, const Vector_nD<T, N>& b)
	{
		return inner_product(
			a.components.begin(), 
			a.components.end(), 
			b.components.begin(), 0.0);
	}
};


template <class T, typename size_t N>
T determinant_nxn(const MatrixX<T>& m)
{
	if (m.cols() != m.rows())
	{
		cout << "Matrix must be square" << endl;
		return 0;
	}

	// We will use this N-vector later, in the dot product operation
	Vector_nD<T, N> a_vector;

	for (size_t i = 0; i < N; i++)
		a_vector.components[i] = m(0, i);

	// We will use these (N - 1) N-vectors later, in the Hodge star operation
	vector<Vector_nD<T, N>> input_vectors;

	for (size_t i = 1; i < N; i++)
	{
		Vector_nD<T, N> b_vector;

		for (size_t j = 0; j < N; j++)
			b_vector.components[j] = m(i, j);

		input_vectors.push_back(b_vector);
	}

	// Compute the Hodge star operation using (N - 1) N-vectors
	Vector_nD<T, N> result = Vector_nD<T, N>::star(input_vectors);

	// Compute the dot product
	T det = Vector_nD<T, N>::dot_product(a_vector, result);

	// These numbers should match
	cout << "Determinant:       " << det << endl;
	cout << "Eigen Determinant: " << m.determinant() << endl << endl;

	return det;
}





int main(int argc, char** argv)
{
	srand(static_cast<unsigned int>(time(0)));

	const size_t N = 10; // Anything larger than 12 takes eons to solve for

	MatrixX<double> m(N, N);

	for (size_t i = 0; i < N; i++)
	{
		for (size_t j = 0; j < N; j++)
		{
			m(i, j) = rand() / static_cast<double>(RAND_MAX);

			if (rand() % 2 == 0)
				m(i, j) = -m(i, j);
		}
	}

	std::chrono::high_resolution_clock::time_point start_time, end_time;
	start_time = std::chrono::high_resolution_clock::now();

	determinant_nxn<double, N>(m);

	end_time = std::chrono::high_resolution_clock::now();
	std::chrono::duration<float, std::milli> elapsed = end_time - start_time;

	cout << "Duration: " << elapsed.count() / 1000.0f << " seconds";

	return 0;
}

Which configuration are you compiling for? I'm assuming Release/optimize, but if for any reason you are doing a debug/unoptimized build, then this could indeed be much slower, since your usage of the cache is suboptimal.

static unordered_map<array<int, (N - 1)>, size_t> cache;

		if (cache.end() != cache.find(perm))
		{
			size_t inversion_count = cache[perm];

What you are writing here is “check if perm is contained in the cache. Then, fetch perm from the cache.”. In unoptimized code or a bad compiler, this will result in two accesses into the container. What you should do is:

if (const auto itr = cache.find(perm); itr != cache.end())
	size_t inversion_count = itr->second;

Similar thing when writing to the cache:

cache[perm] = swap_count;

This is a more minor point, but [perm] actually default-constructs an element, returns the reference and then sets the value. What you do want is

cache.try_emplace(perm, swap_count); 

which constructs swap_count directly.

A decent compiler on release will mostly optimize the difference away, though even clang seems not to be able to fold the [] to a try_emplace (though this is minor). In a debug-build, especially on MSVC with iterator_debug level, the difference will be stark.

Otherwise, it can still make sense why not using a cache is faster. Your main loop here is O(N^2), but for a small N, using the form to /2 the result (potentially even more as it seems). map/unordered_map have a comparatively high overhead compared to just doing brute force on a small N, they are not really cache-friendly. They only save time if the single operation you are trying to do is actually pretty expensive. Having to hash the array also makes it at least O(N) for that part, potentially multiple times, depending on if the STL internally caches the hash in each bucket-entry… Plus a compare every time to make sure that both elements are the same. So, even the few potential inefficient uses of the map aside, I'm not really suprised.

taby said:
nd now that I have it compiling, it is slower than not using a cache lol

Not sure, but looks you could use something else than a map, achieving constant time lookup cost. A matrix of sign bits maybe?

Advertisement

Addition: Since as I said the find - try_emplace pair still uses two lookups into the container unfortunately, one could try to reduce it further to one lookup:

static unordered_map<array<int, (N - 1)>, size_t> cache;

auto [itr, wasAdded] = cache.try_emplace(perm, -1);
if (wasAdded) // placeholder -1 was written, meaning we do have to do the calc
{
   // calculations ...
   
   itr->second = swap_count;
}

This will reduce the second lookup, in any configurations, on all compilers. At least then you'd be sure that, if its slower, the overhead from using the cache does not stem from having to do 2 separate lookup/calculations.

I got rid of the find() call:

#include <Eigen/Dense>
using namespace Eigen;

#include <iostream>
#include <vector>
#include <numeric>
#include <string>
#include <sstream>
#include <algorithm>
#include <array>
#include <map>
#include <unordered_map>
#include <chrono>

using namespace std;


// https://claude.ai/chat/3caf4077-28b5-497f-b704-1b0c336a104d
// https://codereview.stackexchange.com/questions/295232/calculating-the-determinant-of-a-matrix-using-a-purely-analytical-method-that-in


template<class T, size_t N>
class Vector_nD;

template<class T, size_t N>
class hash_
{
public:
	size_t operator()(const array<int, (N - 1)>& key) const
	{
		hash<T> hash_fn;

		size_t r = hash_fn(key[0]);

		for (size_t i = 1; i < N - 1; i++)
			r ^= hash_fn(key[i]);

		return r;
	} 
};

template<class T, size_t N>
class Vector_nD
{
public:
	array<T, N> components;

	T operator[](size_t index) const
	{
		return components[index];
	}

	Vector_nD<T, N>& operator=(const Vector_nD<T, N>& other)
	{
		for(size_t i = 0; i < N; i++)
		components[i] = other.components[i];

		return *this;
	}

	bool operator==(const Vector_nD<T, N>& other) const
	{
		for (size_t i = 0; i < N; i++)
		{
			if (components[i] != other.components[i])
				return false;
		}

		return true;
	}

	Vector_nD<T, N>(const array<T, N>& comps) : components(comps)
	{

	}

	Vector_nD<T, N>(const Vector_nD<T, N>& other) : components(other.components)
	{

	}

	Vector_nD<T, N>(void)
	{
		components.fill(0.0);
	}
	
	static signed char permutation_sign(const array<int, (N - 1)> perm)
	{
		bool sign = true;

		for (size_t i = 0; i < (N - 2); i++)
		{
			for (size_t j = i + 1; j < (N - 1); j++)
			{
				if (perm[i] > perm[j])
				{
					// a swap would be needed here if sorting
					sign = !sign;
				}
			}
		}

		if (sign)
			return 1;
		else
			return -1;
	}

	// Hodge star operator, where k = n - 1
	static Vector_nD star(const vector<Vector_nD<T, N>>& vectors)
	{
		if (vectors.size() != (N - 1))
		{
			cout << "nD operation requires (n - 1) input vectors" << endl;
			return Vector_nD<T, N>();
		}

		array<T, N> result;

		for (size_t i = 0; i < N; i++)
			result[i] = 0.0;

		// These are the indices we'll use for each component calculation
		array<int, (N - 1)> base_indices;

		for (int i = 0; i < (N - 1); i++)
			base_indices[i] = i;

		//unordered_map<array<int, (N - 1)>, signed char, hash_<T, N>> cache;
		map<array<int, (N - 1)>, signed char> cache;

		// Fill cache
		do 
		{
			signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);
			cache[base_indices] = sign;

		} while (next_permutation(
			base_indices.begin(),
			base_indices.end()));


		// Skip k in our calculations - this is equivalent to removing the k-th column
		// For each permutation of the remaining (N - 1) indices
		for (int k = 0; k < N; k++)
		{
			do
			{
				signed char sign = cache[base_indices];
				//signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);

				// Calculate the product for this permutation
				T product = 1.0;
				ostringstream product_oss;

				for (int i = 0; i < (N - 1); i++)
				{
					const int col = base_indices[i];

					// Adjust column index if it's past k
					int actual_col = 0;

					if (col < k)
						actual_col = col;
					else
						actual_col = col + 1;

					product_oss << "v_{" << i << actual_col << "} ";

					product *= vectors[i][actual_col];
				}

				//if (sign == 1)
				//	cout << "x_{" << k << "} += " << product_oss.str() << endl;
				//else
				//	cout << "x_{" << k << "} -= " << product_oss.str() << endl;

				result[k] += sign * product;

			} while(next_permutation(
					base_indices.begin(), 
					base_indices.end()));

			// (N - 1)!
			//cout << x << endl;
		}

		// Flip handedness
		for (size_t k = 0; k < N; k++)
			if (k % 2 == 1)
				result[k] = -result[k];

		cout << endl;

		for (int k = 0; k < N; k++)
			cout << "result[" << k << "] = " << result[k] << endl;

		cout << endl;

		if (N == 3)
		{
			// Demonstrate the traditional cross product too
			double x = vectors[0][0];
			double y = vectors[0][1];
			double z = vectors[0][2];

			double rhs_x = vectors[1][0];
			double rhs_y = vectors[1][1];
			double rhs_z = vectors[1][2];

			double cross_x = y * rhs_z - rhs_y * z;
			double cross_y = z * rhs_x - rhs_z * x;
			double cross_z = x * rhs_y - rhs_x * y;

			cout << cross_x << " " << cross_y << " " << cross_z << endl << endl;
		}

		return Vector_nD(result);
	}

	static T dot_product(const Vector_nD<T, N>& a, const Vector_nD<T, N>& b)
	{
		return inner_product(
			a.components.begin(), 
			a.components.end(), 
			b.components.begin(), 0.0);
	}
};


template <class T, typename size_t N>
T determinant_nxn(const MatrixX<T>& m)
{
	if (m.cols() != m.rows())
	{
		cout << "Matrix must be square" << endl;
		return 0;
	}

	// We will use this N-vector later, in the dot product operation
	Vector_nD<T, N> a_vector;

	for (size_t i = 0; i < N; i++)
		a_vector.components[i] = m(0, i);

	// We will use these (N - 1) N-vectors later, in the Hodge star operation
	vector<Vector_nD<T, N>> input_vectors;

	for (size_t i = 1; i < N; i++)
	{
		Vector_nD<T, N> b_vector;

		for (size_t j = 0; j < N; j++)
			b_vector.components[j] = m(i, j);

		input_vectors.push_back(b_vector);
	}

	// Compute the Hodge star operation using (N - 1) N-vectors
	Vector_nD<T, N> result = Vector_nD<T, N>::star(input_vectors);

	// Compute the dot product
	T det = Vector_nD<T, N>::dot_product(a_vector, result);

	// These numbers should match
	cout << "Determinant:       " << det << endl;
	cout << "Eigen Determinant: " << m.determinant() << endl << endl;

	return det;
}

int main(int argc, char** argv)
{
	srand(static_cast<unsigned int>(time(0)));

	const size_t N = 10; // Anything larger than 12 takes eons to solve for

	MatrixX<double> m(N, N);

	for (size_t i = 0; i < N; i++)
	{
		for (size_t j = 0; j < N; j++)
		{
			m(i, j) = rand() / static_cast<double>(RAND_MAX);

			if (rand() % 2 == 0)
				m(i, j) = -m(i, j);
		}
	}

	std::chrono::high_resolution_clock::time_point start_time, end_time;
	start_time = std::chrono::high_resolution_clock::now();

	determinant_nxn<double, N>(m);

	end_time = std::chrono::high_resolution_clock::now();
	std::chrono::duration<float, std::milli> elapsed = end_time - start_time;

	cout << "Duration: " << elapsed.count() / 1000.0f << " seconds";

	return 0;
}

I tried using a set instead of a map, but no big difference:

#include <Eigen/Dense>
using namespace Eigen;

#include <iostream>
#include <vector>
#include <numeric>
#include <string>
#include <sstream>
#include <algorithm>
#include <array>
#include <chrono>
#include <set>

using namespace std;


// https://claude.ai/chat/3caf4077-28b5-497f-b704-1b0c336a104d
// https://codereview.stackexchange.com/questions/295232/calculating-the-determinant-of-a-matrix-using-a-purely-analytical-method-that-in


template<class T, size_t N>
class Vector_nD
{
public:
	array<T, N> components;

	T operator[](size_t index) const
	{
		return components[index];
	}

	Vector_nD<T, N>& operator=(const Vector_nD<T, N>& other)
	{
		for(size_t i = 0; i < N; i++)
		components[i] = other.components[i];

		return *this;
	}

	bool operator==(const Vector_nD<T, N>& other) const
	{
		for (size_t i = 0; i < N; i++)
		{
			if (components[i] != other.components[i])
				return false;
		}

		return true;
	}


	Vector_nD<T, N>(const array<T, N>& comps) : components(comps)
	{

	}

	Vector_nD<T, N>(const Vector_nD<T, N>& other) : components(other.components)
	{

	}

	Vector_nD<T, N>(void)
	{
		components.fill(0.0);
	}



	static signed char permutation_sign(const array<int, (N - 1)> perm)
	{
		bool sign = true;

		for (size_t i = 0; i < (N - 2); i++)
		{
			for (size_t j = i + 1; j < (N - 1); j++)
			{
				if (perm[i] > perm[j])
				{
					// a swap would be needed here if sorting
					sign = !sign;
				}
			}
		}

		if (sign)
			return 1;
		else
			return -1;
	}



	// Hodge star operator, where k = n - 1
	static Vector_nD star(const vector<Vector_nD<T, N>>& vectors)
	{
		if (vectors.size() != (N - 1))
		{
			cout << "nD operation requires (n - 1) input vectors" << endl;
			return Vector_nD<T, N>();
		}

		array<T, N> result;

		for (size_t i = 0; i < N; i++)
			result[i] = 0.0;

		// These are the indices we'll use for each component calculation
		array<int, (N - 1)> base_indices;

		for (int i = 0; i < (N - 1); i++)
			base_indices[i] = i;

		// Fill cache
		set< array<int, (N - 1)>> cache_set_positive;

		do
		{
			signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);
			
			if (sign == 1)
				cache_set_positive.insert(base_indices);

		} while (next_permutation(
			base_indices.begin(),
			base_indices.end()));


		// Skip k in our calculations - this is equivalent to removing the k-th column
		// For each permutation of the remaining (N - 1) indices
		for (int k = 0; k < N; k++)
		{
			do
			{
				//signed char sign = Vector_nD<T, N>::permutation_sign(base_indices);

				signed char sign = 0;

				auto ci = cache_set_positive.find(base_indices);

				if (ci != cache_set_positive.end())
					sign = 1;
				else
					sign = -1;

				// Calculate the product for this permutation
				T product = 1.0;
				ostringstream product_oss;

				for (int i = 0; i < (N - 1); i++)
				{
					const int col = base_indices[i];

					// Adjust column index if it's past k
					int actual_col = 0;

					if (col < k)
						actual_col = col;
					else
						actual_col = col + 1;

					product_oss << "v_{" << i << actual_col << "} ";

					product *= vectors[i][actual_col];
				}

				//if (sign == 1)
				//	cout << "x_{" << k << "} += " << product_oss.str() << endl;
				//else
				//	cout << "x_{" << k << "} -= " << product_oss.str() << endl;

				result[k] += sign * product;

			} while(next_permutation(
					base_indices.begin(), 
					base_indices.end()));
		}

		// Flip handedness
		for (size_t k = 0; k < N; k++)
			if (k % 2 == 1)
				result[k] = -result[k];

		cout << endl;

		for (int k = 0; k < N; k++)
			cout << "result[" << k << "] = " << result[k] << endl;

		cout << endl;

		if (N == 3)
		{
			// Demonstrate the traditional cross product too
			double x = vectors[0][0];
			double y = vectors[0][1];
			double z = vectors[0][2];

			double rhs_x = vectors[1][0];
			double rhs_y = vectors[1][1];
			double rhs_z = vectors[1][2];

			double cross_x = y * rhs_z - rhs_y * z;
			double cross_y = z * rhs_x - rhs_z * x;
			double cross_z = x * rhs_y - rhs_x * y;

			cout << cross_x << " " << cross_y << " " << cross_z << endl << endl;
		}

		return Vector_nD(result);
	}

	static T dot_product(const Vector_nD<T, N>& a, const Vector_nD<T, N>& b)
	{
		return inner_product(
			a.components.begin(), 
			a.components.end(), 
			b.components.begin(), 0.0);
	}
};


template <class T, typename size_t N>
T determinant_nxn(const MatrixX<T>& m)
{
	if (m.cols() != m.rows())
	{
		cout << "Matrix must be square" << endl;
		return 0;
	}

	// We will use this N-vector later, in the dot product operation
	Vector_nD<T, N> a_vector;

	for (size_t i = 0; i < N; i++)
		a_vector.components[i] = m(0, i);

	// We will use these (N - 1) N-vectors later, in the Hodge star operation
	vector<Vector_nD<T, N>> input_vectors;

	for (size_t i = 1; i < N; i++)
	{
		Vector_nD<T, N> b_vector;

		for (size_t j = 0; j < N; j++)
			b_vector.components[j] = m(i, j);

		input_vectors.push_back(b_vector);
	}

	// Compute the Hodge star operation using (N - 1) N-vectors
	Vector_nD<T, N> result = Vector_nD<T, N>::star(input_vectors);

	// Compute the dot product
	T det = Vector_nD<T, N>::dot_product(a_vector, result);

	// These numbers should match
	cout << "Determinant:       " << det << endl;
	cout << "Eigen Determinant: " << m.determinant() << endl << endl;

	return det;
}





int main(int argc, char** argv)
{
	srand(static_cast<unsigned int>(time(0)));

	const size_t N = 8; // Anything larger than 12 takes eons to solve for

	MatrixX<double> m(N, N);

	for (size_t i = 0; i < N; i++)
	{
		for (size_t j = 0; j < N; j++)
		{
			m(i, j) = rand() / static_cast<double>(RAND_MAX);

			if (rand() % 2 == 0)
				m(i, j) = -m(i, j);
		}
	}

	std::chrono::high_resolution_clock::time_point start_time, end_time;
	start_time = std::chrono::high_resolution_clock::now();

	determinant_nxn<double, N>(m);

	end_time = std::chrono::high_resolution_clock::now();
	std::chrono::duration<float, std::milli> elapsed = end_time - start_time;

	cout << "Duration: " << elapsed.count() / 1000.0f << " seconds";

	return 0;
}

I am convinced that the Levi-Civita function (e.g. permutation sign) foils all attempts to be compressed into some simpler representation or pattern. There’s no shortcut, and if there was, it would be widely known, even by the AIs…. Right?

Advertisement