Sudoku solving with dancing links algorithm.

dlx.hpp

#ifndef DLX_HPP
#define DLX_HPP

#include <vector>
#include <cassert>

namespace dlx
{
	class solver final
	{
		class node
		{
		public:

			//struct for node identification
			//can go further and make node templated, so it can hold any arbitrary data
			struct index_t
			{
				index_t(): row(0), column(0){}

				index_t(std::size_t row, std::size_t column): row(row), column(column){}

				std::size_t row;
				std::size_t column;
			};

			node(bool is_header=false): left(nullptr), right(nullptr), down(nullptr), 
				up(nullptr), is_header_flag(is_header), size_holder(0)
			{
				if(is_header)
				{
					down=this;
					up=this;
				}
			}

			//root is treated like regular header
			node* header() const
			{
				if(is_header() || is_root())
				{
					return const_cast<node*>(this);
				}
				node* current_node=this->up;
				while(!current_node->is_header())
				{
					current_node=current_node->up;
				}
				return current_node;
			}

			inline std::size_t size() const
			{
				assert(is_header());
				return size_holder;
			}

			inline bool is_header() const
			{
				return is_header_flag;
			}

			inline bool is_root() const
			{
				return !up && !down;
			}

			inline auto index() const
			{
				return index_holder;
			}

			node* left;
			node* right;
			node* down;
			node* up;
			bool is_header_flag;
			std::size_t size_holder;
			index_t index_holder;
		};

	public:

		solver(): base(nullptr), size_holder(0)
		{
			init_root();
		}

		//I don't want to consume (&&) problem, or do I?
		solver(std::vector<std::vector<unsigned int>>& problem): 
			base(nullptr), size_holder(0)
		{
			init_root();
			apply_problem(problem);
		}

		void apply_problem(std::vector<std::vector<unsigned int>>& problem)
		{
			//some sanity check needed before applying problem
			problem_sanity_check(problem);
			//place headers first
			for(auto& it: *problem.begin())
			{
				append_header();
			}
			//empty problem omegaLUL
			assert(base->right!=base);
			for(std::size_t row=0; row<problem.size(); row++)
			{
				auto current_header=base->right;
				node* last_node=nullptr;
				for(std::size_t column=0; column<problem[row].size(); column++)
				{
					if(problem[row][column]!=0)
					{
						last_node=append_node(current_header, last_node, node::index_t(row, column));
					}
					current_header=current_header->right;
				}
			}
		}

		inline std::size_t size() const
		{
			return size_holder;
		}

		void clear()
		{
			if(!base)
			{
				return;
			}
			for(auto header=base->right; header!=base;)
			{
				for(auto node=header->down; node!=header;)
				{
					auto temp=node->down;
					delete node;
					node=temp;
				}
				auto temp=header->right;
				delete header;
				header=temp;
			}
			delete base;
			base=nullptr;
			size_holder=0;
		}

		~solver()
		{
			clear();
		}

		inline std::vector<std::vector<std::size_t>> search()
		{
			std::vector<std::vector<node*>> solutions;
			std::vector<node*> solution;
			search(0, solution, solutions);
			return parse_solutions(solutions);
		}

	private:

		inline void problem_sanity_check(std::vector<std::vector<unsigned int>>& problem)
		{
			for(auto& row: problem)
			{
				assert(row.size()==problem.begin()->size());
			}
		}

		std::vector<std::vector<std::size_t>> parse_solutions(std::vector<std::vector<node*>>& solutions)
		{
			std::vector<std::vector<std::size_t>> result;
			for(auto& solution: solutions)
			{
				std::vector<std::size_t> current_result;
				for(auto& element: solution)
				{
					if(element)
					{
						current_result.push_back(element->index().row);
					}
				}
				if(!current_result.empty())
				{
					std::sort(current_result.begin(), current_result.end());
					result.push_back(std::move(current_result));
				}
			}
			return result;
		}

		void search(std::size_t step, std::vector<node*>& solution, 
			std::vector<std::vector<node*>>& solutions)
		{
			if(base->right==base)
			{
				solutions.push_back(solution);
				return;
			}
			auto column=choose_column();
			//check if only headers row left
			if(column->down==column)
			{
				return;
			}
			cover(column);
			for(auto down_node=column->down; down_node!=column; down_node=down_node->down)
			{
				solution.push_back(down_node);
				for(auto right_node=down_node->right; right_node!=down_node; right_node=right_node->right)
				{
					cover(right_node->header());
				}
				search(step+1, solution, solutions);
				auto res=solution.back();
				solution.pop_back();
				for(auto left_node=res->left; left_node!=res; left_node=left_node->left)
				{
					uncover(left_node->header());
				}
			}
			uncover(column);
			return;
		}

		void cover(node* column)
		{
			column->right->left=column->left;
			column->left->right=column->right;
			for(auto down_node=column->down; down_node!=column; down_node=down_node->down)
			{
				for(auto right_node=down_node->right; right_node!=down_node; right_node=right_node->right)
				{
					right_node->down->up=right_node->up;
					right_node->up->down=right_node->down;
					right_node->header()->size_holder--;
				}
			}
		}

		void uncover(node* column)
		{
			for(auto up_node=column->up; up_node!=column; up_node=up_node->up)
			{
				for(auto left_node=up_node->left; left_node!=up_node; left_node=left_node->left)
				{
					left_node->header()->size_holder++;
					left_node->down->up=left_node;
					left_node->up->down=left_node;
				}
			}
			column->right->left=column;
			column->left->right=column;
		}

		//This is the key to performance
		//simple base->right will be too slow
		node* choose_column() const
		{
			auto min_column=base->right;
			for(auto column=base->right; column!=base; column=column->right)
			{
				if(column->size()<min_column->size())
				{
					min_column=column;
				}
			}
			return min_column;
		}

		//The presence of root node is very important too
		//Without root you have to manually handle cases like 1x1, 1x2, 2x2 ...
		void init_root()
		{
			assert(!base);
			base=new node();
			base->left=base;
			base->right=base;
		}

		node* append_header()
		{
			auto current_node=new node(true);
			current_node->left=base->left;
			base->left->right=current_node;
			base->left=current_node;
			current_node->right=base;
			size_holder++;
			return current_node;
		}

		node* append_node(node* header, node* last_node, node::index_t index)
		{
			auto current_node=new node();
			current_node->up=header->up;
			header->up->down=current_node;
			header->up=current_node;
			current_node->down=header;
			if(last_node)
			{
				current_node->left=last_node;
				last_node->right->left=current_node;
				current_node->right=last_node->right;
				last_node->right=current_node;
			}
			else
			{
				current_node->left=current_node;
				current_node->right=current_node;
			}
			current_node->index_holder=index;
			header->size_holder++;
			return current_node;
		}

		node* base;
		std::size_t size_holder;
		std::size_t max_rows;
	};
}

#endif

main.cpp

#if defined(_MSC_VER) && defined(_DEBUG)
	#define _CRTDBG_MAP_ALLOC
	#include <crtdbg.h>
#endif

#include <iostream>
#include <vector>
#include <exception>
#include <algorithm>
#include <cassert>
#include <fstream>
#include <chrono>

#include "dlx.hpp"


auto generate_sudoku_problem()
{
	std::vector<std::vector<unsigned int>> row_column_constraints;
	for(std::size_t it=0; it<9*9; it++)
	{
		for(std::size_t st=1; st<=9; st++)
		{
			std::vector<unsigned int> current(81, 0);
			current[it]=(unsigned int)st;
			row_column_constraints.push_back(std::move(current));
		}
	}
	std::vector<std::vector<unsigned int>> row_number_constraints;
	for(std::size_t it=0; it<9*9; it++)
	{
		for(std::size_t st=0; st<9; st++)
		{
			std::vector<unsigned int> current(81, 0);
			current[(it/9)*9+st]=1;
			row_number_constraints.push_back(std::move(current));
		}
	}
	std::vector<std::vector<unsigned int>> column_number_constraints;
	for(std::size_t it=0; it<9; it++)
	{
		for(std::size_t st=0; st<9*9; st++)
		{
			std::vector<unsigned int> current(81, 0);
			current[st]=1;
			column_number_constraints.push_back(std::move(current));
		}
	}
	std::vector<std::vector<unsigned int>> box_number_constraints;
	for(std::size_t it=0; it<3; it++)
	{
		for(std::size_t st=0; st<3; st++)
		{
			for(std::size_t kt=0; kt<3; kt++)
			{
				for(std::size_t gt=0; gt<9; gt++)
				{
					std::vector<unsigned int> current(81, 0);
					current[it*9*3+gt]=1;
					box_number_constraints.push_back(std::move(current));
				}
			}
			for(std::size_t kt=0; kt<3; kt++)
			{
				for(std::size_t gt=0; gt<9; gt++)
				{
					std::vector<unsigned int> current(81, 0);
					current[it*9*3+1*9+gt]=1;
					box_number_constraints.push_back(std::move(current));
				}
			}
			for(std::size_t kt=0; kt<3; kt++)
			{
				for(std::size_t gt=0; gt<9; gt++)
				{
					std::vector<unsigned int> current(81, 0);
					current[it*9*3+2*9+gt]=1;
					box_number_constraints.push_back(std::move(current));
				}
			}
		}
	}
	std::vector<std::vector<unsigned int>> result;
	for(std::size_t it=0; it<box_number_constraints.size(); it++)
	{
		std::vector<unsigned int> current;
		for(auto& st: row_column_constraints[it])
		{
			current.push_back(st);
		}
		for(auto& st: row_number_constraints[it])
		{
			current.push_back(st);
		}
		for(auto& st: column_number_constraints[it])
		{
			current.push_back(st);
		}
		for(auto& st: box_number_constraints[it])
		{
			current.push_back(st);
		}
		result.push_back(std::move(current));
	}
	return result;
}

std::vector<std::vector<unsigned int>> reduce(std::vector<std::vector<unsigned int>>& problem)
{
	auto matrix=generate_sudoku_problem();
	std::vector<std::size_t> rows_to_remove;
	for(std::size_t it=0; it<problem.size(); it++)
	{
		for(std::size_t st=0; st<problem.begin()->size(); st++)
		{
			if(problem[it][st]==0)
			{
				continue;
			}
			auto row=it;
			auto column=st;
			auto number=problem[it][st];
			auto remove_row=row*9*9+column*9;
			std::size_t rows_removed=0;
			for(std::size_t kt=0; kt<9; kt++)
			{
				if([&matrix, remove_row, number]()
				{
					for(std::size_t it=0; it<81; it++)
					{
						if(matrix[remove_row][it]==number)
						{
							return false;
						}
					}
					return true;
				}())
				{
					rows_to_remove.push_back(remove_row);
				}
				remove_row++;
			}
		}
	}
	for(long long int it=rows_to_remove.size()-1; it>=0; it--)
	{
		matrix.erase(matrix.begin()+rows_to_remove[it]);
	}
	return matrix;
}

void dump_matrix(std::vector<std::vector<unsigned int>>& matrix)
{
	std::ofstream ofs("matrix.log");
	for(auto& it: matrix)
	{
		for(auto& st: it)
		{
			if(st!=0)
			{
				ofs<<st;
			}
			else
			{
				ofs<<" ";
			}
		}
		ofs<<std::endl;
	}
}

void test_case()
{
	std::vector<std::vector<std::vector<unsigned int>>> tests=
	{
		std::vector<std::vector<unsigned int>>(
		{
			{0, 0, 1, 0, 1, 1, 0},
			{1, 0, 0, 1, 0, 0, 1},
			{0, 1, 1, 0, 0, 1, 0},
			{1, 0, 0, 1, 0, 0, 0},
			{0, 1, 0, 0, 0, 0, 1},
			{0, 0, 0, 1, 1, 0, 1}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{1, 0, 0, 1, 0, 0, 1},
			{1, 0, 0, 1, 0, 0, 0},
			{0, 0, 0, 1, 1, 0, 1},
			{0, 0, 1, 0, 1, 1, 0},
			{0, 1, 1, 0, 0, 1, 1},
			{0, 1, 0, 0, 0, 0, 1}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{0},
			{1}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{1},
			{1}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{1},
			{0}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{0, 1},
			{1, 0}
		}),
		std::vector<std::vector<unsigned int>>(
		{
			{0, 0, 0, 0, 1},
			{1, 1, 1, 0, 0},
			{0, 0, 0, 1, 0}
		})
	};
	//0 3 4
	//1 3 5
	//1
	//0 / 1
	//0
	//0 1
	//0 1 2
	for(auto& problem: tests)
	{
		dlx::solver solver;
		solver.apply_problem(problem);
		auto solutions=solver.search();
		for(auto& solution: solutions)
		{
			for(auto& it: solution)
			{
				std::cout<<it<<" ";
			}
			std::cout<<std::endl;
		}
		std::cout<<std::endl;
	}
	return;
}

std::vector<std::vector<unsigned int>> parse_solution(std::vector<std::size_t>& solution, 
	std::vector<std::vector<unsigned int>>& problem)
{
	std::vector<std::vector<unsigned int>> result;
	std::vector<unsigned int> current_row;
	for(std::size_t it=0; it<solution.size(); it++)
	{
		if(it!=0 && it%9==0)
		{
			//I wish I could move here
			result.push_back(current_row);
			current_row.clear();
		}
		//some edgy stuff, what it find returns end()? monkaS
		current_row.push_back(*std::find_if(problem[solution[it]].begin(), problem[solution[it]].end(), [](auto element)
		{
			return element!=0;
		}));
	}
	//can move here
	result.push_back(std::move(current_row));
	return result;
}

int main(int argc, char** argv)
{
	#if defined(_MSC_VER) && defined(_DEBUG)
		_CrtSetDbgFlag(_CRTDBG_ALLOC_MEM_DF|_CRTDBG_LEAK_CHECK_DF);
	#endif

	/*
	std::vector<std::vector<unsigned int>> problem=
	{
		{5, 3, 0, 0, 7, 0, 0, 0, 0},
		{6, 0, 0, 1, 9, 5, 0, 0, 0},
		{0, 9, 8, 0, 0, 0, 0, 6, 0},
		{8, 0, 0, 0, 6, 0, 0, 0, 3},
		{4, 0, 0, 8, 0, 3, 0, 0, 1},
		{7, 0, 0, 0, 2, 0, 0, 0, 6},
		{0, 6, 0, 0, 0, 0, 2, 8, 0},
		{0, 0, 0, 4, 1, 9, 0, 0, 5},
		{0, 0, 0, 0, 8, 0, 0, 7, 9}
	};
	*/
	
	std::vector<std::vector<unsigned int>> problem=
	{
		{8, 0, 0, 0, 0, 0, 0, 0, 0},
		{0, 0, 3, 6, 0, 0, 0, 0, 0},
		{0, 7, 0, 0, 9, 0, 2, 0, 0},
		{0, 5, 0, 0, 0, 7, 0, 0, 0},
		{0, 0, 0, 0, 4, 5, 7, 0, 0},
		{0, 0, 0, 1, 0, 0, 0, 3, 0},
		{0, 0, 1, 0, 0, 0, 0, 6, 8},
		{0, 0, 8, 5, 0, 0, 0, 1, 0},
		{0, 9, 0, 0, 0, 0, 4, 0, 0}
	};
	
	/*
	std::vector<std::vector<unsigned int>> problem=
	{
		{0, 0, 0, 2, 1, 0, 0, 0, 0},
		{0, 0, 7, 3, 0, 0, 0, 0, 0},
		{0, 5, 8, 0, 0, 0, 0, 0, 0},
		{4, 3, 0, 0, 0, 0, 0, 0, 0},
		{2, 0, 0, 0, 0, 0, 0, 0, 8},
		{0, 0, 0, 0, 0, 0, 0, 7, 6},
		{0, 0, 0, 0, 0, 0, 2, 5, 0},
		{0, 0, 0, 0, 0, 7, 3, 0, 0},
		{0, 0, 0, 0, 9, 8, 0, 0, 0},
	};
	*/
	std::cout<<"Initial problem:"<<std::endl;
	for(auto& row: problem)
	{
		for(auto& element: row)
		{
			std::cout<<element<<" ";
		}
		std::cout<<std::endl;
	}
	std::cout<<"Applying problem"<<std::endl;
	auto timer=std::chrono::high_resolution_clock::now();
	auto reduced_problem=reduce(problem);
	auto time_diff=std::chrono::high_resolution_clock::now()-timer;
	auto duration=std::chrono::duration_cast<std::chrono::milliseconds>(time_diff);
	std::cout<<"Time (apply): "<<duration.count()<<" ms"<<std::endl;
	std::cout<<"Solving problem"<<std::endl;
	timer=std::chrono::high_resolution_clock::now();
	dlx::solver solver(reduced_problem);
	auto solutions=solver.search();
	std::cout<<"Solutions found: "<<solutions.size()<<std::endl;
	time_diff=std::chrono::high_resolution_clock::now()-timer;
	duration=std::chrono::duration_cast<std::chrono::milliseconds>(time_diff);
	std::cout<<"Time (solve): "<<duration.count()<<" ms"<<std::endl;
	for(auto& solution: solutions)
	{
		auto solved=parse_solution(solution, reduced_problem);
		for(auto& row: solved)
		{
			for(auto& element: row)
			{
				std::cout<<element<<" ";
			}
			std::cout<<std::endl;
		}
	}
	return 0;
}

xalg.hpp

#ifndef XALG_HPP
#define XALG_HPP

#include <vector>
#include <cassert>
#include <iostream>
#include <algorithm>

namespace xalg
{
	/*
	If A is empty, the problem is solved; terminate successfully.
	Otherwise choose a column, c (deterministically).
	Choose a row, r, such that A[r, c] = 1 (nondeterministically).
	Include r in the partial solution.
	For each j such that A[r, j] = 1,
		delete column j from matrix A;
		for each i such that A[i, j] = 1,
			delete row i from matrix A.
	Repeat this algorithm recursively on the reduced matrix A.
	*/
	bool xalg_core(std::vector<std::vector<unsigned int>>& problem, std::vector<bool> available_rows, 
		std::vector<bool> available_columns, std::vector<std::size_t>& partial_solution, 
		std::vector<std::vector<std::size_t>>& solutions)
	{
		if(std::find(available_columns.begin(), available_columns.end(), true)==available_columns.end())
		{
			std::cout<<"found solution"<<std::endl;
			for(auto& it: partial_solution)
			{
				std::cout<<it<<" ";
			}
			std::cout<<std::endl;
			return true;
		}
		//choose column
		auto count_rows_at_column=[&problem](std::size_t col)
		{
			std::size_t current_count=0;
			for(std::size_t it=0; it<problem.size(); it++)
			{
				if(problem[it][col]!=0)
				{
					current_count++;
				}
			}
			return current_count;
		};
		std::size_t column=[&available_columns]()-> std::size_t
		{
			for(std::size_t it=0; it<available_columns.size(); it++)
			{
				if(available_columns[it])
				{
					return it;
				}
			}
			return 0;
		}();
		std::size_t column_count=count_rows_at_column(column);
		for(std::size_t it=0; it<problem.begin()->size(); it++)
		{
			if(available_columns[it] && count_rows_at_column(it)<column_count)
			{
				column=it;
			}
		}
		std::vector<std::size_t> rows;
		//choose rows, such as problem[row][column]=1
		for(std::size_t it=0; it<problem.size(); it++)
		{
			if(available_rows[it] && problem[it][column]!=0)
			{
				rows.push_back(it);
			}
		}
		for(auto& current_row: rows)
		{
			std::vector<bool> current_available_rows=available_rows;
			std::vector<bool> current_available_columns=available_columns;
			for(std::size_t it=0; it<problem.begin()->size(); it++)
			{
				if(current_available_columns[it] && problem[current_row][it]!=0)
				{
					current_available_columns[it]=false;
					for(std::size_t st=0; st<problem.size(); st++)
					{
						if(current_available_rows[st] && problem[st][it]!=0)
						{
							current_available_rows[st]=false;
						}
					}
				}
			}
			std::vector<std::size_t> current_partial_solution=partial_solution;
			current_partial_solution.push_back(current_row);
			if(xalg_core(problem, std::move(current_available_rows), std::move(current_available_columns), 
				current_partial_solution, solutions))
			{
				std::sort(current_partial_solution.begin(), current_partial_solution.end());
				solutions.push_back(std::move(current_partial_solution));
			}
		}
		return false;
	}

	inline auto xalg(std::vector<std::vector<unsigned int>>& problem)
	{
		#ifdef _DEBUG
		for(auto& it: problem)
		{
			if(it.size()!=problem.begin()->size())
			{
				throw std::runtime_error("problem column size mismatch!");
			}
		}
		#endif
		//shit ton of stack memory monkaS
		std::vector<bool> available_columns(problem.begin()->size(), true);
		std::vector<bool> available_rows(problem.size(), true);
		std::vector<std::size_t> partial_solution;
		std::vector<std::vector<std::size_t>> solutions;
		//actually ignore the return value omegaLUL
		auto res=xalg_core(problem, std::move(available_rows), std::move(available_columns), 
			partial_solution, solutions);
		return solutions;
	}

	inline auto xalg(std::vector<std::vector<unsigned int>>& problem, 
		std::vector<std::size_t>& partial_solution)
	{
		#ifdef _DEBUG
		for(auto& it: problem)
		{
			if(it.size()!=problem.begin()->size())
			{
				throw std::runtime_error("problem column size mismatch!");
			}
		}
		#endif
		//shit ton of stack memory monkaS
		std::vector<bool> available_columns(problem.begin()->size(), true);
		std::vector<bool> available_rows(problem.size(), true);
		std::vector<std::vector<std::size_t>> solutions;
		//actually ignore the return value omegaLUL
		auto res=xalg_core(problem, std::move(available_rows), std::move(available_columns), 
			partial_solution, solutions);
		return solutions;
	}
}

#endif